diff options
Diffstat (limited to 'Eigen/src/PaStiXSupport')
-rw-r--r-- | Eigen/src/PaStiXSupport/PaStiXSupport.h | 394 |
1 files changed, 174 insertions, 220 deletions
diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index f42826208..f687a5e6b 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -35,7 +35,6 @@ namespace Eigen { * * \sa TutorialSparseDirectSolvers */ - template<typename _MatrixType, bool IsStrSym = false> class PastixLU; template<typename _MatrixType, int Options> class PastixLLT; template<typename _MatrixType, int Options> class PastixLDLT; @@ -75,32 +74,34 @@ namespace internal 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) { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } - if (nbrhs == 0) x = NULL; + if (nbrhs == 0) {x = NULL; nbrhs=1;} s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); } 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) { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } - if (nbrhs == 0) x = NULL; + if (nbrhs == 0) {x = NULL; nbrhs=1;} d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); } 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) { + 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<COMPLEX*>(vals), perm, invp, reinterpret_cast<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; + if (nbrhs == 0) {x = NULL; nbrhs=1;} z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); } // Convert the matrix to Fortran-style Numbering template <typename MatrixType> - void EigenToFortranNumbering (MatrixType& mat) + void c_to_fortran_numbering (MatrixType& mat) { if ( !(mat.outerIndexPtr()[0]) ) { @@ -114,7 +115,7 @@ namespace internal // Convert to C-style Numbering template <typename MatrixType> - void EigenToCNumbering (MatrixType& mat) + void fortran_to_c_numbering (MatrixType& mat) { // Check the Numbering if ( mat.outerIndexPtr()[0] == 1 ) @@ -126,38 +127,12 @@ namespace internal --mat.innerIndexPtr()[i]; } } - - // Symmetrize the graph of the input matrix - // In : The Input matrix to symmetrize the pattern - // Out : The output matrix - // StrMatTrans : The structural pattern of the transpose of In; It is - // used to optimize the future symmetrization with the same matrix pattern - // WARNING It is assumed here that successive calls to this routine are done - // with matrices having the same pattern. - template <typename MatrixType> - void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose) - { - eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix"); - if (!hasTranspose) - { //First call to this routine, need to compute the structural pattern of In^T - StrMatTrans = In.transpose(); - // Set the elements of the matrix to zero - for (int i = 0; i < StrMatTrans.rows(); i++) - { - for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it) - it.valueRef() = 0.0; - } - hasTranspose = true; - } - Out = (StrMatTrans + In).eval(); - } - } // This is the base class to interface with PaStiX functions. // Users should not used this class directly. template <class Derived> -class PastixBase +class PastixBase : internal::noncopyable { public: typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; @@ -166,29 +141,19 @@ class PastixBase typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; typedef Matrix<Scalar,Dynamic,1> Vector; + typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix; public: - PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) + PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0) { - m_pastixdata = 0; - m_hasTranspose = false; - PastixInit(); + init(); } ~PastixBase() { - PastixDestroy(); + clean(); } - - // Initialize the Pastix data structure, check the matrix - void PastixInit(); - - // Compute the ordering and the symbolic factorization - Derived& analyzePattern (MatrixType& mat); - - // Compute the numerical factorization - Derived& factorize (MatrixType& mat); /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * @@ -269,7 +234,6 @@ class PastixBase /** Return a reference to a particular index parameter of the DPARM vector * \sa dparm() */ - double& dparm(int idxparam) { return m_dparm(idxparam); @@ -307,17 +271,27 @@ class PastixBase } protected: + + // Initialize the Pastix data structure, check the matrix + void init(); + + // Compute the ordering and the symbolic factorization + void analyzePattern(ColSpMatrix& mat); + + // Compute the numerical factorization + void factorize(ColSpMatrix& mat); + // Free all the data allocated by Pastix - void PastixDestroy() + void clean() { eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); m_iparm(IPARM_START_TASK) = API_TASK_CLEAN; m_iparm(IPARM_END_TASK) = API_TASK_CLEAN; - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(), - m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data()); + internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, + m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); } - Derived& compute (MatrixType& mat); + void compute(ColSpMatrix& mat); int m_initisOk; int m_analysisIsOk; @@ -325,22 +299,12 @@ class PastixBase bool m_isInitialized; mutable ComputationInfo m_info; mutable pastix_data_t *m_pastixdata; // Data structure for pastix - mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input null matrix - mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector - mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix - mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed mutable int m_comm; // The MPI communicator identifier - mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters + 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 int m_ordering; // ordering method to use - mutable int m_amalgamation; // level of amalgamation mutable int m_size; // Size of the matrix - - private: - PastixBase(PastixBase& ) {} - }; /** Initialize the PaStiX data structure. @@ -348,29 +312,29 @@ class PastixBase * \sa iparm() dparm() */ template <class Derived> -void PastixBase<Derived>::PastixInit() +void PastixBase<Derived>::init() { m_size = 0; - m_iparm.resize(IPARM_SIZE); - m_dparm.resize(DPARM_SIZE); + m_iparm.setZero(IPARM_SIZE); + m_dparm.setZero(DPARM_SIZE); m_iparm(IPARM_MODIFY_PARAMETER) = API_NO; - if(m_pastixdata) - { // This trick is used to reset the Pastix internal data between successive - // calls with (structural) different matrices - PastixDestroy(); - m_pastixdata = 0; - m_iparm(IPARM_MODIFY_PARAMETER) = API_YES; - m_hasTranspose = false; - } + pastix(&m_pastixdata, MPI_COMM_WORLD, + 0, 0, 0, 0, + 0, 0, 0, 1, m_iparm.data(), m_dparm.data()); - m_iparm(IPARM_START_TASK) = API_TASK_INIT; - m_iparm(IPARM_END_TASK) = API_TASK_INIT; + m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO; + m_iparm[IPARM_VERBOSE] = 2; + m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH; + m_iparm[IPARM_INCOMPLETE] = API_NO; + m_iparm[IPARM_OOC_LIMIT] = 2000; + m_iparm[IPARM_RHS_MAKING] = API_RHS_B; m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(), - m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data()); - m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; + m_iparm(IPARM_START_TASK) = API_TASK_INIT; + m_iparm(IPARM_END_TASK) = API_TASK_INIT; + internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0, + 0, 0, 0, 0, m_iparm.data(), m_dparm.data()); // Check the returned error if(m_iparm(IPARM_ERROR_NUMBER)) { @@ -384,82 +348,74 @@ void PastixBase<Derived>::PastixInit() } template <class Derived> -Derived& PastixBase<Derived>::compute(MatrixType& mat) +void PastixBase<Derived>::compute(ColSpMatrix& mat) { eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared"); - typedef typename MatrixType::Scalar Scalar; - // Save the size of the current matrix - m_size = mat.rows(); - // Convert the matrix in fortran-style numbering - internal::EigenToFortranNumbering(mat); - analyzePattern(mat); + analyzePattern(mat); factorize(mat); - m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; - if (m_factorizationIsOk) m_isInitialized = true; - //Convert back the matrix -- Is it really necessary here - internal::EigenToCNumbering(mat); - - return derived(); + m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO; + m_isInitialized = m_factorizationIsOk; } template <class Derived> -Derived& PastixBase<Derived>::analyzePattern(MatrixType& mat) -{ - eigen_assert(m_initisOk && "PastixInit should be called first to set the default parameters"); +void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat) +{ + eigen_assert(m_initisOk && "The initialization of PaSTiX failed"); + + // clean previous calls + if(m_size>0) + clean(); + m_size = mat.rows(); m_perm.resize(m_size); m_invp.resize(m_size); - // Convert the matrix in fortran-style numbering - internal::EigenToFortranNumbering(mat); - m_iparm(IPARM_START_TASK) = API_TASK_ORDERING; m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE; - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), - mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data()); + mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); // Check the returned error - if(m_iparm(IPARM_ERROR_NUMBER)) { + if(m_iparm(IPARM_ERROR_NUMBER)) + { m_info = NumericalIssue; m_analysisIsOk = false; } - else { + else + { m_info = Success; m_analysisIsOk = true; } - return derived(); } template <class Derived> -Derived& PastixBase<Derived>::factorize(MatrixType& mat) +void PastixBase<Derived>::factorize(ColSpMatrix& mat) { +// if(&m_cpyMat != &mat) m_cpyMat = 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(); - // Convert the matrix in fortran-style numbering - internal::EigenToFortranNumbering(mat); - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(), - mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data()); + mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data()); // Check the returned error - if(m_iparm(IPARM_ERROR_NUMBER)) { + if(m_iparm(IPARM_ERROR_NUMBER)) + { m_info = NumericalIssue; m_factorizationIsOk = false; m_isInitialized = false; } - else { + else + { m_info = Success; m_factorizationIsOk = true; m_isInitialized = true; } - return derived(); } /* Solve the system */ @@ -475,20 +431,17 @@ bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) co x = b; /* on return, x is overwritten by the computed solution */ for (int i = 0; i < b.cols(); i++){ - m_iparm(IPARM_START_TASK) = API_TASK_SOLVE; - m_iparm(IPARM_END_TASK) = API_TASK_REFINE; - m_iparm(IPARM_RHS_MAKING) = API_RHS_B; - internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(), - m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); + 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, + m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data()); } + // Check the returned error - if(m_iparm(IPARM_ERROR_NUMBER)) { - m_info = NumericalIssue; - return false; - } - else { - return true; - } + m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue; + + return m_iparm(IPARM_ERROR_NUMBER)==0; } /** \ingroup PaStiXSupport_Module @@ -516,14 +469,18 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > public: typedef _MatrixType MatrixType; typedef PastixBase<PastixLU<MatrixType> > Base; - typedef typename MatrixType::Scalar Scalar; - typedef SparseMatrix<Scalar, ColMajor> PaStiXType; + typedef typename Base::ColSpMatrix ColSpMatrix; + typedef typename MatrixType::Index Index; public: - PastixLU():Base() {} + PastixLU() : Base() + { + init(); + } PastixLU(const MatrixType& matrix):Base() { + init(); compute(matrix); } /** Compute the LU supernodal factorization of \p matrix. @@ -533,18 +490,9 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > */ void compute (const MatrixType& matrix) { - // Pastix supports only column-major matrices with a symmetric pattern - Base::PastixInit(); - PaStiXType temp(matrix.rows(), matrix.cols()); - // Symmetrize the graph of the matrix - if (IsStrSym) - temp = matrix; - else - { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose); - } - m_iparm[IPARM_SYM] = API_SYM_NO; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; + m_structureIsUptodate = false; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::compute(temp); } /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. @@ -554,20 +502,9 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > */ void analyzePattern(const MatrixType& matrix) { - - Base::PastixInit(); - /* Pastix supports only column-major matrices with symmetrized patterns */ - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - // Symmetrize the graph of the matrix - if (IsStrSym) - temp = matrix; - else - { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose); - } - - m_iparm(IPARM_SYM) = API_SYM_NO; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; + m_structureIsUptodate = false; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::analyzePattern(temp); } @@ -578,27 +515,48 @@ class PastixLU : public PastixBase< PastixLU<_MatrixType> > */ void factorize(const MatrixType& matrix) { - /* Pastix supports only column-major matrices with symmetrized patterns */ - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - // Symmetrize the graph of the matrix - if (IsStrSym) - temp = matrix; - else - { - internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose); - } - m_iparm(IPARM_SYM) = API_SYM_NO; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::factorize(temp); } protected: + + void init() + { + m_structureIsUptodate = false; + m_iparm(IPARM_SYM) = API_SYM_NO; + m_iparm(IPARM_FACTORIZATION) = API_FACT_LU; + } + + void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) + { + if(IsStrSym) + out = matrix; + else + { + if(!m_structureIsUptodate) + { + // update the transposed structure + m_transposedStructure = matrix.transpose(); + + // Set the elements of the matrix to zero + for (Index j=0; j<m_transposedStructure.outerSize(); ++j) + for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it) + it.valueRef() = 0.0; + + m_structureIsUptodate = true; + } + + out = m_transposedStructure + matrix; + } + internal::c_to_fortran_numbering(out); + } + using Base::m_iparm; using Base::m_dparm; - using Base::m_StrMatTrans; - using Base::m_hasTranspose; - private: - PastixLU(PastixLU& ) {} + ColSpMatrix m_transposedStructure; + bool m_structureIsUptodate; }; /** \ingroup PaStiXSupport_Module @@ -621,15 +579,18 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > public: typedef _MatrixType MatrixType; typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef typename Base::ColSpMatrix ColSpMatrix; public: enum { UpLo = _UpLo }; - PastixLLT():Base() {} + PastixLLT() : Base() + { + init(); + } PastixLLT(const MatrixType& matrix):Base() { + init(); compute(matrix); } @@ -638,13 +599,8 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > */ void compute (const MatrixType& matrix) { - // Pastix supports only lower, column-major matrices - Base::PastixInit(); // This is necessary to let PaStiX initialize its data structure between successive calls to compute - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::compute(temp); } @@ -654,13 +610,8 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > */ void analyzePattern(const MatrixType& matrix) { - Base::PastixInit(); - // Pastix supports only lower, column-major matrices - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::analyzePattern(temp); } /** Compute the LL^T supernodal numerical factorization of \p matrix @@ -668,19 +619,25 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > */ void factorize(const MatrixType& matrix) { - // Pastix supports only lower, column-major matrices - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::factorize(temp); } protected: using Base::m_iparm; - private: - PastixLLT(PastixLLT& ) {} + void init() + { + m_iparm(IPARM_SYM) = API_SYM_YES; + m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT; + } + + void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) + { + // Pastix supports only lower, column-major matrices + out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); + internal::c_to_fortran_numbering(out); + } }; /** \ingroup PaStiXSupport_Module @@ -700,18 +657,21 @@ class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> > template<typename _MatrixType, int _UpLo> class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > { -public: + public: typedef _MatrixType MatrixType; typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef typename Base::ColSpMatrix ColSpMatrix; public: enum { UpLo = _UpLo }; - PastixLDLT():Base() {} + PastixLDLT():Base() + { + init(); + } PastixLDLT(const MatrixType& matrix):Base() { + init(); compute(matrix); } @@ -720,13 +680,8 @@ public: */ void compute (const MatrixType& matrix) { - Base::PastixInit(); - // Pastix supports only lower, column-major matrices - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::compute(temp); } @@ -736,14 +691,8 @@ public: */ void analyzePattern(const MatrixType& matrix) { - Base::PastixInit(); - // Pastix supports only lower, column-major matrices - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::analyzePattern(temp); } /** Compute the LDL^T supernodal numerical factorization of \p matrix @@ -751,21 +700,26 @@ public: */ void factorize(const MatrixType& matrix) { - // Pastix supports only lower, column-major matrices - SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols()); - PermutationMatrix<Dynamic,Dynamic,Index> pnull; - temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull); - - m_iparm(IPARM_SYM) = API_SYM_YES; - m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; + ColSpMatrix temp; + grabMatrix(matrix, temp); Base::factorize(temp); } protected: using Base::m_iparm; - private: - PastixLDLT(PastixLDLT& ) {} + void init() + { + m_iparm(IPARM_SYM) = API_SYM_YES; + m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT; + } + + void grabMatrix(const MatrixType& matrix, ColSpMatrix& out) + { + // Pastix supports only lower, column-major matrices + out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>(); + internal::c_to_fortran_numbering(out); + } }; namespace internal { |