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;
}