Dear all,
Could you please help me with my problem?
I am trying to solve a linear equation. When the matrices are small (200x200), my program works well. But when I try with large matrices (5000x5000), the program runs forever and doesn't give me any result. Here is some code I write:

int np = 5000;
Eigen::MatrixXf U0 = Eigen::MatrixXf::Zero(np, 1);
U0(0, 0) = 1;

_step = 0.125;
Eigen::MatrixXf A = _B - _step * _Q; // _B and _Q are matrices 5000x5000 of float elements.
Eigen::MatrixXf Ut;

int _time = 5;
int ns = ceil(_time / _step); //

for (int i = 1; i <= ns; i++)
{
//Ut = A.colPivHouseholderQr().solve(U0);
Ut = A.householderQr().solve(_B*U0); // my program stays with this line forever.
U0 = Ut;
}


I have 64bit Windows 7 with 16 Gb memory, my IDE is VS2010 (there is only 32bit version). The version of Eigen is 3.2

Solving a linear system of equation is a O(n^3) problem, so in your case the complexity is 5000^3. Nonetheless, it should still solve it in a few seconds. Make sure you compiled your program with optimization ON (-O2 for gcc or clang). If your matrix is actually sparse, i.e., that 'A' is composed of numerous zeros (>90%), then you should use a SparseMatrix with sparse solvers.

Hello ggael,

Could you help me to find out why my program runs slowly?
My IDE's setting: VS2010 x64, optimization is on

Here is an example of my code (I just want to fill sparse matrices and solve linear equations)

Code: Select all

typedef SparseMatrix<double, 0, long int> SpaMat;
typedef SparseLU<SpaMat, COLAMDOrdering<long int>> SpaLUSolver;
typedef SparseQR<SpaMat, COLAMDOrdering<long int>> SpaQRSolver;

Eigen::MatrixXd step1(const SpaMat& _B, const SpaMat& _Q, int _np)
{
    MatrixXd U0 = MatrixXd::Zero(_np, 1);
    U0(0, 0) = 1;

   SpaMat A = _B -  _Q;
   SpaLUSolver dec;
   dec.compute(A);
   MatrixXd Ut = dec.solve(_B * U0);

    return Ut;
}

MatrixXd step2(const SpaMat& _B, const SpaMat& _Q, const MatrixXd& _d)
{ 
    SpaQRSolver dec;   
    dec.compute(_Q);
    MatrixXd p = dec.solve(_B * _d);

    return p;
}

Here is the main function:

Code: Select all

int np = B.size();

SpaMat matrixB(np, np);   
matrixB.reserve(B.size());    //// B is a vector of double, so matrixB is diagonal matrix
   
for (int i = 0; i < B.size(); ++i)
{       
   matrixB.insert(i, i) = B[i];
}   
matrixB.makeCompressed();

SpaMat matrixQ(np, np);
matrixQ.reserve(S.size()); // vectors I and J store indexes, vector S stores non-zero values of sparse matrixQ
   
for (int i = 0; i < S.size(); ++i)
{       
   matrixQ.insert(I[i] - 1, J[i] - 1) = S[i];
}   
matrixQ.makeCompressed();

MatrixXd Ut = step1(matrixB, matrixQ, np);
......................
// create matrix d with size (np,1), all values in d is non-zero
MatrixXd di = step2(matrixB, matrixQ, d);

With the size 7500x7500, the running time is around 8 minutes (5 minutes for the filling process, 3 minutes for step2). I used SuiteSparseQR before and it's faster, as I remember.

And in step2 function, I get different results if I use SparseLU or BiCGSTAB solver (it is strange, I think all results should be the same). For example, the result is wrong with this setting:

Code: Select all

BiCGSTAB<SpaMat, IncompleteLUT<double>> dec;   
dec.preconditioner().setDroptol(10e-15);

With the size 72000x72000, my program takes much time for filling and stops running from step1. Could you tell me how can I make it run with large sparse matrices? I need solve linear equations with much larger sizes.

Have a look at the documentation for the recommended way to assemble a sparse matrix: http://eigen.tuxfamily.org/dox/group__T ... tml#title3

ggael wrote:Have a look at the documentation for the recommended way to assemble a sparse matrix: http://eigen.tuxfamily.org/dox/group__T ... tml#title3

Thanks, I will read it again and more carefully. Maybe I need to do something more with the reserving step.
How about the question with the matrices' size? Is it possible to use Eigen to solve a linear equation with large sparse matrices( For example 100000x100000 double values) ?

You should rather use a basic triplet list to be sure to get reasonable performance.

100000^2 is not that large, but speed mostly depends on the structure and sparsity of the matrix. If you're solving a laplacian-like problem on a 2D domain, then you should get nearly interactive-time performance with direct solvers. If you're working on 3D meshes, then iterative solvers should be faster.

btw, I've just seen that matrixB is a diagonal matrix, then use B.asDiagonal() instead of a sparse matrix! This will be considerably faster!! To compute A=B-Q, first do A = -Q, and then subtract the vector B to the diagonal of A: A.diagonal()-=B (assuming that the diagonal of Q is non zero). If the diagonal of Q is empty, it might be a good idea to add explicit zeros along the diagonal when creating the triplet list.

Thank you so much, I set a wrong size for sparse matrices. It runs much faster now.

Hello ggael,

Could you please tell me the fastest way to detect singular sparse matrix in Eigen? Running SVD takes quite much time for me.
More questions, what solver in Eigen is appropriate with matrix which is very close to singular (1-norm condition number is around 1e+19), and how to optimize that solver? For almost singular matrix (symmetric but not positive definite) with size 7500x7500, the elapsed time is from 5s to 11s.

Try the SarseQR solver, it is rank revealing.