diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index a07b5a9f3..e4a4c3a7b 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -186,7 +186,6 @@ class SparseLU // Triangular solve Map, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); Map< Matrix, 0, OuterStride<> > U (&(X.data()[fsupc]), nsupc, nrhs, OuterStride<>(X.rows()) ); -// Block > U(X, fsupc, 0, nsupc, nrhs); //FIXME TODO Consider more RHS U = A.template triangularView().solve(U); // Matrix-vector product @@ -536,6 +535,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Eliminate the current column info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); + eigen_assert(info==0 && " SINGULAR MATRIX"); if ( info ) { m_info = NumericalIssue; @@ -609,6 +609,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Form the L-segment info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); + eigen_assert(info==0 && " SINGULAR MATRIX"); if ( info ) { std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); VectorBlock u(tempv, 0, segsize); diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h index 42218ba4a..91c795fac 100644 --- a/Eigen/src/SparseLU/SparseLU_pruneL.h +++ b/Eigen/src/SparseLU/SparseLU_pruneL.h @@ -123,7 +123,7 @@ void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, cons if (perm_r(lsub(kmax)) == IND_EMPTY) kmax--; else if ( perm_r(lsub(kmin)) != IND_EMPTY) - kmin--; + kmin++; else { // kmin below pivrow (not yet pivoted), and kmax diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 0bbbb0627..4727cc12b 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp @@ -52,6 +52,7 @@ int main(int argc, char **args) } /* Compute the factorization */ + solver.isSymmetric(true); solver.compute(A); solver._solve(b, x);