00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_LLT_H
00026 #define EIGEN_LLT_H
00027
00028 namespace Eigen {
00029
00030 namespace internal{
00031 template<typename MatrixType, int UpLo> struct LLT_Traits;
00032 }
00033
00061
00062
00063
00064
00065 template<typename _MatrixType, int _UpLo> class LLT
00066 {
00067 public:
00068 typedef _MatrixType MatrixType;
00069 enum {
00070 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00071 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00072 Options = MatrixType::Options,
00073 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00074 };
00075 typedef typename MatrixType::Scalar Scalar;
00076 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00077 typedef typename MatrixType::Index Index;
00078
00079 enum {
00080 PacketSize = internal::packet_traits<Scalar>::size,
00081 AlignmentMask = int(PacketSize)-1,
00082 UpLo = _UpLo
00083 };
00084
00085 typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00086
00093 LLT() : m_matrix(), m_isInitialized(false) {}
00094
00101 LLT(Index size) : m_matrix(size, size),
00102 m_isInitialized(false) {}
00103
00104 LLT(const MatrixType& matrix)
00105 : m_matrix(matrix.rows(), matrix.cols()),
00106 m_isInitialized(false)
00107 {
00108 compute(matrix);
00109 }
00110
00112 inline typename Traits::MatrixU matrixU() const
00113 {
00114 eigen_assert(m_isInitialized && "LLT is not initialized.");
00115 return Traits::getU(m_matrix);
00116 }
00117
00119 inline typename Traits::MatrixL matrixL() const
00120 {
00121 eigen_assert(m_isInitialized && "LLT is not initialized.");
00122 return Traits::getL(m_matrix);
00123 }
00124
00135 template<typename Rhs>
00136 inline const internal::solve_retval<LLT, Rhs>
00137 solve(const MatrixBase<Rhs>& b) const
00138 {
00139 eigen_assert(m_isInitialized && "LLT is not initialized.");
00140 eigen_assert(m_matrix.rows()==b.rows()
00141 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00142 return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00143 }
00144
00145 #ifdef EIGEN2_SUPPORT
00146 template<typename OtherDerived, typename ResultType>
00147 bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00148 {
00149 *result = this->solve(b);
00150 return true;
00151 }
00152
00153 bool isPositiveDefinite() const { return true; }
00154 #endif
00155
00156 template<typename Derived>
00157 void solveInPlace(MatrixBase<Derived> &bAndX) const;
00158
00159 LLT& compute(const MatrixType& matrix);
00160
00165 inline const MatrixType& matrixLLT() const
00166 {
00167 eigen_assert(m_isInitialized && "LLT is not initialized.");
00168 return m_matrix;
00169 }
00170
00171 MatrixType reconstructedMatrix() const;
00172
00173
00179 ComputationInfo info() const
00180 {
00181 eigen_assert(m_isInitialized && "LLT is not initialized.");
00182 return m_info;
00183 }
00184
00185 inline Index rows() const { return m_matrix.rows(); }
00186 inline Index cols() const { return m_matrix.cols(); }
00187
00188 template<typename VectorType>
00189 LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
00190
00191 protected:
00196 MatrixType m_matrix;
00197 bool m_isInitialized;
00198 ComputationInfo m_info;
00199 };
00200
00201 namespace internal {
00202
00203 template<typename Scalar, int UpLo> struct llt_inplace;
00204
00205 template<typename MatrixType, typename VectorType>
00206 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
00207 {
00208 typedef typename MatrixType::Scalar Scalar;
00209 typedef typename MatrixType::RealScalar RealScalar;
00210 typedef typename MatrixType::Index Index;
00211 typedef typename MatrixType::ColXpr ColXpr;
00212 typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00213 typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00214 typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00215 typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00216
00217 int n = mat.cols();
00218 eigen_assert(mat.rows()==n && vec.size()==n);
00219
00220 TempVectorType temp;
00221
00222 if(sigma>0)
00223 {
00224
00225
00226
00227 temp = sqrt(sigma) * vec;
00228
00229 for(int i=0; i<n; ++i)
00230 {
00231 JacobiRotation<Scalar> g;
00232 g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00233
00234 int rs = n-i-1;
00235 if(rs>0)
00236 {
00237 ColXprSegment x(mat.col(i).tail(rs));
00238 TempVecSegment y(temp.tail(rs));
00239 apply_rotation_in_the_plane(x, y, g);
00240 }
00241 }
00242 }
00243 else
00244 {
00245 temp = vec;
00246 RealScalar beta = 1;
00247 for(int j=0; j<n; ++j)
00248 {
00249 RealScalar Ljj = real(mat.coeff(j,j));
00250 RealScalar dj = abs2(Ljj);
00251 Scalar wj = temp.coeff(j);
00252 RealScalar swj2 = sigma*abs2(wj);
00253 RealScalar gamma = dj*beta + swj2;
00254
00255 RealScalar x = dj + swj2/beta;
00256 if (x<=RealScalar(0))
00257 return j;
00258 RealScalar nLjj = sqrt(x);
00259 mat.coeffRef(j,j) = nLjj;
00260 beta += swj2/dj;
00261
00262
00263 Index rs = n-j-1;
00264 if(rs)
00265 {
00266 temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00267 if(gamma != 0)
00268 mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
00269 }
00270 }
00271 }
00272 return -1;
00273 }
00274
00275 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00276 {
00277 typedef typename NumTraits<Scalar>::Real RealScalar;
00278 template<typename MatrixType>
00279 static typename MatrixType::Index unblocked(MatrixType& mat)
00280 {
00281 typedef typename MatrixType::Index Index;
00282
00283 eigen_assert(mat.rows()==mat.cols());
00284 const Index size = mat.rows();
00285 for(Index k = 0; k < size; ++k)
00286 {
00287 Index rs = size-k-1;
00288
00289 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00290 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00291 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00292
00293 RealScalar x = real(mat.coeff(k,k));
00294 if (k>0) x -= A10.squaredNorm();
00295 if (x<=RealScalar(0))
00296 return k;
00297 mat.coeffRef(k,k) = x = sqrt(x);
00298 if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00299 if (rs>0) A21 *= RealScalar(1)/x;
00300 }
00301 return -1;
00302 }
00303
00304 template<typename MatrixType>
00305 static typename MatrixType::Index blocked(MatrixType& m)
00306 {
00307 typedef typename MatrixType::Index Index;
00308 eigen_assert(m.rows()==m.cols());
00309 Index size = m.rows();
00310 if(size<32)
00311 return unblocked(m);
00312
00313 Index blockSize = size/8;
00314 blockSize = (blockSize/16)*16;
00315 blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00316
00317 for (Index k=0; k<size; k+=blockSize)
00318 {
00319
00320
00321
00322
00323 Index bs = (std::min)(blockSize, size-k);
00324 Index rs = size - k - bs;
00325 Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
00326 Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
00327 Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00328
00329 Index ret;
00330 if((ret=unblocked(A11))>=0) return k+ret;
00331 if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00332 if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1);
00333 }
00334 return -1;
00335 }
00336
00337 template<typename MatrixType, typename VectorType>
00338 static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00339 {
00340 return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00341 }
00342 };
00343
00344 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00345 {
00346 typedef typename NumTraits<Scalar>::Real RealScalar;
00347
00348 template<typename MatrixType>
00349 static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00350 {
00351 Transpose<MatrixType> matt(mat);
00352 return llt_inplace<Scalar, Lower>::unblocked(matt);
00353 }
00354 template<typename MatrixType>
00355 static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00356 {
00357 Transpose<MatrixType> matt(mat);
00358 return llt_inplace<Scalar, Lower>::blocked(matt);
00359 }
00360 template<typename MatrixType, typename VectorType>
00361 static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00362 {
00363 Transpose<MatrixType> matt(mat);
00364 return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00365 }
00366 };
00367
00368 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00369 {
00370 typedef const TriangularView<const MatrixType, Lower> MatrixL;
00371 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00372 static inline MatrixL getL(const MatrixType& m) { return m; }
00373 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00374 static bool inplace_decomposition(MatrixType& m)
00375 { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00376 };
00377
00378 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00379 {
00380 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00381 typedef const TriangularView<const MatrixType, Upper> MatrixU;
00382 static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00383 static inline MatrixU getU(const MatrixType& m) { return m; }
00384 static bool inplace_decomposition(MatrixType& m)
00385 { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00386 };
00387
00388 }
00389
00397 template<typename MatrixType, int _UpLo>
00398 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00399 {
00400 eigen_assert(a.rows()==a.cols());
00401 const Index size = a.rows();
00402 m_matrix.resize(size, size);
00403 m_matrix = a;
00404
00405 m_isInitialized = true;
00406 bool ok = Traits::inplace_decomposition(m_matrix);
00407 m_info = ok ? Success : NumericalIssue;
00408
00409 return *this;
00410 }
00411
00417 template<typename _MatrixType, int _UpLo>
00418 template<typename VectorType>
00419 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00420 {
00421 EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00422 eigen_assert(v.size()==m_matrix.cols());
00423 eigen_assert(m_isInitialized);
00424 if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00425 m_info = NumericalIssue;
00426 else
00427 m_info = Success;
00428
00429 return *this;
00430 }
00431
00432 namespace internal {
00433 template<typename _MatrixType, int UpLo, typename Rhs>
00434 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00435 : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00436 {
00437 typedef LLT<_MatrixType,UpLo> LLTType;
00438 EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00439
00440 template<typename Dest> void evalTo(Dest& dst) const
00441 {
00442 dst = rhs();
00443 dec().solveInPlace(dst);
00444 }
00445 };
00446 }
00447
00461 template<typename MatrixType, int _UpLo>
00462 template<typename Derived>
00463 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00464 {
00465 eigen_assert(m_isInitialized && "LLT is not initialized.");
00466 eigen_assert(m_matrix.rows()==bAndX.rows());
00467 matrixL().solveInPlace(bAndX);
00468 matrixU().solveInPlace(bAndX);
00469 }
00470
00474 template<typename MatrixType, int _UpLo>
00475 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00476 {
00477 eigen_assert(m_isInitialized && "LLT is not initialized.");
00478 return matrixL() * matrixL().adjoint().toDenseMatrix();
00479 }
00480
00484 template<typename Derived>
00485 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00486 MatrixBase<Derived>::llt() const
00487 {
00488 return LLT<PlainObject>(derived());
00489 }
00490
00494 template<typename MatrixType, unsigned int UpLo>
00495 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00496 SelfAdjointView<MatrixType, UpLo>::llt() const
00497 {
00498 return LLT<PlainObject,UpLo>(m_matrix);
00499 }
00500
00501 }
00502
00503 #endif // EIGEN_LLT_H