1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) 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_PASTIXSUPPORT_H 11 #define EIGEN_PASTIXSUPPORT_H 12 13 namespace Eigen { 14 15 /** \ingroup PaStiXSupport_Module 16 * \brief Interface to the PaStix solver 17 * 18 * This class is used to solve the linear systems A.X = B via the PaStix library. 19 * The matrix can be either real or complex, symmetric or not. 20 * 21 * \sa TutorialSparseDirectSolvers 22 */ 23 template<typename _MatrixType, bool IsStrSym = false> class PastixLU; 24 template<typename _MatrixType, int Options> class PastixLLT; 25 template<typename _MatrixType, int Options> class PastixLDLT; 26 27 namespace internal 28 { 29 30 template<class Pastix> struct pastix_traits; 31 32 template<typename _MatrixType> 33 struct pastix_traits< PastixLU<_MatrixType> > 34 { 35 typedef _MatrixType MatrixType; 36 typedef typename _MatrixType::Scalar Scalar; 37 typedef typename _MatrixType::RealScalar RealScalar; 38 typedef typename _MatrixType::Index Index; 39 }; 40 41 template<typename _MatrixType, int Options> 42 struct pastix_traits< PastixLLT<_MatrixType,Options> > 43 { 44 typedef _MatrixType MatrixType; 45 typedef typename _MatrixType::Scalar Scalar; 46 typedef typename _MatrixType::RealScalar RealScalar; 47 typedef typename _MatrixType::Index Index; 48 }; 49 50 template<typename _MatrixType, int Options> 51 struct pastix_traits< PastixLDLT<_MatrixType,Options> > 52 { 53 typedef _MatrixType MatrixType; 54 typedef typename _MatrixType::Scalar Scalar; 55 typedef typename _MatrixType::RealScalar RealScalar; 56 typedef typename _MatrixType::Index Index; 57 }; 58 59 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm) 60 { 61 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 62 if (nbrhs == 0) {x = NULL; nbrhs=1;} 63 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 64 } 65 66 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm) 67 { 68 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 69 if (nbrhs == 0) {x = NULL; nbrhs=1;} 70 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 71 } 72 73 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm) 74 { 75 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 76 if (nbrhs == 0) {x = NULL; nbrhs=1;} 77 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm); 78 } 79 80 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm) 81 { 82 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 83 if (nbrhs == 0) {x = NULL; nbrhs=1;} 84 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); 85 } 86 87 // Convert the matrix to Fortran-style Numbering 88 template <typename MatrixType> 89 void c_to_fortran_numbering (MatrixType& mat) 90 { 91 if ( !(mat.outerIndexPtr()[0]) ) 92 { 93 int i; 94 for(i = 0; i <= mat.rows(); ++i) 95 ++mat.outerIndexPtr()[i]; 96 for(i = 0; i < mat.nonZeros(); ++i) 97 ++mat.innerIndexPtr()[i]; 98 } 99 } 100 101 // Convert to C-style Numbering 102 template <typename MatrixType> 103 void fortran_to_c_numbering (MatrixType& mat) 104 { 105 // Check the Numbering 106 if ( mat.outerIndexPtr()[0] == 1 ) 107 { // Convert to C-style numbering 108 int i; 109 for(i = 0; i <= mat.rows(); ++i) 110 --mat.outerIndexPtr()[i]; 111 for(i = 0; i < mat.nonZeros(); ++i) 112 --mat.innerIndexPtr()[i]; 113 } 114 } 115 } 116 117 // This is the base class to interface with PaStiX functions. 118 // Users should not used this class directly. 119 template <class Derived> 120 class PastixBase : internal::noncopyable 121 { 122 public: 123 typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; 124 typedef _MatrixType MatrixType; 125 typedef typename MatrixType::Scalar Scalar; 126 typedef typename MatrixType::RealScalar RealScalar; 127 typedef typename MatrixType::Index Index; 128 typedef Matrix<Scalar,Dynamic,1> Vector; 129 typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix; 130 131 public: 132 133 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0) 134 { 135 init(); 136 } 137 138 ~PastixBase() 139 { 140 clean(); 141 } 142 143 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 144 * 145 * \sa compute() 146 */ 147 template<typename Rhs> 148 inline const internal::solve_retval<PastixBase, Rhs> 149 solve(const MatrixBase<Rhs>& b) const 150 { 151 eigen_assert(m_isInitialized && "Pastix solver is not initialized."); 152 eigen_assert(rows()==b.rows() 153 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 154 return internal::solve_retval<PastixBase, Rhs>(*this, b.derived()); 155 } 156 157 template<typename Rhs,typename Dest> 158 bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; 159 160 /** \internal */ 161 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> 162 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const 163 { 164 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 165 eigen_assert(rows()==b.rows()); 166 167 // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. 168 static const int NbColsAtOnce = 1; 169 int rhsCols = b.cols(); 170 int size = b.rows(); 171 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); 172 for(int k=0; k<rhsCols; k+=NbColsAtOnce) 173 { 174 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); 175 tmp.leftCols(actualCols) = b.middleCols(k,actualCols); 176 tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); 177 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); 178 } 179 } 180 181 Derived& derived() 182 { 183 return *static_cast<Derived*>(this); 184 } 185 const Derived& derived() const 186 { 187 return *static_cast<const Derived*>(this); 188 } 189 190 /** Returns a reference to the integer vector IPARM of PaStiX parameters 191 * to modify the default parameters. 192 * The statistics related to the different phases of factorization and solve are saved here as well 193 * \sa analyzePattern() factorize() 194 */ 195 Array<Index,IPARM_SIZE,1>& iparm() 196 { 197 return m_iparm; 198 } 199 200 /** Return a reference to a particular index parameter of the IPARM vector 201 * \sa iparm() 202 */ 203 204 int& iparm(int idxparam) 205 { 206 return m_iparm(idxparam); 207 } 208 209 /** Returns a reference to the double vector DPARM of PaStiX parameters 210 * The statistics related to the different phases of factorization and solve are saved here as well 211 * \sa analyzePattern() factorize() 212 */ 213 Array<RealScalar,IPARM_SIZE,1>& dparm() 214 { 215 return m_dparm; 216 } 217 218 219 /** Return a reference to a particular index parameter of the DPARM vector 220 * \sa dparm() 221 */ 222 double& dparm(int idxparam) 223 { 224 return m_dparm(idxparam); 225 } 226 227 inline Index cols() const { return m_size; } 228 inline Index rows() const { return m_size; } 229 230 /** \brief Reports whether previous computation was successful. 231 * 232 * \returns \c Success if computation was succesful, 233 * \c NumericalIssue if the PaStiX reports a problem 234 * \c InvalidInput if the input matrix is invalid 235 * 236 * \sa iparm() 237 */ 238 ComputationInfo info() const 239 { 240 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 241 return m_info; 242 } 243 244 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 245 * 246 * \sa compute() 247 */ 248 template<typename Rhs> 249 inline const internal::sparse_solve_retval<PastixBase, Rhs> 250 solve(const SparseMatrixBase<Rhs>& b) const 251 { 252 eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized."); 253 eigen_assert(rows()==b.rows() 254 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 255 return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); 256 } 257 258 protected: 259 260 // Initialize the Pastix data structure, check the matrix 261 void init(); 262 263 // Compute the ordering and the symbolic factorization 264 void analyzePattern(ColSpMatrix& mat); 265 266 // Compute the numerical factorization 267 void factorize(ColSpMatrix& mat); 268 269 // Free all the data allocated by Pastix 270 void clean() 271 { 272 eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 273 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; 274 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; 275 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 276 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 277 } 278 279 void compute(ColSpMatrix& mat); 280 281 int m_initisOk; 282 int m_analysisIsOk; 283 int m_factorizationIsOk; 284 bool m_isInitialized; 285 mutable ComputationInfo m_info; 286 mutable pastix_data_t *m_pastixdata; // Data structure for pastix 287 mutable int m_comm; // The MPI communicator identifier 288 mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters 289 mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters 290 mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector 291 mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector 292 mutable int m_size; // Size of the matrix 293 }; 294 295 /** Initialize the PaStiX data structure. 296 *A first call to this function fills iparm and dparm with the default PaStiX parameters 297 * \sa iparm() dparm() 298 */ 299 template <class Derived> 300 void PastixBase<Derived>::init() 301 { 302 m_size = 0; 303 m_iparm.setZero(IPARM_SIZE); 304 m_dparm.setZero(DPARM_SIZE); 305 306 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; 307 pastix(&m_pastixdata, MPI_COMM_WORLD, 308 0, 0, 0, 0, 309 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); 310 311 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; 312 m_iparm[IPARM_VERBOSE] = 2; 313 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; 314 m_iparm[IPARM_INCOMPLETE] = API_NO; 315 m_iparm[IPARM_OOC_LIMIT] = 2000; 316 m_iparm[IPARM_RHS_MAKING] = API_RHS_B; 317 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 318 319 m_iparm(IPARM_START_TASK) = API_TASK_INIT; 320 m_iparm(IPARM_END_TASK) = API_TASK_INIT; 321 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 322 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); 323 324 // Check the returned error 325 if(m_iparm(IPARM_ERROR_NUMBER)) { 326 m_info = InvalidInput; 327 m_initisOk = false; 328 } 329 else { 330 m_info = Success; 331 m_initisOk = true; 332 } 333 } 334 335 template <class Derived> 336 void PastixBase<Derived>::compute(ColSpMatrix& mat) 337 { 338 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); 339 340 analyzePattern(mat); 341 factorize(mat); 342 343 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 344 m_isInitialized = m_factorizationIsOk; 345 } 346 347 348 template <class Derived> 349 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) 350 { 351 eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); 352 353 // clean previous calls 354 if(m_size>0) 355 clean(); 356 357 m_size = mat.rows(); 358 m_perm.resize(m_size); 359 m_invp.resize(m_size); 360 361 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; 362 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; 363 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 364 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 365 366 // Check the returned error 367 if(m_iparm(IPARM_ERROR_NUMBER)) 368 { 369 m_info = NumericalIssue; 370 m_analysisIsOk = false; 371 } 372 else 373 { 374 m_info = Success; 375 m_analysisIsOk = true; 376 } 377 } 378 379 template <class Derived> 380 void PastixBase<Derived>::factorize(ColSpMatrix& mat) 381 { 382 // if(&m_cpyMat != &mat) m_cpyMat = mat; 383 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); 384 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; 385 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; 386 m_size = mat.rows(); 387 388 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 389 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 390 391 // Check the returned error 392 if(m_iparm(IPARM_ERROR_NUMBER)) 393 { 394 m_info = NumericalIssue; 395 m_factorizationIsOk = false; 396 m_isInitialized = false; 397 } 398 else 399 { 400 m_info = Success; 401 m_factorizationIsOk = true; 402 m_isInitialized = true; 403 } 404 } 405 406 /* Solve the system */ 407 template<typename Base> 408 template<typename Rhs,typename Dest> 409 bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const 410 { 411 eigen_assert(m_isInitialized && "The matrix should be factorized first"); 412 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 413 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 414 int rhs = 1; 415 416 x = b; /* on return, x is overwritten by the computed solution */ 417 418 for (int i = 0; i < b.cols(); i++){ 419 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; 420 m_iparm[IPARM_END_TASK] = API_TASK_REFINE; 421 422 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0, 423 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); 424 } 425 426 // Check the returned error 427 m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; 428 429 return m_iparm(IPARM_ERROR_NUMBER)==0; 430 } 431 432 /** \ingroup PaStiXSupport_Module 433 * \class PastixLU 434 * \brief Sparse direct LU solver based on PaStiX library 435 * 436 * This class is used to solve the linear systems A.X = B with a supernodal LU 437 * factorization in the PaStiX library. The matrix A should be squared and nonsingular 438 * PaStiX requires that the matrix A has a symmetric structural pattern. 439 * This interface can symmetrize the input matrix otherwise. 440 * The vectors or matrices X and B can be either dense or sparse. 441 * 442 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 443 * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false 444 * NOTE : Note that if the analysis and factorization phase are called separately, 445 * the input matrix will be symmetrized at each call, hence it is advised to 446 * symmetrize the matrix in a end-user program and set \p IsStrSym to true 447 * 448 * \sa \ref TutorialSparseDirectSolvers 449 * 450 */ 451 template<typename _MatrixType, bool IsStrSym> 452 class PastixLU : public PastixBase< PastixLU<_MatrixType> > 453 { 454 public: 455 typedef _MatrixType MatrixType; 456 typedef PastixBase<PastixLU<MatrixType> > Base; 457 typedef typename Base::ColSpMatrix ColSpMatrix; 458 typedef typename MatrixType::Index Index; 459 460 public: 461 PastixLU() : Base() 462 { 463 init(); 464 } 465 466 PastixLU(const MatrixType& matrix):Base() 467 { 468 init(); 469 compute(matrix); 470 } 471 /** Compute the LU supernodal factorization of \p matrix. 472 * iparm and dparm can be used to tune the PaStiX parameters. 473 * see the PaStiX user's manual 474 * \sa analyzePattern() factorize() 475 */ 476 void compute (const MatrixType& matrix) 477 { 478 m_structureIsUptodate = false; 479 ColSpMatrix temp; 480 grabMatrix(matrix, temp); 481 Base::compute(temp); 482 } 483 /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 484 * Several ordering methods can be used at this step. See the PaStiX user's manual. 485 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 486 * \sa factorize() 487 */ 488 void analyzePattern(const MatrixType& matrix) 489 { 490 m_structureIsUptodate = false; 491 ColSpMatrix temp; 492 grabMatrix(matrix, temp); 493 Base::analyzePattern(temp); 494 } 495 496 /** Compute the LU supernodal factorization of \p matrix 497 * WARNING The matrix \p matrix should have the same structural pattern 498 * as the same used in the analysis phase. 499 * \sa analyzePattern() 500 */ 501 void factorize(const MatrixType& matrix) 502 { 503 ColSpMatrix temp; 504 grabMatrix(matrix, temp); 505 Base::factorize(temp); 506 } 507 protected: 508 509 void init() 510 { 511 m_structureIsUptodate = false; 512 m_iparm(IPARM_SYM) = API_SYM_NO; 513 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; 514 } 515 516 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 517 { 518 if(IsStrSym) 519 out = matrix; 520 else 521 { 522 if(!m_structureIsUptodate) 523 { 524 // update the transposed structure 525 m_transposedStructure = matrix.transpose(); 526 527 // Set the elements of the matrix to zero 528 for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 529 for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) 530 it.valueRef() = 0.0; 531 532 m_structureIsUptodate = true; 533 } 534 535 out = m_transposedStructure + matrix; 536 } 537 internal::c_to_fortran_numbering(out); 538 } 539 540 using Base::m_iparm; 541 using Base::m_dparm; 542 543 ColSpMatrix m_transposedStructure; 544 bool m_structureIsUptodate; 545 }; 546 547 /** \ingroup PaStiXSupport_Module 548 * \class PastixLLT 549 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 550 * 551 * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization 552 * available in the PaStiX library. The matrix A should be symmetric and positive definite 553 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 554 * The vectors or matrices X and B can be either dense or sparse 555 * 556 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 557 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 558 * 559 * \sa \ref TutorialSparseDirectSolvers 560 */ 561 template<typename _MatrixType, int _UpLo> 562 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > 563 { 564 public: 565 typedef _MatrixType MatrixType; 566 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; 567 typedef typename Base::ColSpMatrix ColSpMatrix; 568 569 public: 570 enum { UpLo = _UpLo }; 571 PastixLLT() : Base() 572 { 573 init(); 574 } 575 576 PastixLLT(const MatrixType& matrix):Base() 577 { 578 init(); 579 compute(matrix); 580 } 581 582 /** Compute the L factor of the LL^T supernodal factorization of \p matrix 583 * \sa analyzePattern() factorize() 584 */ 585 void compute (const MatrixType& matrix) 586 { 587 ColSpMatrix temp; 588 grabMatrix(matrix, temp); 589 Base::compute(temp); 590 } 591 592 /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern 593 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 594 * \sa factorize() 595 */ 596 void analyzePattern(const MatrixType& matrix) 597 { 598 ColSpMatrix temp; 599 grabMatrix(matrix, temp); 600 Base::analyzePattern(temp); 601 } 602 /** Compute the LL^T supernodal numerical factorization of \p matrix 603 * \sa analyzePattern() 604 */ 605 void factorize(const MatrixType& matrix) 606 { 607 ColSpMatrix temp; 608 grabMatrix(matrix, temp); 609 Base::factorize(temp); 610 } 611 protected: 612 using Base::m_iparm; 613 614 void init() 615 { 616 m_iparm(IPARM_SYM) = API_SYM_YES; 617 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; 618 } 619 620 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 621 { 622 // Pastix supports only lower, column-major matrices 623 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 624 internal::c_to_fortran_numbering(out); 625 } 626 }; 627 628 /** \ingroup PaStiXSupport_Module 629 * \class PastixLDLT 630 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 631 * 632 * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization 633 * available in the PaStiX library. The matrix A should be symmetric and positive definite 634 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 635 * The vectors or matrices X and B can be either dense or sparse 636 * 637 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 638 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 639 * 640 * \sa \ref TutorialSparseDirectSolvers 641 */ 642 template<typename _MatrixType, int _UpLo> 643 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > 644 { 645 public: 646 typedef _MatrixType MatrixType; 647 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 648 typedef typename Base::ColSpMatrix ColSpMatrix; 649 650 public: 651 enum { UpLo = _UpLo }; 652 PastixLDLT():Base() 653 { 654 init(); 655 } 656 657 PastixLDLT(const MatrixType& matrix):Base() 658 { 659 init(); 660 compute(matrix); 661 } 662 663 /** Compute the L and D factors of the LDL^T factorization of \p matrix 664 * \sa analyzePattern() factorize() 665 */ 666 void compute (const MatrixType& matrix) 667 { 668 ColSpMatrix temp; 669 grabMatrix(matrix, temp); 670 Base::compute(temp); 671 } 672 673 /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern 674 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 675 * \sa factorize() 676 */ 677 void analyzePattern(const MatrixType& matrix) 678 { 679 ColSpMatrix temp; 680 grabMatrix(matrix, temp); 681 Base::analyzePattern(temp); 682 } 683 /** Compute the LDL^T supernodal numerical factorization of \p matrix 684 * 685 */ 686 void factorize(const MatrixType& matrix) 687 { 688 ColSpMatrix temp; 689 grabMatrix(matrix, temp); 690 Base::factorize(temp); 691 } 692 693 protected: 694 using Base::m_iparm; 695 696 void init() 697 { 698 m_iparm(IPARM_SYM) = API_SYM_YES; 699 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; 700 } 701 702 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 703 { 704 // Pastix supports only lower, column-major matrices 705 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 706 internal::c_to_fortran_numbering(out); 707 } 708 }; 709 710 namespace internal { 711 712 template<typename _MatrixType, typename Rhs> 713 struct solve_retval<PastixBase<_MatrixType>, Rhs> 714 : solve_retval_base<PastixBase<_MatrixType>, Rhs> 715 { 716 typedef PastixBase<_MatrixType> Dec; 717 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 718 719 template<typename Dest> void evalTo(Dest& dst) const 720 { 721 dec()._solve(rhs(),dst); 722 } 723 }; 724 725 template<typename _MatrixType, typename Rhs> 726 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> 727 : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs> 728 { 729 typedef PastixBase<_MatrixType> Dec; 730 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 731 732 template<typename Dest> void evalTo(Dest& dst) const 733 { 734 dec()._solve_sparse(rhs(),dst); 735 } 736 }; 737 738 } // end namespace internal 739 740 } // end namespace Eigen 741 742 #endif 743
