mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-06 19:10:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix SparseLU::signDeterminant() method, and add a SparseLU::determinant() method.

Browse Source
This commit is contained in:
Gael Guennebaud 2015-02-16 19:09:22 +01:00
parent 8768ff3c31
commit f0b1b1df9b
2 changed files with 54 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -302,7 +302,50 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
Scalar signDeterminant()
{
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
return Scalar(m_detPermR);
// Initialize with the determinant of the row matrix
Index det = 1;
// 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.index() == j)
{
if(it.value()<0)
det = -det;
else if(it.value()==0)
return 0;
break;
}
}
}
return det * m_detPermR * m_detPermC;
}
/** \returns The determinant of the matrix.
*
* \sa absDeterminant(), logAbsDeterminant()
*/
Scalar determinant()
{
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,
// 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.index() == j)
{
det *= it.value();
break;
}
}
}
return det * (m_detPermR * m_detPermC);
}
protected:
@ -337,7 +380,7 @@ class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
internal::perfvalues m_perfv;
RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
Index m_nnzL, m_nnzU; // Nonzeros in L and U factors
Index m_detPermR; // Determinant of the coefficient matrix
Index m_detPermR, m_detPermC; // Determinants of the permutation matrices
private:
// Disable copy constructor
SparseLU (const SparseLU& );
@ -617,7 +660,8 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
}
// Update the determinant of the row permutation matrix
if (pivrow != jj) m_detPermR *= -1;
// FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot.
if (pivrow != jj) m_detPermR = -m_detPermR;
// Prune columns (0:jj-1) using column jj
Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
@ -632,10 +676,13 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
jcol += panel_size; // Move to the next panel
} // end for -- end elimination
m_detPermR = m_perm_r.determinant();
m_detPermC = m_perm_c.determinant();
// Count the number of nonzeros in factors
Base::countnz(n, m_nnzL, m_nnzU, m_glu);
// Apply permutation to the L subscripts
Base::fixupL(n, m_perm_r.indices(), m_glu);
Base::fixupL(n, m_perm_r.indices(), m_glu);
// Create supernode matrix L
m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);

3
test/sparselu.cpp
View File

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