LLT.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 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_LLT_H
00011 #define EIGEN_LLT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal{
00016 template<typename MatrixType, int UpLo> struct LLT_Traits;
00017 }
00018 
00046  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00047   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00048   * the strict lower part does not have to store correct values.
00049   */
00050 template<typename _MatrixType, int _UpLo> class LLT
00051 {
00052   public:
00053     typedef _MatrixType MatrixType;
00054     enum {
00055       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00056       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00057       Options = MatrixType::Options,
00058       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00059     };
00060     typedef typename MatrixType::Scalar Scalar;
00061     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00062     typedef typename MatrixType::Index Index;
00063 
00064     enum {
00065       PacketSize = internal::packet_traits<Scalar>::size,
00066       AlignmentMask = int(PacketSize)-1,
00067       UpLo = _UpLo
00068     };
00069 
00070     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00071 
00078     LLT() : m_matrix(), m_isInitialized(false) {}
00079 
00086     LLT(Index size) : m_matrix(size, size),
00087                     m_isInitialized(false) {}
00088 
00089     LLT(const MatrixType& matrix)
00090       : m_matrix(matrix.rows(), matrix.cols()),
00091         m_isInitialized(false)
00092     {
00093       compute(matrix);
00094     }
00095 
00097     inline typename Traits::MatrixU matrixU() const
00098     {
00099       eigen_assert(m_isInitialized && "LLT is not initialized.");
00100       return Traits::getU(m_matrix);
00101     }
00102 
00104     inline typename Traits::MatrixL matrixL() const
00105     {
00106       eigen_assert(m_isInitialized && "LLT is not initialized.");
00107       return Traits::getL(m_matrix);
00108     }
00109 
00120     template<typename Rhs>
00121     inline const internal::solve_retval<LLT, Rhs>
00122     solve(const MatrixBase<Rhs>& b) const
00123     {
00124       eigen_assert(m_isInitialized && "LLT is not initialized.");
00125       eigen_assert(m_matrix.rows()==b.rows()
00126                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00127       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00128     }
00129 
00130     #ifdef EIGEN2_SUPPORT
00131     template<typename OtherDerived, typename ResultType>
00132     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00133     {
00134       *result = this->solve(b);
00135       return true;
00136     }
00137     
00138     bool isPositiveDefinite() const { return true; }
00139     #endif
00140 
00141     template<typename Derived>
00142     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00143 
00144     LLT& compute(const MatrixType& matrix);
00145 
00150     inline const MatrixType& matrixLLT() const
00151     {
00152       eigen_assert(m_isInitialized && "LLT is not initialized.");
00153       return m_matrix;
00154     }
00155 
00156     MatrixType reconstructedMatrix() const;
00157 
00158 
00164     ComputationInfo info() const
00165     {
00166       eigen_assert(m_isInitialized && "LLT is not initialized.");
00167       return m_info;
00168     }
00169 
00170     inline Index rows() const { return m_matrix.rows(); }
00171     inline Index cols() const { return m_matrix.cols(); }
00172 
00173     template<typename VectorType>
00174     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
00175 
00176   protected:
00181     MatrixType m_matrix;
00182     bool m_isInitialized;
00183     ComputationInfo m_info;
00184 };
00185 
00186 namespace internal {
00187 
00188 template<typename Scalar, int UpLo> struct llt_inplace;
00189 
00190 template<typename MatrixType, typename VectorType>
00191 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
00192 {
00193   typedef typename MatrixType::Scalar Scalar;
00194   typedef typename MatrixType::RealScalar RealScalar;
00195   typedef typename MatrixType::Index Index;
00196   typedef typename MatrixType::ColXpr ColXpr;
00197   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00198   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00199   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00200   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00201 
00202   int n = mat.cols();
00203   eigen_assert(mat.rows()==n && vec.size()==n);
00204 
00205   TempVectorType temp;
00206 
00207   if(sigma>0)
00208   {
00209     // This version is based on Givens rotations.
00210     // It is faster than the other one below, but only works for updates,
00211     // i.e., for sigma > 0
00212     temp = sqrt(sigma) * vec;
00213 
00214     for(int i=0; i<n; ++i)
00215     {
00216       JacobiRotation<Scalar> g;
00217       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00218 
00219       int rs = n-i-1;
00220       if(rs>0)
00221       {
00222         ColXprSegment x(mat.col(i).tail(rs));
00223         TempVecSegment y(temp.tail(rs));
00224         apply_rotation_in_the_plane(x, y, g);
00225       }
00226     }
00227   }
00228   else
00229   {
00230     temp = vec;
00231     RealScalar beta = 1;
00232     for(int j=0; j<n; ++j)
00233     {
00234       RealScalar Ljj = real(mat.coeff(j,j));
00235       RealScalar dj = abs2(Ljj);
00236       Scalar wj = temp.coeff(j);
00237       RealScalar swj2 = sigma*abs2(wj);
00238       RealScalar gamma = dj*beta + swj2;
00239 
00240       RealScalar x = dj + swj2/beta;
00241       if (x<=RealScalar(0))
00242         return j;
00243       RealScalar nLjj = sqrt(x);
00244       mat.coeffRef(j,j) = nLjj;
00245       beta += swj2/dj;
00246 
00247       // Update the terms of L
00248       Index rs = n-j-1;
00249       if(rs)
00250       {
00251         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00252         if(gamma != 0)
00253           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
00254       }
00255     }
00256   }
00257   return -1;
00258 }
00259 
00260 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00261 {
00262   typedef typename NumTraits<Scalar>::Real RealScalar;
00263   template<typename MatrixType>
00264   static typename MatrixType::Index unblocked(MatrixType& mat)
00265   {
00266     typedef typename MatrixType::Index Index;
00267     
00268     eigen_assert(mat.rows()==mat.cols());
00269     const Index size = mat.rows();
00270     for(Index k = 0; k < size; ++k)
00271     {
00272       Index rs = size-k-1; // remaining size
00273 
00274       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00275       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00276       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00277 
00278       RealScalar x = real(mat.coeff(k,k));
00279       if (k>0) x -= A10.squaredNorm();
00280       if (x<=RealScalar(0))
00281         return k;
00282       mat.coeffRef(k,k) = x = sqrt(x);
00283       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00284       if (rs>0) A21 *= RealScalar(1)/x;
00285     }
00286     return -1;
00287   }
00288 
00289   template<typename MatrixType>
00290   static typename MatrixType::Index blocked(MatrixType& m)
00291   {
00292     typedef typename MatrixType::Index Index;
00293     eigen_assert(m.rows()==m.cols());
00294     Index size = m.rows();
00295     if(size<32)
00296       return unblocked(m);
00297 
00298     Index blockSize = size/8;
00299     blockSize = (blockSize/16)*16;
00300     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00301 
00302     for (Index k=0; k<size; k+=blockSize)
00303     {
00304       // partition the matrix:
00305       //       A00 |  -  |  -
00306       // lu  = A10 | A11 |  -
00307       //       A20 | A21 | A22
00308       Index bs = (std::min)(blockSize, size-k);
00309       Index rs = size - k - bs;
00310       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00311       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00312       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00313 
00314       Index ret;
00315       if((ret=unblocked(A11))>=0) return k+ret;
00316       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00317       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00318     }
00319     return -1;
00320   }
00321 
00322   template<typename MatrixType, typename VectorType>
00323   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00324   {
00325     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00326   }
00327 };
00328   
00329 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00330 {
00331   typedef typename NumTraits<Scalar>::Real RealScalar;
00332 
00333   template<typename MatrixType>
00334   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00335   {
00336     Transpose<MatrixType> matt(mat);
00337     return llt_inplace<Scalar, Lower>::unblocked(matt);
00338   }
00339   template<typename MatrixType>
00340   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00341   {
00342     Transpose<MatrixType> matt(mat);
00343     return llt_inplace<Scalar, Lower>::blocked(matt);
00344   }
00345   template<typename MatrixType, typename VectorType>
00346   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00347   {
00348     Transpose<MatrixType> matt(mat);
00349     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00350   }
00351 };
00352 
00353 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00354 {
00355   typedef const TriangularView<const MatrixType, Lower> MatrixL;
00356   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00357   static inline MatrixL getL(const MatrixType& m) { return m; }
00358   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00359   static bool inplace_decomposition(MatrixType& m)
00360   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00361 };
00362 
00363 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00364 {
00365   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00366   typedef const TriangularView<const MatrixType, Upper> MatrixU;
00367   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00368   static inline MatrixU getU(const MatrixType& m) { return m; }
00369   static bool inplace_decomposition(MatrixType& m)
00370   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00371 };
00372 
00373 } // end namespace internal
00374 
00382 template<typename MatrixType, int _UpLo>
00383 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00384 {
00385   eigen_assert(a.rows()==a.cols());
00386   const Index size = a.rows();
00387   m_matrix.resize(size, size);
00388   m_matrix = a;
00389 
00390   m_isInitialized = true;
00391   bool ok = Traits::inplace_decomposition(m_matrix);
00392   m_info = ok ? Success : NumericalIssue;
00393 
00394   return *this;
00395 }
00396 
00402 template<typename _MatrixType, int _UpLo>
00403 template<typename VectorType>
00404 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00405 {
00406   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00407   eigen_assert(v.size()==m_matrix.cols());
00408   eigen_assert(m_isInitialized);
00409   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00410     m_info = NumericalIssue;
00411   else
00412     m_info = Success;
00413 
00414   return *this;
00415 }
00416     
00417 namespace internal {
00418 template<typename _MatrixType, int UpLo, typename Rhs>
00419 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00420   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00421 {
00422   typedef LLT<_MatrixType,UpLo> LLTType;
00423   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00424 
00425   template<typename Dest> void evalTo(Dest& dst) const
00426   {
00427     dst = rhs();
00428     dec().solveInPlace(dst);
00429   }
00430 };
00431 }
00432 
00446 template<typename MatrixType, int _UpLo>
00447 template<typename Derived>
00448 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00449 {
00450   eigen_assert(m_isInitialized && "LLT is not initialized.");
00451   eigen_assert(m_matrix.rows()==bAndX.rows());
00452   matrixL().solveInPlace(bAndX);
00453   matrixU().solveInPlace(bAndX);
00454 }
00455 
00459 template<typename MatrixType, int _UpLo>
00460 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00461 {
00462   eigen_assert(m_isInitialized && "LLT is not initialized.");
00463   return matrixL() * matrixL().adjoint().toDenseMatrix();
00464 }
00465 
00469 template<typename Derived>
00470 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00471 MatrixBase<Derived>::llt() const
00472 {
00473   return LLT<PlainObject>(derived());
00474 }
00475 
00479 template<typename MatrixType, unsigned int UpLo>
00480 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00481 SelfAdjointView<MatrixType, UpLo>::llt() const
00482 {
00483   return LLT<PlainObject,UpLo>(m_matrix);
00484 }
00485 
00486 } // end namespace Eigen
00487 
00488 #endif // EIGEN_LLT_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:11