diff options
author | Gael Guennebaud <g.gael@free.fr> | 2016-08-23 23:15:55 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2016-08-23 23:15:55 +0200 |
commit | 8132a126259051f923556fd2c31b9221e2e165e7 (patch) | |
tree | 36a9500212819485cd3f21d59bd182834923317a /Eigen/src/Cholesky | |
parent | bde9b456dcfe78569fbf06ffa9cbf974a810b98e (diff) |
bug #1268: detect faillure in LDLT and report them through info()
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 29 |
1 files changed, 25 insertions, 4 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 69e176110..795d19dce 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -253,7 +253,7 @@ template<typename _MatrixType, int _UpLo> class LDLT ComputationInfo info() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return Success; + return m_info; } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -281,6 +281,7 @@ template<typename _MatrixType, int _UpLo> class LDLT TmpMatrixType m_temporary; internal::SignMatrix m_sign; bool m_isInitialized; + ComputationInfo m_info; }; namespace internal { @@ -298,6 +299,8 @@ template<> struct ldlt_inplace<Lower> typedef typename TranspositionType::StorageIndex IndexType; eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); + bool found_zero_pivot = false; + bool ret = true; if (size <= 1) { @@ -356,9 +359,27 @@ template<> struct ldlt_inplace<Lower> // we should only make sure that we do not introduce INF or NaN values. // Remark that LAPACK also uses 0 as the cutoff value. RealScalar realAkk = numext::real(mat.coeffRef(k,k)); - if((rs>0) && (abs(realAkk) > RealScalar(0))) + bool pivot_is_valid = (abs(realAkk) > RealScalar(0)); + + if(k==0 && !pivot_is_valid) + { + // The entire diagonal is zero, there is nothing more to do + // except filling the transpositions, and checking whether the matrix is zero. + sign = ZeroSign; + for(Index j = 0; j<size; ++j) + { + transpositions.coeffRef(j) = IndexType(j); + ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).all(); + } + return ret; + } + + if((rs>0) && pivot_is_valid) A21 /= realAkk; + if(found_zero_pivot && pivot_is_valid) ret = false; // factorization failed + else if(!pivot_is_valid) found_zero_pivot = true; + if (sign == PositiveSemiDef) { if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite; } else if (sign == NegativeSemiDef) { @@ -369,7 +390,7 @@ template<> struct ldlt_inplace<Lower> } } - return true; + return ret; } // Reference for the algorithm: Davis and Hager, "Multiple Rank @@ -493,7 +514,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputTyp m_temporary.resize(size); m_sign = internal::ZeroSign; - internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); + m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success : NumericalIssue; m_isInitialized = true; return *this; |