LDLT.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
00006 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00007 // Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com >
00008 //
00009 // This Source Code Form is subject to the terms of the Mozilla
00010 // Public License v. 2.0. If a copy of the MPL was not distributed
00011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00012 
00013 #ifndef EIGEN_LDLT_H
00014 #define EIGEN_LDLT_H
00015 
00016 namespace Eigen { 
00017 
00018 namespace internal {
00019   template<typename MatrixType, int UpLo> struct LDLT_Traits;
00020 
00021   // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
00022   enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
00023 }
00024 
00025 /** \ingroup Cholesky_Module
00026   *
00027   * \class LDLT
00028   *
00029   * \brief Robust Cholesky decomposition of a matrix with pivoting
00030   *
00031   * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
00032   * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
00033   *             The other triangular part won't be read.
00034   *
00035   * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite
00036   * matrix \f$ A \f$ such that \f$ A =  P^TLDL^*P \f$, where P is a permutation matrix, L
00037   * is lower triangular with a unit diagonal and D is a diagonal matrix.
00038   *
00039   * The decomposition uses pivoting to ensure stability, so that L will have
00040   * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
00041   * on D also stabilizes the computation.
00042   *
00043   * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
00044   * decomposition to determine whether a system of equations has a solution.
00045   *
00046   * \sa MatrixBase::ldlt(), class LLT
00047   */
00048 template<typename _MatrixType, int _UpLo> class LDLT
00049 {
00050   public:
00051     typedef _MatrixType MatrixType;
00052     enum {
00053       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00054       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00055       Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
00056       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00057       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00058       UpLo = _UpLo
00059     };
00060     typedef typename MatrixType::Scalar Scalar;
00061     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00062     typedef typename MatrixType::Index Index;
00063     typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1>  TmpMatrixType ;
00064 
00065     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime>  TranspositionType ;
00066     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime>  PermutationType ;
00067 
00068     typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
00069 
00070     /** \brief Default Constructor.
00071       *
00072       * The default constructor is useful in cases in which the user intends to
00073       * perform decompositions via LDLT::compute(const MatrixType&).
00074       */
00075     LDLT() 
00076       : m_matrix(), 
00077         m_transpositions(), 
00078         m_sign(internal::ZeroSign),
00079         m_isInitialized(false) 
00080     {}
00081 
00082     /** \brief Default Constructor with memory preallocation
00083       *
00084       * Like the default constructor but with preallocation of the internal data
00085       * according to the specified problem \a size.
00086       * \sa LDLT()
00087       */
00088     LDLT(Index size)
00089       : m_matrix(size, size),
00090         m_transpositions(size),
00091         m_temporary(size),
00092         m_sign(internal::ZeroSign),
00093         m_isInitialized(false)
00094     {}
00095 
00096     /** \brief Constructor with decomposition
00097       *
00098       * This calculates the decomposition for the input \a matrix.
00099       * \sa LDLT(Index size)
00100       */
00101     LDLT(const MatrixType& matrix)
00102       : m_matrix(matrix.rows(), matrix.cols()),
00103         m_transpositions(matrix.rows()),
00104         m_temporary(matrix.rows()),
00105         m_sign(internal::ZeroSign),
00106         m_isInitialized(false)
00107     {
00108       compute(matrix);
00109     }
00110 
00111     /** Clear any existing decomposition
00112      * \sa rankUpdate(w,sigma)
00113      */
00114     void setZero()
00115     {
00116       m_isInitialized = false;
00117     }
00118 
00119     /** \returns a view of the upper triangular matrix U */
00120     inline typename Traits::MatrixU matrixU () const
00121     {
00122       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00123       return Traits::getU(m_matrix);
00124     }
00125 
00126     /** \returns a view of the lower triangular matrix L */
00127     inline typename Traits::MatrixL matrixL () const
00128     {
00129       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00130       return Traits::getL(m_matrix);
00131     }
00132 
00133     /** \returns the permutation matrix P as a transposition sequence.
00134       */
00135     inline const TranspositionType & transpositionsP () const
00136     {
00137       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00138       return m_transpositions;
00139     }
00140 
00141     /** \returns the coefficients of the diagonal matrix D */
00142     inline Diagonal<const MatrixType> vectorD () const
00143     {
00144       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00145       return m_matrix.diagonal();
00146     }
00147 
00148     /** \returns true if the matrix is positive (semidefinite) */
00149     inline bool isPositive () const
00150     {
00151       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00152       return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
00153     }
00154     
00155     #ifdef EIGEN2_SUPPORT
00156     inline bool isPositiveDefinite() const
00157     {
00158       return isPositive ();
00159     }
00160     #endif
00161 
00162     /** \returns true if the matrix is negative (semidefinite) */
00163     inline bool isNegative (void) const
00164     {
00165       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00166       return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
00167     }
00168 
00169     /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
00170       *
00171       * This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> .
00172       *
00173       * \note_about_checking_solutions
00174       *
00175       * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$
00176       * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, 
00177       * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then
00178       * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the
00179       * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
00180       * computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular.
00181       *
00182       * \sa MatrixBase::ldlt()
00183       */
00184     template<typename Rhs>
00185     inline const internal::solve_retval<LDLT, Rhs>
00186     solve (const MatrixBase<Rhs>& b) const
00187     {
00188       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00189       eigen_assert(m_matrix.rows()==b.rows()
00190                 && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
00191       return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
00192     }
00193 
00194     #ifdef EIGEN2_SUPPORT
00195     template<typename OtherDerived, typename ResultType>
00196     bool solve (const MatrixBase<OtherDerived>& b, ResultType *result) const
00197     {
00198       *result = this->solve (b);
00199       return true;
00200     }
00201     #endif
00202 
00203     template<typename Derived>
00204     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
00205 
00206     LDLT& compute(const MatrixType& matrix);
00207 
00208     template <typename Derived>
00209     LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1);
00210 
00211     /** \returns the internal LDLT decomposition matrix
00212       *
00213       * TODO: document the storage layout
00214       */
00215     inline const MatrixType& matrixLDLT () const
00216     {
00217       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00218       return m_matrix;
00219     }
00220 
00221     MatrixType reconstructedMatrix () const;
00222 
00223     inline Index rows() const { return m_matrix.rows(); }
00224     inline Index cols() const { return m_matrix.cols(); }
00225 
00226     /** \brief Reports whether previous computation was successful.
00227       *
00228       * \returns \c Success if computation was succesful,
00229       *          \c NumericalIssue if the matrix.appears to be negative.
00230       */
00231     ComputationInfo info() const
00232     {
00233       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00234       return Success;
00235     }
00236 
00237   protected:
00238     
00239     static void check_template_parameters()
00240     {
00241       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00242     }
00243 
00244     /** \internal
00245       * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
00246       * The strict upper part is used during the decomposition, the strict lower
00247       * part correspond to the coefficients of L (its diagonal is equal to 1 and
00248       * is not stored), and the diagonal entries correspond to D.
00249       */
00250     MatrixType m_matrix;
00251     TranspositionType m_transpositions;
00252     TmpMatrixType m_temporary;
00253     internal::SignMatrix m_sign;
00254     bool m_isInitialized;
00255 };
00256 
00257 namespace internal {
00258 
00259 template<int UpLo> struct ldlt_inplace;
00260 
00261 template<> struct ldlt_inplace<Lower>
00262 {
00263   template<typename MatrixType, typename TranspositionType, typename Workspace>
00264   static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
00265   {
00266     using std::abs;
00267     typedef typename MatrixType::Scalar Scalar;
00268     typedef typename MatrixType::RealScalar RealScalar;
00269     typedef typename MatrixType::Index Index;
00270     eigen_assert(mat.rows()==mat.cols());
00271     const Index size = mat.rows();
00272 
00273     if (size <= 1)
00274     {
00275       transpositions.setIdentity();
00276       if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef;
00277       else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef;
00278       else sign = ZeroSign;
00279       return true;
00280     }
00281 
00282     for (Index k = 0; k < size; ++k)
00283     {
00284       // Find largest diagonal element
00285       Index index_of_biggest_in_corner;
00286       mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
00287       index_of_biggest_in_corner += k;
00288 
00289       transpositions.coeffRef(k) = index_of_biggest_in_corner;
00290       if(k != index_of_biggest_in_corner)
00291       {
00292         // apply the transposition while taking care to consider only
00293         // the lower triangular part
00294         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
00295         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
00296         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
00297         std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
00298         for(int i=k+1;i<index_of_biggest_in_corner;++i)
00299         {
00300           Scalar tmp = mat.coeffRef(i,k);
00301           mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i));
00302           mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp);
00303         }
00304         if(NumTraits<Scalar>::IsComplex)
00305           mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k));
00306       }
00307 
00308       // partition the matrix:
00309       //       A00 |  -  |  -
00310       // lu  = A10 | A11 |  -
00311       //       A20 | A21 | A22
00312       Index rs = size - k - 1;
00313       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00314       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00315       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00316 
00317       if(k>0)
00318       {
00319         temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
00320         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
00321         if(rs>0)
00322           A21.noalias() -= A20 * temp.head(k);
00323       }
00324       
00325       // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot
00326       // was smaller than the cutoff value. However, soince LDLT is not rank-revealing
00327       // we should only make sure we do not introduce INF or NaN values.
00328       // LAPACK also uses 0 as the cutoff value.
00329       RealScalar realAkk = numext::real(mat.coeffRef(k,k));
00330       if((rs>0) && (abs(realAkk) > RealScalar(0)))
00331         A21 /= realAkk;
00332 
00333       if (sign == PositiveSemiDef) {
00334         if (realAkk < 0) sign = Indefinite;
00335       } else if (sign == NegativeSemiDef) {
00336         if (realAkk > 0) sign = Indefinite;
00337       } else if (sign == ZeroSign) {
00338         if (realAkk > 0) sign = PositiveSemiDef;
00339         else if (realAkk < 0) sign = NegativeSemiDef;
00340       }
00341     }
00342 
00343     return true;
00344   }
00345 
00346   // Reference for the algorithm: Davis and Hager, "Multiple Rank
00347   // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
00348   // Trivial rearrangements of their computations (Timothy E. Holy)
00349   // allow their algorithm to work for rank-1 updates even if the
00350   // original matrix is not of full rank.
00351   // Here only rank-1 updates are implemented, to reduce the
00352   // requirement for intermediate storage and improve accuracy
00353   template<typename MatrixType, typename WDerived>
00354   static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1)
00355   {
00356     using numext::isfinite;
00357     typedef typename MatrixType::Scalar Scalar;
00358     typedef typename MatrixType::RealScalar RealScalar;
00359     typedef typename MatrixType::Index Index;
00360 
00361     const Index size = mat.rows();
00362     eigen_assert(mat.cols() == size && w.size()==size);
00363 
00364     RealScalar alpha = 1;
00365 
00366     // Apply the update
00367     for (Index j = 0; j < size; j++)
00368     {
00369       // Check for termination due to an original decomposition of low-rank
00370       if (!(isfinite)(alpha))
00371         break;
00372 
00373       // Update the diagonal terms
00374       RealScalar dj = numext::real(mat.coeff(j,j));
00375       Scalar wj = w.coeff(j);
00376       RealScalar swj2 = sigma*numext::abs2(wj);
00377       RealScalar gamma = dj*alpha + swj2;
00378 
00379       mat.coeffRef(j,j) += swj2/alpha;
00380       alpha += swj2/dj;
00381 
00382 
00383       // Update the terms of L
00384       Index rs = size-j-1;
00385       w.tail(rs) -= wj * mat.col(j).tail(rs);
00386       if(gamma != 0)
00387         mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
00388     }
00389     return true;
00390   }
00391 
00392   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00393   static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
00394   {
00395     // Apply the permutation to the input w
00396     tmp = transpositions * w;
00397 
00398     return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
00399   }
00400 };
00401 
00402 template<> struct ldlt_inplace<Upper>
00403 {
00404   template<typename MatrixType, typename TranspositionType, typename Workspace>
00405   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
00406   {
00407     Transpose<MatrixType> matt(mat);
00408     return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
00409   }
00410 
00411   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00412   static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1)
00413   {
00414     Transpose<MatrixType> matt(mat);
00415     return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
00416   }
00417 };
00418 
00419 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
00420 {
00421   typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
00422   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
00423   static inline MatrixL getL(const MatrixType& m) { return m; }
00424   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00425 };
00426 
00427 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
00428 {
00429   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
00430   typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
00431   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00432   static inline MatrixU getU(const MatrixType& m) { return m; }
00433 };
00434 
00435 } // end namespace internal
00436 
00437 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
00438   */
00439 template<typename MatrixType, int _UpLo>
00440 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00441 {
00442   check_template_parameters();
00443   
00444   eigen_assert(a.rows()==a.cols());
00445   const Index size = a.rows();
00446 
00447   m_matrix = a;
00448 
00449   m_transpositions.resize(size);
00450   m_isInitialized = false;
00451   m_temporary.resize(size);
00452   m_sign = internal::ZeroSign;
00453 
00454   internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
00455 
00456   m_isInitialized = true;
00457   return *this;
00458 }
00459 
00460 /** Update the LDLT decomposition:  given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
00461  * \param w a vector to be incorporated into the decomposition.
00462  * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1.
00463  * \sa setZero()
00464   */
00465 template<typename MatrixType, int _UpLo>
00466 template<typename Derived>
00467 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename LDLT<MatrixType,_UpLo>::RealScalar& sigma)
00468 {
00469   const Index size = w.rows();
00470   if (m_isInitialized)
00471   {
00472     eigen_assert(m_matrix.rows()==size);
00473   }
00474   else
00475   {    
00476     m_matrix.resize(size,size);
00477     m_matrix.setZero();
00478     m_transpositions.resize(size);
00479     for (Index i = 0; i < size; i++)
00480       m_transpositions.coeffRef(i) = i;
00481     m_temporary.resize(size);
00482     m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
00483     m_isInitialized = true;
00484   }
00485 
00486   internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
00487 
00488   return *this;
00489 }
00490 
00491 namespace internal {
00492 template<typename _MatrixType, int _UpLo, typename Rhs>
00493 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
00494   : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
00495 {
00496   typedef LDLT<_MatrixType,_UpLo> LDLTType;
00497   EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
00498 
00499   template<typename Dest> void evalTo(Dest& dst) const
00500   {
00501     eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
00502     // dst = P b
00503     dst = dec().transpositionsP() * rhs();
00504 
00505     // dst = L^-1 (P b)
00506     dec().matrixL().solveInPlace(dst);
00507 
00508     // dst = D^-1 (L^-1 P b)
00509     // more precisely, use pseudo-inverse of D (see bug 241)
00510     using std::abs;
00511     using std::max;
00512     typedef typename LDLTType::MatrixType MatrixType;
00513     typedef typename LDLTType::RealScalar RealScalar;
00514     const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD());
00515     // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
00516     // as motivated by LAPACK's xGELSS:
00517     // RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
00518     // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest
00519     // diagonal element is not well justified and to numerical issues in some cases.
00520     // Moreover, Lapack's xSYTRS routines use 0 for the tolerance.
00521     RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest();
00522     
00523     for (Index i = 0; i < vectorD.size(); ++i) {
00524       if(abs(vectorD(i)) > tolerance)
00525         dst.row(i) /= vectorD(i);
00526       else
00527         dst.row(i).setZero();
00528     }
00529 
00530     // dst = L^-T (D^-1 L^-1 P b)
00531     dec().matrixU().solveInPlace(dst);
00532 
00533     // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
00534     dst = dec().transpositionsP().transpose() * dst;
00535   }
00536 };
00537 }
00538 
00539 /** \internal use x = ldlt_object.solve(x);
00540   *
00541   * This is the \em in-place version of solve().
00542   *
00543   * \param bAndX represents both the right-hand side matrix b and result x.
00544   *
00545   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
00546   *
00547   * This version avoids a copy when the right hand side matrix b is not
00548   * needed anymore.
00549   *
00550   * \sa LDLT::solve(), MatrixBase::ldlt()
00551   */
00552 template<typename MatrixType,int _UpLo>
00553 template<typename Derived>
00554 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00555 {
00556   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00557   eigen_assert(m_matrix.rows() == bAndX.rows());
00558 
00559   bAndX = this->solve(bAndX);
00560 
00561   return true;
00562 }
00563 
00564 /** \returns the matrix represented by the decomposition,
00565  * i.e., it returns the product: P^T L D L^* P.
00566  * This function is provided for debug purpose. */
00567 template<typename MatrixType, int _UpLo>
00568 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix () const
00569 {
00570   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00571   const Index size = m_matrix.rows();
00572   MatrixType res(size,size);
00573 
00574   // P
00575   res.setIdentity();
00576   res = transpositionsP() * res;
00577   // L^* P
00578   res = matrixU() * res;
00579   // D(L^*P)
00580   res = vectorD().real().asDiagonal() * res;
00581   // L(DL^*P)
00582   res = matrixL() * res;
00583   // P^T (LDL^*P)
00584   res = transpositionsP().transpose() * res;
00585 
00586   return res;
00587 }
00588 
00589 /** \cholesky_module
00590   * \returns the Cholesky decomposition with full pivoting without square root of \c *this
00591   */
00592 template<typename MatrixType, unsigned int UpLo>
00593 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00594 SelfAdjointView<MatrixType, UpLo>::ldlt () const
00595 {
00596   return LDLT<PlainObject,UpLo>(m_matrix);
00597 }
00598 
00599 /** \cholesky_module
00600   * \returns the Cholesky decomposition with full pivoting without square root of \c *this
00601   */
00602 template<typename Derived>
00603 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
00604 MatrixBase<Derived>::ldlt () const
00605 {
00606   return LDLT<PlainObject>(derived());
00607 }
00608 
00609 } // end namespace Eigen
00610 
00611 #endif // EIGEN_LDLT_H