Re: [eigen] banded matrices in Eigen

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

In practice this is a bit more tricky. Some months ago Kolja and I came up with the attached example for matrix-free solving. You can easily adapt it to your case. In the future, I hope we'll be able to make this process much simpler since in theory you should only have to provide an operator*.

gael

On Thu, Feb 13, 2014 at 1:36 PM, Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:

On 13.02.2014 12:26, Laura Flores wrote:

If I could create my own class, derived from one of the existing ones
(possibly the Diagonal one?), I could add some methods that suit the
problem I am working on. Therefore, my question would be whether there are
any suggestions
for that.

Looking at the current implementation of ConjugateGradient, it would theoretically be sufficient to write a class which efficiently implements:
  VectorXd res = rhs - mat * x;
and
  tmp.noalias() = mat * p;

Basically, you need to make
  YourClass::operator*(const VectorType& v);
return something similar to GeneralProduct.

The same could be done for the preconditioner and
  YourPreconditioner::solve(...)

I'm afraid neither is entirely trivial, nor guaranteed to be future compatible.

Christoph


--
----------------------------------------------
Dipl.-Inf., Dipl.-Math. Christoph Hertzberg
Cartesium 0.049
Universität Bremen
Enrique-Schmidt-Straße 5
28359 Bremen

Tel: +49 (421) 218-64252
----------------------------------------------

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

class MatrixReplacement;
template<typename Rhs> class MatrixReplacement_ProductReturnType;

class MatrixReplacement : public Eigen::internal::traits<Eigen::SparseMatrix<double> > {
public:
  typedef double Scalar;
  typedef double RealScalar;  

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

  template<typename Rhs>
  MatrixReplacement_ProductReturnType<Rhs> operator*(const Eigen::MatrixBase<Rhs>& x) const {
    return MatrixReplacement_ProductReturnType<Rhs>(*this, x.derived());
  }
  
  // ConjugateGradient calls selfadjointView to work on the lower of upper half of the matrix,
  // or on the entire matrix if Eigen::Lower|Eigen::Upper is specified.
  // In the future, if Eigen::Lower|Eigen::Upper is specified, then this function should not be called at all.
  template<int UpLo>
  const MatrixReplacement& selfadjointView() const
  {
    assert(UpLo==(Eigen::Lower|Eigen::Upper) && "considering only the lower or upper triangular half is not supported");
    return *this;
  }

};

// The proxy class representing the product of a MatrixReplacement with a MatrixBase<>
template<typename Rhs>
class MatrixReplacement_ProductReturnType : public Eigen::ReturnByValue<MatrixReplacement_ProductReturnType<Rhs> > {
public:
  typedef MatrixReplacement::Index Index;
  
  // The ctor store references to the matrix and right-hand-side object (usually a vector).
  MatrixReplacement_ProductReturnType(const MatrixReplacement& matrix, const Rhs& rhs)
    : m_matrix(matrix), m_rhs(rhs)
  {}
  
  Index rows() const { return m_matrix.rows(); }
  Index cols() const { return m_rhs.cols(); }

  // This function is automatically called by Eigen. It must evaluate the product of matrix * rhs into y.
  template<typename Dest>
  void evalTo(Dest& y) const
  {
    y.setZero(4);

    y(0) += 2 * m_rhs(0); y(1) += 1 * m_rhs(0);
    y(0) += 1 * m_rhs(1); y(1) += 2 * m_rhs(1); y(2) += 1 * m_rhs(1);
    y(1) += 1 * m_rhs(2); y(2) += 2 * m_rhs(2); y(3) += 1 * m_rhs(2);
    y(2) += 1 * m_rhs(3); y(3) += 2 * m_rhs(3);
  }

protected:
  const MatrixReplacement& m_matrix;
  typename Rhs::Nested m_rhs;
};

namespace Eigen {
namespace internal {
  template <typename Rhs>
  struct traits<MatrixReplacement_ProductReturnType<Rhs> > {
    // The equivalent plain objet type of the product. This type is used if the product needs to be evaluated into a temporary.
    typedef Eigen::Matrix<typename Rhs::Scalar, Eigen::Dynamic, Rhs::ColsAtCompileTime> ReturnType;
  };
}
}


/*****/


namespace Eigen { 

template <typename _Scalar>
class JacobiPreconditioner
{
    typedef _Scalar Scalar;
    typedef Matrix<Scalar,Dynamic,1> Vector;
    typedef typename Vector::Index Index;

  public:
    // this typedef is only to export the scalar type and compile-time dimensions to solve_retval
    typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;

    JacobiPreconditioner() : m_isInitialized(false) {}

    void setInvDiag(const VectorXd &invdiag) {
      m_invdiag=invdiag;
      m_isInitialized=true;
    }

    Index rows() const { return m_invdiag.size(); }
    Index cols() const { return m_invdiag.size(); }
    
    template<typename MatType>
    JacobiPreconditioner& analyzePattern(const MatType& ) { return *this; }
    
    template<typename MatType>
    JacobiPreconditioner& factorize(const MatType& mat) { return *this; }
    
    template<typename MatType>
    JacobiPreconditioner& compute(const MatType& mat) { return *this; }

    template<typename Rhs, typename Dest>
    void _solve(const Rhs& b, Dest& x) const
    {
      x = m_invdiag.array() * b.array() ;
    }

    template<typename Rhs> inline const internal::solve_retval<JacobiPreconditioner, Rhs>
    solve(const MatrixBase<Rhs>& b) const
    {
      eigen_assert(m_isInitialized && "JacobiPreconditioner is not initialized.");
      eigen_assert(m_invdiag.size()==b.rows()
                && "JacobiPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
      return internal::solve_retval<JacobiPreconditioner, Rhs>(*this, b.derived());
    }

  protected:
    Vector m_invdiag;
    bool m_isInitialized;
};

namespace internal {

template<typename _MatrixType, typename Rhs>
struct solve_retval<JacobiPreconditioner<_MatrixType>, Rhs>
  : solve_retval_base<JacobiPreconditioner<_MatrixType>, Rhs>
{
  typedef JacobiPreconditioner<_MatrixType> Dec;
  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)

  template<typename Dest> void evalTo(Dest& dst) const
  {
    dec()._solve(rhs(),dst);
  }
};

}
}


/*****/


int main()
{
   Eigen::MatrixXd A(4, 4);
   Eigen::VectorXd b(4), x;

   A << 2, 1, 0, 0, 1, 2, 1, 0, 0, 1, 2, 1, 0, 0, 1, 2;
   b << 1, 1, 1, 1;

   // solve Ax = b using CG with matrix-free version:
   Eigen::ConjugateGradient < MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::JacobiPreconditioner<double> > cg;
   
   MatrixReplacement M;

   Eigen::VectorXd invdiag(4);
   invdiag << 1./3., 1./4., 1./4., 1./3.;
   
   cg.preconditioner().setInvDiag(invdiag);

   cg.compute(M);
   x = cg.solve(b);

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

}