diff options
Diffstat (limited to 'Eigen/src/PardisoSupport/PardisoSupport.h')
-rw-r--r-- | Eigen/src/PardisoSupport/PardisoSupport.h | 98 |
1 files changed, 19 insertions, 79 deletions
diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h index b6571069e..054af6635 100644 --- a/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -96,10 +96,17 @@ namespace internal } template<class Derived> -class PardisoImpl : internal::noncopyable +class PardisoImpl : public SparseSolveBase<PardisoImpl<Derived> { + protected: + typedef SparseSolveBase<PardisoImpl<Derived> Base; + using Base::derived; + using Base::m_isInitialized; + typedef internal::pardiso_traits<Derived> Traits; public: + using base::_solve_impl; + typedef typename Traits::MatrixType MatrixType; typedef typename Traits::Scalar Scalar; typedef typename Traits::RealScalar RealScalar; @@ -118,7 +125,7 @@ class PardisoImpl : internal::noncopyable eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type"); m_iparm.setZero(); m_msglvl = 0; // No output - m_initialized = false; + m_isInitialized = false; } ~PardisoImpl() @@ -136,7 +143,7 @@ class PardisoImpl : internal::noncopyable */ ComputationInfo info() const { - eigen_assert(m_initialized && "Decomposition is not initialized."); + eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; } @@ -165,51 +172,14 @@ class PardisoImpl : internal::noncopyable Derived& factorize(const MatrixType& matrix); Derived& compute(const MatrixType& matrix); - - /** \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<PardisoImpl, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_initialized && "Pardiso solver is not initialized."); - eigen_assert(rows()==b.rows() - && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived()); - } - - /** \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<PardisoImpl, Rhs> - solve(const SparseMatrixBase<Rhs>& b) const - { - eigen_assert(m_initialized && "Pardiso solver is not initialized."); - eigen_assert(rows()==b.rows() - && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived()); - } - - Derived& derived() - { - return *static_cast<Derived*>(this); - } - const Derived& derived() const - { - return *static_cast<const Derived*>(this); - } template<typename BDerived, typename XDerived> - bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; + bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; protected: void pardisoRelease() { - if(m_initialized) // Factorization ran at least once + if(m_isInitialized) // Factorization ran at least once { internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0, m_iparm.data(), m_msglvl, 0, 0); @@ -270,7 +240,7 @@ class PardisoImpl : internal::noncopyable mutable SparseMatrixType m_matrix; ComputationInfo m_info; - bool m_initialized, m_analysisIsOk, m_factorizationIsOk; + bool m_analysisIsOk, m_factorizationIsOk; Index m_type, m_msglvl; mutable void *m_pt[64]; mutable ParameterType m_iparm; @@ -298,7 +268,7 @@ Derived& PardisoImpl<Derived>::compute(const MatrixType& a) manageErrorCode(error); m_analysisIsOk = true; m_factorizationIsOk = true; - m_initialized = true; + m_isInitialized = true; return derived(); } @@ -321,7 +291,7 @@ Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) manageErrorCode(error); m_analysisIsOk = true; m_factorizationIsOk = false; - m_initialized = true; + m_isInitialized = true; return derived(); } @@ -345,7 +315,7 @@ Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) template<class Base> template<typename BDerived,typename XDerived> -bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const +bool PardisoImpl<Base>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const { if(m_iparm[0] == 0) // Factorization was not computed return false; @@ -421,7 +391,7 @@ class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> > pardisoInit(Base::ScalarIsComplex ? 13 : 11); } - PardisoLU(const MatrixType& matrix) + explicit PardisoLU(const MatrixType& matrix) : Base() { pardisoInit(Base::ScalarIsComplex ? 13 : 11); @@ -472,7 +442,7 @@ class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> > pardisoInit(Base::ScalarIsComplex ? 4 : 2); } - PardisoLLT(const MatrixType& matrix) + explicit PardisoLLT(const MatrixType& matrix) : Base() { pardisoInit(Base::ScalarIsComplex ? 4 : 2); @@ -530,7 +500,7 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); } - PardisoLDLT(const MatrixType& matrix) + explicit PardisoLDLT(const MatrixType& matrix) : Base() { pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2); @@ -546,36 +516,6 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > } }; -namespace internal { - -template<typename _Derived, typename Rhs> -struct solve_retval<PardisoImpl<_Derived>, Rhs> - : solve_retval_base<PardisoImpl<_Derived>, Rhs> -{ - typedef PardisoImpl<_Derived> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -template<typename Derived, typename Rhs> -struct sparse_solve_retval<PardisoImpl<Derived>, Rhs> - : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs> -{ - typedef PardisoImpl<Derived> 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 // EIGEN_PARDISOSUPPORT_H |