diff options
Diffstat (limited to 'unsupported/Eigen/src/IterativeSolvers')
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/DGMRES.h | 39 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 41 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h | 34 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h | 39 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/MINRES.h | 46 |
5 files changed, 29 insertions, 170 deletions
diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index 9fcc8a8d9..0e1b7d977 100644 --- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -108,6 +108,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > using Base::m_isInitialized; using Base::m_tolerance; public: + using Base::_solve_impl; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -138,25 +139,9 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > ~DGMRES() {} - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template<typename Rhs,typename Guess> - inline const internal::solve_retval_with_guess<DGMRES, Rhs, Guess> - solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "DGMRES is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "DGMRES::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - <DGMRES, Rhs, Guess>(*this, b.derived(), x0); - } - /** \internal */ template<typename Rhs,typename Dest> - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { bool failed = false; for(int j=0; j<b.cols(); ++j) @@ -175,10 +160,10 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > /** \internal */ template<typename Rhs,typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const { x = b; - _solveWithGuess(b,x); + _solve_with_guess_impl(b,x.derived()); } /** * Get the restart value @@ -522,21 +507,5 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, return 0; } -namespace internal { - - template<typename _MatrixType, typename _Preconditioner, typename Rhs> -struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs> - : solve_retval_base<DGMRES<_MatrixType, _Preconditioner>, Rhs> -{ - typedef DGMRES<_MatrixType, _Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; -} // end namespace internal - } // end namespace Eigen #endif diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 67498705b..cd15ce0bf 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -281,6 +281,7 @@ private: int m_restart; public: + using Base::_solve_impl; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -315,25 +316,9 @@ public: */ void set_restart(const int restart) { m_restart=restart; } - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template<typename Rhs,typename Guess> - inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess> - solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "GMRES is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "GMRES::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - <GMRES, Rhs, Guess>(*this, b.derived(), x0); - } - /** \internal */ template<typename Rhs,typename Dest> - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { bool failed = false; for(int j=0; j<b.cols(); ++j) @@ -353,35 +338,17 @@ public: /** \internal */ template<typename Rhs,typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const { x = b; if(x.squaredNorm() == 0) return; // Check Zero right hand side - _solveWithGuess(b,x); + _solve_with_guess_impl(b,x.derived()); } protected: }; - -namespace internal { - - template<typename _MatrixType, typename _Preconditioner, typename Rhs> -struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs> - : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs> -{ - typedef GMRES<_MatrixType, _Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_GMRES_H diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h index 661c1f2e0..dd43de6b3 100644 --- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h +++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h @@ -27,8 +27,11 @@ namespace Eigen { */ template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> > -class IncompleteCholesky : internal::noncopyable +class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > { + protected: + typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base; + using Base::m_isInitialized; public: typedef SparseMatrix<Scalar,ColMajor> MatrixType; typedef _OrderingType OrderingType; @@ -89,7 +92,7 @@ class IncompleteCholesky : internal::noncopyable } template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { eigen_assert(m_factorizationIsOk && "factorize() should be called first"); if (m_perm.rows() == b.rows()) @@ -103,22 +106,13 @@ class IncompleteCholesky : internal::noncopyable x = m_perm * x; x = m_scal.asDiagonal() * x; } - template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed"); - eigen_assert(m_isInitialized && "IncompleteLLT is not initialized."); - eigen_assert(cols()==b.rows() - && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived()); - } + protected: SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC ScalarType m_scal; // The vector for scaling the matrix Scalar m_shift; //The initial shift parameter bool m_analysisIsOk; bool m_factorizationIsOk; - bool m_isInitialized; ComputationInfo m_info; PermutationType m_perm; @@ -256,22 +250,6 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const Idx listCol[rowIdx(jk)].push_back(col); } } -namespace internal { - -template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs> -struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> - : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs> -{ - typedef IncompleteCholesky<_Scalar, _UpLo, OrderingType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal } // end namespace Eigen diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h index 67e780181..7d08c3515 100644 --- a/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h +++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h @@ -13,8 +13,12 @@ namespace Eigen { template <typename _Scalar> -class IncompleteLU +class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> > { + protected: + typedef SparseSolverBase<IncompleteLU<_Scalar> > Base; + using Base::m_isInitialized; + typedef _Scalar Scalar; typedef Matrix<Scalar,Dynamic,1> Vector; typedef typename Vector::Index Index; @@ -23,10 +27,10 @@ class IncompleteLU public: typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; - IncompleteLU() : m_isInitialized(false) {} + IncompleteLU() {} template<typename MatrixType> - IncompleteLU(const MatrixType& mat) : m_isInitialized(false) + IncompleteLU(const MatrixType& mat) { compute(mat); } @@ -71,43 +75,16 @@ class IncompleteLU } template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { x = m_lu.template triangularView<UnitLower>().solve(b); x = m_lu.template triangularView<Upper>().solve(x); } - template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "IncompleteLU is not initialized."); - eigen_assert(cols()==b.rows() - && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived()); - } - protected: FactorType m_lu; - bool m_isInitialized; -}; - -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<IncompleteLU<_MatrixType>, Rhs> - : solve_retval_base<IncompleteLU<_MatrixType>, Rhs> -{ - typedef IncompleteLU<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } }; -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_INCOMPLETE_LU_H diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h index 98f9ecc17..aaf42c78a 100644 --- a/unsupported/Eigen/src/IterativeSolvers/MINRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/MINRES.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu> -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -217,6 +217,7 @@ namespace Eigen { using Base::m_info; using Base::m_isInitialized; public: + using Base::_solve_impl; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -244,26 +245,10 @@ namespace Eigen { /** Destructor. */ ~MINRES(){} - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template<typename Rhs,typename Guess> - inline const internal::solve_retval_with_guess<MINRES, Rhs, Guess> - solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "MINRES is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "MINRES::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - <MINRES, Rhs, Guess>(*this, b.derived(), x0); - } - + /** \internal */ template<typename Rhs,typename Dest> - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; @@ -284,33 +269,16 @@ namespace Eigen { /** \internal */ template<typename Rhs,typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const { x.setZero(); - _solveWithGuess(b,x); + _solve_with_guess_impl(b,x.derived()); } protected: }; - - namespace internal { - - template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs> - struct solve_retval<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs> - : solve_retval_base<MINRES<_MatrixType,_UpLo,_Preconditioner>, Rhs> - { - typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } - }; - - } // end namespace internal - + } // end namespace Eigen #endif // EIGEN_MINRES_H |