Appcrash on running UmfPackLU<> with Eigen library

Question

I am compiling and trying to run the UMfPackLU<SparseMatrix<>> routine from Eigen 3.2.9 and UMFPACK v4.5 libraries with TDM-GCC 5.1.0 on Win64 platform. But I am getting Appcrash with exception code c0000005.

What I need to implement is the following:

     _ _        _ _
A = | P |, B = | R |, where P and Q are sparse and Z is 0 with 3 cols
    | Q |      | Z |
    |_ _|      |_ _|

X = A\B;

What I am doing (excerpt only) is the following:

#define num_t double
...
SparseMatrix<num_t,RowMajor> A(P.rows()+Q.rows(), P.cols());
A.topRows(P.rows()) = P;
A.bottomRows(Q.rows()) = Q;
Matrix<num_t, Dynamic, 3> B(P.rows()+Q.rows(), 3);
B.topLeftCorner(P.rows(), 3) = R;
B.bottomLeftCorner(Q.rows(), 3) = S;

UmfPackLU<SparseMatrix<num_t>> solver(A.transpose()*A);
auto AtB = A.transpose()*B;
X.col(0) = solver.solve(AtB.col(0)); // @@@ segmentation error here @@@
X.col(1) = solver.solve(AtB.col(1));
X.col(2) = solver.solve(AtB.col(2));

Note the SparseMatrix<> is in RowMajor format.

On debugging with gdb: I get Program received signal SIGSEGV, Segmentation fault. at the line marked as above`.

Instead of UmfPackLU<SparseMatrix<>>, solving with SimplicialLLT<SparseMatrix<>>, SimplicialLDLT<SparseMatrix<>> or CholmodDecomposition<SparseMatrix<>> is working correctly.

Thanks in advance for any help.

Answer 1

Please feel free to correct/enhance/establish my answer.

What I have found is that, I need to instantiate an explicit ColMajor matrix AtA before feeding it to the solver (RowMajor is not working) as follows:

SparseMatrix<num_t, ColMajor> AtA = A.transpose()*A;
UmfPackLU<SparseMatrix<num_t>> solver(AtA);

Is this a requirement due to Eigen's lazy evaluation implementation for calling an external routine?

Answer 2

This is a shortcoming in Eigen 3.2.9 that has been fixed a while ago in the 3.3 branch. It's now fixed in the 3.2 branch too (changeset 1e7d97fea51d).

You can workaround the issue by calling compute(...) instead of the constructor:

UmfPackLU<SparseMatrix<num_t>> solver;
solver.compute(A.transpose()*A);