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
00026
00027 #ifndef EIGEN_LDLT_H
00028 #define EIGEN_LDLT_H
00029
00030 namespace internal {
00031 template<typename MatrixType, int UpLo> struct LDLT_Traits;
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
00059 template<typename _MatrixType, int _UpLo> class LDLT
00060 {
00061 public:
00062 typedef _MatrixType MatrixType;
00063 enum {
00064 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00065 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00066 Options = MatrixType::Options & ~RowMajorBit,
00067 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00068 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00069 UpLo = _UpLo
00070 };
00071 typedef typename MatrixType::Scalar Scalar;
00072 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00073 typedef typename MatrixType::Index Index;
00074 typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
00075
00076 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00077 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00078
00079 typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
00080
00081
00082
00083
00084
00085
00086 LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
00087
00088
00089
00090
00091
00092
00093
00094 LDLT(Index size)
00095 : m_matrix(size, size),
00096 m_transpositions(size),
00097 m_temporary(size),
00098 m_isInitialized(false)
00099 {}
00100
00101 LDLT(const MatrixType& matrix)
00102 : m_matrix(matrix.rows(), matrix.cols()),
00103 m_transpositions(matrix.rows()),
00104 m_temporary(matrix.rows()),
00105 m_isInitialized(false)
00106 {
00107 compute(matrix);
00108 }
00109
00110
00111 inline typename Traits::MatrixU matrixU() const
00112 {
00113 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00114 return Traits::getU(m_matrix);
00115 }
00116
00117
00118 inline typename Traits::MatrixL matrixL() const
00119 {
00120 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00121 return Traits::getL(m_matrix);
00122 }
00123
00124
00125
00126 inline const TranspositionType& transpositionsP() const
00127 {
00128 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00129 return m_transpositions;
00130 }
00131
00132
00133 inline Diagonal<const MatrixType> vectorD(void) const
00134 {
00135 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00136 return m_matrix.diagonal();
00137 }
00138
00139
00140 inline bool isPositive(void) const
00141 {
00142 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00143 return m_sign == 1;
00144 }
00145
00146
00147 inline bool isNegative(void) const
00148 {
00149 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00150 return m_sign == -1;
00151 }
00152
00153
00154
00155
00156
00157
00158
00159 template<typename Rhs>
00160 inline const internal::solve_retval<LDLT, Rhs>
00161 solve(const MatrixBase<Rhs>& b) const
00162 {
00163 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00164 eigen_assert(m_matrix.rows()==b.rows()
00165 && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
00166 return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
00167 }
00168
00169 template<typename Derived>
00170 bool solveInPlace(MatrixBase<Derived> &bAndX) const;
00171
00172 LDLT& compute(const MatrixType& matrix);
00173
00174
00175
00176
00177
00178 inline const MatrixType& matrixLDLT() const
00179 {
00180 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00181 return m_matrix;
00182 }
00183
00184 MatrixType reconstructedMatrix() const;
00185
00186 inline Index rows() const { return m_matrix.rows(); }
00187 inline Index cols() const { return m_matrix.cols(); }
00188
00189 protected:
00190
00191
00192
00193
00194
00195
00196
00197 MatrixType m_matrix;
00198 TranspositionType m_transpositions;
00199 TmpMatrixType m_temporary;
00200 int m_sign;
00201 bool m_isInitialized;
00202 };
00203
00204 namespace internal {
00205
00206 template<int UpLo> struct ldlt_inplace;
00207
00208 template<> struct ldlt_inplace<Lower>
00209 {
00210 template<typename MatrixType, typename TranspositionType, typename Workspace>
00211 static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00212 {
00213 typedef typename MatrixType::Scalar Scalar;
00214 typedef typename MatrixType::RealScalar RealScalar;
00215 typedef typename MatrixType::Index Index;
00216 eigen_assert(mat.rows()==mat.cols());
00217 const Index size = mat.rows();
00218
00219 if (size <= 1)
00220 {
00221 transpositions.setIdentity();
00222 if(sign)
00223 *sign = real(mat.coeff(0,0))>0 ? 1:-1;
00224 return true;
00225 }
00226
00227 RealScalar cutoff = 0, biggest_in_corner;
00228
00229 for (Index k = 0; k < size; ++k)
00230 {
00231
00232 Index index_of_biggest_in_corner;
00233 biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
00234 index_of_biggest_in_corner += k;
00235
00236 if(k == 0)
00237 {
00238
00239
00240
00241 cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
00242
00243 if(sign)
00244 *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
00245 }
00246
00247
00248 if(biggest_in_corner < cutoff)
00249 {
00250 for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
00251 break;
00252 }
00253
00254 transpositions.coeffRef(k) = index_of_biggest_in_corner;
00255 if(k != index_of_biggest_in_corner)
00256 {
00257
00258
00259 Index s = size-index_of_biggest_in_corner-1;
00260 mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
00261 mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
00262 std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
00263 for(int i=k+1;i<index_of_biggest_in_corner;++i)
00264 {
00265 Scalar tmp = mat.coeffRef(i,k);
00266 mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
00267 mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
00268 }
00269 if(NumTraits<Scalar>::IsComplex)
00270 mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
00271 }
00272
00273
00274
00275
00276
00277 Index rs = size - k - 1;
00278 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00279 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00280 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00281
00282 if(k>0)
00283 {
00284 temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
00285 mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
00286 if(rs>0)
00287 A21.noalias() -= A20 * temp.head(k);
00288 }
00289 if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
00290 A21 /= mat.coeffRef(k,k);
00291 }
00292
00293 return true;
00294 }
00295 };
00296
00297 template<> struct ldlt_inplace<Upper>
00298 {
00299 template<typename MatrixType, typename TranspositionType, typename Workspace>
00300 static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00301 {
00302 Transpose<MatrixType> matt(mat);
00303 return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
00304 }
00305 };
00306
00307 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
00308 {
00309 typedef TriangularView<MatrixType, UnitLower> MatrixL;
00310 typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
00311 inline static MatrixL getL(const MatrixType& m) { return m; }
00312 inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00313 };
00314
00315 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
00316 {
00317 typedef TriangularView<typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
00318 typedef TriangularView<MatrixType, UnitUpper> MatrixU;
00319 inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00320 inline static MatrixU getU(const MatrixType& m) { return m; }
00321 };
00322
00323 }
00324
00325
00326
00327 template<typename MatrixType, int _UpLo>
00328 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00329 {
00330 eigen_assert(a.rows()==a.cols());
00331 const Index size = a.rows();
00332
00333 m_matrix = a;
00334
00335 m_transpositions.resize(size);
00336 m_isInitialized = false;
00337 m_temporary.resize(size);
00338
00339 internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
00340
00341 m_isInitialized = true;
00342 return *this;
00343 }
00344
00345 namespace internal {
00346 template<typename _MatrixType, int _UpLo, typename Rhs>
00347 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
00348 : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
00349 {
00350 typedef LDLT<_MatrixType,_UpLo> LDLTType;
00351 EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
00352
00353 template<typename Dest> void evalTo(Dest& dst) const
00354 {
00355 eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
00356
00357 dst = dec().transpositionsP() * rhs();
00358
00359
00360 dec().matrixL().solveInPlace(dst);
00361
00362
00363 dst = dec().vectorD().asDiagonal().inverse() * dst;
00364
00365
00366 dec().matrixU().solveInPlace(dst);
00367
00368
00369 dst = dec().transpositionsP().transpose() * dst;
00370 }
00371 };
00372 }
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 template<typename MatrixType,int _UpLo>
00388 template<typename Derived>
00389 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00390 {
00391 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00392 const Index size = m_matrix.rows();
00393 eigen_assert(size == bAndX.rows());
00394
00395 bAndX = this->solve(bAndX);
00396
00397 return true;
00398 }
00399
00400
00401
00402
00403 template<typename MatrixType, int _UpLo>
00404 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
00405 {
00406 eigen_assert(m_isInitialized && "LDLT is not initialized.");
00407 const Index size = m_matrix.rows();
00408 MatrixType res(size,size);
00409
00410
00411 res.setIdentity();
00412 res = transpositionsP() * res;
00413
00414 res = matrixU() * res;
00415
00416 res = vectorD().asDiagonal() * res;
00417
00418 res = matrixL() * res;
00419
00420 res = transpositionsP().transpose() * res;
00421
00422 return res;
00423 }
00424
00425
00426
00427
00428 template<typename MatrixType, unsigned int UpLo>
00429 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00430 SelfAdjointView<MatrixType, UpLo>::ldlt() const
00431 {
00432 return LDLT<PlainObject,UpLo>(m_matrix);
00433 }
00434
00435
00436
00437
00438 template<typename Derived>
00439 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
00440 MatrixBase<Derived>::ldlt() const
00441 {
00442 return LDLT<PlainObject>(derived());
00443 }
00444
00445 #endif // EIGEN_LDLT_H