Hi everyone,
I have a sparse lower triangular matrix A. I need to compute the product A_transpose * A. I realize that the result is not as sparse as the original lower triangular (but is still relatively sparse), so I would still like to store it in a sparse matrix.

The issue is that this works fine for problem sizes up to around 150,000 (i.e. 150,000 x 150,000 sparse matrix). But once I cross 200,000 rows, the sparse matrix product gives me a bad_alloc exception. How do I avoid this issue?

Code: Select all

SparseMatrix<Real> B = A.transpose() * A;

Both A and B are column major. Some sparsity information:

Success case: For a 161,826 rows. A has 4,691,214 non-zeros (0.0179% fill in). The product B has 9,220,662 non-zeros (0.0352% fill in).

Failure case: For a 215,768 x 215,768 matrix, A has 8303988 nonzeros 0.01178% filled.

Thanks a lot!

right, the problem seems that the current heuristic to preallocate memory of the result assumes a random structure for both operands and then the sparsity of the result is largely underestimated ! (0.3% fill in in your case). I've just committed a more optimistic one.

ggael wrote:right, the problem seems that the current heuristic to preallocate memory of the result assumes a random structure for both operands and then the sparsity of the result is largely underestimated ! (0.3% fill in in your case). I've just committed a more optimistic one.

Thanks! I'm new to the Eigen development process - I'm using the 3.0.4 release. What would I need to do to use your newly committed patch?

you need the "unstable source code from the development branch" as explained there:
http://eigen.tuxfamily.org/index.php?ti ... e#Download

and the respective doc is there:

http://eigen.tuxfamily.org/dox-devel/Tu ... parse.html