Bugzilla – Attachment 513 Details for
Bug 926
Patch: Improved GMRES Implementation
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
[patch]
Improved GMRES Implementation
gmres_patch.txt (text/plain), 25.68 KB, created by
ettersi
on 2015-01-03 16:56:27 UTC
(
hide
)
Description:
Improved GMRES Implementation
Filename:
MIME Type:
Creator:
ettersi
Created:
2015-01-03 16:56:27 UTC
Size:
25.68 KB
patch
obsolete
># HG changeset patch ># User Simon Etter <ettersi@student.ethz.ch> ># Date 1420299812 -3600 ># Sat Jan 03 16:43:32 2015 +0100 ># Branch 3.2 ># Node ID cece072adcab1febf4fdc1d48a6f6e990cb04a02 ># Parent 33204d4b8dd628b4bb7ed9d5ae4fd8bca91e72aa >Improved GMRES implementation. > >Improvements include: > - Termination based on relative error quantity (like e.g. ConjugateGradient) > - Document termination criterion > - Properly set info() and error() fields after call to solve() > - Rename set_restart() -> setRestart() and get_restart() -> restart() to follow the Eigen naming convention > - More verbous in-code documentation > - Use x = 0 as default initial guess as documented > - Avoid excessive memory use if user specifies large restart parameter > >diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h >--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h >+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h >@@ -1,383 +1,358 @@ > // This file is part of Eigen, a lightweight C++ template library > // for linear algebra. > // > // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> > // Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de> >+// Copyright (C) 2014 Simon Etter > // > // 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 > // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. > > #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. >- * >- */ >+/** \internal Low-level Generalized Minimal Residual algorithm >+ * \param mat The matrix A >+ * \param rhs The right hand side vector b >+ * \param x On input and initial solution, on output the computed solution. >+ * \param precond A preconditioner being able to efficiently solve for an >+ * approximation of Ax=b (regardless of b) >+ * \param iters On input the maximal number of iteration, on output the number of performed iterations. >+ * \param restart Number of iterations before restarting. >+ * \param tol_error On input the tolerance error, on output an estimation of the relative error. >+ * >+ * 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) { >+void gmres( >+ const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond, >+ int &iters, const int &restart, typename Dest::RealScalar & tol_error >+) { >+ 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; // Full Matrix Type > >- using std::sqrt; >- using std::abs; >+ const int m = mat.rows(); >+ const int maxIters = iters; >+ iters = 0; > >- typedef typename Dest::RealScalar RealScalar; >- typedef typename Dest::Scalar Scalar; >- typedef Matrix < Scalar, Dynamic, 1 > VectorType; >- typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType; >+ // Quick exit for zero right-hand side >+ if (rhs == VectorType::Zero(m)) { >+ x.setZero(); >+ tol_error = 0; >+ return; >+ } >+ >+ RealScalar rhsNorm = precond.solve(rhs).norm(); >+ RealScalar threshold = tol_error*rhsNorm; > >- RealScalar tol = tol_error; >- const int maxIters = iters; >- iters = 0; >+ // Quick exit if tolerance already satisfied >+ VectorType residual = precond.solve(rhs - mat*x); // Initial (preconditioned) residual >+ RealScalar residualNorm = residual.norm(); >+ if (residualNorm < threshold) { >+ tol_error = residualNorm/rhsNorm; >+ return; >+ } > >- const int m = mat.rows(); >+ /* >+ * Algorithm outline: >+ * - A := mat, M := precond >+ * - Generate Arnoldi matrices H' (upper Hessenberg), V (orthonormal columns) such that M^{-1} A V.leftCols(k) = V H' >+ * V is stored implicitly as a sequence of Householder reflectors V = V_{0} ... V_{k} >+ * - Determine orthogonal G such that H := G^* H' is upper triangular >+ * G is stored implicitly as a sequence of Jacobi (= Givens) rotations G = G_{0} ... G_{k-1} >+ * - Approximate solution is given by x = V.leftCols(k) y with y least-squares solution of H y = G^* V^* (b - A x0) >+ */ >+ bool stop = false; >+ while (!stop) { >+ VectorType w = VectorType::Zero(restart + 1); // w := G^* V^* b >+ FMatrixType H = FMatrixType::Zero(m, restart + 1); // Upper triangular part: H, >+ // Strictly lower triangular part: Householder vectors for V >+ VectorType tau = VectorType::Zero(restart + 1); // Householder coefficients >+ std::vector<JacobiRotation<Scalar> > G(restart); // Givens rotations for G > >- VectorType p0 = rhs - mat*x; >- VectorType r0 = precond.solve(p0); >- >- // is initial guess already good enough? >- if(abs(r0.norm()) < tol) { >- return true; >- } >+ // Generate first Householder reflector V_{0} >+ VectorType r = precond.solve(rhs - mat*x); >+ RealScalar beta; >+ VectorType e(m-1); >+ r.makeHouseholder(e, tau.coeffRef(0), beta); >+ H.col(0).tail(m-1) = e; >+ w(0)=(Scalar) beta; > >- VectorType w = VectorType::Zero(restart + 1); >+ Scalar workspace; > >- FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix >- VectorType tau = VectorType::Zero(restart + 1); >- std::vector < JacobiRotation < Scalar > > G(restart); >+ int kmax = std::min(m, maxIters - iters); >+ stop = kmax <= restart; >+ kmax = std::min(kmax, restart); >+ int k = 1; >+ for (; k <= kmax; ++k, ++iters) { >+ VectorType v = VectorType::Unit(m, k-1); > >- // 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; >+ // Apply Householder reflectors V_{0} ... V_{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); >+ } >+ // At this point, v is the last Arnoldi vector V.col(k-1) > >- for (int k = 1; k <= restart; ++k) { >+ // Apply matrix M^{-1} A to v >+ v = precond.solve(mat*v); > >- ++iters; >+ // Apply Householder reflectors V_{k-1} ... V_{0} 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); >+ } >+ if (k < m) { >+ // Generate new Householder reflector V_{k} ... >+ RealScalar beta; >+ VectorType e(m-k-1); >+ v.tail(m-k).makeHouseholder(e, tau.coeffRef(k), beta); >+ H.col(k).tail(m-k-1) = e; >+ // ... and apply it to v >+ v.coeffRef(k) = beta; >+ v.tail(m-k-1).setZero(); >+ } >+ // At this point, v = (H').col(k-1) > >- VectorType v = VectorType::Unit(m, k - 1), workspace(m); >+ // Apply Givens rotations G_{0}^* ... G_{k-2}^* to v >+ for (int i = 0; i < k - 1; ++i) { >+ v.applyOnTheLeft(i, i+1, G[i].adjoint()); >+ } >+ if (k < m && v(k) != Scalar(0)) { >+ // Generate new Givens rotation G_{k-1} ... >+ G[k-1].makeGivens(v(k-1), v(k)); >+ // ... and apply it to v and w >+ v.applyOnTheLeft(k-1, k, G[k-1].adjoint()); >+ w.applyOnTheLeft(k-1, k, G[k-1].adjoint()); >+ } > >- // 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()); >- } >+ // At this point, v = H.col(k-1) >+ H.col(k-1).head(k) = v.head(k); > >- // apply matrix M to v: v = mat * v; >- VectorType t=mat*v; >- v=precond.solve(t); >+ // Check for convergence >+ if (abs(w(k)) < threshold) { >+ stop = true; >+ break; >+ } >+ } >+ k = std::min(k, kmax); > >- // 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()); >- } >+ // Solve least squares problem >+ VectorType y = w.head(k); >+ H.topLeftCorner(k,k).template triangularView<Eigen::Upper>().solveInPlace(y); > >- if (v.tail(m - k).norm() != 0.0) { >+ // Use Horner-like scheme to calculate solution vector >+ VectorType dx = y(k - 1) * VectorType::Unit(m, k - 1); >+ // Apply Householder reflection V_{k-1} to dx >+ dx.tail(m-k+1).applyHouseholderOnTheLeft(H.col(k-1).tail(m-k), tau.coeffRef(k-1), &workspace); >+ for (int i = k - 2; i >= 0; --i) { >+ dx += y(i) * VectorType::Unit(m, i); >+ // Apply Householder reflection V_{i} to dx >+ dx.tail(m-i).applyHouseholderOnTheLeft(H.col(i).tail(m-i-1), tau.coeffRef(i), &workspace); >+ } > >- 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); >- >- if (stop || k == restart) { >- >- // solve upper triangular system >- VectorType y = w.head(k); >- H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); >- >- // use Horner-like scheme to calculate solution vector >- VectorType x_new = y(k - 1) * 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()); >- >- 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()); >- } >- >- x += x_new; >- >- if (stop) { >- return true; >- } else { >- k=0; >- >- // 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); >- >- // generate first Householder vector >- RealScalar beta; >- r0.makeHouseholder(e, tau.coeffRef(0), beta); >- w(0)=(Scalar) beta; >- H.bottomLeftCorner(m - 1, 1) = e; >- >- } >- >- } >- >- >- >- } >- >- return false; >- >+ // Update solution >+ x += dx; >+ tol_error = abs(w(k)); >+ } >+ tol_error /= rhsNorm; > } > >-} >+} > >-template< typename _MatrixType, >- typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > >+template< typename _MatrixType, typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > > class GMRES; > > namespace internal { > > template< typename _MatrixType, typename _Preconditioner> > struct traits<GMRES<_MatrixType,_Preconditioner> > > { >- typedef _MatrixType MatrixType; >- typedef _Preconditioner Preconditioner; >+ typedef _MatrixType MatrixType; >+ typedef _Preconditioner Preconditioner; > }; > > } > > /** \ingroup IterativeLinearSolvers_Module >- * \brief A GMRES solver for sparse square problems >- * >- * This class allows to solve for A.x = b sparse linear problems using a generalized minimal >- * residual method. 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 _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner >- * >- * 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; >- * VectorXd x(n), b(n); >- * SparseMatrix<double> A(n,n); >- * // fill A and b >- * GMRES<SparseMatrix<double> > solver(A); >- * x = solver.solve(b); >- * std::cout << "#iterations: " << solver.iterations() << std::endl; >- * std::cout << "estimated error: " << solver.error() << std::endl; >- * // 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. Here is a step by >- * step execution example starting with a random guess and printing the evolution >- * of the estimated error: >- * * \code >- * x = VectorXd::Random(n); >- * solver.setMaxIterations(1); >- * int i = 0; >- * do { >- * x = solver.solveWithGuess(b,x); >- * std::cout << i << " : " << solver.error() << std::endl; >- * ++i; >- * } while (solver.info()!=Success && i<100); >- * \endcode >- * Note that such a step by step excution is slightly slower. >- * >- * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner >- */ >+ * \brief A GMRES solver for sparse square problems >+ * >+ * This class allows to solve for A.x = b sparse linear problems using a generalized minimal >+ * residual method. 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 _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner >+ * >+ * The algorithm is terminated once the maximal number of iterations is reached or the error >+ * \f$ \| M^{-1}.(b - A.x) \| / \| M^{-1} b \| \$ is below a given tolerance. 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; >+ * VectorXd x(n), b(n); >+ * SparseMatrix<double> A(n,n); >+ * // fill A and b >+ * GMRES<SparseMatrix<double> > solver(A); >+ * x = solver.solve(b); >+ * std::cout << "#iterations: " << solver.iterations() << std::endl; >+ * std::cout << "estimated error: " << solver.error() << std::endl; >+ * // 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. Here is a step by >+ * step execution example starting with a random guess and printing the evolution >+ * of the estimated error: >+ * * \code >+ * x = VectorXd::Random(n); >+ * solver.setMaxIterations(1); >+ * int i = 0; >+ * do { >+ * x = solver.solveWithGuess(b,x); >+ * std::cout << i << " : " << solver.error() << std::endl; >+ * ++i; >+ * } while (solver.info()!=Success && i<100); >+ * \endcode >+ * Note that such a step by step excution is slightly slower. >+ * >+ * \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::m_error; >- using Base::m_iterations; >- using Base::m_info; >- using Base::m_isInitialized; >- >+ typedef IterativeSolverBase<GMRES> Base; >+ using Base::mp_matrix; >+ using Base::m_error; >+ using Base::m_iterations; >+ using Base::m_info; >+ using Base::m_isInitialized; >+ > private: >- int m_restart; >- >+ int m_restart; >+ > public: >- typedef _MatrixType MatrixType; >- typedef typename MatrixType::Scalar Scalar; >- typedef typename MatrixType::Index Index; >- typedef typename MatrixType::RealScalar RealScalar; >- typedef _Preconditioner Preconditioner; >+ typedef _MatrixType MatrixType; >+ typedef typename MatrixType::Scalar Scalar; >+ typedef typename MatrixType::Index Index; >+ typedef typename MatrixType::RealScalar RealScalar; >+ typedef _Preconditioner Preconditioner; > > public: > >- /** Default constructor. */ >- GMRES() : Base(), m_restart(30) {} >+ /** Default constructor. */ >+ 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. >+ /** 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) {} >+ >+ ~GMRES() {} >+ >+ /** Get the number of iterations after that a restart is performed. > */ >- GMRES(const MatrixType& A) : Base(A), m_restart(30) {} >+ int restart() { return m_restart; } > >- ~GMRES() {} >- >- /** Get the number of iterations after that a restart is performed. >- */ >- int 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 >+ /** Set the number of iterations after that a restart is performed. >+ * \param restart number of iterations for a restarti, default is 30. >+ */ >+ void setRestart(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 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 >- { >- bool failed = false; >- for(int j=0; j<b.cols(); ++j) >+ } >+ >+ /** \internal */ >+ template<typename Rhs,typename Dest> >+ void _solveWithGuess(const Rhs& b, Dest& x) const >+ { >+ for(int j=0; j<b.cols(); ++j) >+ { >+ m_iterations = Base::maxIterations(); >+ m_error = Base::m_tolerance; >+ int restart = std::min(m_restart, m_iterations); >+ typename Dest::ColXpr xj(x,j); >+ internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, restart, m_error); >+ } >+ m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; >+ m_isInitialized = true; >+ } >+ >+ /** \internal */ >+ template<typename Rhs,typename Dest> >+ void _solve(const Rhs& b, Dest& x) const > { >- 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)) >- failed = true; >+ x.setZero(); >+ _solveWithGuess(b,x); > } >- m_info = failed ? NumericalIssue >- : m_error <= Base::m_tolerance ? Success >- : NoConvergence; >- m_isInitialized = true; >- } >- >- /** \internal */ >- template<typename Rhs,typename Dest> >- void _solve(const Rhs& b, Dest& x) const >- { >- x = b; >- if(x.squaredNorm() == 0) return; // Check Zero right hand side >- _solveWithGuess(b,x); >- } >- >-protected: > > }; > > > namespace internal { > >- template<typename _MatrixType, typename _Preconditioner, typename Rhs> >+template<typename _MatrixType, typename _Preconditioner, typename Rhs> > struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs> >- : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs> >+: solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs> > { >- typedef GMRES<_MatrixType, _Preconditioner> Dec; >- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) >+ typedef GMRES<_MatrixType, _Preconditioner> Dec; >+ EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) > >- template<typename Dest> void evalTo(Dest& dst) const >- { >- dec()._solve(rhs(),dst); >- } >+ template<typename Dest> void evalTo(Dest& dst) const >+ { >+ dec()._solve(rhs(),dst); >+ } > }; > > } // end namespace internal > > } // end namespace Eigen > > #endif // EIGEN_GMRES_H
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Diff
View Attachment As Raw
Actions:
View
|
Diff
Attachments on
bug 926
: 513