![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_PASTIXSUPPORT_H 00011 #define EIGEN_PASTIXSUPPORT_H 00012 00013 namespace Eigen { 00014 00015 #if defined(DCOMPLEX) 00016 #define PASTIX_COMPLEX COMPLEX 00017 #define PASTIX_DCOMPLEX DCOMPLEX 00018 #else 00019 #define PASTIX_COMPLEX std::complex<float> 00020 #define PASTIX_DCOMPLEX std::complex<double> 00021 #endif 00022 00031 template<typename _MatrixType, bool IsStrSym = false> class PastixLU; 00032 template<typename _MatrixType, int Options> class PastixLLT; 00033 template<typename _MatrixType, int Options> class PastixLDLT; 00034 00035 namespace internal 00036 { 00037 00038 template<class Pastix> struct pastix_traits; 00039 00040 template<typename _MatrixType> 00041 struct pastix_traits< PastixLU<_MatrixType> > 00042 { 00043 typedef _MatrixType MatrixType; 00044 typedef typename _MatrixType::Scalar Scalar; 00045 typedef typename _MatrixType::RealScalar RealScalar; 00046 typedef typename _MatrixType::StorageIndex StorageIndex; 00047 }; 00048 00049 template<typename _MatrixType, int Options> 00050 struct pastix_traits< PastixLLT<_MatrixType,Options> > 00051 { 00052 typedef _MatrixType MatrixType; 00053 typedef typename _MatrixType::Scalar Scalar; 00054 typedef typename _MatrixType::RealScalar RealScalar; 00055 typedef typename _MatrixType::StorageIndex StorageIndex; 00056 }; 00057 00058 template<typename _MatrixType, int Options> 00059 struct pastix_traits< PastixLDLT<_MatrixType,Options> > 00060 { 00061 typedef _MatrixType MatrixType; 00062 typedef typename _MatrixType::Scalar Scalar; 00063 typedef typename _MatrixType::RealScalar RealScalar; 00064 typedef typename _MatrixType::StorageIndex StorageIndex; 00065 }; 00066 00067 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) 00068 { 00069 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 00070 if (nbrhs == 0) {x = NULL; nbrhs=1;} 00071 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 00072 } 00073 00074 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) 00075 { 00076 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 00077 if (nbrhs == 0) {x = NULL; nbrhs=1;} 00078 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 00079 } 00080 00081 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) 00082 { 00083 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 00084 if (nbrhs == 0) {x = NULL; nbrhs=1;} 00085 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); 00086 } 00087 00088 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) 00089 { 00090 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } 00091 if (nbrhs == 0) {x = NULL; nbrhs=1;} 00092 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); 00093 } 00094 00095 // Convert the matrix to Fortran-style Numbering 00096 template <typename MatrixType> 00097 void c_to_fortran_numbering (MatrixType& mat) 00098 { 00099 if ( !(mat.outerIndexPtr()[0]) ) 00100 { 00101 int i; 00102 for(i = 0; i <= mat.rows(); ++i) 00103 ++mat.outerIndexPtr()[i]; 00104 for(i = 0; i < mat.nonZeros(); ++i) 00105 ++mat.innerIndexPtr()[i]; 00106 } 00107 } 00108 00109 // Convert to C-style Numbering 00110 template <typename MatrixType> 00111 void fortran_to_c_numbering (MatrixType& mat) 00112 { 00113 // Check the Numbering 00114 if ( mat.outerIndexPtr()[0] == 1 ) 00115 { // Convert to C-style numbering 00116 int i; 00117 for(i = 0; i <= mat.rows(); ++i) 00118 --mat.outerIndexPtr()[i]; 00119 for(i = 0; i < mat.nonZeros(); ++i) 00120 --mat.innerIndexPtr()[i]; 00121 } 00122 } 00123 } 00124 00125 // This is the base class to interface with PaStiX functions. 00126 // Users should not used this class directly. 00127 template <class Derived> 00128 class PastixBase : public SparseSolverBase<Derived> 00129 { 00130 protected: 00131 typedef SparseSolverBase<Derived> Base; 00132 using Base::derived; 00133 using Base::m_isInitialized; 00134 public: 00135 using Base::_solve_impl; 00136 00137 typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; 00138 typedef _MatrixType MatrixType; 00139 typedef typename MatrixType::Scalar Scalar; 00140 typedef typename MatrixType::RealScalar RealScalar; 00141 typedef typename MatrixType::StorageIndex StorageIndex; 00142 typedef Matrix<Scalar,Dynamic,1> Vector; 00143 typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix; 00144 enum { 00145 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00146 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00147 }; 00148 00149 public: 00150 00151 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0) 00152 { 00153 init(); 00154 } 00155 00156 ~PastixBase() 00157 { 00158 clean(); 00159 } 00160 00161 template<typename Rhs,typename Dest> 00162 bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; 00163 00169 Array<StorageIndex,IPARM_SIZE,1>& iparm() 00170 { 00171 return m_iparm; 00172 } 00173 00178 int& iparm(int idxparam) 00179 { 00180 return m_iparm(idxparam); 00181 } 00182 00187 Array<double,DPARM_SIZE,1>& dparm() 00188 { 00189 return m_dparm; 00190 } 00191 00192 00196 double& dparm(int idxparam) 00197 { 00198 return m_dparm(idxparam); 00199 } 00200 00201 inline Index cols() const { return m_size; } 00202 inline Index rows() const { return m_size; } 00203 00212 ComputationInfo info() const 00213 { 00214 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00215 return m_info; 00216 } 00217 00218 protected: 00219 00220 // Initialize the Pastix data structure, check the matrix 00221 void init(); 00222 00223 // Compute the ordering and the symbolic factorization 00224 void analyzePattern(ColSpMatrix& mat); 00225 00226 // Compute the numerical factorization 00227 void factorize(ColSpMatrix& mat); 00228 00229 // Free all the data allocated by Pastix 00230 void clean() 00231 { 00232 eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 00233 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; 00234 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; 00235 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 00236 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 00237 } 00238 00239 void compute(ColSpMatrix& mat); 00240 00241 int m_initisOk; 00242 int m_analysisIsOk; 00243 int m_factorizationIsOk; 00244 mutable ComputationInfo m_info; 00245 mutable pastix_data_t *m_pastixdata; // Data structure for pastix 00246 mutable int m_comm; // The MPI communicator identifier 00247 mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters 00248 mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters 00249 mutable Matrix<StorageIndex,Dynamic,1> m_perm; // Permutation vector 00250 mutable Matrix<StorageIndex,Dynamic,1> m_invp; // Inverse permutation vector 00251 mutable int m_size; // Size of the matrix 00252 }; 00253 00258 template <class Derived> 00259 void PastixBase<Derived>::init() 00260 { 00261 m_size = 0; 00262 m_iparm.setZero(IPARM_SIZE); 00263 m_dparm.setZero(DPARM_SIZE); 00264 00265 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; 00266 pastix(&m_pastixdata, MPI_COMM_WORLD, 00267 0, 0, 0, 0, 00268 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); 00269 00270 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; 00271 m_iparm[IPARM_VERBOSE] = API_VERBOSE_NOT; 00272 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; 00273 m_iparm[IPARM_INCOMPLETE] = API_NO; 00274 m_iparm[IPARM_OOC_LIMIT] = 2000; 00275 m_iparm[IPARM_RHS_MAKING] = API_RHS_B; 00276 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 00277 00278 m_iparm(IPARM_START_TASK) = API_TASK_INIT; 00279 m_iparm(IPARM_END_TASK) = API_TASK_INIT; 00280 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, 00281 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); 00282 00283 // Check the returned error 00284 if(m_iparm(IPARM_ERROR_NUMBER)) { 00285 m_info = InvalidInput; 00286 m_initisOk = false; 00287 } 00288 else { 00289 m_info = Success; 00290 m_initisOk = true; 00291 } 00292 } 00293 00294 template <class Derived> 00295 void PastixBase<Derived>::compute(ColSpMatrix& mat) 00296 { 00297 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); 00298 00299 analyzePattern(mat); 00300 factorize(mat); 00301 00302 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; 00303 } 00304 00305 00306 template <class Derived> 00307 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) 00308 { 00309 eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); 00310 00311 // clean previous calls 00312 if(m_size>0) 00313 clean(); 00314 00315 m_size = internal::convert_index<int>(mat.rows()); 00316 m_perm.resize(m_size); 00317 m_invp.resize(m_size); 00318 00319 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; 00320 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; 00321 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 00322 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 00323 00324 // Check the returned error 00325 if(m_iparm(IPARM_ERROR_NUMBER)) 00326 { 00327 m_info = NumericalIssue; 00328 m_analysisIsOk = false; 00329 } 00330 else 00331 { 00332 m_info = Success; 00333 m_analysisIsOk = true; 00334 } 00335 } 00336 00337 template <class Derived> 00338 void PastixBase<Derived>::factorize(ColSpMatrix& mat) 00339 { 00340 // if(&m_cpyMat != &mat) m_cpyMat = mat; 00341 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); 00342 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; 00343 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; 00344 m_size = internal::convert_index<int>(mat.rows()); 00345 00346 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), 00347 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); 00348 00349 // Check the returned error 00350 if(m_iparm(IPARM_ERROR_NUMBER)) 00351 { 00352 m_info = NumericalIssue; 00353 m_factorizationIsOk = false; 00354 m_isInitialized = false; 00355 } 00356 else 00357 { 00358 m_info = Success; 00359 m_factorizationIsOk = true; 00360 m_isInitialized = true; 00361 } 00362 } 00363 00364 /* Solve the system */ 00365 template<typename Base> 00366 template<typename Rhs,typename Dest> 00367 bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const 00368 { 00369 eigen_assert(m_isInitialized && "The matrix should be factorized first"); 00370 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, 00371 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 00372 int rhs = 1; 00373 00374 x = b; /* on return, x is overwritten by the computed solution */ 00375 00376 for (int i = 0; i < b.cols(); i++){ 00377 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; 00378 m_iparm[IPARM_END_TASK] = API_TASK_REFINE; 00379 00380 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0, 00381 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); 00382 } 00383 00384 // Check the returned error 00385 m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; 00386 00387 return m_iparm(IPARM_ERROR_NUMBER)==0; 00388 } 00389 00411 template<typename _MatrixType, bool IsStrSym> 00412 class PastixLU : public PastixBase< PastixLU<_MatrixType> > 00413 { 00414 public: 00415 typedef _MatrixType MatrixType; 00416 typedef PastixBase<PastixLU<MatrixType> > Base; 00417 typedef typename Base::ColSpMatrix ColSpMatrix; 00418 typedef typename MatrixType::StorageIndex StorageIndex; 00419 00420 public: 00421 PastixLU() : Base() 00422 { 00423 init(); 00424 } 00425 00426 explicit PastixLU(const MatrixType& matrix):Base() 00427 { 00428 init(); 00429 compute(matrix); 00430 } 00436 void compute (const MatrixType& matrix) 00437 { 00438 m_structureIsUptodate = false; 00439 ColSpMatrix temp; 00440 grabMatrix(matrix, temp); 00441 Base::compute(temp); 00442 } 00448 void analyzePattern(const MatrixType& matrix) 00449 { 00450 m_structureIsUptodate = false; 00451 ColSpMatrix temp; 00452 grabMatrix(matrix, temp); 00453 Base::analyzePattern(temp); 00454 } 00455 00461 void factorize(const MatrixType& matrix) 00462 { 00463 ColSpMatrix temp; 00464 grabMatrix(matrix, temp); 00465 Base::factorize(temp); 00466 } 00467 protected: 00468 00469 void init() 00470 { 00471 m_structureIsUptodate = false; 00472 m_iparm(IPARM_SYM) = API_SYM_NO; 00473 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; 00474 } 00475 00476 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 00477 { 00478 if(IsStrSym) 00479 out = matrix; 00480 else 00481 { 00482 if(!m_structureIsUptodate) 00483 { 00484 // update the transposed structure 00485 m_transposedStructure = matrix.transpose(); 00486 00487 // Set the elements of the matrix to zero 00488 for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 00489 for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) 00490 it.valueRef() = 0.0; 00491 00492 m_structureIsUptodate = true; 00493 } 00494 00495 out = m_transposedStructure + matrix; 00496 } 00497 internal::c_to_fortran_numbering(out); 00498 } 00499 00500 using Base::m_iparm; 00501 using Base::m_dparm; 00502 00503 ColSpMatrix m_transposedStructure; 00504 bool m_structureIsUptodate; 00505 }; 00506 00523 template<typename _MatrixType, int _UpLo> 00524 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > 00525 { 00526 public: 00527 typedef _MatrixType MatrixType; 00528 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; 00529 typedef typename Base::ColSpMatrix ColSpMatrix; 00530 00531 public: 00532 enum { UpLo = _UpLo }; 00533 PastixLLT() : Base() 00534 { 00535 init(); 00536 } 00537 00538 explicit PastixLLT(const MatrixType& matrix):Base() 00539 { 00540 init(); 00541 compute(matrix); 00542 } 00543 00547 void compute (const MatrixType& matrix) 00548 { 00549 ColSpMatrix temp; 00550 grabMatrix(matrix, temp); 00551 Base::compute(temp); 00552 } 00553 00558 void analyzePattern(const MatrixType& matrix) 00559 { 00560 ColSpMatrix temp; 00561 grabMatrix(matrix, temp); 00562 Base::analyzePattern(temp); 00563 } 00567 void factorize(const MatrixType& matrix) 00568 { 00569 ColSpMatrix temp; 00570 grabMatrix(matrix, temp); 00571 Base::factorize(temp); 00572 } 00573 protected: 00574 using Base::m_iparm; 00575 00576 void init() 00577 { 00578 m_iparm(IPARM_SYM) = API_SYM_YES; 00579 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; 00580 } 00581 00582 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 00583 { 00584 out.resize(matrix.rows(), matrix.cols()); 00585 // Pastix supports only lower, column-major matrices 00586 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 00587 internal::c_to_fortran_numbering(out); 00588 } 00589 }; 00590 00607 template<typename _MatrixType, int _UpLo> 00608 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > 00609 { 00610 public: 00611 typedef _MatrixType MatrixType; 00612 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 00613 typedef typename Base::ColSpMatrix ColSpMatrix; 00614 00615 public: 00616 enum { UpLo = _UpLo }; 00617 PastixLDLT():Base() 00618 { 00619 init(); 00620 } 00621 00622 explicit PastixLDLT(const MatrixType& matrix):Base() 00623 { 00624 init(); 00625 compute(matrix); 00626 } 00627 00631 void compute (const MatrixType& matrix) 00632 { 00633 ColSpMatrix temp; 00634 grabMatrix(matrix, temp); 00635 Base::compute(temp); 00636 } 00637 00642 void analyzePattern(const MatrixType& matrix) 00643 { 00644 ColSpMatrix temp; 00645 grabMatrix(matrix, temp); 00646 Base::analyzePattern(temp); 00647 } 00651 void factorize(const MatrixType& matrix) 00652 { 00653 ColSpMatrix temp; 00654 grabMatrix(matrix, temp); 00655 Base::factorize(temp); 00656 } 00657 00658 protected: 00659 using Base::m_iparm; 00660 00661 void init() 00662 { 00663 m_iparm(IPARM_SYM) = API_SYM_YES; 00664 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; 00665 } 00666 00667 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) 00668 { 00669 // Pastix supports only lower, column-major matrices 00670 out.resize(matrix.rows(), matrix.cols()); 00671 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); 00672 internal::c_to_fortran_numbering(out); 00673 } 00674 }; 00675 00676 } // end namespace Eigen 00677 00678 #endif