diff options
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 5 | ||||
| -rw-r--r-- | Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h | 4 | ||||
| -rw-r--r-- | test/sparse_solver.h | 33 | ||||
| -rw-r--r-- | test/sparselu.cpp | 3 |
4 files changed, 40 insertions, 5 deletions
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index b02796293..d72d7f150 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -249,14 +249,13 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); // Initialize with the determinant of the row matrix Scalar det = Scalar(1.); - //Note that the diagonal blocks of U are stored in supernodes, + // Note that the diagonal blocks of U are stored in supernodes, // which are available in the L part :) for (Index j = 0; j < this->cols(); ++j) { for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) { - if(it.row() < j) continue; - if(it.row() == j) + if(it.index() == j) { det *= abs(it.value()); break; diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index 6b0b83e0b..e8ee35a94 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator m_idval(mat.colIndexPtr()[outer]), m_startidval(m_idval), m_endidval(mat.colIndexPtr()[outer+1]), - m_idrow(mat.rowIndexPtr()[outer]), - m_endidrow(mat.rowIndexPtr()[outer+1]) + m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]), + m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1]) {} inline InnerIterator& operator++() { diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 8c8d7f939..b10eaebcc 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -133,7 +133,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& Scalar refDet = dA.determinant(); VERIFY_IS_APPROX(refDet,solver.determinant()); } +template<typename Solver, typename DenseMat> +void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) +{ + using std::abs; + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n"; + return; + } + Scalar refDet = abs(dA.determinant()); + VERIFY_IS_APPROX(refDet,solver.absDeterminant()); +} template<typename Solver, typename DenseMat> int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) @@ -333,3 +349,20 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver) check_sparse_determinant(solver, A, dA); } } + +template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver) +{ + typedef typename Solver::MatrixType Mat; + typedef typename Mat::Scalar Scalar; + typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; + + // generate the problem + Mat A; + DenseMatrix dA; + generate_sparse_square_problem(solver, A, dA, 30); + A.makeCompressed(); + for (int i = 0; i < g_repeat; i++) { + check_sparse_abs_determinant(solver, A, dA); + } +} + diff --git a/test/sparselu.cpp b/test/sparselu.cpp index 37980defc..52371cb12 100644 --- a/test/sparselu.cpp +++ b/test/sparselu.cpp @@ -44,6 +44,9 @@ template<typename T> void test_sparselu_T() check_sparse_square_solving(sparselu_colamd); check_sparse_square_solving(sparselu_amd); check_sparse_square_solving(sparselu_natural); + + check_sparse_square_abs_determinant(sparselu_colamd); + check_sparse_square_abs_determinant(sparselu_amd); } void test_sparselu() |