diff options
Diffstat (limited to 'unsupported/Eigen/src/IterativeSolvers/GMRES.h')
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 409 |
1 files changed, 191 insertions, 218 deletions
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 7ba13afd2..5a82b0df6 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -11,193 +11,197 @@ #ifndef EIGEN_GMRES_H #define EIGEN_GMRES_H -namespace Eigen { +namespace Eigen { namespace internal { /** - * Generalized Minimal Residual Algorithm based on the - * Arnoldi algorithm implemented with Householder reflections. - * - * Parameters: - * \param mat matrix of linear system of equations - * \param Rhs right hand side vector of linear system of equations - * \param x on input: initial guess, on output: solution - * \param precond preconditioner used - * \param iters on input: maximum number of iterations to perform - * on output: number of iterations performed - * \param restart number of iterations for a restart - * \param tol_error on input: residual tolerance - * on output: residuum achieved - * - * \sa IterativeMethods::bicgstab() - * - * - * For references, please see: - * - * Saad, Y. and Schultz, M. H. - * GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. - * SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869. - * - * Saad, Y. - * Iterative Methods for Sparse Linear Systems. - * Society for Industrial and Applied Mathematics, Philadelphia, 2003. - * - * Walker, H. F. - * Implementations of the GMRES method. - * Comput.Phys.Comm. 53, 1989, pp. 311 - 320. - * - * Walker, H. F. - * Implementation of the GMRES Method using Householder Transformations. - * SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163. - * - */ +* Generalized Minimal Residual Algorithm based on the +* Arnoldi algorithm implemented with Householder reflections. +* +* Parameters: +* \param mat matrix of linear system of equations +* \param Rhs right hand side vector of linear system of equations +* \param x on input: initial guess, on output: solution +* \param precond preconditioner used +* \param iters on input: maximum number of iterations to perform +* on output: number of iterations performed +* \param restart number of iterations for a restart +* \param tol_error on input: relative residual tolerance +* on output: residuum achieved +* +* \sa IterativeMethods::bicgstab() +* +* +* For references, please see: +* +* Saad, Y. and Schultz, M. H. +* GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. +* SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869. +* +* Saad, Y. +* Iterative Methods for Sparse Linear Systems. +* Society for Industrial and Applied Mathematics, Philadelphia, 2003. +* +* Walker, H. F. +* Implementations of the GMRES method. +* Comput.Phys.Comm. 53, 1989, pp. 311 - 320. +* +* Walker, H. F. +* Implementation of the GMRES Method using Householder Transformations. +* SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163. +* +*/ template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond, - int &iters, const int &restart, typename Dest::RealScalar & tol_error) { + Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) { - using std::sqrt; - using std::abs; + using std::sqrt; + using std::abs; - typedef typename Dest::RealScalar RealScalar; - typedef typename Dest::Scalar Scalar; - typedef Matrix < Scalar, Dynamic, 1 > VectorType; - typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType; + typedef typename Dest::RealScalar RealScalar; + typedef typename Dest::Scalar Scalar; + typedef Matrix < Scalar, Dynamic, 1 > VectorType; + typedef Matrix < Scalar, Dynamic, Dynamic, ColMajor> FMatrixType; - RealScalar tol = tol_error; - const int maxIters = iters; - iters = 0; + RealScalar tol = tol_error; + const Index maxIters = iters; + iters = 0; - const int m = mat.rows(); + const Index m = mat.rows(); - VectorType p0 = rhs - mat*x; - VectorType r0 = precond.solve(p0); - - // is initial guess already good enough? - if(abs(r0.norm()) < tol) { - return true; - } + // residual and preconditioned residual + VectorType p0 = rhs - mat*x; + VectorType r0 = precond.solve(p0); - VectorType w = VectorType::Zero(restart + 1); + const RealScalar r0Norm = r0.norm(); - FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix - VectorType tau = VectorType::Zero(restart + 1); - std::vector < JacobiRotation < Scalar > > G(restart); - - // generate first Householder vector - VectorType e(m-1); - RealScalar beta; - r0.makeHouseholder(e, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; - H.bottomLeftCorner(m - 1, 1) = e; - - for (int k = 1; k <= restart; ++k) { - - ++iters; - - VectorType v = VectorType::Unit(m, k - 1), workspace(m); - - // apply Householder reflections H_{1} ... H_{k-1} to v - for (int i = k - 1; i >= 0; --i) { - v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } - - // apply matrix M to v: v = mat * v; - VectorType t=mat*v; - v=precond.solve(t); - - // apply Householder reflections H_{k-1} ... H_{1} to v - for (int i = 0; i < k; ++i) { - v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } - - if (v.tail(m - k).norm() != 0.0) { - - if (k <= restart) { - - // generate new Householder vector - VectorType e(m - k - 1); - RealScalar beta; - v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta); - H.col(k).tail(m - k - 1) = e; - - // apply Householder reflection H_{k} to v - v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data()); - - } - } - - if (k > 1) { - for (int i = 0; i < k - 1; ++i) { - // apply old Givens rotations to v - v.applyOnTheLeft(i, i + 1, G[i].adjoint()); - } - } - - if (k<m && v(k) != (Scalar) 0) { - // determine next Givens rotation - G[k - 1].makeGivens(v(k - 1), v(k)); - - // apply Givens rotation to v and w - v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); - w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); - - } - - // insert coefficients into upper matrix triangle - H.col(k - 1).head(k) = v.head(k); - - bool stop=(k==m || abs(w(k)) < tol || iters == maxIters); + // is initial guess already good enough? + if(r0Norm == 0) + { + tol_error = 0; + return true; + } - if (stop || k == restart) { + // storage for Hessenberg matrix and Householder data + FMatrixType H = FMatrixType::Zero(m, restart + 1); + VectorType w = VectorType::Zero(restart + 1); + VectorType tau = VectorType::Zero(restart + 1); - // solve upper triangular system - VectorType y = w.head(k); - H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); + // storage for Jacobi rotations + std::vector < JacobiRotation < Scalar > > G(restart); + + // storage for temporaries + VectorType t(m), v(m), workspace(m), x_new(m); + + // generate first Householder vector + Ref<VectorType> H0_tail = H.col(0).tail(m - 1); + RealScalar beta; + r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); + w(0) = Scalar(beta); + + for (Index k = 1; k <= restart; ++k) + { + ++iters; - // use Horner-like scheme to calculate solution vector - VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); + v = VectorType::Unit(m, k - 1); - // apply Householder reflection H_{k} to x_new - x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); + // apply Householder reflections H_{1} ... H_{k-1} to v + // TODO: use a HouseholderSequence + for (Index i = k - 1; i >= 0; --i) { + v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } - for (int i = k - 2; i >= 0; --i) { - x_new += y(i) * VectorType::Unit(m, i); - // apply Householder reflection H_{i} to x_new - x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } + // apply matrix M to v: v = mat * v; + t.noalias() = mat * v; + v = precond.solve(t); - x += x_new; + // apply Householder reflections H_{k-1} ... H_{1} to v + // TODO: use a HouseholderSequence + for (Index i = 0; i < k; ++i) { + v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } - if (stop) { - return true; - } else { - k=0; + if (v.tail(m - k).norm() != 0.0) + { + if (k <= restart) + { + // generate new Householder vector + Ref<VectorType> Hk_tail = H.col(k).tail(m - k - 1); + v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta); + + // apply Householder reflection H_{k} to v + v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data()); + } + } - // reset data for a restart r0 = rhs - mat * x; - VectorType p0=mat*x; - VectorType p1=precond.solve(p0); - r0 = rhs - p1; -// r0_sqnorm = r0.squaredNorm(); - w = VectorType::Zero(restart + 1); - H = FMatrixType::Zero(m, restart + 1); - tau = VectorType::Zero(restart + 1); + if (k > 1) + { + for (Index i = 0; i < k - 1; ++i) + { + // apply old Givens rotations to v + v.applyOnTheLeft(i, i + 1, G[i].adjoint()); + } + } - // generate first Householder vector - RealScalar beta; - r0.makeHouseholder(e, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; - H.bottomLeftCorner(m - 1, 1) = e; + if (k<m && v(k) != (Scalar) 0) + { + // determine next Givens rotation + G[k - 1].makeGivens(v(k - 1), v(k)); - } + // apply Givens rotation to v and w + v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); + w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); + } - } + // insert coefficients into upper matrix triangle + H.col(k-1).head(k) = v.head(k); + tol_error = abs(w(k)) / r0Norm; + bool stop = (k==m || tol_error < tol || iters == maxIters); + if (stop || k == restart) + { + // solve upper triangular system + Ref<VectorType> y = w.head(k); + H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y); + + // use Horner-like scheme to calculate solution vector + x_new.setZero(); + for (Index i = k - 1; i >= 0; --i) + { + x_new(i) += y(i); + // apply Householder reflection H_{i} to x_new + x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } + + x += x_new; + + if(stop) + { + return true; + } + else + { + k=0; + + // reset data for restart + p0.noalias() = rhs - mat*x; + r0 = precond.solve(p0); + + // clear Hessenberg matrix and Householder data + H.setZero(); + w.setZero(); + tau.setZero(); + + // generate first Householder vector + r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); + w(0) = Scalar(beta); + } + } + } - } - - return false; + return false; } @@ -230,7 +234,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> > * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations * and NumTraits<Scalar>::epsilon() for the tolerance. - * + * * This class can be used as the direct solver classes. Here is a typical usage example: * \code * int n = 10000; @@ -244,29 +248,31 @@ struct traits<GMRES<_MatrixType,_Preconditioner> > * // update b, and solve again * x = solver.solve(b); * \endcode - * + * * By default the iterations start with x=0 as an initial guess of the solution. * One can control the start using the solveWithGuess() method. * + * GMRES can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. + * * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename _MatrixType, typename _Preconditioner> class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> > { typedef IterativeSolverBase<GMRES> Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; using Base::m_isInitialized; - + private: - int m_restart; - + Index m_restart; + 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; @@ -276,95 +282,62 @@ public: GMRES() : Base(), m_restart(30) {} /** Initialize the solver with matrix \a A for further \c Ax=b solving. - * + * * This constructor is a shortcut for the default constructor followed * by a call to compute(). - * + * * \warning this class stores a reference to the matrix A as well as some * precomputed values that depend on it. Therefore, if \a A is changed * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - GMRES(const MatrixType& A) : Base(A), m_restart(30) {} + template<typename MatrixDerived> + explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {} ~GMRES() {} - + /** Get the number of iterations after that a restart is performed. */ - int get_restart() { return m_restart; } - + Index get_restart() { return m_restart; } + /** Set the number of iterations after that a restart is performed. * \param restart number of iterations for a restarti, default is 30. */ - 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); - } - + void set_restart(const Index restart) { m_restart=restart; } + /** \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) + for(Index j=0; j<b.cols(); ++j) { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - + typename Dest::ColXpr xj(x,j); - if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) + if(!internal::gmres(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) failed = true; } m_info = failed ? NumericalIssue - : m_error <= Base::m_tolerance ? Success - : NoConvergence; + : m_error <= Base::m_tolerance ? Success + : NoConvergence; m_isInitialized = true; } /** \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 |