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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_LLT_H
00026 #define EIGEN_LLT_H
00027 
00028 namespace internal{
00029 template<typename MatrixType, int UpLo> struct LLT_Traits;
00030 }
00031 
00054  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00055   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00056   * the strict lower part does not have to store correct values.
00057   */
00058 template<typename _MatrixType, int _UpLo> class LLT
00059 {
00060   public:
00061     typedef _MatrixType MatrixType;
00062     enum {
00063       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065       Options = MatrixType::Options,
00066       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00067     };
00068     typedef typename MatrixType::Scalar Scalar;
00069     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00070     typedef typename MatrixType::Index Index;
00071 
00072     enum {
00073       PacketSize = internal::packet_traits<Scalar>::size,
00074       AlignmentMask = int(PacketSize)-1,
00075       UpLo = _UpLo
00076     };
00077 
00078     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00079 
00086     LLT() : m_matrix(), m_isInitialized(false) {}
00087 
00094     LLT(Index size) : m_matrix(size, size),
00095                     m_isInitialized(false) {}
00096 
00097     LLT(const MatrixType& matrix)
00098       : m_matrix(matrix.rows(), matrix.cols()),
00099         m_isInitialized(false)
00100     {
00101       compute(matrix);
00102     }
00103 
00105     inline typename Traits::MatrixU matrixU() const
00106     {
00107       eigen_assert(m_isInitialized && "LLT is not initialized.");
00108       return Traits::getU(m_matrix);
00109     }
00110 
00112     inline typename Traits::MatrixL matrixL() const
00113     {
00114       eigen_assert(m_isInitialized && "LLT is not initialized.");
00115       return Traits::getL(m_matrix);
00116     }
00117 
00128     template<typename Rhs>
00129     inline const internal::solve_retval<LLT, Rhs>
00130     solve(const MatrixBase<Rhs>& b) const
00131     {
00132       eigen_assert(m_isInitialized && "LLT is not initialized.");
00133       eigen_assert(m_matrix.rows()==b.rows()
00134                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00135       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00136     }
00137 
00138     #ifdef EIGEN2_SUPPORT
00139     template<typename OtherDerived, typename ResultType>
00140     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00141     {
00142       *result = this->solve(b);
00143       return true;
00144     }
00145     
00146     bool isPositiveDefinite() const { return true; }
00147     #endif
00148 
00149     template<typename Derived>
00150     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00151 
00152     LLT& compute(const MatrixType& matrix);
00153 
00158     inline const MatrixType& matrixLLT() const
00159     {
00160       eigen_assert(m_isInitialized && "LLT is not initialized.");
00161       return m_matrix;
00162     }
00163 
00164     MatrixType reconstructedMatrix() const;
00165 
00166 
00172     ComputationInfo info() const
00173     {
00174       eigen_assert(m_isInitialized && "LLT is not initialized.");
00175       return m_info;
00176     }
00177 
00178     inline Index rows() const { return m_matrix.rows(); }
00179     inline Index cols() const { return m_matrix.cols(); }
00180 
00181   protected:
00186     MatrixType m_matrix;
00187     bool m_isInitialized;
00188     ComputationInfo m_info;
00189 };
00190 
00191 namespace internal {
00192 
00193 template<int UpLo> struct llt_inplace;
00194 
00195 template<> struct llt_inplace<Lower>
00196 {
00197   template<typename MatrixType>
00198   static typename MatrixType::Index unblocked(MatrixType& mat)
00199   {
00200     typedef typename MatrixType::Index Index;
00201     typedef typename MatrixType::Scalar Scalar;
00202     typedef typename MatrixType::RealScalar RealScalar;
00203     
00204     eigen_assert(mat.rows()==mat.cols());
00205     const Index size = mat.rows();
00206     for(Index k = 0; k < size; ++k)
00207     {
00208       Index rs = size-k-1; // remaining size
00209 
00210       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00211       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00212       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00213 
00214       RealScalar x = real(mat.coeff(k,k));
00215       if (k>0) x -= A10.squaredNorm();
00216       if (x<=RealScalar(0))
00217         return k;
00218       mat.coeffRef(k,k) = x = sqrt(x);
00219       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00220       if (rs>0) A21 *= RealScalar(1)/x;
00221     }
00222     return -1;
00223   }
00224 
00225   template<typename MatrixType>
00226   static typename MatrixType::Index blocked(MatrixType& m)
00227   {
00228     typedef typename MatrixType::Index Index;
00229     eigen_assert(m.rows()==m.cols());
00230     Index size = m.rows();
00231     if(size<32)
00232       return unblocked(m);
00233 
00234     Index blockSize = size/8;
00235     blockSize = (blockSize/16)*16;
00236     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00237 
00238     for (Index k=0; k<size; k+=blockSize)
00239     {
00240       // partition the matrix:
00241       //       A00 |  -  |  -
00242       // lu  = A10 | A11 |  -
00243       //       A20 | A21 | A22
00244       Index bs = (std::min)(blockSize, size-k);
00245       Index rs = size - k - bs;
00246       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00247       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00248       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00249 
00250       Index ret;
00251       if((ret=unblocked(A11))>=0) return k+ret;
00252       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00253       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00254     }
00255     return -1;
00256   }
00257 };
00258 
00259 template<> struct llt_inplace<Upper>
00260 {
00261   template<typename MatrixType>
00262   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00263   {
00264     Transpose<MatrixType> matt(mat);
00265     return llt_inplace<Lower>::unblocked(matt);
00266   }
00267   template<typename MatrixType>
00268   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00269   {
00270     Transpose<MatrixType> matt(mat);
00271     return llt_inplace<Lower>::blocked(matt);
00272   }
00273 };
00274 
00275 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00276 {
00277   typedef TriangularView<MatrixType, Lower> MatrixL;
00278   typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU;
00279   inline static MatrixL getL(const MatrixType& m) { return m; }
00280   inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00281   static bool inplace_decomposition(MatrixType& m)
00282   { return llt_inplace<Lower>::blocked(m)==-1; }
00283 };
00284 
00285 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00286 {
00287   typedef TriangularView<typename MatrixType::AdjointReturnType, Lower> MatrixL;
00288   typedef TriangularView<MatrixType, Upper> MatrixU;
00289   inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00290   inline static MatrixU getU(const MatrixType& m) { return m; }
00291   static bool inplace_decomposition(MatrixType& m)
00292   { return llt_inplace<Upper>::blocked(m)==-1; }
00293 };
00294 
00295 } // end namespace internal
00296 
00302 template<typename MatrixType, int _UpLo>
00303 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00304 {
00305   assert(a.rows()==a.cols());
00306   const Index size = a.rows();
00307   m_matrix.resize(size, size);
00308   m_matrix = a;
00309 
00310   m_isInitialized = true;
00311   bool ok = Traits::inplace_decomposition(m_matrix);
00312   m_info = ok ? Success : NumericalIssue;
00313 
00314   return *this;
00315 }
00316 
00317 namespace internal {
00318 template<typename _MatrixType, int UpLo, typename Rhs>
00319 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00320   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00321 {
00322   typedef LLT<_MatrixType,UpLo> LLTType;
00323   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00324 
00325   template<typename Dest> void evalTo(Dest& dst) const
00326   {
00327     dst = rhs();
00328     dec().solveInPlace(dst);
00329   }
00330 };
00331 }
00332 
00346 template<typename MatrixType, int _UpLo>
00347 template<typename Derived>
00348 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00349 {
00350   eigen_assert(m_isInitialized && "LLT is not initialized.");
00351   eigen_assert(m_matrix.rows()==bAndX.rows());
00352   matrixL().solveInPlace(bAndX);
00353   matrixU().solveInPlace(bAndX);
00354 }
00355 
00359 template<typename MatrixType, int _UpLo>
00360 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00361 {
00362   eigen_assert(m_isInitialized && "LLT is not initialized.");
00363   return matrixL() * matrixL().adjoint().toDenseMatrix();
00364 }
00365 
00369 template<typename Derived>
00370 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00371 MatrixBase<Derived>::llt() const
00372 {
00373   return LLT<PlainObject>(derived());
00374 }
00375 
00379 template<typename MatrixType, unsigned int UpLo>
00380 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00381 SelfAdjointView<MatrixType, UpLo>::llt() const
00382 {
00383   return LLT<PlainObject,UpLo>(m_matrix);
00384 }
00385 
00386 #endif // EIGEN_LLT_H


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