00001
00002
00003
00004
00005
00006
00007
00008
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
00047
00048
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
00211
00212
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
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;
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
00307
00308
00309
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);
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 }
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 }
00489
00490 #endif // EIGEN_LLT_H