![]() |
Eigen
3.3.3
|
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 00048 /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) 00049 * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, 00050 * the strict lower part does not have to store correct values. 00051 */ 00052 template<typename _MatrixType, int _UpLo> class LLT 00053 { 00054 public: 00055 typedef _MatrixType MatrixType; 00056 enum { 00057 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00058 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00059 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00060 }; 00061 typedef typename MatrixType::Scalar Scalar; 00062 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00063 typedef Eigen::Index Index; 00064 typedef typename MatrixType::StorageIndex StorageIndex; 00065 00066 enum { 00067 PacketSize = internal::packet_traits<Scalar>::size, 00068 AlignmentMask = int(PacketSize)-1, 00069 UpLo = _UpLo 00070 }; 00071 00072 typedef internal::LLT_Traits<MatrixType,UpLo> Traits; 00073 00080 LLT() : m_matrix(), m_isInitialized(false) {} 00081 00088 explicit LLT(Index size) : m_matrix(size, size), 00089 m_isInitialized(false) {} 00090 00091 template<typename InputType> 00092 explicit LLT(const EigenBase<InputType>& matrix) 00093 : m_matrix(matrix.rows(), matrix.cols()), 00094 m_isInitialized(false) 00095 { 00096 compute(matrix.derived()); 00097 } 00098 00106 template<typename InputType> 00107 explicit LLT(EigenBase<InputType>& matrix) 00108 : m_matrix(matrix.derived()), 00109 m_isInitialized(false) 00110 { 00111 compute(matrix.derived()); 00112 } 00113 00115 inline typename Traits::MatrixU matrixU() const 00116 { 00117 eigen_assert(m_isInitialized && "LLT is not initialized."); 00118 return Traits::getU(m_matrix); 00119 } 00120 00122 inline typename Traits::MatrixL matrixL() const 00123 { 00124 eigen_assert(m_isInitialized && "LLT is not initialized."); 00125 return Traits::getL(m_matrix); 00126 } 00127 00138 template<typename Rhs> 00139 inline const Solve<LLT, Rhs> 00140 solve(const MatrixBase<Rhs>& b) const 00141 { 00142 eigen_assert(m_isInitialized && "LLT is not initialized."); 00143 eigen_assert(m_matrix.rows()==b.rows() 00144 && "LLT::solve(): invalid number of rows of the right hand side matrix b"); 00145 return Solve<LLT, Rhs>(*this, b.derived()); 00146 } 00147 00148 template<typename Derived> 00149 void solveInPlace(MatrixBase<Derived> &bAndX) const; 00150 00151 template<typename InputType> 00152 LLT& compute(const EigenBase<InputType>& matrix); 00153 00157 RealScalar rcond() const 00158 { 00159 eigen_assert(m_isInitialized && "LLT is not initialized."); 00160 eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative"); 00161 return internal::rcond_estimate_helper(m_l1_norm, *this); 00162 } 00163 00168 inline const MatrixType& matrixLLT() const 00169 { 00170 eigen_assert(m_isInitialized && "LLT is not initialized."); 00171 return m_matrix; 00172 } 00173 00174 MatrixType reconstructedMatrix() const; 00175 00176 00182 ComputationInfo info() const 00183 { 00184 eigen_assert(m_isInitialized && "LLT is not initialized."); 00185 return m_info; 00186 } 00187 00193 const LLT& adjoint() const { return *this; }; 00194 00195 inline Index rows() const { return m_matrix.rows(); } 00196 inline Index cols() const { return m_matrix.cols(); } 00197 00198 template<typename VectorType> 00199 LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); 00200 00201 #ifndef EIGEN_PARSED_BY_DOXYGEN 00202 template<typename RhsType, typename DstType> 00203 EIGEN_DEVICE_FUNC 00204 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00205 #endif 00206 00207 protected: 00208 00209 static void check_template_parameters() 00210 { 00211 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00212 } 00213 00218 MatrixType m_matrix; 00219 RealScalar m_l1_norm; 00220 bool m_isInitialized; 00221 ComputationInfo m_info; 00222 }; 00223 00224 namespace internal { 00225 00226 template<typename Scalar, int UpLo> struct llt_inplace; 00227 00228 template<typename MatrixType, typename VectorType> 00229 static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) 00230 { 00231 using std::sqrt; 00232 typedef typename MatrixType::Scalar Scalar; 00233 typedef typename MatrixType::RealScalar RealScalar; 00234 typedef typename MatrixType::ColXpr ColXpr; 00235 typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; 00236 typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; 00237 typedef Matrix<Scalar,Dynamic,1> TempVectorType; 00238 typedef typename TempVectorType::SegmentReturnType TempVecSegment; 00239 00240 Index n = mat.cols(); 00241 eigen_assert(mat.rows()==n && vec.size()==n); 00242 00243 TempVectorType temp; 00244 00245 if(sigma>0) 00246 { 00247 // This version is based on Givens rotations. 00248 // It is faster than the other one below, but only works for updates, 00249 // i.e., for sigma > 0 00250 temp = sqrt(sigma) * vec; 00251 00252 for(Index i=0; i<n; ++i) 00253 { 00254 JacobiRotation<Scalar> g; 00255 g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); 00256 00257 Index rs = n-i-1; 00258 if(rs>0) 00259 { 00260 ColXprSegment x(mat.col(i).tail(rs)); 00261 TempVecSegment y(temp.tail(rs)); 00262 apply_rotation_in_the_plane(x, y, g); 00263 } 00264 } 00265 } 00266 else 00267 { 00268 temp = vec; 00269 RealScalar beta = 1; 00270 for(Index j=0; j<n; ++j) 00271 { 00272 RealScalar Ljj = numext::real(mat.coeff(j,j)); 00273 RealScalar dj = numext::abs2(Ljj); 00274 Scalar wj = temp.coeff(j); 00275 RealScalar swj2 = sigma*numext::abs2(wj); 00276 RealScalar gamma = dj*beta + swj2; 00277 00278 RealScalar x = dj + swj2/beta; 00279 if (x<=RealScalar(0)) 00280 return j; 00281 RealScalar nLjj = sqrt(x); 00282 mat.coeffRef(j,j) = nLjj; 00283 beta += swj2/dj; 00284 00285 // Update the terms of L 00286 Index rs = n-j-1; 00287 if(rs) 00288 { 00289 temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs); 00290 if(gamma != 0) 00291 mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs); 00292 } 00293 } 00294 } 00295 return -1; 00296 } 00297 00298 template<typename Scalar> struct llt_inplace<Scalar, Lower> 00299 { 00300 typedef typename NumTraits<Scalar>::Real RealScalar; 00301 template<typename MatrixType> 00302 static Index unblocked(MatrixType& mat) 00303 { 00304 using std::sqrt; 00305 00306 eigen_assert(mat.rows()==mat.cols()); 00307 const Index size = mat.rows(); 00308 for(Index k = 0; k < size; ++k) 00309 { 00310 Index rs = size-k-1; // remaining size 00311 00312 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 00313 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 00314 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 00315 00316 RealScalar x = numext::real(mat.coeff(k,k)); 00317 if (k>0) x -= A10.squaredNorm(); 00318 if (x<=RealScalar(0)) 00319 return k; 00320 mat.coeffRef(k,k) = x = sqrt(x); 00321 if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); 00322 if (rs>0) A21 /= x; 00323 } 00324 return -1; 00325 } 00326 00327 template<typename MatrixType> 00328 static Index blocked(MatrixType& m) 00329 { 00330 eigen_assert(m.rows()==m.cols()); 00331 Index size = m.rows(); 00332 if(size<32) 00333 return unblocked(m); 00334 00335 Index blockSize = size/8; 00336 blockSize = (blockSize/16)*16; 00337 blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128)); 00338 00339 for (Index k=0; k<size; k+=blockSize) 00340 { 00341 // partition the matrix: 00342 // A00 | - | - 00343 // lu = A10 | A11 | - 00344 // A20 | A21 | A22 00345 Index bs = (std::min)(blockSize, size-k); 00346 Index rs = size - k - bs; 00347 Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs); 00348 Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs); 00349 Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs); 00350 00351 Index ret; 00352 if((ret=unblocked(A11))>=0) return k+ret; 00353 if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); 00354 if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck 00355 } 00356 return -1; 00357 } 00358 00359 template<typename MatrixType, typename VectorType> 00360 static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 00361 { 00362 return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); 00363 } 00364 }; 00365 00366 template<typename Scalar> struct llt_inplace<Scalar, Upper> 00367 { 00368 typedef typename NumTraits<Scalar>::Real RealScalar; 00369 00370 template<typename MatrixType> 00371 static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat) 00372 { 00373 Transpose<MatrixType> matt(mat); 00374 return llt_inplace<Scalar, Lower>::unblocked(matt); 00375 } 00376 template<typename MatrixType> 00377 static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat) 00378 { 00379 Transpose<MatrixType> matt(mat); 00380 return llt_inplace<Scalar, Lower>::blocked(matt); 00381 } 00382 template<typename MatrixType, typename VectorType> 00383 static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 00384 { 00385 Transpose<MatrixType> matt(mat); 00386 return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma); 00387 } 00388 }; 00389 00390 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> 00391 { 00392 typedef const TriangularView<const MatrixType, Lower> MatrixL; 00393 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU; 00394 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 00395 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 00396 static bool inplace_decomposition(MatrixType& m) 00397 { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; } 00398 }; 00399 00400 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> 00401 { 00402 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL; 00403 typedef const TriangularView<const MatrixType, Upper> MatrixU; 00404 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); } 00405 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); } 00406 static bool inplace_decomposition(MatrixType& m) 00407 { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; } 00408 }; 00409 00410 } // end namespace internal 00411 00419 template<typename MatrixType, int _UpLo> 00420 template<typename InputType> 00421 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a) 00422 { 00423 check_template_parameters(); 00424 00425 eigen_assert(a.rows()==a.cols()); 00426 const Index size = a.rows(); 00427 m_matrix.resize(size, size); 00428 m_matrix = a.derived(); 00429 00430 // Compute matrix L1 norm = max abs column sum. 00431 m_l1_norm = RealScalar(0); 00432 // TODO move this code to SelfAdjointView 00433 for (Index col = 0; col < size; ++col) { 00434 RealScalar abs_col_sum; 00435 if (_UpLo == Lower) 00436 abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>(); 00437 else 00438 abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>(); 00439 if (abs_col_sum > m_l1_norm) 00440 m_l1_norm = abs_col_sum; 00441 } 00442 00443 m_isInitialized = true; 00444 bool ok = Traits::inplace_decomposition(m_matrix); 00445 m_info = ok ? Success : NumericalIssue; 00446 00447 return *this; 00448 } 00449 00455 template<typename _MatrixType, int _UpLo> 00456 template<typename VectorType> 00457 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) 00458 { 00459 EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); 00460 eigen_assert(v.size()==m_matrix.cols()); 00461 eigen_assert(m_isInitialized); 00462 if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0) 00463 m_info = NumericalIssue; 00464 else 00465 m_info = Success; 00466 00467 return *this; 00468 } 00469 00470 #ifndef EIGEN_PARSED_BY_DOXYGEN 00471 template<typename _MatrixType,int _UpLo> 00472 template<typename RhsType, typename DstType> 00473 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const 00474 { 00475 dst = rhs; 00476 solveInPlace(dst); 00477 } 00478 #endif 00479 00490 template<typename MatrixType, int _UpLo> 00491 template<typename Derived> 00492 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 00493 { 00494 eigen_assert(m_isInitialized && "LLT is not initialized."); 00495 eigen_assert(m_matrix.rows()==bAndX.rows()); 00496 matrixL().solveInPlace(bAndX); 00497 matrixU().solveInPlace(bAndX); 00498 } 00499 00503 template<typename MatrixType, int _UpLo> 00504 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const 00505 { 00506 eigen_assert(m_isInitialized && "LLT is not initialized."); 00507 return matrixL() * matrixL().adjoint().toDenseMatrix(); 00508 } 00509 00514 template<typename Derived> 00515 inline const LLT<typename MatrixBase<Derived>::PlainObject> 00516 MatrixBase<Derived>::llt() const 00517 { 00518 return LLT<PlainObject>(derived()); 00519 } 00520 00525 template<typename MatrixType, unsigned int UpLo> 00526 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 00527 SelfAdjointView<MatrixType, UpLo>::llt() const 00528 { 00529 return LLT<PlainObject,UpLo>(m_matrix); 00530 } 00531 00532 } // end namespace Eigen 00533 00534 #endif // EIGEN_LLT_H