function [x, error, iter, flag] = gmres( A, x, b, M, restrt, max_it, tol )
% -- Iterative template routine --
% Univ. of Tennessee and Oak Ridge National Laboratory
% October 1, 1993
% Details of this algorithm are described in "Templates for the
% Solution of Linear Systems: Building Blocks for Iterative
% Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
% Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
% 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
% [x, error, iter, flag] = gmres( A, x, b, M, restrt, max_it, tol )
%
% gmres.m solves the linear system Ax=b
% using the Generalized Minimal residual ( GMRESm ) method with restarts .
%
% input A REAL nonsymmetric positive definite matrix
% x REAL initial guess vector
% b REAL right hand side vector
% M REAL preconditioner matrix
% restrt INTEGER number of iterations between restarts
% max_it INTEGER maximum number of iterations
% tol REAL error tolerance
%
% output x REAL solution vector
% error REAL error norm
% iter INTEGER number of iterations performed
% flag INTEGER: 0 = solution found to tolerance
% 1 = no convergence given max_it
iter = 0; % initialization
flag = 0;
bnrm2 = norm( b );
if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end
r = M \ ( b-A*x );
error = norm( r ) / bnrm2;
if ( error < tol ) return, end
[n,n] = size(A); % initialize workspace
m = restrt;
V(1:n,1:m+1) = zeros(n,m+1);
H(1:m+1,1:m) = zeros(m+1,m);
cs(1:m) = zeros(m,1);
sn(1:m) = zeros(m,1);
e1 = zeros(n,1);
e1(1) = 1.0;
for iter = 1:max_it, % begin iteration
r = M \ ( b-A*x );
V(:,1) = r / norm( r );
s = norm( r )*e1;
for i = 1:m, % construct orthonormal
w = M \ (A*V(:,i)); % basis using Gram-Schmidt
for k = 1:i,
H(k,i)= w'*V(:,k);
w = w - H(k,i)*V(:,k);
end
H(i+1,i) = norm( w );
V(:,i+1) = w / H(i+1,i);
for k = 1:i-1, % apply Givens rotation
temp = cs(k)*H(k,i) + sn(k)*H(k+1,i);
H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
H(k,i) = temp;
end
[cs(i),sn(i)] = rotmat( H(i,i), H(i+1,i) ); % form i-th rotation matrix
temp = cs(i)*s(i); % approximate residual norm
s(i+1) = -sn(i)*s(i);
s(i) = temp;
H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
H(i+1,i) = 0.0;
error = abs(s(i+1)) / bnrm2;
if ( error <= tol ), % update approximation
y = H(1:i,1:i) \ s(1:i); % and exit
x = x + V(:,1:i)*y;
break;
end
end
if ( error <= tol ), break, end
y = H(1:m,1:m) \ s(1:m);
x = x + V(:,1:m)*y; % update approximation
r = M \ ( b-A*x ) % compute residual
s(i+1) = norm(r);
error = s(i+1) / bnrm2; % check convergence
if ( error <= tol ), break, end;
end
if ( error > tol ) flag = 1; end; % converged
% END of gmres.m
.