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

bug #1076: fix scaling in IncompleteCholesky, improve doc, add read-only access to the different factors, remove debugging code.

Browse Source
This commit is contained in:
Gael Guennebaud 2015-10-06 13:25:45 +02:00
parent f25bdc707f
commit 752a0e5339
1 changed files with 25 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
View File

@ -24,6 +24,11 @@ namespace Eigen {
* matrix. It is advised to give a row-oriented sparse matrix * matrix. It is advised to give a row-oriented sparse matrix
* \tparam _UpLo The triangular part of the matrix to reference. * \tparam _UpLo The triangular part of the matrix to reference.
* \tparam _OrderingType * \tparam _OrderingType
*
* It performs the following incomplete factorization: \f$ S P A P' S \approx L L' \f$
* where L is a lower triangular factor, S if a diagonal scaling matrix, and P is a
* fill-in reducing permutation as computed of the ordering method.
*
*/ */
template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> > template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
@ -86,6 +91,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
if(pinv.size()>0) m_perm = pinv.inverse(); if(pinv.size()>0) m_perm = pinv.inverse();
else m_perm.resize(0); else m_perm.resize(0);
m_analysisIsOk = true; m_analysisIsOk = true;
m_isInitialized = true;
} }
template<typename MatrixType> template<typename MatrixType>
@ -110,9 +116,17 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up
x = m_scale.asDiagonal() * x; x = m_scale.asDiagonal() * x;
if (m_perm.rows() == b.rows()) if (m_perm.rows() == b.rows())
x = m_perm.inverse() * x; x = m_perm.inverse() * x;
} }
/** \returns the sparse lower triangular factor L */
const FactorType& matrixL() const { return m_L; }
/** \returns a vector representing the scaling factor S */
const VectorRx& scalingS() const { return m_scale; }
/** \returns the fill-in reducing permutation P (can be empty for a natural ordering) */
const PermutationType permutationP() const { return m_perm; }
protected: protected:
FactorType m_L; // The lower part stored in CSC FactorType m_L; // The lower part stored in CSC
VectorRx m_scale; // The vector for scaling the matrix VectorRx m_scale; // The vector for scaling the matrix
@ -177,12 +191,20 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
m_scale = m_scale.cwiseSqrt().cwiseSqrt(); m_scale = m_scale.cwiseSqrt().cwiseSqrt();
for (Index j = 0; j < n; ++j)
if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
m_scale(j) = RealScalar(1)/m_scale(j);
else
m_scale(j) = 1;
// FIXME disable scaling if not needed, i.e., if it is roughtly uniform? (this will make solve() faster)
// Scale and compute the shift for the matrix // Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest(); RealScalar mindiag = NumTraits<RealScalar>::highest();
for (Index j = 0; j < n; j++) for (Index j = 0; j < n; j++)
{ {
for (Index k = colPtr[j]; k < colPtr[j+1]; k++) for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
vals[k] /= (m_scale(j)*m_scale(rowIdx[k])); vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored"); eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag); mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
} }
@ -240,7 +262,6 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
// Scale the current column // Scale the current column
if(numext::real(diag) <= 0) if(numext::real(diag) <= 0)
{ {
std::cerr << "\nNegative diagonal during Incomplete factorization at position " << j << " (value = " << diag << ")\n";
m_info = NumericalIssue; m_info = NumericalIssue;
return; return;
} }
@ -276,7 +297,6 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
} }
m_factorizationIsOk = true; m_factorizationIsOk = true;
m_isInitialized = true;
m_info = Success; m_info = Success;
} }