ggael wrote:For dense symmetric matrices, now Eigen has quite good performances.

But here you are interested in sparse matrices that is very different. So far we don't provide any support for that. So you have to use another lib to do the actual decomposition. BLAS/LAPACK won't help you here because they are for dense objects only. Nevertheless you can still assemble your SparseMatrix using Eigen's API, and then send it the other library. Eigen uses a standard storage format so that should not be a problem.

A popular lib for that is SLEPc: http://www.grycap.upv.es/slepc.

Code: Select all

#include <arpack++/arssym.h>
#include <vector>
#include <eigen2/Eigen/Eigen>

class TMatrixForArpackpp {
private:
    Eigen::SparseMatrix< double > &M;
public:
    TMatrixForArpackpp( Eigen::SparseMatrix< double > &_M ) :
        M(_M) {
        ;
    };
    void MultMv( double *v, double *w ) {
        Eigen::Map< Eigen::VectorXd > v_map = Eigen::VectorXd::Map( v, M.rows() );
        Eigen::Map< Eigen::VectorXd > w_map = Eigen::VectorXd::Map( w, M.cols() );
        w_map = M.marked<Eigen::SelfAdjoint|Eigen::LowerTriangular>()*v_map;
    }
};


void Lanczos_symm( Eigen::SparseMatrix< double > &M, std::vector< double > &eigenValues, std::vector< Eigen::VectorXd > &eigenVectors ) {
    int n = M.cols();
    TMatrixForArpackpp M_arp( M );
    ARSymStdEig< double, TMatrixForArpackpp > prob( n, n-1, &M_arp, &TMatrixForArpackpp::MultMv, "SM" ); //"SM" to get smallest eigenvalues
    prob.FindEigenvectors();
    //prob.FindEigenvalues();
    int m = prob.ConvergedEigenvalues();
    eigenValues.resize(m);
    eigenVectors.resize(m);
    eigenValues.assign( prob.RawEigenvalues(), prob.RawEigenvalues()+m );
    for( int i=0; i<m; i++ )
        eigenVectors[i] = Eigen::VectorXd::Map( prob.RawEigenvector(i), n);
    /*eigenValues = *( prob.StlEigenvalues() );
    eigenVectors = *( prob.StlEigenvectors() );*/
}


Code: Select all

-lblas -lgfortran -lgfortranbegin -lnsl -larpack -larpack++