GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 03:39:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix SparseLU::absDeterminant and add respective unit test

Browse Source
This commit is contained in:
Gael Guennebaud 2014-10-17 16:52:56 +02:00
parent a13bc22204
commit a370b1f2e2
4 changed files with 40 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/src/SparseLU/SparseLU.h
View File

@ -249,14 +249,13 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
// Initialize with the determinant of the row matrix // Initialize with the determinant of the row matrix
Scalar det = Scalar(1.); 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 :) // which are available in the L part :)
for (Index j = 0; j < this->cols(); ++j) for (Index j = 0; j < this->cols(); ++j)
{ {
for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
{ {
if(it.row() < j) continue; if(it.index() == j)
if(it.row() == j)
{ {
det *= abs(it.value()); det *= abs(it.value());
break; break;

4
Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
View File

@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
m_idval(mat.colIndexPtr()[outer]), m_idval(mat.colIndexPtr()[outer]),
m_startidval(m_idval), m_startidval(m_idval),
m_endidval(mat.colIndexPtr()[outer+1]), m_endidval(mat.colIndexPtr()[outer+1]),
m_idrow(mat.rowIndexPtr()[outer]), m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
m_endidrow(mat.rowIndexPtr()[outer+1]) m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
{} {}
inline InnerIterator& operator++() inline InnerIterator& operator++()
{ {

33
test/sparse_solver.h
View File

@ -133,7 +133,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType&
Scalar refDet = dA.determinant(); Scalar refDet = dA.determinant();
VERIFY_IS_APPROX(refDet,solver.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> template<typename Solver, typename DenseMat>
int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) 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); 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);
}
}

3
test/sparselu.cpp
View File

@ -44,6 +44,9 @@ template<typename T> void test_sparselu_T()
check_sparse_square_solving(sparselu_colamd); check_sparse_square_solving(sparselu_colamd);
check_sparse_square_solving(sparselu_amd); check_sparse_square_solving(sparselu_amd);
check_sparse_square_solving(sparselu_natural); check_sparse_square_solving(sparselu_natural);
check_sparse_square_abs_determinant(sparselu_colamd);
check_sparse_square_abs_determinant(sparselu_amd);
} }
void test_sparselu() void test_sparselu()