From 449ec69bab227bae4bc41959c4835973c382f5f3 Mon Sep 17 00:00:00 2001
From: luo_longshan <luo_longshan@iapcm.ac.cn>
Date: Tue, 2 Jan 2024 01:12:59 +0800
Subject: [PATCH] add 1D Poission Eq. matrix-free example

---
 doc/examples/matrixfree.cpp | 145 ++++++++++++++++++++++++++++++++++++
 1 file changed, 145 insertions(+)
 create mode 100644 doc/examples/matrixfree.cpp

diff --git a/doc/examples/matrixfree.cpp b/doc/examples/matrixfree.cpp
new file mode 100644
index 000000000..4a5b8b3fb
--- /dev/null
+++ b/doc/examples/matrixfree.cpp
@@ -0,0 +1,145 @@
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/IterativeSolvers>
+
+class MatrixReplacement;
+using Eigen::SparseMatrix;
+
+namespace Eigen {
+namespace internal {
+  // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
+  template<>
+  struct traits<MatrixReplacement> :  public Eigen::internal::traits<Eigen::SparseMatrix<double> >
+  {};
+}
+}
+
+// Example of a matrix-free wrapper from a user type to Eigen's compatible type
+// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix.
+class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
+public:
+  // Required typedefs, constants, and method:
+  typedef double Scalar;
+  typedef double RealScalar;
+  typedef int StorageIndex;
+  enum {
+    ColsAtCompileTime = Eigen::Dynamic,
+    MaxColsAtCompileTime = Eigen::Dynamic,
+    IsRowMajor = false
+  };
+
+  Index rows() const { return _n; }
+  Index cols() const { return _n; }
+
+  template<typename Rhs>
+  Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
+    return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
+  }
+
+  // Custom API:
+  MatrixReplacement(int n){ _n = n; }
+
+private:
+  int _n;
+};
+
+
+// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
+namespace Eigen {
+namespace internal {
+
+  template<typename Rhs>
+  struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
+  : generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
+  {
+    typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
+
+    template<typename Dest>
+    static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
+    {
+      // This method should implement "dst += alpha * lhs * rhs" inplace,
+      // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
+      assert(alpha==Scalar(1) && "scaling is not implemented");
+      EIGEN_ONLY_USED_FOR_DEBUG(alpha);
+
+      // LLS:
+      // the matrix comes from 1D Poisson Eq. -u_{xx} = f
+      // FD scheme is: 2 U_{i} - U_{i-1}-U_{i+1} = f_i, i=0,...,n-1
+      // BC condition: U_{-1}=0, U_{n}=0
+      int n = lhs.cols();
+      Scalar *dst_ptr = dst.data();
+      const Scalar *rhs_ptr = rhs.data();
+
+      for(int i=0;i<n;i++) {
+        Scalar v=rhs_ptr[i];
+        Scalar vp,vm;
+        if(i==n-1) {
+          vp = 0.0;
+        } else {
+          vp = rhs_ptr[i+1];
+        }
+        if(i==0) {
+          vm = 0.0;
+        } else {
+          vm = rhs_ptr[i-1];
+        }
+        dst_ptr[i]=2.0*v-vp-vm;
+      }
+    }
+  };
+
+}
+}
+
+int main()
+{
+  int n = 4;
+
+  MatrixReplacement A(n);
+
+  Eigen::VectorXd b(n), x;
+  b.setOnes();
+
+  // Solve Ax = b using various iterative solver with matrix-free version:
+  {
+    Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
+    cg.compute(A);
+    x = cg.solve(b);
+    std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
+    for(int i=0;i<n;i++) std::cout << "x["<<i<<"] = "<<x[i]<<"\n";
+  }
+
+  {
+    Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg;
+    bicg.compute(A);
+    x = bicg.solve(b);
+    std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
+    for(int i=0;i<n;i++) std::cout << "x["<<i<<"] = "<<x[i]<<"\n";
+  }
+
+  {
+    Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
+    gmres.compute(A);
+    x = gmres.solve(b);
+    std::cout << "GMRES:    #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
+    for(int i=0;i<n;i++) std::cout << "x["<<i<<"] = "<<x[i]<<"\n";
+  }
+
+  {
+    Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
+    gmres.compute(A);
+    x = gmres.solve(b);
+    std::cout << "DGMRES:   #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
+    for(int i=0;i<n;i++) std::cout << "x["<<i<<"] = "<<x[i]<<"\n";
+  }
+
+  {
+    Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres;
+    minres.compute(A);
+    x = minres.solve(b);
+    std::cout << "MINRES:   #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
+    for(int i=0;i<n;i++) std::cout << "x["<<i<<"] = "<<x[i]<<"\n";
+  }
+}