1 #include <iostream>
2 #include <Eigen/Core>
3 #include <Eigen/Dense>
4 #include <Eigen/IterativeLinearSolvers>
5 #include <unsupported/Eigen/IterativeSolvers>
6
7 class MatrixReplacement;
8 using Eigen::SparseMatrix;
9
10 namespace Eigen {
11 namespace 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.
21 class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
22 public:
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
rows() const33 Index rows() const { return mp_mat->rows(); }
cols() const34 Index cols() const { return mp_mat->cols(); }
35
36 template<typename Rhs>
operator *(const Eigen::MatrixBase<Rhs> & x) const37 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:
MatrixReplacement()42 MatrixReplacement() : mp_mat(0) {}
43
attachMyMatrix(const SparseMatrix<double> & mat)44 void attachMyMatrix(const SparseMatrix<double> &mat) {
45 mp_mat = &mat;
46 }
my_matrix() const47 const SparseMatrix<double> my_matrix() const { return *mp_mat; }
48
49 private:
50 const SparseMatrix<double> *mp_mat;
51 };
52
53
54 // Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
55 namespace Eigen {
56 namespace 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>
scaleAndAddToEigen::internal::generic_product_impl65 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
main()81 int 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