Iterative solvers such as ConjugateGradient and BiCGSTAB can be used in a matrix free context. To this end, user must provide a wrapper class inheriting EigenBase<> and implementing the following methods:
Index
rows()
and Index
cols()
: returns number of rows and columns respectively
operator*
with your type and an Eigen dense column vector (its actual implementation goes in a specialization of the internal::generic_product_impl class)
Eigen::internal::traits<>
must also be specialized for the wrapper type.
Here is a complete example wrapping an Eigen::SparseMatrix:
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
#include <unsupported/Eigen/IterativeSolvers>
class MatrixReplacement;
template<>
struct traits<MatrixReplacement> : public Eigen::internal::traits<Eigen::SparseMatrix<double> >
{};
}
}
public:
typedef double Scalar;
typedef double RealScalar;
typedef int StorageIndex;
enum {
IsRowMajor = false
};
Index rows()
const {
return mp_mat->rows(); }
Index cols()
const {
return mp_mat->cols(); }
template<typename Rhs>
}
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;
};
template<typename Rhs>
struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct>
: 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)
{
assert(alpha==Scalar(1) && "scaling is not implemented");
EIGEN_ONLY_USED_FOR_DEBUG(alpha);
for(
Index i=0; i<lhs.cols(); ++i)
dst += rhs(i) * lhs.my_matrix().col(i);
}
};
}
}
int main()
{
int n = 10;
S = S.transpose()*S;
MatrixReplacement A;
A.attachMyMatrix(S);
{
std::cout <<
"CG: #iterations: " << cg.
iterations() <<
", estimated error: " << cg.
error() << std::endl;
}
{
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;
}
}
Output:
CG: #iterations: 20, estimated error: 8.86333e-14
BiCGSTAB: #iterations: 20, estimated error: 2.10809e-15
GMRES: #iterations: 10, estimated error: 0
DGMRES: #iterations: 20, estimated error: 1.10455e-28
MINRES: #iterations: 20, estimated error: 2.94473e-14