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