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);
}
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