Function [x , normrn] = gmres(A,b,maxiter)
function [x , normrn] = gmres(A,b,maxiter)
% solves A x = b using gmres
% input: A - m by m matrix
% b - m by 1 vector
% output: x - approximate solution
% normrn - norm(b-A*x) at each step of algorithm
% Also
% plot normrn/ norm(b) versus n (step number)
% Example:
% m = 200;
% A = eigmat(2,1/2,m);
% xtrue = randn(m,1);
% b = A*xtrue;
% [x , normrn] = gmres(A,b,20);
% CAUTION: Matlab has a built-in function with the same
% name (gmres). If you want to use the new version
% of gmres make sure your active folder contains the
% new version of gmres in file gmres.m.
Q = [];
H = 0;
m = length(A);
normb = norm(b);
normrn=normb;
Q(:,1) = b / normb;
for n = 1:maxiter
% Arnoldi step calculations (Algorithm 33.1 for Trefethen and Bau)
v = A*Q(:,n);
for j = 1:n
H(j,n) = Q(:,j)'* v;
v = v - H(j,n)*Q(:,j);
end
Hn = H(1:n,1:n);
H(n+1,n) = norm(v);
if H(n+1,n) == 0, break, end % breakdown so stop
Q(:,n+1) = v / H(n+1,n);
e1 = [1;zeros(n,1)];
y = H \ (normb*e1); % This can be done much more quickly
% using Givens rotations.
% For simplicity we just use Matlab's \
normrn = [normrn,norm(H*y-normb*e1)]; % remember residual norm
end
x = Q(:,1:n)*y;
clf
semilogy(0:n,normrn/normb,'--o'),shg
xlabel('step number of gmres algorithm')
ylabel(' norm( b - A * xn ) / norm(b)')
title('convergence of residual in gmres')
grid
% file example35p1.m
% example 35.1 from Trefethen + Bau
m=200;
maxiter = 20;
A=eigmat(2,.5,m);
lam=eig(A);
figure(1)
clf
plot(lam,'x');
axis square
title('eigenvalues of A')
ylabel('imag part')
xlabel('real part')
grid
xtrue=randn(m,1);
b = A*xtrue;
figure(2)
[x , normrn] = gmres(A,b,maxiter);
[pic]
[pic]
% file example35p2.m
% example 35.2 from Trefethen + Bau
m=200;
maxiter = 20;
th = (0:m-1)*pi / (m-1);
d=(-2 +2* sin(th))+sqrt(-1)*cos(th);
A=eigmat(2,.5,m)+diag(d);
lam=eig(A);
figure(1)
clf
plot(lam,'x');
axis square
title('eigenvalues of A')
ylabel('imag part')
xlabel('real part')
grid
xtrue=randn(m,1);
b = A*xtrue;
figure(2)
[x , normrn] = gmres(A,b,maxiter);
[pic] [pic]
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.