34 template<
typename MatrixType,
int UpLo>
struct LDLT_Traits;
60 template<
typename _MatrixType,
int _UpLo>
class LDLT
72 typedef typename MatrixType::Scalar
Scalar;
74 typedef typename MatrixType::Index
Index;
80 typedef internal::LDLT_Traits<MatrixType,UpLo>
Traits;
125 inline typename Traits::MatrixU
matrixU()
const
132 inline typename Traits::MatrixL
matrixL()
const
160 #ifdef EIGEN2_SUPPORT
161 inline bool isPositiveDefinite()
const
189 template<
typename Rhs>
190 inline const internal::solve_retval<LDLT, Rhs>
195 &&
"LDLT::solve(): invalid number of rows of the right hand side matrix b");
196 return internal::solve_retval<LDLT, Rhs>(*
this, b.derived());
199 #ifdef EIGEN2_SUPPORT
200 template<
typename OtherDerived,
typename ResultType>
203 *result = this->
solve(b);
208 template<
typename Derived>
213 template <
typename Derived>
259 template<
int UpLo>
struct ldlt_inplace;
261 template<>
struct ldlt_inplace<Lower>
263 template<
typename MatrixType,
typename TranspositionType,
typename Workspace>
264 static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp,
int* sign=0)
266 typedef typename MatrixType::Scalar Scalar;
267 typedef typename MatrixType::RealScalar RealScalar;
268 typedef typename MatrixType::Index Index;
270 const Index size = mat.rows();
274 transpositions.setIdentity();
276 *sign =
real(mat.coeff(0,0))>0 ? 1:-1;
280 RealScalar cutoff(0), biggest_in_corner;
282 for (Index k = 0; k < size; ++k)
285 Index index_of_biggest_in_corner;
286 biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
287 index_of_biggest_in_corner += k;
297 *sign =
real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
301 if(biggest_in_corner < cutoff)
303 for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
307 transpositions.coeffRef(k) = index_of_biggest_in_corner;
308 if(k != index_of_biggest_in_corner)
312 Index s = size-index_of_biggest_in_corner-1;
313 mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
314 mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
315 std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
316 for(
int i=k+1;i<index_of_biggest_in_corner;++i)
318 Scalar tmp = mat.coeffRef(i,k);
319 mat.coeffRef(i,k) =
conj(mat.coeffRef(index_of_biggest_in_corner,i));
320 mat.coeffRef(index_of_biggest_in_corner,i) =
conj(tmp);
323 mat.coeffRef(index_of_biggest_in_corner,k) =
conj(mat.coeff(index_of_biggest_in_corner,k));
330 Index rs = size - k - 1;
337 temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
338 mat.
coeffRef(k,k) -= (A10 * temp.head(k)).value();
340 A21.noalias() -= A20 * temp.head(k);
342 if((rs>0) && (
abs(mat.coeffRef(k,k)) > cutoff))
356 template<
typename MatrixType,
typename WDerived>
357 static bool updateInPlace(MatrixType& mat,
MatrixBase<WDerived>& w,
typename MatrixType::RealScalar sigma=1)
360 typedef typename MatrixType::Scalar Scalar;
361 typedef typename MatrixType::RealScalar RealScalar;
362 typedef typename MatrixType::Index Index;
364 const Index size = mat.rows();
367 RealScalar alpha = 1;
370 for (Index j = 0; j < size; j++)
377 RealScalar dj =
real(mat.coeff(j,j));
378 Scalar wj = w.coeff(j);
379 RealScalar swj2 = sigma*
abs2(wj);
380 RealScalar gamma = dj*alpha + swj2;
382 mat.coeffRef(j,j) += swj2/alpha;
388 w.
tail(rs) -= wj * mat.col(j).tail(rs);
390 mat.col(j).tail(rs) += (sigma*
conj(wj)/gamma)*w.
tail(rs);
395 template<
typename MatrixType,
typename TranspositionType,
typename Workspace,
typename WType>
396 static bool update(MatrixType& mat,
const TranspositionType& transpositions, Workspace& tmp,
const WType& w,
typename MatrixType::RealScalar sigma=1)
399 tmp = transpositions * w;
401 return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
405 template<>
struct ldlt_inplace<
Upper>
407 template<
typename MatrixType,
typename TranspositionType,
typename Workspace>
408 static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp,
int* sign=0)
410 Transpose<MatrixType> matt(mat);
411 return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
414 template<
typename MatrixType,
typename TranspositionType,
typename Workspace,
typename WType>
415 static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w,
typename MatrixType::RealScalar sigma=1)
417 Transpose<MatrixType> matt(mat);
418 return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
422 template<
typename MatrixType>
struct LDLT_Traits<MatrixType,Lower>
424 typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
425 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
426 static inline MatrixL getL(
const MatrixType& m) {
return m; }
427 static inline MatrixU getU(
const MatrixType& m) {
return m.adjoint(); }
430 template<
typename MatrixType>
struct LDLT_Traits<MatrixType,
Upper>
432 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
433 typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
434 static inline MatrixL getL(
const MatrixType& m) {
return m.adjoint(); }
435 static inline MatrixU getU(
const MatrixType& m) {
return m; }
442 template<
typename MatrixType,
int _UpLo>
446 const Index size = a.rows();
450 m_transpositions.resize(size);
451 m_isInitialized =
false;
452 m_temporary.resize(size);
454 internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
456 m_isInitialized =
true;
465 template<
typename MatrixType,
int _UpLo>
466 template<
typename Derived>
469 const Index size = w.rows();
476 m_matrix.resize(size,size);
478 m_transpositions.resize(size);
479 for (
Index i = 0; i < size; i++)
480 m_transpositions.coeffRef(i) = i;
481 m_temporary.resize(size);
483 m_isInitialized =
true;
486 internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
492 template<
typename _MatrixType,
int _UpLo,
typename Rhs>
493 struct solve_retval<
LDLT<_MatrixType,_UpLo>, Rhs>
494 : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
499 template<typename Dest>
void evalTo(Dest& dst)
const
501 eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
503 dst = dec().transpositionsP() * rhs();
506 dec().matrixL().solveInPlace(dst);
512 typedef typename LDLTType::MatrixType MatrixType;
513 typedef typename LDLTType::Scalar Scalar;
514 typedef typename LDLTType::RealScalar RealScalar;
518 for (Index i = 0; i < vectorD.size(); ++i) {
519 if(
abs(vectorD(i)) > tolerance)
520 dst.row(i) /= vectorD(i);
522 dst.row(i).setZero();
526 dec().matrixU().solveInPlace(dst);
529 dst = dec().transpositionsP().transpose() * dst;
547 template<
typename MatrixType,
int _UpLo>
548 template<
typename Derived>
551 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
552 const Index size = m_matrix.rows();
555 bAndX = this->solve(bAndX);
563 template<
typename MatrixType,
int _UpLo>
566 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
567 const Index size = m_matrix.rows();
568 MatrixType res(size,size);
572 res = transpositionsP() * res;
574 res = matrixU() * res;
576 res = vectorD().asDiagonal() * res;
578 res = matrixL() * res;
580 res = transpositionsP().transpose() * res;
588 template<
typename MatrixType,
unsigned int UpLo>
598 template<
typename Derived>
607 #endif // EIGEN_LDLT_H