5
votes

I want to solve a linear algebraic equation Ax = b using Eigen solvers. In my case, A is a complex sparse matrix(26410*26410), b is a real vector (26410*1).

I use mex file in MATLAB to map the sparse matrix A and vector b to Eigen accepted format. The reason why I use Eigen solver is to hope it would be faster than solving directly in MATLAB using x = A\b.

However, after tried LDLT, SparseLU, CG and BiCGSTAB, I found the results are not very satisfying:

LDLT takes 1.462s with norm(A*x - b)/norm(b) = 331; SparseLU takes 37.994s with 1.5193e-4; BiCGSTAB takes 95.217s with 4.5977e-4; On the contrast, directly use x = A\b in MATLAB consumes 13.992s with norm of the error 2.606e-5.

I know it is a little stupid and also time consuming to map the sparse matrix A and vector b in MATLAB workspace to Eigen. But I am wondering whether the results I got are the best results which Eigen can give? Anyone can give me some pointers? Should I try some other linear equation solvers? Thanks a lot in advance! The following is the main part of codes.

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

    //input vars

    //temp var
    size_t nrows;

    //output vars
    //double *x;

    //GetData
    /* check inputs 
    ...*/

    //"mxArray2SCM" is a sub-function for map the complex sparse matrix in Eigen
    SpCMat A = mxArray2SCM(prhs[0]);
    //SpMat A = mxArray2SM(prhs[0]);

    //"mxArray2ECV" is a sub-function for map the real vector in Eigen
    Eigen::VectorXcd b = mxArray2ECV(prhs[1]);
    //Eigen::VectorXd b = mxArray2EV(prhs[1]);   
    nrows = b.size();
    //Computation
    Eigen::VectorXcd x(nrows);
    //SparseLU<SparseMatrix<CD> > solver;
    BiCGSTAB<SparseMatrix<CD>,IncompleteLUT<CD> > BiCG;
    //BiCG.preconditioner().setDroptol(0.001);
    BiCG.compute(A);
    if(BiCG.info()!=Success){
        //decomposition failed
        return;
    }
    x = BiCG.solve(b);

    //Output results
    plhs[0] = ECV2mxArray(x);   
}
1
Do you have compiler optimizations enabled?Baum mit Augen
Yes. I use "mex -O -largeArrayDims xxx.cpp" command to compile. By the way, version of MatLAB is R2016b and compiler is VS2015 Professional.@Baum mit AugenL.Ferron
Why would you guess that Eigen would be faster than Matlab? Matlab relies on highly optimized and multi-threaded external libraries such as UmfPack and MKL. Eigen can also make use of them. For instance, if you use Eigen::UmfpackLU and link to UmfPack, you will probably get the same performance at MatLab.ggael
Since I have to do many iterations for my optimization, and in each iteration a linear algebra equation is solved, then I expect to find a way to reduce the computation time for each iteration. Thus, I tried Eigen to expect it could be faster. It seems that I was wrong. As you said, UmfPack maybe can provide the same performance as MatLab. Could you give me some tips for using that in Eigen? Thanks a lot! @ggaelL.Ferron
Based on what the matrix you describe, Matlab’s A\b will use a sparse LU solver. It’s not surprising LDLT fails (provides non-epsilon error) because Eigen’s LDLT is for square-Hermitian matrixes only.Ahmed Fasih

1 Answers

0
votes

Have you considered using PetSc for Krylov solvers or SLEPc to compute eigenvalues?

Make sure you analyze the eigenspectrum before using a specific Krylov solver (CG works only for symmetric positive definite matrices).

PETSc has quite a few solvers that you can try out based on your eigenspectrum.

You can check Y. Saad's book on how these solvers work.

If your matrix is not symmetric positive definite GMRES is a good option.