Eigen Sparse Matrix example

Introduction

We have already studied that a sparse matrix is said to be a matrix in which many or most elements have zero value. Here we will see an example of how we can implement it and completely wrap a sparse matrix. Iterative solvers like ConjugateGradient and BiCGSTAB can be used in a matrix-free context. To this end, the user must provide a wrapper class inheriting EigenBase<> and implementing the given methods:

• Index rows() and Index cols(): They return the number of rows and columns, respectively.
• operator* with your datatype and an Eigen dense column vector.

Example

User-defined files

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
#include <unsupported/Eigen/IterativeSolvers>
class MatReplacement;
using Eigen::SparseMatrix;
namespace Eigen {
	namespace internal {
 		// MatrixReplacement is similar to a SparseMatrix, so let us inherit its traits:
 		template<>
 		struct traits<MatReplacement> :  public Eigen::internal::traits<Eigen::SparseM
 		atrix<double> >
 		{};
	}
}

Class Matrix replacement

This method shows an example of a matrix-free wrapper from a user-type to Eigen's compatible type.

class MatReplacement : public Eigen::EigenBase<Mat_Replacement> {
public:
 	// We are required for constants, typedefs, and method:
 	typedef double Scalar;
 	typedef double RealScalar;
 	typedef int StorageIndex;
 	enum {
   		ColsAtCompileTime = Eigen::Dynamic,
   		MaximumColsAtCompileTime = Eigen::Dynamic,
   		IsRowMajor = false
 		};
 	Index rows() const { return mat_rep->rows(); }
 	Index cols() const { return mat_rep->cols(); }
 	template<typename Rhs>
 	Eigen::Product<MatReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Ei
 	gen::MatrixBase<Rhs>& soln) const 
 	{
  		return Eigen::Product<MatReplacement,Rhs,Eigen::AliasFreeProduct>(*this, s
  		oln.derived());
 	}
 	
 	// Custom API:
 	MatReplacement() : mat_rep(0) {}
 	void attachMyMatrix(const SparseMatrix<double> &mat) 
 	{
  		mat_rep = &mat;
 	}
 	const SparseMatrix<double> matrix_with_me() const { return *mat_rep; } 
private:
	const SparseMatrix<double> *mat_rep;
};

Method for replacement

Here, the Implementation of Matrix Replacement * Eigen::DenseVector is done through a specialization of internal::generic_product_impl as given below:

namespace Eigen {
namespace internal {
	template<typename Rhs>
 	struct generic_product_impl<MatReplacement, Rhs, SparseShape, DenseShape, G
 	emvProduct> // GEMV stands for matrix-vector
 	: generic_product_impl_base<MatReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
 	{
   		typedef typename Product<MatReplacement,Rhs>::Scalar Scalar;
   		template<typename Dest>
   		static void scaleAndAddTo(Dest& dst, const MatReplacement& lhs, const Rhs
   		& rhs, const Scalar& alpha)
   		{
     		// This method implements "dst += alpha * lhs * rhs" inplace,
     		// but we dont need to worry as for iterative solvers, alpha is equal 
     		to 1 always     
	 		assert(alpha==Scalar(1) && "scaling isn't implemented");
     		EIGEN_ONLY_USED_FOR_DEBUG(alpha);
     		for(Index ind=0; ind<lhs.cols(); ++ind)
       			dst = dst+ rhs(ind) * lhs.matrix_with_me().col(ind);
   		}
 	 };
  }
}

Main method

int main()
{
 	int num = 10;
 	Eigen::SparseMatrix<double> Sp = Eigen::MatrixXd::Random(num,num).sparseVie
 	w(0.5,1);
 	Sp = Sp.transpose()*Sp;
 	MatReplacement mat1;
 	mat1.attachMyMatrix(Sp);
 	Eigen::VectorXd vec(num), soln;
 	vec.setRandom();
 	// Solve mat1*soln = vec with the help of various iterative solver with matrix-free v
 	ersion:
 	{
   		Eigen::ConjugateGradient<MatReplacement, Eigen::Lower|Eigen::Upper, Eige
   		n::IdentityPreconditioner> cg;
   		cg.compute(mat1);
   		soln= cg.solve(vec);
   		std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error is: " 
   		<< cg.error() << std::endl;
 	}
 	{
   		Eigen::BiCGSTAB<MatReplacement, Eigen::IdentityPreconditioner> bicg;
   		bicg.compute(mat1);
   		soln= bicg.solve(vec);
   		std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated erro
   		r: " << bicg.error() << std::endl;
 	}
 	{
   		Eigen::GMRES<MatReplacement, Eigen::IdentityPreconditioner> gmres;
   		gmres.compute(mat1);
   		soln= gmres.solve(vec);
   		std::cout << "GMRES:    #iterations: " << gmres.iterations() << ", estimated erro
   		r: " << gmres.error() << std::endl;
 	}
 	{
   		Eigen::DGMRES<MatReplacement, Eigen::IdentityPreconditioner> gmres;
   		gmres.compute(mat1);
   		soln= gmres.solve(vec);
   		std::cout << "DGMRES:   #iterations: " << gmres.iterations() << ", estimated err
   		or: " << gmres.error() << std::endl;
 	}
 	{
   		Eigen::MINRES<MatReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPr
   		econditioner> minres;
   		minres.compute(mat1);
   		soln= minres.solve(vec);
   		std::cout << "MINRES:   #iterations: " << minres.iterations() << ", estimated err
   		or: " << minres.error() << std::endl;
 	}
}