Hello. I am new to Eigen2, but am already very impressed, nice work all. I have a question involving Sparse Matrices, and the Cholesky Decomposition SparseLLT. I searched the archives and documentation, and even browsed the source code a bit but am stymied. (I admit I am not too familiar with template programming however). Hopefully in the future I may be able to help out some.

Here is an example piece of code I am struggling with. I am attempting a Cholesky Decomposition of a known SPD matrix, using SparseMatrix and SparseLLT. The code works well up to dimension 2999 X 2999 matrices.

Code: Select all

#include
#include
#include
#include
#include

using namespace std;
using namespace Eigen;

int main(int argc, char** argv)
{
   if (argc  m1(rows, cols);
    
   m1.startFill();
    m1.fill(0, 0) = (float)(2.0);
    m1.fill(0, 1) = (float)(-1.0);
   for (int ir=1; ir, RowMajor > CholFact;
    CholFact.compute(m1);
}


And compiling with gcc version 4.3.2 on a 64 bit Linux machine x86_64, Intel chip.
g++ -I/path/to/eigen2 test.cpp -o test

The code executes and returns meaningful results up to dimension 2999. At 3000 the following error is reported upon execution:

test: /home/tirons/src/eigen2/Eigen/src/Sparse/SparseMatrix.h:165: typename Eigen::ei_traits >::Scalar& Eigen::SparseMatrix::fill(int, int) [with _Scalar = double, int _Flags = 2048]: Assertion `m_data.index(m_data.size()-1)<inner && "wrong sorted insertion"' failed.
Aborted

Which suggests to me that I am filling the matrix wrong? If I don't fill the off-diagonal entries I do not get this error. I have tried with DynamicSparseMatrix, and also fillRand. But I get the same error. The matrix vector multiplication works for these larger sized matrices. Am I doing something wrong here? The error is triggered upon the call to CholFact.compute(m1). Any help is appreciated.

Regards,
Trevor

Hi,
The only reference here is Gael, and he is away for a week for a conference, but since you are having suddenly an error with a special size, 3000, this suggests that you are doing something that depends on a magic value. Looking at SparseMatrix.h, I can see a magic value:

at line 139 (in trunk):

Code: Select all

    inline void startFill(int reserveSize = 1000)


So i would suggest that you try passing explicitly a larger value to startFill...

I have no idea whether this is supposed to be required or if you hit a bug in Eigen -- I know nothing about sparse matrices.

Thanks for the reply. I tried several values for the argument to startFill but am still seeing this error. Hopefully Gael will have some insight when he returns. If I figure it out before then, I'll post. Cheers.

Your suggestion helped me at least track down the problem. On line 178 in SparseLLT.h the line

Code: Select all

m_matrix.fill(it.index(), j) = it.value() * y;


Commenting out this line allows the program to execute, but the result is obviously incorrect.

Alternatively changing the fills to fillrand lets the program execute. However, the result is wrong *above 3000*. Below 3000, it is fine. I'm not sure if the algorithm breaks down at 3000, or if there is something else going, on however.

the problem is that currently the default sparse LLT algorithm assumes the input matrix is column major and that only the lower triangular part is filled.

I just added an assert.

anyway, the default sparse LLT is very basic and so it is pretty bad in term of both performance and filling. I highly recommend you to use one of our backends (Cholmod or Taucs).

kinetic wrote:Your suggestion helped me at least track down the problem. On line 178 in SparseLLT.h the line

Code: Select all

m_matrix.fill(it.index(), j) = it.value() * y;


Commenting out this line allows the program to execute, but the result is obviously incorrect.

Alternatively changing the fills to fillrand lets the program execute. However, the result is wrong *above 3000*. Below 3000, it is fine. I'm not sure if the algorithm breaks down at 3000, or if there is something else going, on however.

Thanks. In this case not linking to an external library is more important than performance. But in the future I'll use the bindings. Nice library. It is working for me now.