16 #include "./InternalHeaderCheck.h"
21 template<
typename MatrixType_,
int UpLo_>
struct traits<LDLT<MatrixType_, UpLo_> >
24 typedef MatrixXpr XprKind;
25 typedef SolverStorage StorageKind;
26 typedef int StorageIndex;
30 template<
typename MatrixType,
int UpLo>
struct LDLT_Traits;
33 enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
61 template<
typename MatrixType_,
int UpLo_>
class LDLT
65 typedef MatrixType_ MatrixType;
69 EIGEN_GENERIC_PUBLIC_INTERFACE(
LDLT)
71 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
72 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
80 typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
90 m_sign(internal::ZeroSign),
91 m_isInitialized(false)
102 m_transpositions(
size),
104 m_sign(internal::ZeroSign),
105 m_isInitialized(false)
114 template<
typename InputType>
116 : m_matrix(matrix.rows(), matrix.cols()),
117 m_transpositions(matrix.rows()),
118 m_temporary(matrix.rows()),
119 m_sign(internal::ZeroSign),
120 m_isInitialized(false)
131 template<
typename InputType>
134 m_transpositions(matrix.rows()),
135 m_temporary(matrix.rows()),
136 m_sign(internal::ZeroSign),
137 m_isInitialized(false)
147 m_isInitialized =
false;
151 inline typename Traits::MatrixU
matrixU()
const
153 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
154 return Traits::getU(m_matrix);
158 inline typename Traits::MatrixL
matrixL()
const
160 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
161 return Traits::getL(m_matrix);
168 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
169 return m_transpositions;
175 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
176 return m_matrix.diagonal();
182 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
183 return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
189 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
190 return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
193 #ifdef EIGEN_PARSED_BY_DOXYGEN
209 template<
typename Rhs>
214 template<
typename Derived>
217 template<
typename InputType>
225 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
226 return internal::rcond_estimate_helper(m_l1_norm, *
this);
229 template <
typename Derived>
238 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
251 EIGEN_DEVICE_FUNC
inline EIGEN_CONSTEXPR
Index rows() const EIGEN_NOEXCEPT {
return m_matrix.rows(); }
252 EIGEN_DEVICE_FUNC
inline EIGEN_CONSTEXPR
Index cols() const EIGEN_NOEXCEPT {
return m_matrix.cols(); }
261 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
265 #ifndef EIGEN_PARSED_BY_DOXYGEN
266 template<
typename RhsType,
typename DstType>
267 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
269 template<
bool Conjugate,
typename RhsType,
typename DstType>
270 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
275 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
284 RealScalar m_l1_norm;
285 TranspositionType m_transpositions;
286 TmpMatrixType m_temporary;
287 internal::SignMatrix m_sign;
288 bool m_isInitialized;
294 template<
int UpLo>
struct ldlt_inplace;
296 template<>
struct ldlt_inplace<
Lower>
298 template<
typename MatrixType,
typename TranspositionType,
typename Workspace>
299 static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix&
sign)
302 typedef typename MatrixType::Scalar Scalar;
303 typedef typename MatrixType::RealScalar RealScalar;
304 typedef typename TranspositionType::StorageIndex IndexType;
305 eigen_assert(mat.rows()==mat.cols());
306 const Index size = mat.rows();
307 bool found_zero_pivot =
false;
312 transpositions.setIdentity();
313 if(size==0)
sign = ZeroSign;
314 else if (numext::real(mat.coeff(0,0)) >
static_cast<RealScalar
>(0) )
sign = PositiveSemiDef;
315 else if (numext::real(mat.coeff(0,0)) <
static_cast<RealScalar
>(0))
sign = NegativeSemiDef;
316 else sign = ZeroSign;
320 for (
Index k = 0; k < size; ++k)
323 Index index_of_biggest_in_corner;
324 mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
325 index_of_biggest_in_corner += k;
327 transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner);
328 if(k != index_of_biggest_in_corner)
332 Index s = size-index_of_biggest_in_corner-1;
333 mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
334 mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
335 std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
336 for(
Index i=k+1;i<index_of_biggest_in_corner;++i)
338 Scalar tmp = mat.coeffRef(i,k);
339 mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i));
340 mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp);
342 if(NumTraits<Scalar>::IsComplex)
343 mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k));
350 Index rs = size - k - 1;
357 temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
358 mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
360 A21.noalias() -= A20 * temp.head(k);
367 RealScalar realAkk = numext::real(mat.coeffRef(k,k));
368 bool pivot_is_valid = (
abs(realAkk) > RealScalar(0));
370 if(k==0 && !pivot_is_valid)
375 for(
Index j = 0; j<size; ++j)
377 transpositions.coeffRef(j) = IndexType(j);
378 ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).
all();
383 if((rs>0) && pivot_is_valid)
386 ret = ret && (A21.array()==Scalar(0)).
all();
388 if(found_zero_pivot && pivot_is_valid) ret =
false;
389 else if(!pivot_is_valid) found_zero_pivot =
true;
391 if (
sign == PositiveSemiDef) {
392 if (realAkk <
static_cast<RealScalar
>(0))
sign = Indefinite;
393 }
else if (
sign == NegativeSemiDef) {
394 if (realAkk >
static_cast<RealScalar
>(0))
sign = Indefinite;
395 }
else if (
sign == ZeroSign) {
396 if (realAkk >
static_cast<RealScalar
>(0))
sign = PositiveSemiDef;
397 else if (realAkk <
static_cast<RealScalar
>(0))
sign = NegativeSemiDef;
411 template<
typename MatrixType,
typename WDerived>
412 static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w,
const typename MatrixType::RealScalar& sigma=1)
414 using numext::isfinite;
415 typedef typename MatrixType::Scalar Scalar;
416 typedef typename MatrixType::RealScalar RealScalar;
418 const Index size = mat.rows();
419 eigen_assert(mat.cols() == size && w.size()==size);
421 RealScalar alpha = 1;
424 for (
Index j = 0; j < size; j++)
431 RealScalar dj = numext::real(mat.coeff(j,j));
432 Scalar wj = w.coeff(j);
433 RealScalar swj2 = sigma*numext::abs2(wj);
434 RealScalar gamma = dj*alpha + swj2;
436 mat.coeffRef(j,j) += swj2/alpha;
442 w.tail(rs) -= wj * mat.col(j).tail(rs);
443 if(!numext::is_exactly_zero(gamma))
444 mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
449 template<
typename MatrixType,
typename TranspositionType,
typename Workspace,
typename WType>
450 static bool update(MatrixType& mat,
const TranspositionType& transpositions, Workspace& tmp,
const WType& w,
const typename MatrixType::RealScalar& sigma=1)
453 tmp = transpositions * w;
455 return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
459 template<>
struct ldlt_inplace<
Upper>
461 template<
typename MatrixType,
typename TranspositionType,
typename Workspace>
462 static EIGEN_STRONG_INLINE
bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix&
sign)
464 Transpose<MatrixType> matt(mat);
465 return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp,
sign);
468 template<
typename MatrixType,
typename TranspositionType,
typename Workspace,
typename WType>
469 static EIGEN_STRONG_INLINE
bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w,
const typename MatrixType::RealScalar& sigma=1)
471 Transpose<MatrixType> matt(mat);
472 return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
476 template<
typename MatrixType>
struct LDLT_Traits<MatrixType,
Lower>
478 typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
479 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
480 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m); }
481 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m.adjoint()); }
484 template<
typename MatrixType>
struct LDLT_Traits<MatrixType,
Upper>
486 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
487 typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
488 static inline MatrixL getL(
const MatrixType& m) {
return MatrixL(m.adjoint()); }
489 static inline MatrixU getU(
const MatrixType& m) {
return MatrixU(m); }
496 template<
typename MatrixType,
int UpLo_>
497 template<
typename InputType>
506 m_l1_norm = RealScalar(0);
508 for (
Index col = 0; col < size; ++col) {
509 RealScalar abs_col_sum;
511 abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
513 abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
514 if (abs_col_sum > m_l1_norm)
515 m_l1_norm = abs_col_sum;
518 m_transpositions.resize(size);
519 m_isInitialized =
false;
520 m_temporary.resize(size);
521 m_sign = internal::ZeroSign;
523 m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ?
Success :
NumericalIssue;
525 m_isInitialized =
true;
534 template<
typename MatrixType,
int UpLo_>
535 template<
typename Derived>
538 typedef typename TranspositionType::StorageIndex IndexType;
542 eigen_assert(m_matrix.rows()==size);
546 m_matrix.resize(size,size);
548 m_transpositions.resize(size);
549 for (
Index i = 0; i < size; i++)
550 m_transpositions.coeffRef(i) = IndexType(i);
551 m_temporary.resize(size);
552 m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
553 m_isInitialized =
true;
556 internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
561 #ifndef EIGEN_PARSED_BY_DOXYGEN
562 template<
typename MatrixType_,
int UpLo_>
563 template<
typename RhsType,
typename DstType>
566 _solve_impl_transposed<true>(rhs, dst);
569 template<
typename MatrixType_,
int UpLo_>
570 template<
bool Conjugate,
typename RhsType,
typename DstType>
571 void LDLT<MatrixType_,UpLo_>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
574 dst = m_transpositions * rhs;
578 matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
584 const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD());
592 RealScalar tolerance = (std::numeric_limits<RealScalar>::min)();
593 for (
Index i = 0; i < vecD.size(); ++i)
595 if(
abs(vecD(i)) > tolerance)
596 dst.row(i) /= vecD(i);
598 dst.row(i).setZero();
603 matrixL().transpose().template conjugateIf<Conjugate>().solveInPlace(dst);
607 dst = m_transpositions.transpose() * dst;
624 template<
typename MatrixType,
int UpLo_>
625 template<
typename Derived>
626 bool LDLT<MatrixType,UpLo_>::solveInPlace(MatrixBase<Derived> &bAndX)
const
628 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
629 eigen_assert(m_matrix.rows() == bAndX.rows());
631 bAndX = this->solve(bAndX);
639 template<
typename MatrixType,
int UpLo_>
642 eigen_assert(m_isInitialized &&
"LDLT is not initialized.");
643 const Index size = m_matrix.rows();
644 MatrixType res(size,size);
648 res = transpositionsP() * res;
650 res = matrixU() * res;
652 res = vectorD().real().asDiagonal() * res;
654 res = matrixL() * res;
656 res = transpositionsP().transpose() * res;
665 template<
typename MatrixType,
unsigned int UpLo>
676 template<
typename Derived>
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:67
Robust Cholesky decomposition of a matrix with pivoting.
Definition: LDLT.h:63
const TranspositionType & transpositionsP() const
Definition: LDLT.h:166
Traits::MatrixL matrixL() const
Definition: LDLT.h:158
Diagonal< const MatrixType > vectorD() const
Definition: LDLT.h:173
const Solve< LDLT, Rhs > solve(const MatrixBase< Rhs > &b) const
const MatrixType & matrixLDLT() const
Definition: LDLT.h:236
void setZero()
Definition: LDLT.h:145
RealScalar rcond() const
Definition: LDLT.h:223
bool isNegative(void) const
Definition: LDLT.h:187
LDLT()
Default Constructor.
Definition: LDLT.h:87
bool isPositive() const
Definition: LDLT.h:180
LDLT(const EigenBase< InputType > &matrix)
Constructor with decomposition.
Definition: LDLT.h:115
LDLT(Index size)
Default Constructor with memory preallocation.
Definition: LDLT.h:100
LDLT(EigenBase< InputType > &matrix)
Constructs a LDLT factorization from a given matrix.
Definition: LDLT.h:132
MatrixType reconstructedMatrix() const
Definition: LDLT.h:640
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LDLT.h:259
Traits::MatrixU matrixU() const
Definition: LDLT.h:151
const LDLT & adjoint() const
Definition: LDLT.h:249
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const LDLT< PlainObject > ldlt() const
Definition: LDLT.h:678
const LDLT< PlainObject, UpLo > ldlt() const
Definition: LDLT.h:667
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
LDLT< MatrixType_, UpLo_ > & derived()
Definition: EigenBase.h:48
static const Eigen::internal::all_t all
Definition: IndexedViewHelper.h:189
ComputationInfo
Definition: Constants.h:442
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sign_op< typename Derived::Scalar >, const Derived > sign(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_isfinite_op< typename Derived::Scalar >, const Derived > isfinite(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
EIGEN_CONSTEXPR Index size() const EIGEN_NOEXCEPT
Definition: EigenBase.h:69