CholmodSupport.h (25441B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 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 10 #ifndef EIGEN_CHOLMODSUPPORT_H 11 #define EIGEN_CHOLMODSUPPORT_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 template<typename Scalar> struct cholmod_configure_matrix; 18 19 template<> struct cholmod_configure_matrix<double> { 20 template<typename CholmodType> 21 static void run(CholmodType& mat) { 22 mat.xtype = CHOLMOD_REAL; 23 mat.dtype = CHOLMOD_DOUBLE; 24 } 25 }; 26 27 template<> struct cholmod_configure_matrix<std::complex<double> > { 28 template<typename CholmodType> 29 static void run(CholmodType& mat) { 30 mat.xtype = CHOLMOD_COMPLEX; 31 mat.dtype = CHOLMOD_DOUBLE; 32 } 33 }; 34 35 // Other scalar types are not yet supported by Cholmod 36 // template<> struct cholmod_configure_matrix<float> { 37 // template<typename CholmodType> 38 // static void run(CholmodType& mat) { 39 // mat.xtype = CHOLMOD_REAL; 40 // mat.dtype = CHOLMOD_SINGLE; 41 // } 42 // }; 43 // 44 // template<> struct cholmod_configure_matrix<std::complex<float> > { 45 // template<typename CholmodType> 46 // static void run(CholmodType& mat) { 47 // mat.xtype = CHOLMOD_COMPLEX; 48 // mat.dtype = CHOLMOD_SINGLE; 49 // } 50 // }; 51 52 } // namespace internal 53 54 /** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object. 55 * Note that the data are shared. 56 */ 57 template<typename _Scalar, int _Options, typename _StorageIndex> 58 cholmod_sparse viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_StorageIndex> > mat) 59 { 60 cholmod_sparse res; 61 res.nzmax = mat.nonZeros(); 62 res.nrow = mat.rows(); 63 res.ncol = mat.cols(); 64 res.p = mat.outerIndexPtr(); 65 res.i = mat.innerIndexPtr(); 66 res.x = mat.valuePtr(); 67 res.z = 0; 68 res.sorted = 1; 69 if(mat.isCompressed()) 70 { 71 res.packed = 1; 72 res.nz = 0; 73 } 74 else 75 { 76 res.packed = 0; 77 res.nz = mat.innerNonZeroPtr(); 78 } 79 80 res.dtype = 0; 81 res.stype = -1; 82 83 if (internal::is_same<_StorageIndex,int>::value) 84 { 85 res.itype = CHOLMOD_INT; 86 } 87 else if (internal::is_same<_StorageIndex,SuiteSparse_long>::value) 88 { 89 res.itype = CHOLMOD_LONG; 90 } 91 else 92 { 93 eigen_assert(false && "Index type not supported yet"); 94 } 95 96 // setup res.xtype 97 internal::cholmod_configure_matrix<_Scalar>::run(res); 98 99 res.stype = 0; 100 101 return res; 102 } 103 104 template<typename _Scalar, int _Options, typename _Index> 105 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat) 106 { 107 cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived())); 108 return res; 109 } 110 111 template<typename _Scalar, int _Options, typename _Index> 112 const cholmod_sparse viewAsCholmod(const SparseVector<_Scalar,_Options,_Index>& mat) 113 { 114 cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived())); 115 return res; 116 } 117 118 /** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix. 119 * The data are not copied but shared. */ 120 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo> 121 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat) 122 { 123 cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.matrix().const_cast_derived())); 124 125 if(UpLo==Upper) res.stype = 1; 126 if(UpLo==Lower) res.stype = -1; 127 // swap stype for rowmajor matrices (only works for real matrices) 128 EIGEN_STATIC_ASSERT((_Options & RowMajorBit) == 0 || NumTraits<_Scalar>::IsComplex == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 129 if(_Options & RowMajorBit) res.stype *=-1; 130 131 return res; 132 } 133 134 /** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix. 135 * The data are not copied but shared. */ 136 template<typename Derived> 137 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat) 138 { 139 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 140 typedef typename Derived::Scalar Scalar; 141 142 cholmod_dense res; 143 res.nrow = mat.rows(); 144 res.ncol = mat.cols(); 145 res.nzmax = res.nrow * res.ncol; 146 res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride(); 147 res.x = (void*)(mat.derived().data()); 148 res.z = 0; 149 150 internal::cholmod_configure_matrix<Scalar>::run(res); 151 152 return res; 153 } 154 155 /** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix. 156 * The data are not copied but shared. */ 157 template<typename Scalar, int Flags, typename StorageIndex> 158 MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm) 159 { 160 return MappedSparseMatrix<Scalar,Flags,StorageIndex> 161 (cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol], 162 static_cast<StorageIndex*>(cm.p), static_cast<StorageIndex*>(cm.i),static_cast<Scalar*>(cm.x) ); 163 } 164 165 namespace internal { 166 167 // template specializations for int and long that call the correct cholmod method 168 169 #define EIGEN_CHOLMOD_SPECIALIZE0(ret, name) \ 170 template<typename _StorageIndex> inline ret cm_ ## name (cholmod_common &Common) { return cholmod_ ## name (&Common); } \ 171 template<> inline ret cm_ ## name<SuiteSparse_long> (cholmod_common &Common) { return cholmod_l_ ## name (&Common); } 172 173 #define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1) \ 174 template<typename _StorageIndex> inline ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_ ## name (&a1, &Common); } \ 175 template<> inline ret cm_ ## name<SuiteSparse_long> (t1& a1, cholmod_common &Common) { return cholmod_l_ ## name (&a1, &Common); } 176 177 EIGEN_CHOLMOD_SPECIALIZE0(int, start) 178 EIGEN_CHOLMOD_SPECIALIZE0(int, finish) 179 180 EIGEN_CHOLMOD_SPECIALIZE1(int, free_factor, cholmod_factor*, L) 181 EIGEN_CHOLMOD_SPECIALIZE1(int, free_dense, cholmod_dense*, X) 182 EIGEN_CHOLMOD_SPECIALIZE1(int, free_sparse, cholmod_sparse*, A) 183 184 EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A) 185 186 template<typename _StorageIndex> inline cholmod_dense* cm_solve (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_solve (sys, &L, &B, &Common); } 187 template<> inline cholmod_dense* cm_solve<SuiteSparse_long> (int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) { return cholmod_l_solve (sys, &L, &B, &Common); } 188 189 template<typename _StorageIndex> inline cholmod_sparse* cm_spsolve (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_spsolve (sys, &L, &B, &Common); } 190 template<> inline cholmod_sparse* cm_spsolve<SuiteSparse_long> (int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) { return cholmod_l_spsolve (sys, &L, &B, &Common); } 191 192 template<typename _StorageIndex> 193 inline int cm_factorize_p (cholmod_sparse* A, double beta[2], _StorageIndex* fset, std::size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_factorize_p (A, beta, fset, fsize, L, &Common); } 194 template<> 195 inline int cm_factorize_p<SuiteSparse_long> (cholmod_sparse* A, double beta[2], SuiteSparse_long* fset, std::size_t fsize, cholmod_factor* L, cholmod_common &Common) { return cholmod_l_factorize_p (A, beta, fset, fsize, L, &Common); } 196 197 #undef EIGEN_CHOLMOD_SPECIALIZE0 198 #undef EIGEN_CHOLMOD_SPECIALIZE1 199 200 } // namespace internal 201 202 203 enum CholmodMode { 204 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt 205 }; 206 207 208 /** \ingroup CholmodSupport_Module 209 * \class CholmodBase 210 * \brief The base class for the direct Cholesky factorization of Cholmod 211 * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT 212 */ 213 template<typename _MatrixType, int _UpLo, typename Derived> 214 class CholmodBase : public SparseSolverBase<Derived> 215 { 216 protected: 217 typedef SparseSolverBase<Derived> Base; 218 using Base::derived; 219 using Base::m_isInitialized; 220 public: 221 typedef _MatrixType MatrixType; 222 enum { UpLo = _UpLo }; 223 typedef typename MatrixType::Scalar Scalar; 224 typedef typename MatrixType::RealScalar RealScalar; 225 typedef MatrixType CholMatrixType; 226 typedef typename MatrixType::StorageIndex StorageIndex; 227 enum { 228 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 229 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 230 }; 231 232 public: 233 234 CholmodBase() 235 : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false) 236 { 237 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); 238 m_shiftOffset[0] = m_shiftOffset[1] = 0.0; 239 internal::cm_start<StorageIndex>(m_cholmod); 240 } 241 242 explicit CholmodBase(const MatrixType& matrix) 243 : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false) 244 { 245 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY); 246 m_shiftOffset[0] = m_shiftOffset[1] = 0.0; 247 internal::cm_start<StorageIndex>(m_cholmod); 248 compute(matrix); 249 } 250 251 ~CholmodBase() 252 { 253 if(m_cholmodFactor) 254 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod); 255 internal::cm_finish<StorageIndex>(m_cholmod); 256 } 257 258 inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); } 259 inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); } 260 261 /** \brief Reports whether previous computation was successful. 262 * 263 * \returns \c Success if computation was successful, 264 * \c NumericalIssue if the matrix.appears to be negative. 265 */ 266 ComputationInfo info() const 267 { 268 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 269 return m_info; 270 } 271 272 /** Computes the sparse Cholesky decomposition of \a matrix */ 273 Derived& compute(const MatrixType& matrix) 274 { 275 analyzePattern(matrix); 276 factorize(matrix); 277 return derived(); 278 } 279 280 /** Performs a symbolic decomposition on the sparsity pattern of \a matrix. 281 * 282 * This function is particularly useful when solving for several problems having the same structure. 283 * 284 * \sa factorize() 285 */ 286 void analyzePattern(const MatrixType& matrix) 287 { 288 if(m_cholmodFactor) 289 { 290 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod); 291 m_cholmodFactor = 0; 292 } 293 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>()); 294 m_cholmodFactor = internal::cm_analyze<StorageIndex>(A, m_cholmod); 295 296 this->m_isInitialized = true; 297 this->m_info = Success; 298 m_analysisIsOk = true; 299 m_factorizationIsOk = false; 300 } 301 302 /** Performs a numeric decomposition of \a matrix 303 * 304 * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been performed. 305 * 306 * \sa analyzePattern() 307 */ 308 void factorize(const MatrixType& matrix) 309 { 310 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 311 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>()); 312 internal::cm_factorize_p<StorageIndex>(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod); 313 314 // If the factorization failed, minor is the column at which it did. On success minor == n. 315 this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue); 316 m_factorizationIsOk = true; 317 } 318 319 /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations. 320 * See the Cholmod user guide for details. */ 321 cholmod_common& cholmod() { return m_cholmod; } 322 323 #ifndef EIGEN_PARSED_BY_DOXYGEN 324 /** \internal */ 325 template<typename Rhs,typename Dest> 326 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 327 { 328 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 329 const Index size = m_cholmodFactor->n; 330 EIGEN_UNUSED_VARIABLE(size); 331 eigen_assert(size==b.rows()); 332 333 // Cholmod needs column-major storage without inner-stride, which corresponds to the default behavior of Ref. 334 Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b.derived()); 335 336 cholmod_dense b_cd = viewAsCholmod(b_ref); 337 cholmod_dense* x_cd = internal::cm_solve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod); 338 if(!x_cd) 339 { 340 this->m_info = NumericalIssue; 341 return; 342 } 343 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) 344 // NOTE Actually, the copy can be avoided by calling cholmod_solve2 instead of cholmod_solve 345 dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols()); 346 internal::cm_free_dense<StorageIndex>(x_cd, m_cholmod); 347 } 348 349 /** \internal */ 350 template<typename RhsDerived, typename DestDerived> 351 void _solve_impl(const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest) const 352 { 353 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 354 const Index size = m_cholmodFactor->n; 355 EIGEN_UNUSED_VARIABLE(size); 356 eigen_assert(size==b.rows()); 357 358 // note: cs stands for Cholmod Sparse 359 Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived()); 360 cholmod_sparse b_cs = viewAsCholmod(b_ref); 361 cholmod_sparse* x_cs = internal::cm_spsolve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod); 362 if(!x_cs) 363 { 364 this->m_info = NumericalIssue; 365 return; 366 } 367 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) 368 // NOTE cholmod_spsolve in fact just calls the dense solver for blocks of 4 columns at a time (similar to Eigen's sparse solver) 369 dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs); 370 internal::cm_free_sparse<StorageIndex>(x_cs, m_cholmod); 371 } 372 #endif // EIGEN_PARSED_BY_DOXYGEN 373 374 375 /** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization. 376 * 377 * During the numerical factorization, an offset term is added to the diagonal coefficients:\n 378 * \c d_ii = \a offset + \c d_ii 379 * 380 * The default is \a offset=0. 381 * 382 * \returns a reference to \c *this. 383 */ 384 Derived& setShift(const RealScalar& offset) 385 { 386 m_shiftOffset[0] = double(offset); 387 return derived(); 388 } 389 390 /** \returns the determinant of the underlying matrix from the current factorization */ 391 Scalar determinant() const 392 { 393 using std::exp; 394 return exp(logDeterminant()); 395 } 396 397 /** \returns the log determinant of the underlying matrix from the current factorization */ 398 Scalar logDeterminant() const 399 { 400 using std::log; 401 using numext::real; 402 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 403 404 RealScalar logDet = 0; 405 Scalar *x = static_cast<Scalar*>(m_cholmodFactor->x); 406 if (m_cholmodFactor->is_super) 407 { 408 // Supernodal factorization stored as a packed list of dense column-major blocs, 409 // as described by the following structure: 410 411 // super[k] == index of the first column of the j-th super node 412 StorageIndex *super = static_cast<StorageIndex*>(m_cholmodFactor->super); 413 // pi[k] == offset to the description of row indices 414 StorageIndex *pi = static_cast<StorageIndex*>(m_cholmodFactor->pi); 415 // px[k] == offset to the respective dense block 416 StorageIndex *px = static_cast<StorageIndex*>(m_cholmodFactor->px); 417 418 Index nb_super_nodes = m_cholmodFactor->nsuper; 419 for (Index k=0; k < nb_super_nodes; ++k) 420 { 421 StorageIndex ncols = super[k + 1] - super[k]; 422 StorageIndex nrows = pi[k + 1] - pi[k]; 423 424 Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows+1)); 425 logDet += sk.real().log().sum(); 426 } 427 } 428 else 429 { 430 // Simplicial factorization stored as standard CSC matrix. 431 StorageIndex *p = static_cast<StorageIndex*>(m_cholmodFactor->p); 432 Index size = m_cholmodFactor->n; 433 for (Index k=0; k<size; ++k) 434 logDet += log(real( x[p[k]] )); 435 } 436 if (m_cholmodFactor->is_ll) 437 logDet *= 2.0; 438 return logDet; 439 }; 440 441 template<typename Stream> 442 void dumpMemory(Stream& /*s*/) 443 {} 444 445 protected: 446 mutable cholmod_common m_cholmod; 447 cholmod_factor* m_cholmodFactor; 448 double m_shiftOffset[2]; 449 mutable ComputationInfo m_info; 450 int m_factorizationIsOk; 451 int m_analysisIsOk; 452 }; 453 454 /** \ingroup CholmodSupport_Module 455 * \class CholmodSimplicialLLT 456 * \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod 457 * 458 * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization 459 * using the Cholmod library. 460 * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical interest. 461 * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices 462 * X and B can be either dense or sparse. 463 * 464 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 465 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 466 * or Upper. Default is Lower. 467 * 468 * \implsparsesolverconcept 469 * 470 * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. 471 * 472 * \warning Only double precision real and complex scalar types are supported by Cholmod. 473 * 474 * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLLT 475 */ 476 template<typename _MatrixType, int _UpLo = Lower> 477 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> > 478 { 479 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base; 480 using Base::m_cholmod; 481 482 public: 483 484 typedef _MatrixType MatrixType; 485 486 CholmodSimplicialLLT() : Base() { init(); } 487 488 CholmodSimplicialLLT(const MatrixType& matrix) : Base() 489 { 490 init(); 491 this->compute(matrix); 492 } 493 494 ~CholmodSimplicialLLT() {} 495 protected: 496 void init() 497 { 498 m_cholmod.final_asis = 0; 499 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 500 m_cholmod.final_ll = 1; 501 } 502 }; 503 504 505 /** \ingroup CholmodSupport_Module 506 * \class CholmodSimplicialLDLT 507 * \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod 508 * 509 * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization 510 * using the Cholmod library. 511 * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical interest. 512 * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices 513 * X and B can be either dense or sparse. 514 * 515 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 516 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 517 * or Upper. Default is Lower. 518 * 519 * \implsparsesolverconcept 520 * 521 * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. 522 * 523 * \warning Only double precision real and complex scalar types are supported by Cholmod. 524 * 525 * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLDLT 526 */ 527 template<typename _MatrixType, int _UpLo = Lower> 528 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> > 529 { 530 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base; 531 using Base::m_cholmod; 532 533 public: 534 535 typedef _MatrixType MatrixType; 536 537 CholmodSimplicialLDLT() : Base() { init(); } 538 539 CholmodSimplicialLDLT(const MatrixType& matrix) : Base() 540 { 541 init(); 542 this->compute(matrix); 543 } 544 545 ~CholmodSimplicialLDLT() {} 546 protected: 547 void init() 548 { 549 m_cholmod.final_asis = 1; 550 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 551 } 552 }; 553 554 /** \ingroup CholmodSupport_Module 555 * \class CholmodSupernodalLLT 556 * \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod 557 * 558 * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization 559 * using the Cholmod library. 560 * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM. 561 * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices 562 * X and B can be either dense or sparse. 563 * 564 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 565 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 566 * or Upper. Default is Lower. 567 * 568 * \implsparsesolverconcept 569 * 570 * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. 571 * 572 * \warning Only double precision real and complex scalar types are supported by Cholmod. 573 * 574 * \sa \ref TutorialSparseSolverConcept 575 */ 576 template<typename _MatrixType, int _UpLo = Lower> 577 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> > 578 { 579 typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base; 580 using Base::m_cholmod; 581 582 public: 583 584 typedef _MatrixType MatrixType; 585 586 CholmodSupernodalLLT() : Base() { init(); } 587 588 CholmodSupernodalLLT(const MatrixType& matrix) : Base() 589 { 590 init(); 591 this->compute(matrix); 592 } 593 594 ~CholmodSupernodalLLT() {} 595 protected: 596 void init() 597 { 598 m_cholmod.final_asis = 1; 599 m_cholmod.supernodal = CHOLMOD_SUPERNODAL; 600 } 601 }; 602 603 /** \ingroup CholmodSupport_Module 604 * \class CholmodDecomposition 605 * \brief A general Cholesky factorization and solver based on Cholmod 606 * 607 * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization 608 * using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices 609 * X and B can be either dense or sparse. 610 * 611 * This variant permits to change the underlying Cholesky method at runtime. 612 * On the other hand, it does not provide access to the result of the factorization. 613 * The default is to let Cholmod automatically choose between a simplicial and supernodal factorization. 614 * 615 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 616 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 617 * or Upper. Default is Lower. 618 * 619 * \implsparsesolverconcept 620 * 621 * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. 622 * 623 * \warning Only double precision real and complex scalar types are supported by Cholmod. 624 * 625 * \sa \ref TutorialSparseSolverConcept 626 */ 627 template<typename _MatrixType, int _UpLo = Lower> 628 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> > 629 { 630 typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base; 631 using Base::m_cholmod; 632 633 public: 634 635 typedef _MatrixType MatrixType; 636 637 CholmodDecomposition() : Base() { init(); } 638 639 CholmodDecomposition(const MatrixType& matrix) : Base() 640 { 641 init(); 642 this->compute(matrix); 643 } 644 645 ~CholmodDecomposition() {} 646 647 void setMode(CholmodMode mode) 648 { 649 switch(mode) 650 { 651 case CholmodAuto: 652 m_cholmod.final_asis = 1; 653 m_cholmod.supernodal = CHOLMOD_AUTO; 654 break; 655 case CholmodSimplicialLLt: 656 m_cholmod.final_asis = 0; 657 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 658 m_cholmod.final_ll = 1; 659 break; 660 case CholmodSupernodalLLt: 661 m_cholmod.final_asis = 1; 662 m_cholmod.supernodal = CHOLMOD_SUPERNODAL; 663 break; 664 case CholmodLDLt: 665 m_cholmod.final_asis = 1; 666 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 667 break; 668 default: 669 break; 670 } 671 } 672 protected: 673 void init() 674 { 675 m_cholmod.final_asis = 1; 676 m_cholmod.supernodal = CHOLMOD_AUTO; 677 } 678 }; 679 680 } // end namespace Eigen 681 682 #endif // EIGEN_CHOLMODSUPPORT_H