diff options
Diffstat (limited to 'Eigen/src')
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 55 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/LLT.h | 2 |
2 files changed, 47 insertions, 10 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 7c8e1eb04..8cfc256bb 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -62,14 +62,21 @@ template<typename _MatrixType> class LDLT typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<int, 1, MatrixType::RowsAtCompileTime> IntRowVectorType; - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via LDLT::compute(const MatrixType&). - */ + /** \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via LDLT::compute(const MatrixType&). + */ LDLT() : m_matrix(), m_p(), m_transpositions(), m_isInitialized(false) {} + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa LDLT() + */ + LDLT(int size) : m_matrix(size,size), m_p(size), m_transpositions(size), m_isInitialized(false) {} + LDLT(const MatrixType& matrix) : m_matrix(matrix.rows(), matrix.cols()), m_p(matrix.rows()), @@ -148,6 +155,8 @@ template<typename _MatrixType> class LDLT return m_matrix; } + const MatrixType reconstructedMatrix() const; + inline int rows() const { return m_matrix.rows(); } inline int cols() const { return m_matrix.cols(); } @@ -175,6 +184,10 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) m_matrix = a; + m_p.resize(size); + m_transpositions.resize(size); + m_isInitialized = false; + if (size <= 1) { m_p.setZero(); m_transpositions.setZero(); @@ -228,8 +241,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) continue; } - RealScalar Djj = ei_real(m_matrix.coeff(j,j) - m_matrix.row(j).head(j) - .dot(m_matrix.col(j).head(j))); + RealScalar Djj = ei_real(m_matrix.coeff(j,j) - m_matrix.row(j).head(j).dot(m_matrix.col(j).head(j))); m_matrix.coeffRef(j,j) = Djj; int endSize = size - j - 1; @@ -238,7 +250,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) * m_matrix.col(j).head(j).conjugate(); m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate() - - _temporary.tail(endSize).transpose(); + - _temporary.tail(endSize).transpose(); if(ei_abs(Djj) > cutoff) { @@ -308,6 +320,31 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const return true; } +/** \returns the matrix represented by the decomposition, + * i.e., it returns the product: P^T L D L^* P. + * This function is provided for debug purpose. */ +template<typename MatrixType> +const MatrixType LDLT<MatrixType>::reconstructedMatrix() const +{ + ei_assert(m_isInitialized && "LDLT is not initialized."); + const int size = m_matrix.rows(); + MatrixType res(size,size); + res.setIdentity(); + + // PI + for(int i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i)); + // L^* P + res = matrixL().adjoint() * res; + // D(L^*P) + res = vectorD().asDiagonal() * res; + // L(DL^*P) + res = matrixL() * res; + // P^T (LDL^*P) + for (int i = size-1; i >= 0; --i) res.row(m_transpositions.coeff(i)).swap(res.row(i)); + + return res; +} + /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this */ diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 474b82406..96e1e5f73 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -117,7 +117,7 @@ template<typename _MatrixType, int _UpLo> class LLT && "LLT::solve(): invalid number of rows of the right hand side matrix b"); return ei_solve_retval<LLT, Rhs>(*this, b.derived()); } - + template<typename Derived> bool solveInPlace(MatrixBase<Derived> &bAndX) const; |