1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
#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 mp_mat->rows(); }
Index cols() const { return mp_mat->cols(); }
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() : mp_mat(0) {}
void attachMyMatrix(const SparseMatrix<double> &mat) {
mp_mat = &mat;
}
const SparseMatrix<double> my_matrix() const { return *mp_mat; }
private:
const SparseMatrix<double> *mp_mat;
};
// 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");
// Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
// but let's do something fancier (and less efficient):
for(Index i=0; i<lhs.cols(); ++i)
dst += rhs(i) * lhs.my_matrix().col(i);
}
};
}
}
int main()
{
int n = 10;
Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
S = S.transpose()*S;
MatrixReplacement A;
A.attachMyMatrix(S);
Eigen::VectorXd b(n), x;
b.setRandom();
// 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;
}
{
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;
}
{
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;
}
{
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;
}
{
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;
}
}
|
