Homepage > Man Pages > Category > File Formats
Homepage > Man Pages > Name > G

gmres

man page of gmres

gmres: generalized minimum residual method

NAME

gmres - generalized minimum residual method
SYNOPSIS
template <class Operator, class Vector, class Preconditioner, class Matrix, class Real, class Int> int gmres (const Operator &A, Vector &x, const Vector &b, const Preconditioner &M, Matrix &H, Int m, Int &max_iter, Real &tol);
EXAMPLE
The simplest call to gmres has the folling form: int m = 6; dns H(m+1,m+1); int status = gmres(a, x, b, EYE, H, m, 100, 1e-7);

DESCRIPTION

gmres solves the unsymmetric linear system Ax = b using the generalized minimum residual method. The return value indicates convergence within max_iter (input) iterations (0), or no convergence within max_iter iterations (1). Upon successful return, output arguments have the following values: x approximate solution to Ax = b max_iter the number of iterations performed before the tolerance was reached tol the residual after the final iteration In addition, M specifies a preconditioner, H specifies a matrix to hold the coefficients of the upper Hessenberg matrix constructed by the gmres iterations, m specifies the number of iterations for each restart. gmres requires two matrices as input, A and H. The matrix A, which will typically be a sparse matrix) corresponds to the matrix in the linear system Ax=b. The matrix H, which will be typically a dense matrix, corresponds to the upper Hessenberg matrix H that is constructed during the gmres iterations. Within gmres, H is used in a different way than A, so its class must supply different functionality. That is, A is only accessed though its matrix-vector and transpose- matrix-vector multiplication functions. On the other hand, gmres solves a dense upper triangular linear system of equations on H. Therefore, the class to which H belongs must provide H(i,j) operator for element acess.

NOTE

It is important to remember that we use the convention that indices are 0-based. That is H(0,0) is the first component of the matrix H. Also, the type of the matrix must be compatible with the type of single vector entry. That is, operations such as H(i,j)*x(j) must be able to be carried out. gmres is an iterative template routine. gmres follows the algorithm described on p. 20 in @quotation Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition, R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst, SIAM, 1994, ftp.netlib.org/templates/templates.ps. @end quotation The present implementation is inspired from IML++ 1.2 iterative method library, //math.nist.gov/iml++.

IMPLEMENTATION

template < class Matrix, class Vector, class Int > void Update(Vector &x, Int k, Matrix &h, Vector &s, Vector v[]) { Vector y(s); // Backsolve: for (Int i = k; i >= 0; i--) { y(i) /= h(i,i); for (Int j = i - 1; j >= 0; j--) y(j) -= h(j,i) * y(i); } for (Int j = 0; j <= k; j++) x += v[j] * y(j); } #ifdef TO_CLEAN template < class Real > Real abs(Real x) { return (x > Real(0) ? x : -x); } #endif // TO_CLEAN template < class Operator, class Vector, class Preconditioner, class Matrix, class Real, class Int > int gmres(const Operator &A, Vector &x, const Vector &b, const Preconditioner &M, Matrix &H, const Int &m, Int &max_iter, Real &tol) { Real resid; Int i, j = 1, k; Vector s(m+1), cs(m+1), sn(m+1), w; Real normb = norm(M.solve(b)); Vector r = M.solve(b - A * x); Real beta = norm(r); if (normb == Real(0)) normb = 1; if ((resid = norm(r) / normb) <= tol) { tol = resid; max_iter = 0; return 0; } Vector *v = new Vector[m+1]; while (j <= max_iter) { v[0] = r * (1.0 / beta); // ??? r / beta s = 0.0; s(0) = beta; for (i = 0; i < m && j <= max_iter; i++, j++) { w = M.solve(A * v[i]); for (k = 0; k <= i; k++) { H(k, i) = dot(w, v[k]); w -= H(k, i) * v[k]; } H(i+1, i) = norm(w); v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i) for (k = 0; k < i; k++) ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k)); GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i)); if ((resid = abs(s(i+1)) / normb) < tol) { Update(x, i, H, s, v); tol = resid; max_iter = j; delete [] v; return 0; } } Update(x, m - 1, H, s, v); r = M.solve(b - A * x); beta = norm(r); if ((resid = beta / normb) < tol) { tol = resid; max_iter = j; delete [] v; return 0; } } tol = resid; delete [] v; return 1; } template<class Real> void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) { if (dy == Real(0)) { cs = 1.0; sn = 0.0; } else if (abs(dy) > abs(dx)) { Real temp = dx / dy; sn = 1.0 / ::sqrt( 1.0 + temp*temp ); cs = temp * sn; } else { Real temp = dy / dx; cs = 1.0 / ::sqrt( 1.0 + temp*temp ); sn = temp * cs; } } template<class Real> void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) { Real temp = cs * dx + sn * dy; dy = -sn * dx + cs * dy; dx = temp; } GMRES(5)
 
 
 

Copyright © 2011–2018 by topics-of-interest.com . All rights reserved. Hosted by all-inkl.
Contact · Imprint · Privacy

Page generated in 88.55ms.

dfptutorial.com | autoresponder.name | www.daelim-forum.com