diff options
Diffstat (limited to 'eigen/Eigen/src/PaStiXSupport')
| -rw-r--r-- | eigen/Eigen/src/PaStiXSupport/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h | 151 |
2 files changed, 50 insertions, 107 deletions
diff --git a/eigen/Eigen/src/PaStiXSupport/CMakeLists.txt b/eigen/Eigen/src/PaStiXSupport/CMakeLists.txt deleted file mode 100644 index 28c657e..0000000 --- a/eigen/Eigen/src/PaStiXSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_PastixSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_PastixSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/PaStiXSupport COMPONENT Devel - ) diff --git a/eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h b/eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h index 20acc02..d2ebfd7 100644 --- a/eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -10,6 +10,8 @@ #ifndef EIGEN_PASTIXSUPPORT_H #define EIGEN_PASTIXSUPPORT_H +namespace Eigen { + #if defined(DCOMPLEX) #define PASTIX_COMPLEX COMPLEX #define PASTIX_DCOMPLEX DCOMPLEX @@ -18,8 +20,6 @@ #define PASTIX_DCOMPLEX std::complex<double> #endif -namespace Eigen { - /** \ingroup PaStiXSupport_Module * \brief Interface to the PaStix solver * @@ -43,7 +43,7 @@ namespace internal typedef _MatrixType MatrixType; typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; - typedef typename _MatrixType::Index Index; + typedef typename _MatrixType::StorageIndex StorageIndex; }; template<typename _MatrixType, int Options> @@ -52,7 +52,7 @@ namespace internal typedef _MatrixType MatrixType; typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; - typedef typename _MatrixType::Index Index; + typedef typename _MatrixType::StorageIndex StorageIndex; }; template<typename _MatrixType, int Options> @@ -61,7 +61,7 @@ namespace internal typedef _MatrixType MatrixType; typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; - typedef typename _MatrixType::Index Index; + typedef typename _MatrixType::StorageIndex StorageIndex; }; 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) @@ -82,14 +82,14 @@ namespace internal { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); + c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); } 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) { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); + z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); } // Convert the matrix to Fortran-style Numbering @@ -125,20 +125,30 @@ namespace internal // This is the base class to interface with PaStiX functions. // Users should not used this class directly. template <class Derived> -class PastixBase : internal::noncopyable +class PastixBase : public SparseSolverBase<Derived> { + protected: + typedef SparseSolverBase<Derived> Base; + using Base::derived; + using Base::m_isInitialized; public: + using Base::_solve_impl; + typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; + typedef typename MatrixType::StorageIndex StorageIndex; typedef Matrix<Scalar,Dynamic,1> Vector; typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; public: - PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0) + PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0) { init(); } @@ -147,39 +157,16 @@ class PastixBase : internal::noncopyable { clean(); } - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. - * - * \sa compute() - */ - template<typename Rhs> - inline const internal::solve_retval<PastixBase, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "Pastix solver is not initialized."); - eigen_assert(rows()==b.rows() - && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<PastixBase, Rhs>(*this, b.derived()); - } template<typename Rhs,typename Dest> - bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; + bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; - Derived& derived() - { - return *static_cast<Derived*>(this); - } - const Derived& derived() const - { - return *static_cast<const Derived*>(this); - } - /** Returns a reference to the integer vector IPARM of PaStiX parameters * to modify the default parameters. * The statistics related to the different phases of factorization and solve are saved here as well * \sa analyzePattern() factorize() */ - Array<Index,IPARM_SIZE,1>& iparm() + Array<StorageIndex,IPARM_SIZE,1>& iparm() { return m_iparm; } @@ -197,7 +184,7 @@ class PastixBase : internal::noncopyable * The statistics related to the different phases of factorization and solve are saved here as well * \sa analyzePattern() factorize() */ - Array<RealScalar,IPARM_SIZE,1>& dparm() + Array<double,DPARM_SIZE,1>& dparm() { return m_dparm; } @@ -228,20 +215,6 @@ class PastixBase : internal::noncopyable return m_info; } - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. - * - * \sa compute() - */ - template<typename Rhs> - inline const internal::sparse_solve_retval<PastixBase, Rhs> - solve(const SparseMatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized."); - eigen_assert(rows()==b.rows() - && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); - } - protected: // Initialize the Pastix data structure, check the matrix @@ -268,14 +241,13 @@ class PastixBase : internal::noncopyable int m_initisOk; int m_analysisIsOk; int m_factorizationIsOk; - bool m_isInitialized; mutable ComputationInfo m_info; mutable pastix_data_t *m_pastixdata; // Data structure for pastix mutable int m_comm; // The MPI communicator identifier - mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters - mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters - mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector - mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector + mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters + mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters + mutable Matrix<StorageIndex,Dynamic,1> m_perm; // Permutation vector + mutable Matrix<StorageIndex,Dynamic,1> m_invp; // Inverse permutation vector mutable int m_size; // Size of the matrix }; @@ -296,7 +268,7 @@ void PastixBase<Derived>::init() 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; - m_iparm[IPARM_VERBOSE] = 2; + m_iparm[IPARM_VERBOSE] = API_VERBOSE_NOT; m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; m_iparm[IPARM_INCOMPLETE] = API_NO; m_iparm[IPARM_OOC_LIMIT] = 2000; @@ -328,7 +300,6 @@ void PastixBase<Derived>::compute(ColSpMatrix& mat) factorize(mat); m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; - m_isInitialized = m_factorizationIsOk; } @@ -341,7 +312,7 @@ void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) if(m_size>0) clean(); - m_size = mat.rows(); + m_size = internal::convert_index<int>(mat.rows()); m_perm.resize(m_size); m_invp.resize(m_size); @@ -370,7 +341,7 @@ void PastixBase<Derived>::factorize(ColSpMatrix& mat) eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase"); m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT; m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT; - m_size = mat.rows(); + m_size = internal::convert_index<int>(mat.rows()); internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); @@ -393,7 +364,7 @@ void PastixBase<Derived>::factorize(ColSpMatrix& mat) /* Solve the system */ template<typename Base> template<typename Rhs,typename Dest> -bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const +bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const { eigen_assert(m_isInitialized && "The matrix should be factorized first"); EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, @@ -406,7 +377,7 @@ bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) co m_iparm[IPARM_START_TASK] = API_TASK_SOLVE; m_iparm[IPARM_END_TASK] = API_TASK_REFINE; - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0, + internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0, m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); } @@ -431,8 +402,10 @@ bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) co * NOTE : Note that if the analysis and factorization phase are called separately, * the input matrix will be symmetrized at each call, hence it is advised to * symmetrize the matrix in a end-user program and set \p IsStrSym to true - * - * \sa \ref TutorialSparseDirectSolvers + * + * \implsparsesolverconcept + * + * \sa \ref TutorialSparseSolverConcept, class SparseLU * */ template<typename _MatrixType, bool IsStrSym> @@ -442,7 +415,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > typedef _MatrixType MatrixType; typedef PastixBase<PastixLU<MatrixType> > Base; typedef typename Base::ColSpMatrix ColSpMatrix; - typedef typename MatrixType::Index Index; + typedef typename MatrixType::StorageIndex StorageIndex; public: PastixLU() : Base() @@ -450,7 +423,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > init(); } - PastixLU(const MatrixType& matrix):Base() + explicit PastixLU(const MatrixType& matrix):Base() { init(); compute(matrix); @@ -542,8 +515,10 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > * * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX - * - * \sa \ref TutorialSparseDirectSolvers + * + * \implsparsesolverconcept + * + * \sa \ref TutorialSparseSolverConcept, class SimplicialLLT */ template<typename _MatrixType, int _UpLo> class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > @@ -560,7 +535,7 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > init(); } - PastixLLT(const MatrixType& matrix):Base() + explicit PastixLLT(const MatrixType& matrix):Base() { init(); compute(matrix); @@ -606,6 +581,7 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) { + out.resize(matrix.rows(), matrix.cols()); // Pastix supports only lower, column-major matrices out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); internal::c_to_fortran_numbering(out); @@ -623,8 +599,10 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > * * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX - * - * \sa \ref TutorialSparseDirectSolvers + * + * \implsparsesolverconcept + * + * \sa \ref TutorialSparseSolverConcept, class SimplicialLDLT */ template<typename _MatrixType, int _UpLo> class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > @@ -641,7 +619,7 @@ class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > init(); } - PastixLDLT(const MatrixType& matrix):Base() + explicit PastixLDLT(const MatrixType& matrix):Base() { init(); compute(matrix); @@ -689,41 +667,12 @@ class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) { // Pastix supports only lower, column-major matrices + out.resize(matrix.rows(), matrix.cols()); out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); internal::c_to_fortran_numbering(out); } }; -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<PastixBase<_MatrixType>, Rhs> - : solve_retval_base<PastixBase<_MatrixType>, Rhs> -{ - typedef PastixBase<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -template<typename _MatrixType, typename Rhs> -struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> - : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs> -{ - typedef PastixBase<_MatrixType> Dec; - EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - this->defaultEvalTo(dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif |