From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- .../Eigen/src/IterativeSolvers/MINRES.h | 84 ++++++++-------------- 1 file changed, 31 insertions(+), 53 deletions(-) (limited to 'eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h') diff --git a/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h index 670f274..256990c 100644 --- a/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h +++ b/eigen/unsupported/Eigen/src/IterativeSolvers/MINRES.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Giacomo Po -// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2011-2014 Gael Guennebaud // // 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 @@ -29,7 +29,7 @@ namespace Eigen { template EIGEN_DONT_INLINE void minres(const MatrixType& mat, const Rhs& rhs, Dest& x, - const Preconditioner& precond, int& iters, + const Preconditioner& precond, Index& iters, typename Dest::RealScalar& tol_error) { using std::sqrt; @@ -48,8 +48,8 @@ namespace Eigen { } // initialize - const int maxIters(iters); // initialize maxIters to iters - const int N(mat.cols()); // the size of the matrix + const Index maxIters(iters); // initialize maxIters to iters + const Index N(mat.cols()); // the size of the matrix const RealScalar threshold2(tol_error*tol_error*rhsNorm2); // convergence threshold (compared to residualNorm2) // Initialize preconditioned Lanczos @@ -144,7 +144,6 @@ namespace Eigen { template< typename _MatrixType, int _UpLo=Lower, typename _Preconditioner = IdentityPreconditioner> -// typename _Preconditioner = IdentityPreconditioner > // preconditioner must be positive definite class MINRES; namespace internal { @@ -166,8 +165,8 @@ namespace Eigen { * The vectors x and b can be either dense or sparse. * * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. - * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower - * or Upper. Default is Lower. + * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower, + * Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower. * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner * * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() @@ -192,6 +191,8 @@ namespace Eigen { * By default the iterations start with x=0 as an initial guess of the solution. * One can control the start using the solveWithGuess() method. * + * MINRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. + * * \sa class ConjugateGradient, BiCGSTAB, SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename _MatrixType, int _UpLo, typename _Preconditioner> @@ -199,15 +200,15 @@ namespace Eigen { { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; 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; typedef typename MatrixType::RealScalar RealScalar; typedef _Preconditioner Preconditioner; @@ -233,42 +234,36 @@ 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 - inline const internal::solve_retval_with_guess - solveWithGuess(const MatrixBase& 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 - (*this, b.derived(), x0); - } - + /** \internal */ template - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { + typedef typename Base::MatrixWrapper MatrixWrapper; + typedef typename Base::ActualMatrixType ActualMatrixType; + enum { + TransposeInput = (!MatrixWrapper::MatrixFree) + && (UpLo==(Lower|Upper)) + && (!MatrixType::IsRowMajor) + && (!NumTraits::IsComplex) + }; + typedef typename internal::conditional, ActualMatrixType const&>::type RowMajorWrapper; + EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY); typedef typename internal::conditional - >::type MatrixWrapperType; - + RowMajorWrapper, + typename MatrixWrapper::template ConstSelfAdjointViewReturnType::Type + >::type SelfAdjointWrapper; + m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - + RowMajorWrapper row_mat(matrix()); for(int j=0; j - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, MatrixBase &x) const { x.setZero(); - _solveWithGuess(b,x); + _solve_with_guess_impl(b,x.derived()); } protected: }; - - namespace internal { - - template - struct solve_retval, Rhs> - : solve_retval_base, Rhs> - { - typedef MINRES<_MatrixType,_UpLo,_Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } - }; - - } // end namespace internal - + } // end namespace Eigen #endif // EIGEN_MINRES_H -- cgit v1.2.3