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 Derived& derived() 161 { 162 return *static_cast<Derived*>(this); 163 } 164 const Derived& derived() const 165 { 166 return *static_cast<const Derived*>(this); 167 } 168 169 /** Returns a reference to the integer vector IPARM of PaStiX parameters 170 * to modify the default parameters. 171 * The statistics related to the different phases of factorization and solve are saved here as well 172 * \sa analyzePattern() factorize() 173 */ 174 Array<Index,IPARM_SIZE,1>& iparm() 175 { 176 return m_iparm; 177 } 178 179 /** Return a reference to a particular index parameter of the IPARM vector 180 * \sa iparm() 181 */ 182 183 int& iparm(int idxparam) 184 { 185 return m_iparm(idxparam); 186 } 187 188 /** Returns a reference to the double vector DPARM of PaStiX parameters 189 * The statistics related to the different phases of factorization and solve are saved here as well 190 * \sa analyzePattern() factorize() 191 */ 192 Array<RealScalar,IPARM_SIZE,1>& dparm() 193 { 194 return m_dparm; 195 } 196 197 198 /** Return a reference to a particular index parameter of the DPARM vector 199 * \sa dparm() 200 */ 201 double& dparm(int idxparam) 202 { 203 return m_dparm(idxparam); 204 } 205 206 inline Index cols() const { return m_size; } 207 inline Index rows() const { return m_size; } 208 209 /** \brief Reports whether previous computation was successful. 210 * 211 * \returns \c Success if computation was succesful, 212 * \c NumericalIssue if the PaStiX reports a problem 213 * \c InvalidInput if the input matrix is invalid 214 * 215 * \sa iparm() 216 */ 217 ComputationInfo info() const 218 { 219 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 220 return m_info; 221 } 222 223 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 224 * 225 * \sa compute() 226 */ 227 template<typename Rhs> 228 inline const internal::sparse_solve_retval<PastixBase, Rhs> 229 solve(const SparseMatrixBase<Rhs>& b) const 230 { 231 eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized."); 232 eigen_assert(rows()==b.rows() 233 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); 234 return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); 235 } 236 237 protected: 238 239 // Initialize the Pastix data structure, check the matrix 240 void init(); 241 242 // Compute the ordering and the symbolic factorization 243 void analyzePattern(ColSpMatrix& mat); 244 245 // Compute the numerical factorization 246 void factorize(ColSpMatrix& mat); 247 248 // Free all the data allocated by Pastix 249 void clean() 250 { 251 eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 252 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; 253 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; 254 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 255 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 256 } 257 258 void compute(ColSpMatrix& mat); 259 260 int m_initisOk; 261 int m_analysisIsOk; 262 int m_factorizationIsOk; 263 bool m_isInitialized; 264 mutable ComputationInfo m_info; 265 mutable pastix_data_t *m_pastixdata; // Data structure for pastix 266 mutable int m_comm; // The MPI communicator identifier 267 mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters 268 mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters 269 mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector 270 mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector 271 mutable int m_size; // Size of the matrix 272 }; 273 274 /** Initialize the PaStiX data structure. 275 *A first call to this function fills iparm and dparm with the default PaStiX parameters 276 * \sa iparm() dparm() 277 */ 278 template <class Derived> 279 void PastixBase<Derived>::init() 280 { 281 m_size = 0; 282 m_iparm.setZero(IPARM_SIZE); 283 m_dparm.setZero(DPARM_SIZE); 284 285 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; 286 pastix(&m_pastixdata, MPI_COMM_WORLD, 287 0, 0, 0, 0, 288 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); 289 290 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; 291 m_iparm[IPARM_VERBOSE] = 2; 292 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; 293 m_iparm[IPARM_INCOMPLETE] = API_NO; 294 m_iparm[IPARM_OOC_LIMIT] = 2000; 295 m_iparm[IPARM_RHS_MAKING] = API_RHS_B; 296 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 297 298 m_iparm(IPARM_START_TASK) = API_TASK_INIT; 299 m_iparm(IPARM_END_TASK) = API_TASK_INIT; 300 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 301 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); 302 303 // Check the returned error 304 if(m_iparm(IPARM_ERROR_NUMBER)) { 305 m_info = InvalidInput; 306 m_initisOk = false; 307 } 308 else { 309 m_info = Success; 310 m_initisOk = true; 311 } 312 } 313 314 template <class Derived> 315 void PastixBase<Derived>::compute(ColSpMatrix& mat) 316 { 317 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); 318 319 analyzePattern(mat); 320 factorize(mat); 321 322 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 323 m_isInitialized = m_factorizationIsOk; 324 } 325 326 327 template <class Derived> 328 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) 329 { 330 eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); 331 332 // clean previous calls 333 if(m_size>0) 334 clean(); 335 336 m_size = mat.rows(); 337 m_perm.resize(m_size); 338 m_invp.resize(m_size); 339 340 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; 341 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; 342 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 343 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 344 345 // Check the returned error 346 if(m_iparm(IPARM_ERROR_NUMBER)) 347 { 348 m_info = NumericalIssue; 349 m_analysisIsOk = false; 350 } 351 else 352 { 353 m_info = Success; 354 m_analysisIsOk = true; 355 } 356 } 357 358 template <class Derived> 359 void PastixBase<Derived>::factorize(ColSpMatrix& mat) 360 { 361 // if(&m_cpyMat != &mat) m_cpyMat = mat; 362 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); 363 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; 364 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; 365 m_size = mat.rows(); 366 367 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 368 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 369 370 // Check the returned error 371 if(m_iparm(IPARM_ERROR_NUMBER)) 372 { 373 m_info = NumericalIssue; 374 m_factorizationIsOk = false; 375 m_isInitialized = false; 376 } 377 else 378 { 379 m_info = Success; 380 m_factorizationIsOk = true; 381 m_isInitialized = true; 382 } 383 } 384 385 /* Solve the system */ 386 template<typename Base> 387 template<typename Rhs,typename Dest> 388 bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const 389 { 390 eigen_assert(m_isInitialized && "The matrix should be factorized first"); 391 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 392 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 393 int rhs = 1; 394 395 x = b; /* on return, x is overwritten by the computed solution */ 396 397 for (int i = 0; i < b.cols(); i++){ 398 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; 399 m_iparm[IPARM_END_TASK] = API_TASK_REFINE; 400 401 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0, 402 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); 403 } 404 405 // Check the returned error 406 m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; 407 408 return m_iparm(IPARM_ERROR_NUMBER)==0; 409 } 410 411 /** \ingroup PaStiXSupport_Module 412 * \class PastixLU 413 * \brief Sparse direct LU solver based on PaStiX library 414 * 415 * This class is used to solve the linear systems A.X = B with a supernodal LU 416 * factorization in the PaStiX library. The matrix A should be squared and nonsingular 417 * PaStiX requires that the matrix A has a symmetric structural pattern. 418 * This interface can symmetrize the input matrix otherwise. 419 * The vectors or matrices X and B can be either dense or sparse. 420 * 421 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 422 * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false 423 * NOTE : Note that if the analysis and factorization phase are called separately, 424 * the input matrix will be symmetrized at each call, hence it is advised to 425 * symmetrize the matrix in a end-user program and set \p IsStrSym to true 426 * 427 * \sa \ref TutorialSparseDirectSolvers 428 * 429 */ 430 template<typename _MatrixType, bool IsStrSym> 431 class PastixLU : public PastixBase< PastixLU<_MatrixType> > 432 { 433 public: 434 typedef _MatrixType MatrixType; 435 typedef PastixBase<PastixLU<MatrixType> > Base; 436 typedef typename Base::ColSpMatrix ColSpMatrix; 437 typedef typename MatrixType::Index Index; 438 439 public: 440 PastixLU() : Base() 441 { 442 init(); 443 } 444 445 PastixLU(const MatrixType& matrix):Base() 446 { 447 init(); 448 compute(matrix); 449 } 450 /** Compute the LU supernodal factorization of \p matrix. 451 * iparm and dparm can be used to tune the PaStiX parameters. 452 * see the PaStiX user's manual 453 * \sa analyzePattern() factorize() 454 */ 455 void compute (const MatrixType& matrix) 456 { 457 m_structureIsUptodate = false; 458 ColSpMatrix temp; 459 grabMatrix(matrix, temp); 460 Base::compute(temp); 461 } 462 /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 463 * Several ordering methods can be used at this step. See the PaStiX user's manual. 464 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 465 * \sa factorize() 466 */ 467 void analyzePattern(const MatrixType& matrix) 468 { 469 m_structureIsUptodate = false; 470 ColSpMatrix temp; 471 grabMatrix(matrix, temp); 472 Base::analyzePattern(temp); 473 } 474 475 /** Compute the LU supernodal factorization of \p matrix 476 * WARNING The matrix \p matrix should have the same structural pattern 477 * as the same used in the analysis phase. 478 * \sa analyzePattern() 479 */ 480 void factorize(const MatrixType& matrix) 481 { 482 ColSpMatrix temp; 483 grabMatrix(matrix, temp); 484 Base::factorize(temp); 485 } 486 protected: 487 488 void init() 489 { 490 m_structureIsUptodate = false; 491 m_iparm(IPARM_SYM) = API_SYM_NO; 492 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; 493 } 494 495 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 496 { 497 if(IsStrSym) 498 out = matrix; 499 else 500 { 501 if(!m_structureIsUptodate) 502 { 503 // update the transposed structure 504 m_transposedStructure = matrix.transpose(); 505 506 // Set the elements of the matrix to zero 507 for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 508 for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) 509 it.valueRef() = 0.0; 510 511 m_structureIsUptodate = true; 512 } 513 514 out = m_transposedStructure + matrix; 515 } 516 internal::c_to_fortran_numbering(out); 517 } 518 519 using Base::m_iparm; 520 using Base::m_dparm; 521 522 ColSpMatrix m_transposedStructure; 523 bool m_structureIsUptodate; 524 }; 525 526 /** \ingroup PaStiXSupport_Module 527 * \class PastixLLT 528 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 529 * 530 * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization 531 * available in the PaStiX library. The matrix A should be symmetric and positive definite 532 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 533 * The vectors or matrices X and B can be either dense or sparse 534 * 535 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 536 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 537 * 538 * \sa \ref TutorialSparseDirectSolvers 539 */ 540 template<typename _MatrixType, int _UpLo> 541 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > 542 { 543 public: 544 typedef _MatrixType MatrixType; 545 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; 546 typedef typename Base::ColSpMatrix ColSpMatrix; 547 548 public: 549 enum { UpLo = _UpLo }; 550 PastixLLT() : Base() 551 { 552 init(); 553 } 554 555 PastixLLT(const MatrixType& matrix):Base() 556 { 557 init(); 558 compute(matrix); 559 } 560 561 /** Compute the L factor of the LL^T supernodal factorization of \p matrix 562 * \sa analyzePattern() factorize() 563 */ 564 void compute (const MatrixType& matrix) 565 { 566 ColSpMatrix temp; 567 grabMatrix(matrix, temp); 568 Base::compute(temp); 569 } 570 571 /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern 572 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 573 * \sa factorize() 574 */ 575 void analyzePattern(const MatrixType& matrix) 576 { 577 ColSpMatrix temp; 578 grabMatrix(matrix, temp); 579 Base::analyzePattern(temp); 580 } 581 /** Compute the LL^T supernodal numerical factorization of \p matrix 582 * \sa analyzePattern() 583 */ 584 void factorize(const MatrixType& matrix) 585 { 586 ColSpMatrix temp; 587 grabMatrix(matrix, temp); 588 Base::factorize(temp); 589 } 590 protected: 591 using Base::m_iparm; 592 593 void init() 594 { 595 m_iparm(IPARM_SYM) = API_SYM_YES; 596 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; 597 } 598 599 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 600 { 601 // Pastix supports only lower, column-major matrices 602 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 603 internal::c_to_fortran_numbering(out); 604 } 605 }; 606 607 /** \ingroup PaStiXSupport_Module 608 * \class PastixLDLT 609 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library 610 * 611 * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization 612 * available in the PaStiX library. The matrix A should be symmetric and positive definite 613 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX 614 * The vectors or matrices X and B can be either dense or sparse 615 * 616 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 617 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX 618 * 619 * \sa \ref TutorialSparseDirectSolvers 620 */ 621 template<typename _MatrixType, int _UpLo> 622 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > 623 { 624 public: 625 typedef _MatrixType MatrixType; 626 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 627 typedef typename Base::ColSpMatrix ColSpMatrix; 628 629 public: 630 enum { UpLo = _UpLo }; 631 PastixLDLT():Base() 632 { 633 init(); 634 } 635 636 PastixLDLT(const MatrixType& matrix):Base() 637 { 638 init(); 639 compute(matrix); 640 } 641 642 /** Compute the L and D factors of the LDL^T factorization of \p matrix 643 * \sa analyzePattern() factorize() 644 */ 645 void compute (const MatrixType& matrix) 646 { 647 ColSpMatrix temp; 648 grabMatrix(matrix, temp); 649 Base::compute(temp); 650 } 651 652 /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern 653 * The result of this operation can be used with successive matrices having the same pattern as \p matrix 654 * \sa factorize() 655 */ 656 void analyzePattern(const MatrixType& matrix) 657 { 658 ColSpMatrix temp; 659 grabMatrix(matrix, temp); 660 Base::analyzePattern(temp); 661 } 662 /** Compute the LDL^T supernodal numerical factorization of \p matrix 663 * 664 */ 665 void factorize(const MatrixType& matrix) 666 { 667 ColSpMatrix temp; 668 grabMatrix(matrix, temp); 669 Base::factorize(temp); 670 } 671 672 protected: 673 using Base::m_iparm; 674 675 void init() 676 { 677 m_iparm(IPARM_SYM) = API_SYM_YES; 678 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; 679 } 680 681 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 682 { 683 // Pastix supports only lower, column-major matrices 684 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 685 internal::c_to_fortran_numbering(out); 686 } 687 }; 688 689 namespace internal { 690 691 template<typename _MatrixType, typename Rhs> 692 struct solve_retval<PastixBase<_MatrixType>, Rhs> 693 : solve_retval_base<PastixBase<_MatrixType>, Rhs> 694 { 695 typedef PastixBase<_MatrixType> Dec; 696 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 697 698 template<typename Dest> void evalTo(Dest& dst) const 699 { 700 dec()._solve(rhs(),dst); 701 } 702 }; 703 704 template<typename _MatrixType, typename Rhs> 705 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> 706 : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs> 707 { 708 typedef PastixBase<_MatrixType> Dec; 709 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 710 711 template<typename Dest> void evalTo(Dest& dst) const 712 { 713 this->defaultEvalTo(dst); 714 } 715 }; 716 717 } // end namespace internal 718 719 } // end namespace Eigen 720 721 #endif 722
