diff options
author | 2011-01-17 11:09:03 +0100 | |
---|---|---|
committer | 2011-01-17 11:09:03 +0100 | |
commit | ef3e690a0c410e7309cb7a0ff60154943375db03 (patch) | |
tree | 3d42e1b60a1d1814ce734463375d2a5894efdc55 /Eigen/src/Cholesky/LLT.h | |
parent | 8b6c1caa3e41d5659e35b1852b791a0f9bd880d2 (diff) |
return the index of the first non positive diagonal entry (more useful than simply true or false)
Diffstat (limited to 'Eigen/src/Cholesky/LLT.h')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 24 |
1 files changed, 13 insertions, 11 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ceb532055..53a0543fe 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -184,11 +184,12 @@ template<int UpLo> struct llt_inplace; template<> struct llt_inplace<Lower> { template<typename MatrixType> - static bool unblocked(MatrixType& mat) + static typename MatrixType::Index unblocked(MatrixType& mat) { + typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; + eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); for(Index k = 0; k < size; ++k) @@ -202,16 +203,16 @@ template<> struct llt_inplace<Lower> RealScalar x = real(mat.coeff(k,k)); if (k>0) x -= A10.squaredNorm(); if (x<=RealScalar(0)) - return false; + return k; mat.coeffRef(k,k) = x = sqrt(x); if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); if (rs>0) A21 *= RealScalar(1)/x; } - return true; + return -1; } template<typename MatrixType> - static bool blocked(MatrixType& m) + static typename MatrixType::Index blocked(MatrixType& m) { typedef typename MatrixType::Index Index; eigen_assert(m.rows()==m.cols()); @@ -235,24 +236,25 @@ template<> struct llt_inplace<Lower> Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs); Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs); - if(!unblocked(A11)) return false; + Index ret; + if((ret=unblocked(A11))>=0) return k+ret; if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck } - return true; + return -1; } }; template<> struct llt_inplace<Upper> { template<typename MatrixType> - static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat) + static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); return llt_inplace<Lower>::unblocked(matt); } template<typename MatrixType> - static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat) + static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat) { Transpose<MatrixType> matt(mat); return llt_inplace<Lower>::blocked(matt); @@ -266,7 +268,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Lower>::blocked(m); } + { return llt_inplace<Lower>::blocked(m)==-1; } }; template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> @@ -276,7 +278,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) - { return llt_inplace<Upper>::blocked(m); } + { return llt_inplace<Upper>::blocked(m)==-1; } }; } // end namespace internal |