summaryrefslogtreecommitdiffhomepage
path: root/eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h
diff options
context:
space:
mode:
Diffstat (limited to 'eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h')
-rw-r--r--eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h406
1 files changed, 189 insertions, 217 deletions
diff --git a/eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index ea5deb2..5a82b0d 100644
--- a/eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/eigen/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,10 +282,10 @@ 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
@@ -289,83 +295,49 @@ public:
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
generated by cgit v1.2.3 (git 2.39.1) at 2025-02-23 11:33:38 +0000