[eigen] Does eigen store zeros in sparse matrices ?

[Brad Bell <bradbell@xxxxxxxxxx>]

I have an application where some elements of result may, or may not, be 
zero. In addition, I have to inform the an optimizer of which elements 
may be no-zero. Therefore, I want to do an Eigen sparse matrix 
calculation and know which results may be non-zero.

It appears that Eigen (version 3.2.7) matrix multiply, but not its 
Cholesky solve, works for this application. For example, if A is a 
sparse identity matrix, and x contains a component with the value zero, 
then y =  A * x will also contain the zero element.  On the other hand, 
if I use a Choleksy factorization to solve the equation x = A * y, the 
result y will not contain the zero element. I am including source code 
that demonstrates this below:


# 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 = 5;
    sparse_matrix A(n, n), x(n, 1), y(n, 1);
    for(int i = 0; i < n; i++)
    {   A.insert(i, i) = 1.0;         // n by n identity matrix
        x.insert(i, 0) = double(i);   // transpose of (0, 1, ..., n-1)
    }
    y = A * x;
    std::cout << "Test if Eignes's sparse matrix multiply stores zeros\n";
    for(int j = 0; j < y.outerSize(); j++)
    {   for(column_itr itr(y, j); itr; ++itr)
        {   int row      = itr.row();
            int col      = itr.col();
            double value = itr.value();
            std::cout << "y(" << row << "," << col << ")=" << value << 
"\n";
        }
    }
    sparse_cholesky C;  // begin sparse Cholesky solve x = A * y
    C.analyzePattern(A);
    C.factorize(A);
    y = C.solve(x);     // end sparse Cholesky solve
    std::cout << "Test if Eignes's sparse Cholesky solve stores zeros\n";
    for(int j = 0; j < y.outerSize(); j++)
    {   for(column_itr itr(y, j); itr; ++itr)
        {   int row      = itr.row();
            int col      = itr.col();
            double value = itr.value();
            std::cout << "y(" << row << "," << col << ")=" << value << 
"\n";
        }
    }
    return 0;
}