LDLT.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 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 //
00008 // Eigen is free software; you can redistribute it and/or
00009 // modify it under the terms of the GNU Lesser General Public
00010 // License as published by the Free Software Foundation; either
00011 // version 3 of the License, or (at your option) any later version.
00012 //
00013 // Alternatively, you can redistribute it and/or
00014 // modify it under the terms of the GNU General Public License as
00015 // published by the Free Software Foundation; either version 2 of
00016 // the License, or (at your option) any later version.
00017 //
00018 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00019 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00020 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00021 // GNU General Public License for more details.
00022 //
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License and a copy of the GNU General Public License along with
00025 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 #ifndef EIGEN_LDLT_H
00028 #define EIGEN_LDLT_H
00029 
00030 namespace internal {
00031 template<typename MatrixType, int UpLo> struct LDLT_Traits;
00032 }
00033 
00055  /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE
00056   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00057   * the strict lower part does not have to store correct values.
00058   */
00059 template<typename _MatrixType, int _UpLo> class LDLT
00060 {
00061   public:
00062     typedef _MatrixType MatrixType;
00063     enum {
00064       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00065       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00066       Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
00067       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00068       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00069       UpLo = _UpLo
00070     };
00071     typedef typename MatrixType::Scalar Scalar;
00072     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00073     typedef typename MatrixType::Index Index;
00074     typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
00075 
00076     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00077     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00078 
00079     typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
00080 
00086     LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
00087 
00094     LDLT(Index size)
00095       : m_matrix(size, size),
00096         m_transpositions(size),
00097         m_temporary(size),
00098         m_isInitialized(false)
00099     {}
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_isInitialized(false)
00106     {
00107       compute(matrix);
00108     }
00109 
00111     inline typename Traits::MatrixU matrixU() const
00112     {
00113       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00114       return Traits::getU(m_matrix);
00115     }
00116 
00118     inline typename Traits::MatrixL matrixL() const
00119     {
00120       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00121       return Traits::getL(m_matrix);
00122     }
00123 
00126     inline const TranspositionType& transpositionsP() const
00127     {
00128       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00129       return m_transpositions;
00130     }
00131 
00133     inline Diagonal<const MatrixType> vectorD(void) const
00134     {
00135       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00136       return m_matrix.diagonal();
00137     }
00138 
00140     inline bool isPositive(void) const
00141     {
00142       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00143       return m_sign == 1;
00144     }
00145     
00146     #ifdef EIGEN2_SUPPORT
00147     inline bool isPositiveDefinite() const
00148     {
00149       return isPositive();
00150     }
00151     #endif
00152 
00154     inline bool isNegative(void) const
00155     {
00156       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00157       return m_sign == -1;
00158     }
00159 
00166     template<typename Rhs>
00167     inline const internal::solve_retval<LDLT, Rhs>
00168     solve(const MatrixBase<Rhs>& b) const
00169     {
00170       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00171       eigen_assert(m_matrix.rows()==b.rows()
00172                 && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
00173       return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
00174     }
00175 
00176     #ifdef EIGEN2_SUPPORT
00177     template<typename OtherDerived, typename ResultType>
00178     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00179     {
00180       *result = this->solve(b);
00181       return true;
00182     }
00183     #endif
00184 
00185     template<typename Derived>
00186     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
00187 
00188     LDLT& compute(const MatrixType& matrix);
00189 
00194     inline const MatrixType& matrixLDLT() const
00195     {
00196       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00197       return m_matrix;
00198     }
00199 
00200     MatrixType reconstructedMatrix() const;
00201 
00202     inline Index rows() const { return m_matrix.rows(); }
00203     inline Index cols() const { return m_matrix.cols(); }
00204 
00205   protected:
00206 
00213     MatrixType m_matrix;
00214     TranspositionType m_transpositions;
00215     TmpMatrixType m_temporary;
00216     int m_sign;
00217     bool m_isInitialized;
00218 };
00219 
00220 namespace internal {
00221 
00222 template<int UpLo> struct ldlt_inplace;
00223 
00224 template<> struct ldlt_inplace<Lower>
00225 {
00226   template<typename MatrixType, typename TranspositionType, typename Workspace>
00227   static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00228   {
00229     typedef typename MatrixType::Scalar Scalar;
00230     typedef typename MatrixType::RealScalar RealScalar;
00231     typedef typename MatrixType::Index Index;
00232     eigen_assert(mat.rows()==mat.cols());
00233     const Index size = mat.rows();
00234 
00235     if (size <= 1)
00236     {
00237       transpositions.setIdentity();
00238       if(sign)
00239         *sign = real(mat.coeff(0,0))>0 ? 1:-1;
00240       return true;
00241     }
00242 
00243     RealScalar cutoff = 0, biggest_in_corner;
00244 
00245     for (Index k = 0; k < size; ++k)
00246     {
00247       // Find largest diagonal element
00248       Index index_of_biggest_in_corner;
00249       biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
00250       index_of_biggest_in_corner += k;
00251 
00252       if(k == 0)
00253       {
00254         // The biggest overall is the point of reference to which further diagonals
00255         // are compared; if any diagonal is negligible compared
00256         // to the largest overall, the algorithm bails.
00257         cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
00258 
00259         if(sign)
00260           *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
00261       }
00262 
00263       // Finish early if the matrix is not full rank.
00264       if(biggest_in_corner < cutoff)
00265       {
00266         for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
00267         break;
00268       }
00269 
00270       transpositions.coeffRef(k) = index_of_biggest_in_corner;
00271       if(k != index_of_biggest_in_corner)
00272       {
00273         // apply the transposition while taking care to consider only
00274         // the lower triangular part
00275         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
00276         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
00277         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
00278         std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
00279         for(int i=k+1;i<index_of_biggest_in_corner;++i)
00280         {
00281           Scalar tmp = mat.coeffRef(i,k);
00282           mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
00283           mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
00284         }
00285         if(NumTraits<Scalar>::IsComplex)
00286           mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
00287       }
00288 
00289       // partition the matrix:
00290       //       A00 |  -  |  -
00291       // lu  = A10 | A11 |  -
00292       //       A20 | A21 | A22
00293       Index rs = size - k - 1;
00294       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00295       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00296       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00297 
00298       if(k>0)
00299       {
00300         temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
00301         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
00302         if(rs>0)
00303           A21.noalias() -= A20 * temp.head(k);
00304       }
00305       if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
00306         A21 /= mat.coeffRef(k,k);
00307     }
00308 
00309     return true;
00310   }
00311 };
00312 
00313 template<> struct ldlt_inplace<Upper>
00314 {
00315   template<typename MatrixType, typename TranspositionType, typename Workspace>
00316   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00317   {
00318     Transpose<MatrixType> matt(mat);
00319     return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
00320   }
00321 };
00322 
00323 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
00324 {
00325   typedef TriangularView<MatrixType, UnitLower> MatrixL;
00326   typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
00327   inline static MatrixL getL(const MatrixType& m) { return m; }
00328   inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00329 };
00330 
00331 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
00332 {
00333   typedef TriangularView<typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
00334   typedef TriangularView<MatrixType, UnitUpper> MatrixU;
00335   inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00336   inline static MatrixU getU(const MatrixType& m) { return m; }
00337 };
00338 
00339 } // end namespace internal
00340 
00343 template<typename MatrixType, int _UpLo>
00344 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00345 {
00346   eigen_assert(a.rows()==a.cols());
00347   const Index size = a.rows();
00348 
00349   m_matrix = a;
00350 
00351   m_transpositions.resize(size);
00352   m_isInitialized = false;
00353   m_temporary.resize(size);
00354 
00355   internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
00356 
00357   m_isInitialized = true;
00358   return *this;
00359 }
00360 
00361 namespace internal {
00362 template<typename _MatrixType, int _UpLo, typename Rhs>
00363 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
00364   : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
00365 {
00366   typedef LDLT<_MatrixType,_UpLo> LDLTType;
00367   EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
00368 
00369   template<typename Dest> void evalTo(Dest& dst) const
00370   {
00371     eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
00372     // dst = P b
00373     dst = dec().transpositionsP() * rhs();
00374 
00375     // dst = L^-1 (P b)
00376     dec().matrixL().solveInPlace(dst);
00377 
00378     // dst = D^-1 (L^-1 P b)
00379     dst = dec().vectorD().asDiagonal().inverse() * dst;
00380 
00381     // dst = L^-T (D^-1 L^-1 P b)
00382     dec().matrixU().solveInPlace(dst);
00383 
00384     // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
00385     dst = dec().transpositionsP().transpose() * dst;
00386   }
00387 };
00388 }
00389 
00403 template<typename MatrixType,int _UpLo>
00404 template<typename Derived>
00405 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00406 {
00407   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00408   const Index size = m_matrix.rows();
00409   eigen_assert(size == bAndX.rows());
00410 
00411   bAndX = this->solve(bAndX);
00412 
00413   return true;
00414 }
00415 
00419 template<typename MatrixType, int _UpLo>
00420 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
00421 {
00422   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00423   const Index size = m_matrix.rows();
00424   MatrixType res(size,size);
00425 
00426   // P
00427   res.setIdentity();
00428   res = transpositionsP() * res;
00429   // L^* P
00430   res = matrixU() * res;
00431   // D(L^*P)
00432   res = vectorD().asDiagonal() * res;
00433   // L(DL^*P)
00434   res = matrixL() * res;
00435   // P^T (LDL^*P)
00436   res = transpositionsP().transpose() * res;
00437 
00438   return res;
00439 }
00440 
00444 template<typename MatrixType, unsigned int UpLo>
00445 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00446 SelfAdjointView<MatrixType, UpLo>::ldlt() const
00447 {
00448   return LDLT<PlainObject,UpLo>(m_matrix);
00449 }
00450 
00454 template<typename Derived>
00455 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
00456 MatrixBase<Derived>::ldlt() const
00457 {
00458   return LDLT<PlainObject>(derived());
00459 }
00460 
00461 #endif // EIGEN_LDLT_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:56