From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- eigen/Eigen/src/PaStiXSupport/CMakeLists.txt | 6 - eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h | 151 +++++++++----------------- 2 files changed, 50 insertions(+), 107 deletions(-) delete mode 100644 eigen/Eigen/src/PaStiXSupport/CMakeLists.txt (limited to 'eigen/Eigen/src/PaStiXSupport') 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 #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 @@ -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 @@ -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(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); + c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); } void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex *vals, int *perm, int * invp, std::complex *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(vals), perm, invp, reinterpret_cast(x), nbrhs, iparm, dparm); + z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast(vals), perm, invp, reinterpret_cast(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 PastixBase : internal::noncopyable +class PastixBase : public SparseSolverBase { + protected: + typedef SparseSolverBase Base; + using Base::derived; + using Base::m_isInitialized; public: + using Base::_solve_impl; + typedef typename internal::pastix_traits::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 Vector; typedef SparseMatrix 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 - inline const internal::solve_retval - solve(const MatrixBase& 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(*this, b.derived()); - } template - bool _solve (const MatrixBase &b, MatrixBase &x) const; + bool _solve_impl(const MatrixBase &b, MatrixBase &x) const; - Derived& derived() - { - return *static_cast(this); - } - const Derived& derived() const - { - return *static_cast(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& iparm() + Array& 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& dparm() + Array& 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 - inline const internal::sparse_solve_retval - solve(const SparseMatrixBase& 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(*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 m_iparm; // integer vector for the input parameters - mutable Matrix m_dparm; // Scalar vector for the input parameters - mutable Matrix m_perm; // Permutation vector - mutable Matrix m_invp; // Inverse permutation vector + mutable Array m_iparm; // integer vector for the input parameters + mutable Array m_dparm; // Scalar vector for the input parameters + mutable Matrix m_perm; // Permutation vector + mutable Matrix m_invp; // Inverse permutation vector mutable int m_size; // Size of the matrix }; @@ -296,7 +268,7 @@ void PastixBase::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::compute(ColSpMatrix& mat) factorize(mat); m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; - m_isInitialized = m_factorizationIsOk; } @@ -341,7 +312,7 @@ void PastixBase::analyzePattern(ColSpMatrix& mat) if(m_size>0) clean(); - m_size = mat.rows(); + m_size = internal::convert_index(mat.rows()); m_perm.resize(m_size); m_invp.resize(m_size); @@ -370,7 +341,7 @@ void PastixBase::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(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::factorize(ColSpMatrix& mat) /* Solve the system */ template template -bool PastixBase::_solve (const MatrixBase &b, MatrixBase &x) const +bool PastixBase::_solve_impl(const MatrixBase &b, MatrixBase &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::_solve (const MatrixBase &b, MatrixBase &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(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::_solve (const MatrixBase &b, MatrixBase &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 @@ -442,7 +415,7 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > typedef _MatrixType MatrixType; typedef PastixBase > 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 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() = matrix.template selfadjointView(); 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 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() = matrix.template selfadjointView(); internal::c_to_fortran_numbering(out); } }; -namespace internal { - -template -struct solve_retval, Rhs> - : solve_retval_base, Rhs> -{ - typedef PastixBase<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -template -struct sparse_solve_retval, Rhs> - : sparse_solve_retval_base, Rhs> -{ - typedef PastixBase<_MatrixType> Dec; - EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) - - template void evalTo(Dest& dst) const - { - this->defaultEvalTo(dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif -- cgit v1.2.3