Despite what ggael says here, implementation of implicit matrices is slightly more tedious than simply defining three functions. Here is a simple example that compiles on GCC 4.8.2 using Eigen 3.2.2.

Code: Select all

#include <iostream>
#include <complex>
#include <Eigen/Core>
#include <Eigen/IterativeLinearSolvers>


// Base class for all implicit matrices
template<
    typename _Scalar,
    int _RowsAtCompileTime = Eigen::Dynamic,
    int _ColsAtCompileTime = Eigen::Dynamic,
    int _Options = 0,
    int _MaxRowsAtCompileTime = _RowsAtCompileTime,
    int _MaxColsAtCompileTime = _ColsAtCompileTime
>
struct ImplicitMatrixBase {
    typedef _Scalar Scalar;
    typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;

    static constexpr int RowsAtCompileTime = _RowsAtCompileTime;
    static constexpr int ColsAtCompileTime = _ColsAtCompileTime;
    static constexpr int MaxRowsAtCompileTime = _MaxRowsAtCompileTime;
    static constexpr int MaxColsAtCompileTime = _MaxColsAtCompileTime;
    typedef int Index;
};


// A concrete implicit matrix
template<class Scalar>
struct Laplacian : ImplicitMatrixBase<Scalar> {
    typedef ImplicitMatrixBase<Scalar> Base;
    typedef typename Base::Index Index;
    typedef Eigen::Matrix<Scalar,Base::ColsAtCompileTime,1> Vector;

    Laplacian(Index n) : n_(n) {}

    Index rows() const { return n_; }
    Index cols() const { return n_; }

    template<class D>
    Vector operator*(const Eigen::MatrixBase<D>& x) const {
        Vector y(x.size());
        y[0] = 2*x[0] - x[1];
        for (int i = 1; i < x.size()-1; ++i) y[i] = -x[i-1] + 2*x[i] - x[i+1];
        y[x.size()-1] = -x[x.size() -2] + 2*x[x.size()-1];
        return y;
    }

private:
    Index n_;
};

int main() {
    using namespace Eigen;
    unsigned n = 100;

    Laplacian<double> lapl(n);
    VectorXd x = VectorXd::Random(n);
    VectorXd b = VectorXd::Random(n);

    BiCGSTAB<Laplacian<double>, IdentityPreconditioner> solver(lapl);
    solver.setTolerance(1e-6);
    x = solver.solveWithGuess(b,x);

    std::cout << "#iterations: " << solver.iterations() << std::endl;
    std::cout << "estimated error: " << solver.error() << std::endl;

    return 0;
}


Thank you for sharing this example.

You can use Eigen::NumTraits<Scalar>::Real to define the RealScalar type.

I posted exactly for comments like this one!
Updated my code accordingly.