Hi,

I've implemented a simple incomplete LU solver for sparse matrices and have a problem solving the resulting L*U*x=b system. If LU is a sparse matrix and x,y,b dense vectors, then the code

Code: Select all

y = LU.triangularView<Eigen::UnitLower>().solve(b);
x = LU.triangularView<Eigen::Upper>().solve(y);


yields the runtime error

Code: Select all

Assertion failed: (it && it.index() == i), function run, file /.../Eigen/src/Sparse/TriangularSolver.h, line 96.


When I use a dense matrix DLU = LU.toDense(); instead, the code runs just fine.

I'm using the newest packed dev version of Eigen with GCC 4.3.4.

I just noticed that

Code: Select all

Eigen::SparseVector<double> w;
w.norm();


yields a compilation error:

Code: Select all

/.../eigen/Eigen/src/Sparse/SparseDot.h: In member function 'typename Eigen::NumTraits<typename Eigen::ei_traits<MatrixType>::Scalar>::Real Eigen::SparseMatrixBase<Derived>::squaredNorm() const [with Derived = Eigen::SparseVector<double, 0>]':
/.../eigen/Eigen/src/Sparse/SparseDot.h:94:   instantiated from 'typename Eigen::NumTraits<typename Eigen::ei_traits<MatrixType>::Scalar>::Real Eigen::SparseMatrixBase<Derived>::norm() const [with Derived = Eigen::SparseVector<double, 0>]'
../directsparsesolver.hpp:54:   instantiated from here
/.../eigen/Eigen/src/Sparse/SparseDot.h:87: error: 'const class Eigen::SparseMatrixBase<Eigen::SparseVector<double, 0> >' has no member named 'cwise'


while

Code: Select all

w.cwiseProduct(w);

compiles fine.

Hi,

w.norm() has been just fixed.

Regarding your issue with the triangular solver, the pb is that for efficiency reasons it assumes the matrix is triangular, i.e., only the respective triangular part is stored. Otherwise, for instance for a column major lower triangular matrix, we would have to manually skip all the entries until find the first one holding to the lower part, and for a upper triangular matrix, check every time we have not crossed the diagonal... Well actually, we could add such tests, if the matrix is already triangular, they should have no cost.

Actually, for sparse matrices, there is almost no memory advantage in packing two triangular matrices into a single one. So, anyway, I would suggest you to split the LU matrix into a L and U one.