GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 02:39:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

Correct GMRES:

Browse Source 
* Fix error in calculation of residual at restart.
* Use relative residual as stopping criterion.
* Improve documentation.
This commit is contained in:
Kolja Brix 2014-08-02 18:39:15 +02:00
parent e51da9c3a8
commit 953ec08089
1 changed files with 80 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
unsupported/Eigen/src/IterativeSolvers/GMRES.h
View File

@ -11,7 +11,7 @@
#ifndef EIGEN_GMRES_H #ifndef EIGEN_GMRES_H
#define EIGEN_GMRES_H #define EIGEN_GMRES_H
namespace Eigen { namespace Eigen {
namespace internal { namespace internal {
@ -27,11 +27,11 @@ namespace internal {
* \param iters on input: maximum number of iterations to perform * \param iters on input: maximum number of iterations to perform
* on output: number of iterations performed * on output: number of iterations performed
* \param restart number of iterations for a restart * \param restart number of iterations for a restart
* \param tol_error on input: residual tolerance * \param tol_error on input: relative residual tolerance
* on output: residuum achieved * on output: residuum achieved
* *
* \sa IterativeMethods::bicgstab() * \sa IterativeMethods::bicgstab()
* *
* *
* For references, please see: * For references, please see:
* *
@ -70,18 +70,24 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
const int m = mat.rows(); const int m = mat.rows();
VectorType p0 = rhs - mat*x; // residual and preconditioned residual
const VectorType p0 = rhs - mat*x;
VectorType r0 = precond.solve(p0); VectorType r0 = precond.solve(p0);
const RealScalar r0Norm = r0.norm();
// is initial guess already good enough? // is initial guess already good enough?
if(abs(r0.norm()) < tol) { if(r0Norm == 0) {
return true; tol_error=0;
return true;
} }
// storage for Hessenberg matrix and Householder data
FMatrixType H = FMatrixType::Zero(m, restart + 1);
VectorType w = VectorType::Zero(restart + 1); VectorType w = VectorType::Zero(restart + 1);
FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix
VectorType tau = VectorType::Zero(restart + 1); VectorType tau = VectorType::Zero(restart + 1);
// storage for Jacobi rotations
std::vector < JacobiRotation < Scalar > > G(restart); std::vector < JacobiRotation < Scalar > > G(restart);
// generate first Householder vector // generate first Householder vector
@ -112,11 +118,10 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
} }
if (v.tail(m - k).norm() != 0.0) { if (v.tail(m - k).norm() != 0.0) {
if (k <= restart) { if (k <= restart) {
// generate new Householder vector // generate new Householder vector
VectorType e(m - k - 1); VectorType e(m - k - 1);
RealScalar beta; RealScalar beta;
v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta); v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
H.col(k).tail(m - k - 1) = e; H.col(k).tail(m - k - 1) = e;
@ -125,78 +130,77 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition
v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data()); v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
} }
} }
if (k > 1) { if (k > 1) {
for (int i = 0; i < k - 1; ++i) { for (int i = 0; i < k - 1; ++i) {
// apply old Givens rotations to v // apply old Givens rotations to v
v.applyOnTheLeft(i, i + 1, G[i].adjoint()); v.applyOnTheLeft(i, i + 1, G[i].adjoint());
} }
} }
if (k<m && v(k) != (Scalar) 0) { 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 // determine next Givens rotation
v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint()); G[k - 1].makeGivens(v(k - 1), v(k));
w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
} // 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 // insert coefficients into upper matrix triangle
H.col(k - 1).head(k) = v.head(k); H.col(k - 1).head(k) = v.head(k);
bool stop=(k==m || abs(w(k)) < tol || iters == maxIters); bool stop=(k==m || abs(w(k)) < tol * r0Norm || iters == maxIters);
if (stop || k == restart) { if (stop || k == restart) {
// solve upper triangular system // solve upper triangular system
VectorType y = w.head(k); VectorType y = w.head(k);
H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y); H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
// use Horner-like scheme to calculate solution vector // use Horner-like scheme to calculate solution vector
VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
// apply Householder reflection H_{k} to x_new // 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()); 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) { for (int i = k - 2; i >= 0; --i) {
x_new += y(i) * VectorType::Unit(m, i); x_new += y(i) * VectorType::Unit(m, i);
// apply Householder reflection H_{i} to x_new // 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_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
} }
x += x_new; x += x_new;
if (stop) { if (stop) {
return true; return true;
} else { } else {
k=0; k=0;
// reset data for a restart r0 = rhs - mat * x; // reset data for restart
VectorType p0=mat*x; const VectorType p0 = rhs - mat*x;
VectorType p1=precond.solve(p0); r0 = 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 // clear Hessenberg matrix and Householder data
RealScalar beta; H = FMatrixType::Zero(m, restart + 1);
r0.makeHouseholder(e, tau.coeffRef(0), beta); w = VectorType::Zero(restart + 1);
w(0)=(Scalar) beta; tau = VectorType::Zero(restart + 1);
H.bottomLeftCorner(m - 1, 1) = e;
} // 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; 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() * 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 setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
* and NumTraits<Scalar>::epsilon() for the tolerance. * and NumTraits<Scalar>::epsilon() for the tolerance.
* *
* This class can be used as the direct solver classes. Here is a typical usage example: * This class can be used as the direct solver classes. Here is a typical usage example:
* \code * \code
* int n = 10000; * int n = 10000;
@ -244,7 +248,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> >
* // update b, and solve again * // update b, and solve again
* x = solver.solve(b); * x = solver.solve(b);
* \endcode * \endcode
* *
* By default the iterations start with x=0 as an initial guess of the solution. * 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 * 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 * step execution example starting with a random guess and printing the evolution
@ -260,7 +264,7 @@ struct traits<GMRES<_MatrixType,_Preconditioner> >
* } while (solver.info()!=Success && i<100); * } while (solver.info()!=Success && i<100);
* \endcode * \endcode
* Note that such a step by step excution is slightly slower. * Note that such a step by step excution is slightly slower.
* *
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/ */
template< typename _MatrixType, typename _Preconditioner> template< typename _MatrixType, typename _Preconditioner>
@ -272,10 +276,10 @@ class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
using Base::m_iterations; using Base::m_iterations;
using Base::m_info; using Base::m_info;
using Base::m_isInitialized; using Base::m_isInitialized;
private: private:
int m_restart; int m_restart;
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
@ -289,10 +293,10 @@ public:
GMRES() : Base(), m_restart(30) {} GMRES() : Base(), m_restart(30) {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving. /** Initialize the solver with matrix \a A for further \c Ax=b solving.
* *
* This constructor is a shortcut for the default constructor followed * This constructor is a shortcut for the default constructor followed
* by a call to compute(). * by a call to compute().
* *
* \warning this class stores a reference to the matrix A as well as some * \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 * precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new * this class becomes invalid. Call compute() to update it with the new
@ -301,16 +305,16 @@ public:
GMRES(const MatrixType& A) : Base(A), m_restart(30) {} GMRES(const MatrixType& A) : Base(A), m_restart(30) {}
~GMRES() {} ~GMRES() {}
/** Get the number of iterations after that a restart is performed. /** Get the number of iterations after that a restart is performed.
*/ */
int get_restart() { return m_restart; } int get_restart() { return m_restart; }
/** Set the number of iterations after that a restart is performed. /** Set the number of iterations after that a restart is performed.
* \param restart number of iterations for a restarti, default is 30. * \param restart number of iterations for a restarti, default is 30.
*/ */
void set_restart(const int restart) { m_restart=restart; } 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 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution. * \a x0 as an initial solution.
* *
@ -326,17 +330,17 @@ public:
return internal::solve_retval_with_guess return internal::solve_retval_with_guess
<GMRES, Rhs, Guess>(*this, b.derived(), x0); <GMRES, Rhs, Guess>(*this, b.derived(), x0);
} }
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solveWithGuess(const Rhs& b, Dest& x) const void _solveWithGuess(const Rhs& b, Dest& x) const
{ {
bool failed = false; bool failed = false;
for(int j=0; j<b.cols(); ++j) for(int j=0; j<b.cols(); ++j)
{ {
m_iterations = Base::maxIterations(); m_iterations = Base::maxIterations();
m_error = Base::m_tolerance; m_error = Base::m_tolerance;
typename Dest::ColXpr xj(x,j); 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(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
failed = true; failed = true;