Please, help us to better know about our user community by answering the following short survey:
Eigen  3.4.99 (git rev 199c5f2b47eb1f8e5a2d20e60f07e97cd95a6ba6)
Matrix-free solvers

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:

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;
namespace Eigen {
namespace internal {
// MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
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> {
// 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>
// Custom API:
MatrixReplacement() : mp_mat(0) {}
void attachMyMatrix(const SparseMatrix<double> &mat) {
mp_mat = &mat;
const SparseMatrix<double> my_matrix() const { return *mp_mat; }
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;
S = S.transpose()*S;
MatrixReplacement A;
Eigen::VectorXd b(n), x;
// Solve Ax = b using various iterative solver with matrix-free version:
x = cg.solve(b);
std::cout << "CG: #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
x = bicg.solve(b);
std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
x = gmres.solve(b);
std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
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;
x = minres.solve(b);
std::cout << "MINRES: #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
A bi conjugate gradient stabilized solver for sparse square problems.
Definition: BiCGSTAB.h:159
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:159
static const RandomReturnType Random()
Definition: Random.h:113
Derived & derived()
Definition: EigenBase.h:46
RealScalar error() const
Definition: IterativeSolverBase.h:305
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:238
Index iterations() const
Definition: IterativeSolverBase.h:296
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Derived & setRandom(Index size)
Definition: Random.h:151
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:75
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:88
Namespace containing all symbols from the Eigen library.
Definition: Core:134
The Index type as used for the API.
Definition: Meta.h:74
const Product< MatrixDerived, PermutationDerived, AliasFreeProduct > operator*(const MatrixBase< MatrixDerived > &matrix, const PermutationBase< PermutationDerived > &permutation)
Definition: PermutationMatrix.h:515
const int Dynamic
Definition: Constants.h:22
Definition: EigenBase.h:30


CG:       #iterations: 16, estimated error: 1.48963e-16
BiCGSTAB: #iterations: 15, estimated error: 1.23108e-16
GMRES:    #iterations: 10, estimated error: 0
DGMRES:   #iterations: 20, estimated error: 2.1154e-16
MINRES:   #iterations: 17, estimated error: 9.30361e-17