mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

return the index of the first non positive diagonal entry (more useful than simply true or false)

Browse Source
This commit is contained in:
Gael Guennebaud 2011-01-17 11:09:03 +01:00
parent 8b6c1caa3e
commit ef3e690a0c
1 changed files with 13 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Eigen/src/Cholesky/LLT.h
View File

@ -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