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 internal{
00029 template<typename MatrixType, int UpLo> struct LLT_Traits;
00030 }
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
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
00080
00081
00082
00083
00084
00085
00086 LLT() : m_matrix(), m_isInitialized(false) {}
00087
00088
00089
00090
00091
00092
00093
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
00104
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
00111
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
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
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 template<typename Derived>
00139 void solveInPlace(MatrixBase<Derived> &bAndX) const;
00140
00141 LLT& compute(const MatrixType& matrix);
00142
00143
00144
00145
00146
00147 inline const MatrixType& matrixLLT() const
00148 {
00149 eigen_assert(m_isInitialized && "LLT is not initialized.");
00150 return m_matrix;
00151 }
00152
00153 MatrixType reconstructedMatrix() const;
00154
00155
00156
00157
00158
00159
00160
00161 ComputationInfo info() const
00162 {
00163 eigen_assert(m_isInitialized && "LLT is not initialized.");
00164 return m_info;
00165 }
00166
00167 inline Index rows() const { return m_matrix.rows(); }
00168 inline Index cols() const { return m_matrix.cols(); }
00169
00170 protected:
00171
00172
00173
00174
00175 MatrixType m_matrix;
00176 bool m_isInitialized;
00177 ComputationInfo m_info;
00178 };
00179
00180 namespace internal {
00181
00182 template<int UpLo> struct llt_inplace;
00183
00184 template<> struct llt_inplace<Lower>
00185 {
00186 template<typename MatrixType>
00187 static bool unblocked(MatrixType& mat)
00188 {
00189 typedef typename MatrixType::Scalar Scalar;
00190 typedef typename MatrixType::RealScalar RealScalar;
00191 typedef typename MatrixType::Index Index;
00192 eigen_assert(mat.rows()==mat.cols());
00193 const Index size = mat.rows();
00194 for(Index k = 0; k < size; ++k)
00195 {
00196 Index rs = size-k-1;
00197
00198 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00199 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00200 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00201
00202 RealScalar x = real(mat.coeff(k,k));
00203 if (k>0) x -= A10.squaredNorm();
00204 if (x<=RealScalar(0))
00205 return false;
00206 mat.coeffRef(k,k) = x = sqrt(x);
00207 if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00208 if (rs>0) A21 *= RealScalar(1)/x;
00209 }
00210 return true;
00211 }
00212
00213 template<typename MatrixType>
00214 static bool blocked(MatrixType& m)
00215 {
00216 typedef typename MatrixType::Index Index;
00217 eigen_assert(m.rows()==m.cols());
00218 Index size = m.rows();
00219 if(size<32)
00220 return unblocked(m);
00221
00222 Index blockSize = size/8;
00223 blockSize = (blockSize/16)*16;
00224 blockSize = std::min(std::max(blockSize,Index(8)), Index(128));
00225
00226 for (Index k=0; k<size; k+=blockSize)
00227 {
00228
00229
00230
00231
00232 Index bs = std::min(blockSize, size-k);
00233 Index rs = size - k - bs;
00234 Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
00235 Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
00236 Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00237
00238 if(!unblocked(A11)) return false;
00239 if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00240 if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1);
00241 }
00242 return true;
00243 }
00244 };
00245
00246 template<> struct llt_inplace<Upper>
00247 {
00248 template<typename MatrixType>
00249 static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
00250 {
00251 Transpose<MatrixType> matt(mat);
00252 return llt_inplace<Lower>::unblocked(matt);
00253 }
00254 template<typename MatrixType>
00255 static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
00256 {
00257 Transpose<MatrixType> matt(mat);
00258 return llt_inplace<Lower>::blocked(matt);
00259 }
00260 };
00261
00262 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00263 {
00264 typedef TriangularView<MatrixType, Lower> MatrixL;
00265 typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU;
00266 inline static MatrixL getL(const MatrixType& m) { return m; }
00267 inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00268 static bool inplace_decomposition(MatrixType& m)
00269 { return llt_inplace<Lower>::blocked(m); }
00270 };
00271
00272 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00273 {
00274 typedef TriangularView<typename MatrixType::AdjointReturnType, Lower> MatrixL;
00275 typedef TriangularView<MatrixType, Upper> MatrixU;
00276 inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00277 inline static MatrixU getU(const MatrixType& m) { return m; }
00278 static bool inplace_decomposition(MatrixType& m)
00279 { return llt_inplace<Upper>::blocked(m); }
00280 };
00281
00282 }
00283
00284
00285
00286
00287
00288
00289 template<typename MatrixType, int _UpLo>
00290 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00291 {
00292 assert(a.rows()==a.cols());
00293 const Index size = a.rows();
00294 m_matrix.resize(size, size);
00295 m_matrix = a;
00296
00297 m_isInitialized = true;
00298 bool ok = Traits::inplace_decomposition(m_matrix);
00299 m_info = ok ? Success : NumericalIssue;
00300
00301 return *this;
00302 }
00303
00304 namespace internal {
00305 template<typename _MatrixType, int UpLo, typename Rhs>
00306 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00307 : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00308 {
00309 typedef LLT<_MatrixType,UpLo> LLTType;
00310 EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00311
00312 template<typename Dest> void evalTo(Dest& dst) const
00313 {
00314 dst = rhs();
00315 dec().solveInPlace(dst);
00316 }
00317 };
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 template<typename MatrixType, int _UpLo>
00334 template<typename Derived>
00335 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00336 {
00337 eigen_assert(m_isInitialized && "LLT is not initialized.");
00338 eigen_assert(m_matrix.rows()==bAndX.rows());
00339 matrixL().solveInPlace(bAndX);
00340 matrixU().solveInPlace(bAndX);
00341 }
00342
00343
00344
00345
00346 template<typename MatrixType, int _UpLo>
00347 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00348 {
00349 eigen_assert(m_isInitialized && "LLT is not initialized.");
00350 return matrixL() * matrixL().adjoint().toDenseMatrix();
00351 }
00352
00353
00354
00355
00356 template<typename Derived>
00357 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00358 MatrixBase<Derived>::llt() const
00359 {
00360 return LLT<PlainObject>(derived());
00361 }
00362
00363
00364
00365
00366 template<typename MatrixType, unsigned int UpLo>
00367 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00368 SelfAdjointView<MatrixType, UpLo>::llt() const
00369 {
00370 return LLT<PlainObject,UpLo>(m_matrix);
00371 }
00372
00373 #endif // EIGEN_LLT_H