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   using std::sqrt;
00194   typedef typename MatrixType::Scalar Scalar;
00195   typedef typename MatrixType::RealScalar RealScalar;
00196   typedef typename MatrixType::Index Index;
00197   typedef typename MatrixType::ColXpr ColXpr;
00198   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00199   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00200   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00201   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00202 
00203   Index n = mat.cols();
00204   eigen_assert(mat.rows()==n && vec.size()==n);
00205 
00206   TempVectorType temp;
00207 
00208   if(sigma>0)
00209   {
00210     // This version is based on Givens rotations.
00211     // It is faster than the other one below, but only works for updates,
00212     // i.e., for sigma > 0
00213     temp = sqrt(sigma) * vec;
00214 
00215     for(Index i=0; i<n; ++i)
00216     {
00217       JacobiRotation<Scalar> g;
00218       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00219 
00220       Index rs = n-i-1;
00221       if(rs>0)
00222       {
00223         ColXprSegment x(mat.col(i).tail(rs));
00224         TempVecSegment y(temp.tail(rs));
00225         apply_rotation_in_the_plane(x, y, g);
00226       }
00227     }
00228   }
00229   else
00230   {
00231     temp = vec;
00232     RealScalar beta = 1;
00233     for(Index j=0; j<n; ++j)
00234     {
00235       RealScalar Ljj = numext::real(mat.coeff(j,j));
00236       RealScalar dj = numext::abs2(Ljj);
00237       Scalar wj = temp.coeff(j);
00238       RealScalar swj2 = sigma*numext::abs2(wj);
00239       RealScalar gamma = dj*beta + swj2;
00240 
00241       RealScalar x = dj + swj2/beta;
00242       if (x<=RealScalar(0))
00243         return j;
00244       RealScalar nLjj = sqrt(x);
00245       mat.coeffRef(j,j) = nLjj;
00246       beta += swj2/dj;
00247 
00248       // Update the terms of L
00249       Index rs = n-j-1;
00250       if(rs)
00251       {
00252         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00253         if(gamma != 0)
00254           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
00255       }
00256     }
00257   }
00258   return -1;
00259 }
00260 
00261 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00262 {
00263   typedef typename NumTraits<Scalar>::Real RealScalar;
00264   template<typename MatrixType>
00265   static typename MatrixType::Index unblocked(MatrixType& mat)
00266   {
00267     using std::sqrt;
00268     typedef typename MatrixType::Index Index;
00269     
00270     eigen_assert(mat.rows()==mat.cols());
00271     const Index size = mat.rows();
00272     for(Index k = 0; k < size; ++k)
00273     {
00274       Index rs = size-k-1; // remaining size
00275 
00276       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00277       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00278       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00279 
00280       RealScalar x = numext::real(mat.coeff(k,k));
00281       if (k>0) x -= A10.squaredNorm();
00282       if (x<=RealScalar(0))
00283         return k;
00284       mat.coeffRef(k,k) = x = sqrt((double)x);
00285       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00286       if (rs>0) A21 *= RealScalar(1)/x;
00287     }
00288     return -1;
00289   }
00290 
00291   template<typename MatrixType>
00292   static typename MatrixType::Index blocked(MatrixType& m)
00293   {
00294     typedef typename MatrixType::Index Index;
00295     eigen_assert(m.rows()==m.cols());
00296     Index size = m.rows();
00297     if(size<32)
00298       return unblocked(m);
00299 
00300     Index blockSize = size/8;
00301     blockSize = (blockSize/16)*16;
00302     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00303 
00304     for (Index k=0; k<size; k+=blockSize)
00305     {
00306       // partition the matrix:
00307       //       A00 |  -  |  -
00308       // lu  = A10 | A11 |  -
00309       //       A20 | A21 | A22
00310       Index bs = (std::min)(blockSize, size-k);
00311       Index rs = size - k - bs;
00312       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00313       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00314       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00315 
00316       Index ret;
00317       if((ret=unblocked(A11))>=0) return k+ret;
00318       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00319       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00320     }
00321     return -1;
00322   }
00323 
00324   template<typename MatrixType, typename VectorType>
00325   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00326   {
00327     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00328   }
00329 };
00330   
00331 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00332 {
00333   typedef typename NumTraits<Scalar>::Real RealScalar;
00334 
00335   template<typename MatrixType>
00336   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00337   {
00338     Transpose<MatrixType> matt(mat);
00339     return llt_inplace<Scalar, Lower>::unblocked(matt);
00340   }
00341   template<typename MatrixType>
00342   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00343   {
00344     Transpose<MatrixType> matt(mat);
00345     return llt_inplace<Scalar, Lower>::blocked(matt);
00346   }
00347   template<typename MatrixType, typename VectorType>
00348   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00349   {
00350     Transpose<MatrixType> matt(mat);
00351     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00352   }
00353 };
00354 
00355 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00356 {
00357   typedef const TriangularView<const MatrixType, Lower> MatrixL;
00358   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00359   static inline MatrixL getL(const MatrixType& m) { return m; }
00360   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00361   static bool inplace_decomposition(MatrixType& m)
00362   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00363 };
00364 
00365 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00366 {
00367   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00368   typedef const TriangularView<const MatrixType, Upper> MatrixU;
00369   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00370   static inline MatrixU getU(const MatrixType& m) { return m; }
00371   static bool inplace_decomposition(MatrixType& m)
00372   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00373 };
00374 
00375 } // end namespace internal
00376 
00384 template<typename MatrixType, int _UpLo>
00385 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00386 {
00387   eigen_assert(a.rows()==a.cols());
00388   const Index size = a.rows();
00389   m_matrix.resize(size, size);
00390   m_matrix = a;
00391 
00392   m_isInitialized = true;
00393   bool ok = Traits::inplace_decomposition(m_matrix);
00394   m_info = ok ? Success : NumericalIssue;
00395 
00396   return *this;
00397 }
00398 
00404 template<typename _MatrixType, int _UpLo>
00405 template<typename VectorType>
00406 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00407 {
00408   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00409   eigen_assert(v.size()==m_matrix.cols());
00410   eigen_assert(m_isInitialized);
00411   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00412     m_info = NumericalIssue;
00413   else
00414     m_info = Success;
00415 
00416   return *this;
00417 }
00418     
00419 namespace internal {
00420 template<typename _MatrixType, int UpLo, typename Rhs>
00421 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00422   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00423 {
00424   typedef LLT<_MatrixType,UpLo> LLTType;
00425   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00426 
00427   template<typename Dest> void evalTo(Dest& dst) const
00428   {
00429     dst = rhs();
00430     dec().solveInPlace(dst);
00431   }
00432 };
00433 }
00434 
00448 template<typename MatrixType, int _UpLo>
00449 template<typename Derived>
00450 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00451 {
00452   eigen_assert(m_isInitialized && "LLT is not initialized.");
00453   eigen_assert(m_matrix.rows()==bAndX.rows());
00454   matrixL().solveInPlace(bAndX);
00455   matrixU().solveInPlace(bAndX);
00456 }
00457 
00461 template<typename MatrixType, int _UpLo>
00462 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00463 {
00464   eigen_assert(m_isInitialized && "LLT is not initialized.");
00465   return matrixL() * matrixL().adjoint().toDenseMatrix();
00466 }
00467 
00471 template<typename Derived>
00472 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00473 MatrixBase<Derived>::llt() const
00474 {
00475   return LLT<PlainObject>(derived());
00476 }
00477 
00481 template<typename MatrixType, unsigned int UpLo>
00482 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00483 SelfAdjointView<MatrixType, UpLo>::llt() const
00484 {
00485   return LLT<PlainObject,UpLo>(m_matrix);
00486 }
00487 
00488 } // end namespace Eigen
00489 
00490 #endif // EIGEN_LLT_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:59:05