Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 4 | // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 5 | // |
| 6 | // This Source Code Form is subject to the terms of the Mozilla |
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 9 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 10 | #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H |
| 11 | #define EIGEN_SIMPLICIAL_CHOLESKY_H |
| 12 | |
| 13 | namespace Eigen { |
| 14 | |
| 15 | enum SimplicialCholeskyMode { |
| 16 | SimplicialCholeskyLLT, |
| 17 | SimplicialCholeskyLDLT |
| 18 | }; |
| 19 | |
| 20 | /** \ingroup SparseCholesky_Module |
| 21 | * \brief A direct sparse Cholesky factorizations |
| 22 | * |
| 23 | * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are |
| 24 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where |
| 25 | * X and B can be either dense or sparse. |
| 26 | * |
| 27 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| 28 | * such that the factorized matrix is P A P^-1. |
| 29 | * |
| 30 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| 31 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower |
| 32 | * or Upper. Default is Lower. |
| 33 | * |
| 34 | */ |
| 35 | template<typename Derived> |
| 36 | class SimplicialCholeskyBase : internal::noncopyable |
| 37 | { |
| 38 | public: |
| 39 | typedef typename internal::traits<Derived>::MatrixType MatrixType; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 40 | typedef typename internal::traits<Derived>::OrderingType OrderingType; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 41 | enum { UpLo = internal::traits<Derived>::UpLo }; |
| 42 | typedef typename MatrixType::Scalar Scalar; |
| 43 | typedef typename MatrixType::RealScalar RealScalar; |
| 44 | typedef typename MatrixType::Index Index; |
| 45 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; |
| 46 | typedef Matrix<Scalar,Dynamic,1> VectorType; |
| 47 | |
| 48 | public: |
| 49 | |
| 50 | /** Default constructor */ |
| 51 | SimplicialCholeskyBase() |
| 52 | : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) |
| 53 | {} |
| 54 | |
| 55 | SimplicialCholeskyBase(const MatrixType& matrix) |
| 56 | : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) |
| 57 | { |
| 58 | derived().compute(matrix); |
| 59 | } |
| 60 | |
| 61 | ~SimplicialCholeskyBase() |
| 62 | { |
| 63 | } |
| 64 | |
| 65 | Derived& derived() { return *static_cast<Derived*>(this); } |
| 66 | const Derived& derived() const { return *static_cast<const Derived*>(this); } |
| 67 | |
| 68 | inline Index cols() const { return m_matrix.cols(); } |
| 69 | inline Index rows() const { return m_matrix.rows(); } |
| 70 | |
| 71 | /** \brief Reports whether previous computation was successful. |
| 72 | * |
| 73 | * \returns \c Success if computation was succesful, |
| 74 | * \c NumericalIssue if the matrix.appears to be negative. |
| 75 | */ |
| 76 | ComputationInfo info() const |
| 77 | { |
| 78 | eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| 79 | return m_info; |
| 80 | } |
| 81 | |
| 82 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. |
| 83 | * |
| 84 | * \sa compute() |
| 85 | */ |
| 86 | template<typename Rhs> |
| 87 | inline const internal::solve_retval<SimplicialCholeskyBase, Rhs> |
| 88 | solve(const MatrixBase<Rhs>& b) const |
| 89 | { |
| 90 | eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); |
| 91 | eigen_assert(rows()==b.rows() |
| 92 | && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b"); |
| 93 | return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); |
| 94 | } |
| 95 | |
| 96 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. |
| 97 | * |
| 98 | * \sa compute() |
| 99 | */ |
| 100 | template<typename Rhs> |
| 101 | inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs> |
| 102 | solve(const SparseMatrixBase<Rhs>& b) const |
| 103 | { |
| 104 | eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); |
| 105 | eigen_assert(rows()==b.rows() |
| 106 | && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); |
| 107 | return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); |
| 108 | } |
| 109 | |
| 110 | /** \returns the permutation P |
| 111 | * \sa permutationPinv() */ |
| 112 | const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const |
| 113 | { return m_P; } |
| 114 | |
| 115 | /** \returns the inverse P^-1 of the permutation P |
| 116 | * \sa permutationP() */ |
| 117 | const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const |
| 118 | { return m_Pinv; } |
| 119 | |
| 120 | /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization. |
| 121 | * |
| 122 | * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n |
| 123 | * \c d_ii = \a offset + \a scale * \c d_ii |
| 124 | * |
| 125 | * The default is the identity transformation with \a offset=0, and \a scale=1. |
| 126 | * |
| 127 | * \returns a reference to \c *this. |
| 128 | */ |
| 129 | Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) |
| 130 | { |
| 131 | m_shiftOffset = offset; |
| 132 | m_shiftScale = scale; |
| 133 | return derived(); |
| 134 | } |
| 135 | |
| 136 | #ifndef EIGEN_PARSED_BY_DOXYGEN |
| 137 | /** \internal */ |
| 138 | template<typename Stream> |
| 139 | void dumpMemory(Stream& s) |
| 140 | { |
| 141 | int total = 0; |
| 142 | s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; |
| 143 | s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; |
| 144 | s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| 145 | s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| 146 | s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| 147 | s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| 148 | s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; |
| 149 | } |
| 150 | |
| 151 | /** \internal */ |
| 152 | template<typename Rhs,typename Dest> |
| 153 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const |
| 154 | { |
| 155 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); |
| 156 | eigen_assert(m_matrix.rows()==b.rows()); |
| 157 | |
| 158 | if(m_info!=Success) |
| 159 | return; |
| 160 | |
| 161 | if(m_P.size()>0) |
| 162 | dest = m_P * b; |
| 163 | else |
| 164 | dest = b; |
| 165 | |
| 166 | if(m_matrix.nonZeros()>0) // otherwise L==I |
| 167 | derived().matrixL().solveInPlace(dest); |
| 168 | |
| 169 | if(m_diag.size()>0) |
| 170 | dest = m_diag.asDiagonal().inverse() * dest; |
| 171 | |
| 172 | if (m_matrix.nonZeros()>0) // otherwise U==I |
| 173 | derived().matrixU().solveInPlace(dest); |
| 174 | |
| 175 | if(m_P.size()>0) |
| 176 | dest = m_Pinv * dest; |
| 177 | } |
| 178 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 179 | #endif // EIGEN_PARSED_BY_DOXYGEN |
| 180 | |
| 181 | protected: |
| 182 | |
| 183 | /** Computes the sparse Cholesky decomposition of \a matrix */ |
| 184 | template<bool DoLDLT> |
| 185 | void compute(const MatrixType& matrix) |
| 186 | { |
| 187 | eigen_assert(matrix.rows()==matrix.cols()); |
| 188 | Index size = matrix.cols(); |
| 189 | CholMatrixType ap(size,size); |
| 190 | ordering(matrix, ap); |
| 191 | analyzePattern_preordered(ap, DoLDLT); |
| 192 | factorize_preordered<DoLDLT>(ap); |
| 193 | } |
| 194 | |
| 195 | template<bool DoLDLT> |
| 196 | void factorize(const MatrixType& a) |
| 197 | { |
| 198 | eigen_assert(a.rows()==a.cols()); |
| 199 | int size = a.cols(); |
| 200 | CholMatrixType ap(size,size); |
| 201 | ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); |
| 202 | factorize_preordered<DoLDLT>(ap); |
| 203 | } |
| 204 | |
| 205 | template<bool DoLDLT> |
| 206 | void factorize_preordered(const CholMatrixType& a); |
| 207 | |
| 208 | void analyzePattern(const MatrixType& a, bool doLDLT) |
| 209 | { |
| 210 | eigen_assert(a.rows()==a.cols()); |
| 211 | int size = a.cols(); |
| 212 | CholMatrixType ap(size,size); |
| 213 | ordering(a, ap); |
| 214 | analyzePattern_preordered(ap,doLDLT); |
| 215 | } |
| 216 | void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); |
| 217 | |
| 218 | void ordering(const MatrixType& a, CholMatrixType& ap); |
| 219 | |
| 220 | /** keeps off-diagonal entries; drops diagonal entries */ |
| 221 | struct keep_diag { |
| 222 | inline bool operator() (const Index& row, const Index& col, const Scalar&) const |
| 223 | { |
| 224 | return row!=col; |
| 225 | } |
| 226 | }; |
| 227 | |
| 228 | mutable ComputationInfo m_info; |
| 229 | bool m_isInitialized; |
| 230 | bool m_factorizationIsOk; |
| 231 | bool m_analysisIsOk; |
| 232 | |
| 233 | CholMatrixType m_matrix; |
| 234 | VectorType m_diag; // the diagonal coefficients (LDLT mode) |
| 235 | VectorXi m_parent; // elimination tree |
| 236 | VectorXi m_nonZerosPerCol; |
| 237 | PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation |
| 238 | PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation |
| 239 | |
| 240 | RealScalar m_shiftOffset; |
| 241 | RealScalar m_shiftScale; |
| 242 | }; |
| 243 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 244 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT; |
| 245 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT; |
| 246 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 247 | |
| 248 | namespace internal { |
| 249 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 250 | template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 251 | { |
| 252 | typedef _MatrixType MatrixType; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 253 | typedef _Ordering OrderingType; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 254 | enum { UpLo = _UpLo }; |
| 255 | typedef typename MatrixType::Scalar Scalar; |
| 256 | typedef typename MatrixType::Index Index; |
| 257 | typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; |
| 258 | typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL; |
| 259 | typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; |
| 260 | static inline MatrixL getL(const MatrixType& m) { return m; } |
| 261 | static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } |
| 262 | }; |
| 263 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 264 | template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 265 | { |
| 266 | typedef _MatrixType MatrixType; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 267 | typedef _Ordering OrderingType; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 268 | enum { UpLo = _UpLo }; |
| 269 | typedef typename MatrixType::Scalar Scalar; |
| 270 | typedef typename MatrixType::Index Index; |
| 271 | typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; |
| 272 | typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL; |
| 273 | typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; |
| 274 | static inline MatrixL getL(const MatrixType& m) { return m; } |
| 275 | static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } |
| 276 | }; |
| 277 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 278 | template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 279 | { |
| 280 | typedef _MatrixType MatrixType; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 281 | typedef _Ordering OrderingType; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 282 | enum { UpLo = _UpLo }; |
| 283 | }; |
| 284 | |
| 285 | } |
| 286 | |
| 287 | /** \ingroup SparseCholesky_Module |
| 288 | * \class SimplicialLLT |
| 289 | * \brief A direct sparse LLT Cholesky factorizations |
| 290 | * |
| 291 | * This class provides a LL^T Cholesky factorizations of sparse matrices that are |
| 292 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where |
| 293 | * X and B can be either dense or sparse. |
| 294 | * |
| 295 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| 296 | * such that the factorized matrix is P A P^-1. |
| 297 | * |
| 298 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| 299 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower |
| 300 | * or Upper. Default is Lower. |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 301 | * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 302 | * |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 303 | * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 304 | */ |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 305 | template<typename _MatrixType, int _UpLo, typename _Ordering> |
| 306 | class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 307 | { |
| 308 | public: |
| 309 | typedef _MatrixType MatrixType; |
| 310 | enum { UpLo = _UpLo }; |
| 311 | typedef SimplicialCholeskyBase<SimplicialLLT> Base; |
| 312 | typedef typename MatrixType::Scalar Scalar; |
| 313 | typedef typename MatrixType::RealScalar RealScalar; |
| 314 | typedef typename MatrixType::Index Index; |
| 315 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; |
| 316 | typedef Matrix<Scalar,Dynamic,1> VectorType; |
| 317 | typedef internal::traits<SimplicialLLT> Traits; |
| 318 | typedef typename Traits::MatrixL MatrixL; |
| 319 | typedef typename Traits::MatrixU MatrixU; |
| 320 | public: |
| 321 | /** Default constructor */ |
| 322 | SimplicialLLT() : Base() {} |
| 323 | /** Constructs and performs the LLT factorization of \a matrix */ |
| 324 | SimplicialLLT(const MatrixType& matrix) |
| 325 | : Base(matrix) {} |
| 326 | |
| 327 | /** \returns an expression of the factor L */ |
| 328 | inline const MatrixL matrixL() const { |
| 329 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); |
| 330 | return Traits::getL(Base::m_matrix); |
| 331 | } |
| 332 | |
| 333 | /** \returns an expression of the factor U (= L^*) */ |
| 334 | inline const MatrixU matrixU() const { |
| 335 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); |
| 336 | return Traits::getU(Base::m_matrix); |
| 337 | } |
| 338 | |
| 339 | /** Computes the sparse Cholesky decomposition of \a matrix */ |
| 340 | SimplicialLLT& compute(const MatrixType& matrix) |
| 341 | { |
| 342 | Base::template compute<false>(matrix); |
| 343 | return *this; |
| 344 | } |
| 345 | |
| 346 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 347 | * |
| 348 | * This function is particularly useful when solving for several problems having the same structure. |
| 349 | * |
| 350 | * \sa factorize() |
| 351 | */ |
| 352 | void analyzePattern(const MatrixType& a) |
| 353 | { |
| 354 | Base::analyzePattern(a, false); |
| 355 | } |
| 356 | |
| 357 | /** Performs a numeric decomposition of \a matrix |
| 358 | * |
| 359 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| 360 | * |
| 361 | * \sa analyzePattern() |
| 362 | */ |
| 363 | void factorize(const MatrixType& a) |
| 364 | { |
| 365 | Base::template factorize<false>(a); |
| 366 | } |
| 367 | |
| 368 | /** \returns the determinant of the underlying matrix from the current factorization */ |
| 369 | Scalar determinant() const |
| 370 | { |
| 371 | Scalar detL = Base::m_matrix.diagonal().prod(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 372 | return numext::abs2(detL); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 373 | } |
| 374 | }; |
| 375 | |
| 376 | /** \ingroup SparseCholesky_Module |
| 377 | * \class SimplicialLDLT |
| 378 | * \brief A direct sparse LDLT Cholesky factorizations without square root. |
| 379 | * |
| 380 | * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are |
| 381 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where |
| 382 | * X and B can be either dense or sparse. |
| 383 | * |
| 384 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| 385 | * such that the factorized matrix is P A P^-1. |
| 386 | * |
| 387 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| 388 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower |
| 389 | * or Upper. Default is Lower. |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 390 | * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 391 | * |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 392 | * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 393 | */ |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 394 | template<typename _MatrixType, int _UpLo, typename _Ordering> |
| 395 | class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 396 | { |
| 397 | public: |
| 398 | typedef _MatrixType MatrixType; |
| 399 | enum { UpLo = _UpLo }; |
| 400 | typedef SimplicialCholeskyBase<SimplicialLDLT> Base; |
| 401 | typedef typename MatrixType::Scalar Scalar; |
| 402 | typedef typename MatrixType::RealScalar RealScalar; |
| 403 | typedef typename MatrixType::Index Index; |
| 404 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; |
| 405 | typedef Matrix<Scalar,Dynamic,1> VectorType; |
| 406 | typedef internal::traits<SimplicialLDLT> Traits; |
| 407 | typedef typename Traits::MatrixL MatrixL; |
| 408 | typedef typename Traits::MatrixU MatrixU; |
| 409 | public: |
| 410 | /** Default constructor */ |
| 411 | SimplicialLDLT() : Base() {} |
| 412 | |
| 413 | /** Constructs and performs the LLT factorization of \a matrix */ |
| 414 | SimplicialLDLT(const MatrixType& matrix) |
| 415 | : Base(matrix) {} |
| 416 | |
| 417 | /** \returns a vector expression of the diagonal D */ |
| 418 | inline const VectorType vectorD() const { |
| 419 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| 420 | return Base::m_diag; |
| 421 | } |
| 422 | /** \returns an expression of the factor L */ |
| 423 | inline const MatrixL matrixL() const { |
| 424 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| 425 | return Traits::getL(Base::m_matrix); |
| 426 | } |
| 427 | |
| 428 | /** \returns an expression of the factor U (= L^*) */ |
| 429 | inline const MatrixU matrixU() const { |
| 430 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| 431 | return Traits::getU(Base::m_matrix); |
| 432 | } |
| 433 | |
| 434 | /** Computes the sparse Cholesky decomposition of \a matrix */ |
| 435 | SimplicialLDLT& compute(const MatrixType& matrix) |
| 436 | { |
| 437 | Base::template compute<true>(matrix); |
| 438 | return *this; |
| 439 | } |
| 440 | |
| 441 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 442 | * |
| 443 | * This function is particularly useful when solving for several problems having the same structure. |
| 444 | * |
| 445 | * \sa factorize() |
| 446 | */ |
| 447 | void analyzePattern(const MatrixType& a) |
| 448 | { |
| 449 | Base::analyzePattern(a, true); |
| 450 | } |
| 451 | |
| 452 | /** Performs a numeric decomposition of \a matrix |
| 453 | * |
| 454 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| 455 | * |
| 456 | * \sa analyzePattern() |
| 457 | */ |
| 458 | void factorize(const MatrixType& a) |
| 459 | { |
| 460 | Base::template factorize<true>(a); |
| 461 | } |
| 462 | |
| 463 | /** \returns the determinant of the underlying matrix from the current factorization */ |
| 464 | Scalar determinant() const |
| 465 | { |
| 466 | return Base::m_diag.prod(); |
| 467 | } |
| 468 | }; |
| 469 | |
| 470 | /** \deprecated use SimplicialLDLT or class SimplicialLLT |
| 471 | * \ingroup SparseCholesky_Module |
| 472 | * \class SimplicialCholesky |
| 473 | * |
| 474 | * \sa class SimplicialLDLT, class SimplicialLLT |
| 475 | */ |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 476 | template<typename _MatrixType, int _UpLo, typename _Ordering> |
| 477 | class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 478 | { |
| 479 | public: |
| 480 | typedef _MatrixType MatrixType; |
| 481 | enum { UpLo = _UpLo }; |
| 482 | typedef SimplicialCholeskyBase<SimplicialCholesky> Base; |
| 483 | typedef typename MatrixType::Scalar Scalar; |
| 484 | typedef typename MatrixType::RealScalar RealScalar; |
| 485 | typedef typename MatrixType::Index Index; |
| 486 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; |
| 487 | typedef Matrix<Scalar,Dynamic,1> VectorType; |
| 488 | typedef internal::traits<SimplicialCholesky> Traits; |
| 489 | typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; |
| 490 | typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; |
| 491 | public: |
| 492 | SimplicialCholesky() : Base(), m_LDLT(true) {} |
| 493 | |
| 494 | SimplicialCholesky(const MatrixType& matrix) |
| 495 | : Base(), m_LDLT(true) |
| 496 | { |
| 497 | compute(matrix); |
| 498 | } |
| 499 | |
| 500 | SimplicialCholesky& setMode(SimplicialCholeskyMode mode) |
| 501 | { |
| 502 | switch(mode) |
| 503 | { |
| 504 | case SimplicialCholeskyLLT: |
| 505 | m_LDLT = false; |
| 506 | break; |
| 507 | case SimplicialCholeskyLDLT: |
| 508 | m_LDLT = true; |
| 509 | break; |
| 510 | default: |
| 511 | break; |
| 512 | } |
| 513 | |
| 514 | return *this; |
| 515 | } |
| 516 | |
| 517 | inline const VectorType vectorD() const { |
| 518 | eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); |
| 519 | return Base::m_diag; |
| 520 | } |
| 521 | inline const CholMatrixType rawMatrix() const { |
| 522 | eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); |
| 523 | return Base::m_matrix; |
| 524 | } |
| 525 | |
| 526 | /** Computes the sparse Cholesky decomposition of \a matrix */ |
| 527 | SimplicialCholesky& compute(const MatrixType& matrix) |
| 528 | { |
| 529 | if(m_LDLT) |
| 530 | Base::template compute<true>(matrix); |
| 531 | else |
| 532 | Base::template compute<false>(matrix); |
| 533 | return *this; |
| 534 | } |
| 535 | |
| 536 | /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| 537 | * |
| 538 | * This function is particularly useful when solving for several problems having the same structure. |
| 539 | * |
| 540 | * \sa factorize() |
| 541 | */ |
| 542 | void analyzePattern(const MatrixType& a) |
| 543 | { |
| 544 | Base::analyzePattern(a, m_LDLT); |
| 545 | } |
| 546 | |
| 547 | /** Performs a numeric decomposition of \a matrix |
| 548 | * |
| 549 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| 550 | * |
| 551 | * \sa analyzePattern() |
| 552 | */ |
| 553 | void factorize(const MatrixType& a) |
| 554 | { |
| 555 | if(m_LDLT) |
| 556 | Base::template factorize<true>(a); |
| 557 | else |
| 558 | Base::template factorize<false>(a); |
| 559 | } |
| 560 | |
| 561 | /** \internal */ |
| 562 | template<typename Rhs,typename Dest> |
| 563 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const |
| 564 | { |
| 565 | eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); |
| 566 | eigen_assert(Base::m_matrix.rows()==b.rows()); |
| 567 | |
| 568 | if(Base::m_info!=Success) |
| 569 | return; |
| 570 | |
| 571 | if(Base::m_P.size()>0) |
| 572 | dest = Base::m_P * b; |
| 573 | else |
| 574 | dest = b; |
| 575 | |
| 576 | if(Base::m_matrix.nonZeros()>0) // otherwise L==I |
| 577 | { |
| 578 | if(m_LDLT) |
| 579 | LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); |
| 580 | else |
| 581 | LLTTraits::getL(Base::m_matrix).solveInPlace(dest); |
| 582 | } |
| 583 | |
| 584 | if(Base::m_diag.size()>0) |
| 585 | dest = Base::m_diag.asDiagonal().inverse() * dest; |
| 586 | |
| 587 | if (Base::m_matrix.nonZeros()>0) // otherwise I==I |
| 588 | { |
| 589 | if(m_LDLT) |
| 590 | LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); |
| 591 | else |
| 592 | LLTTraits::getU(Base::m_matrix).solveInPlace(dest); |
| 593 | } |
| 594 | |
| 595 | if(Base::m_P.size()>0) |
| 596 | dest = Base::m_Pinv * dest; |
| 597 | } |
| 598 | |
| 599 | Scalar determinant() const |
| 600 | { |
| 601 | if(m_LDLT) |
| 602 | { |
| 603 | return Base::m_diag.prod(); |
| 604 | } |
| 605 | else |
| 606 | { |
| 607 | Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 608 | return numext::abs2(detL); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 609 | } |
| 610 | } |
| 611 | |
| 612 | protected: |
| 613 | bool m_LDLT; |
| 614 | }; |
| 615 | |
| 616 | template<typename Derived> |
| 617 | void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap) |
| 618 | { |
| 619 | eigen_assert(a.rows()==a.cols()); |
| 620 | const Index size = a.rows(); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 621 | // Note that amd compute the inverse permutation |
| 622 | { |
| 623 | CholMatrixType C; |
| 624 | C = a.template selfadjointView<UpLo>(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 625 | |
| 626 | OrderingType ordering; |
| 627 | ordering(C,m_Pinv); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 628 | } |
| 629 | |
| 630 | if(m_Pinv.size()>0) |
| 631 | m_P = m_Pinv.inverse(); |
| 632 | else |
| 633 | m_P.resize(0); |
| 634 | |
| 635 | ap.resize(size,size); |
| 636 | ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); |
| 637 | } |
| 638 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 639 | namespace internal { |
| 640 | |
| 641 | template<typename Derived, typename Rhs> |
| 642 | struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs> |
| 643 | : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> |
| 644 | { |
| 645 | typedef SimplicialCholeskyBase<Derived> Dec; |
| 646 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) |
| 647 | |
| 648 | template<typename Dest> void evalTo(Dest& dst) const |
| 649 | { |
| 650 | dec().derived()._solve(rhs(),dst); |
| 651 | } |
| 652 | }; |
| 653 | |
| 654 | template<typename Derived, typename Rhs> |
| 655 | struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs> |
| 656 | : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> |
| 657 | { |
| 658 | typedef SimplicialCholeskyBase<Derived> Dec; |
| 659 | EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) |
| 660 | |
| 661 | template<typename Dest> void evalTo(Dest& dst) const |
| 662 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 663 | this->defaultEvalTo(dst); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 664 | } |
| 665 | }; |
| 666 | |
| 667 | } // end namespace internal |
| 668 | |
| 669 | } // end namespace Eigen |
| 670 | |
| 671 | #endif // EIGEN_SIMPLICIAL_CHOLESKY_H |