Re: [eigen] Sparse Cholesky

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

You are using SimplicialLDLT which performs P*A*P^T = L * D * L^T with D a diagonal matrix, and L a lower triangular matrix with ones on the diagonal. Use SimplicialLLT if you want the P*A*P^T = L  * L^T factorization with square roots.

gael.

On Sun, Oct 1, 2017 at 3:07 PM, Brad Bell <bradbell@xxxxxxxxxx> wrote:

I do am trying to get a copy of the lower triangular matrix corresponding to a sparse Cholesky factorization. It seems to me that, in the case of a 1 by 1 matrix, this should be a 1 by 1 matrix with entry equal to the square root of the original matrix. But the following program is printing 1 instead of 2 for its result. What am I doing wrong ?

# include <iostream>
# include <Eigen/Sparse>
int main(void)
{   typedef Eigen::SparseMatrix<double, Eigen::ColMajor> sparse_matrix;
    typedef Eigen::SimplicialLDLT<sparse_matrix, Eigen::Lower> sparse_cholesky;
    typedef sparse_matrix::InnerIterator                       column_itr;
    //
    int n = 1;
    sparse_matrix H(n, n);
    H.insert(0, 0) = 4.0;
    //
    sparse_cholesky C;
    C.compute( H );
    sparse_matrix L = C.matrixL();
    //
    for(int j = 0; j < L.outerSize(); j++)
    {   for(column_itr itr(L, j); itr; ++itr)
        {   std::cout << "L(" << itr.row() << "," << itr.col() << ") = ";
            std::cout << itr.value() << "\n";
        }
    }
    return 0;
}