1#include <iostream>
2#include <Eigen/Core>
3#include <Eigen/Dense>
4#include <Eigen/IterativeLinearSolvers>
5#include <unsupported/Eigen/IterativeSolvers>
6
7class MatrixReplacement;
8using Eigen::SparseMatrix;
9
10namespace Eigen {
11namespace internal {
12  // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
13  template<>
14  struct traits<MatrixReplacement> :  public Eigen::internal::traits<Eigen::SparseMatrix<double> >
15  {};
16}
17}
18
19// Example of a matrix-free wrapper from a user type to Eigen's compatible type
20// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix.
21class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
22public:
23  // Required typedefs, constants, and method:
24  typedef double Scalar;
25  typedef double RealScalar;
26  typedef int StorageIndex;
27  enum {
28    ColsAtCompileTime = Eigen::Dynamic,
29    MaxColsAtCompileTime = Eigen::Dynamic,
30    IsRowMajor = false
31  };
32
33  Index rows() const { return mp_mat->rows(); }
34  Index cols() const { return mp_mat->cols(); }
35
36  template<typename Rhs>
37  Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
38    return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
39  }
40
41  // Custom API:
42  MatrixReplacement() : mp_mat(0) {}
43
44  void attachMyMatrix(const SparseMatrix<double> &mat) {
45    mp_mat = &mat;
46  }
47  const SparseMatrix<double> my_matrix() const { return *mp_mat; }
48
49private:
50  const SparseMatrix<double> *mp_mat;
51};
52
53
54// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
55namespace Eigen {
56namespace internal {
57
58  template<typename Rhs>
59  struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
60  : generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
61  {
62    typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
63
64    template<typename Dest>
65    static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
66    {
67      // This method should implement "dst += alpha * lhs * rhs" inplace,
68      // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
69      assert(alpha==Scalar(1) && "scaling is not implemented");
70
71      // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
72      // but let's do something fancier (and less efficient):
73      for(Index i=0; i<lhs.cols(); ++i)
74        dst += rhs(i) * lhs.my_matrix().col(i);
75    }
76  };
77
78}
79}
80
81int main()
82{
83  int n = 10;
84  Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
85  S = S.transpose()*S;
86
87  MatrixReplacement A;
88  A.attachMyMatrix(S);
89
90  Eigen::VectorXd b(n), x;
91  b.setRandom();
92
93  // Solve Ax = b using various iterative solver with matrix-free version:
94  {
95    Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
96    cg.compute(A);
97    x = cg.solve(b);
98    std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
99  }
100
101  {
102    Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg;
103    bicg.compute(A);
104    x = bicg.solve(b);
105    std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
106  }
107
108  {
109    Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
110    gmres.compute(A);
111    x = gmres.solve(b);
112    std::cout << "GMRES:    #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
113  }
114
115  {
116    Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
117    gmres.compute(A);
118    x = gmres.solve(b);
119    std::cout << "DGMRES:   #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
120  }
121
122  {
123    Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres;
124    minres.compute(A);
125    x = minres.solve(b);
126    std::cout << "MINRES:   #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
127  }
128}
129