% ................... gmres.m ......................
%
% The following notations differs slightly from those used in the
% code presented in Section 11.8 of the book "Concise Numerical Mathematics".
% The notations used here are as close as possible to
% those used in the Section 11.6 and 11.7 of the book.
% Several comments have been added.
clear;
% Some initializations
N = 64;
x = zeros(N,1); res = zeros(N,1);
u = zeros(2,1); y = zeros(N,1);
Q = zeros(N,N); R = zeros(N,N);
w = zeros(2,N);
% End of some initializations
% Vectors are initialized here in the
% form zeros(.,1).
% Start definition of a test matrix A and a right-hand side b.
A = zeros(N,N);
for i = 1:N
for j = 1:N
if (j == i-1) A(i,i-1)= 1; end
end
end
A(1,N)= 1;
b = zeros(N,1); xstar = zeros(N,1);
b(1) = 1; xstar(N) = 1;
y = zeros(N,1); y(1) = norm(b);
% End definition of a test matrix A and a right-hand side b.
%
% The vector xstar is the solution to Ax=b. The vector
% appears in Theorem 11.28 in the book (denoted by c_n there)
% and in fact is the right-hand
% side in the minimization problem with the Hessenberg matrix $ H_n $.
% Actually, y = c_n coincides with the first unit vector up
% to a constant factor (see (11.33) in the book).
h = zeros(N,1); % At the n-th iteration step this vector is used to store the
% n-th column of the Hessenberg matrix H_n.
go_ahead = 1;
n = 1;
Q(:,1) = b/norm(b); % The matrix Q contains the orthonormal system of
myeps = 0.000001; % Arnoldi process.
% Begin of iterations; n is the iteration number.
while (go_ahead == 1)
v = A*Q(:,n); % n-th step of Arnoldi process:
z = v; % Computation of (n+1)-th orthogonal vector Q(:,n+1),
% and also computation of the n-th column of the
% Hessenberg matrix which in fact
% is saved in the vector h here.
for k = 1:n
h(k) = Q(:,k)'*v;
z = z - h(k)*Q(:,k);
end
h(n+1) = norm(z);
if ( (abs(h(n+1)) <= myeps) | (n==N) ) % This is the termination case
go_ahead = 0; % for the Arnoldi process.
else
Q(:,n+1) = z/h(n+1);
end % End of n-th step of Arnoldi process.
% Next the factorization H_n = T_n S_n
% of the Hessenberg matrix H_n starts,
% with T_n being orthogonal, and S_n is a generalized
% upper triangualar (n+1)-by-n matrix.
% Its triangular n-by-n submatrix is denoted by R_n in the book (and by $ R $
% here in this code), and the entries in the last row
% of S_n all vanish.
% The first n-1 columns of matrix R_n are already given
% from the previous transformations, and we
% start the application of previous orthogonal transformation
% as a first step for the computation of n-th column of R_n.
for k = 1:n-1 %
u = h(k:k+1); % This is the procedure described in (11.37) in the book,
v = w(:,k); % and see Remark 11.29 for its technical realization.
h(k:k+1) = u - 2*(v'*u)*v;
end
R(1:n-1,n) = h(1:n-1);
% End of application of previous orthogonal transformations.
% Here, the two-dimensional vector w(:,k) corresponds to the vector
% w^[k] representing the Householder transformation
% W^[k] = I - 2w^[k](w^[k])^transpose.
% We recall that the details are given in Section 11.7.3 in the book,
% especially the considerations after (11.37) and Remark 11.29 there.
% Next computation of new orthogonal transformation
% W^[n] = I - 2w^[n](w^[n])^transpose in R^(2x2) (see also Remark 11.29 in the book)
% to complete the computation of n-th column of R_n.
% Also the rest of the n-th column of the upper triangular matrix R
% used in the orthogonal factorization of the Hessenberg matrix H_n
% is computed according to the procedure decribed in Section 11.7.3
% of the book.
if (go_ahead == 0) % (this is the termination case)
R(n,n) = h(n);
else % (no termination case)
u = h(n:n+1);
if (abs(u(1)) <= myeps) sigma = -norm(u);
else sigma = -norm(u)*u(1)/abs(u(1));
end
u(1) = u(1) - sigma; u(2) = u(2);
w(:,n) = u/norm(u);
R(n,n) = sigma;
R(n+1,n) = 0;
%
% Next we consider the computation
% of the first n entries of the vector (T_n)^transpose c_n,
% with the result being saved in the variable y here.
% The vector c_n here again is the vector
% appearing in Theorem 11.28 in the book.
%
% These computations can be realized by modifying
% only two entries of the vector T_(n-1)^transpose c_(n-1)
% (which is available in the variable y)
% via a Householder transformation.
u = y(n:n+1); v = w(:,n);
y(n:n+1) = u - 2*(v'*u)*v;
end
% Start computation of the vector z (coincides with the vector z_n in Theorem 11.28
% of the book)
% and the n-th GMRES iteration x_n = Q_n z_n.
%
% In the sequel we recall some basic facts about this procedure;
% for simplicity the considerations are restricted to the non-termination case:
% The unknown vector z_n satisfies R_n z_n = y, where the vector y coincides
% with the first n entries of the vector (T_n)^transpose c_n.
% Here c_n again coincides
% with the first unit vector up to a constant factor (see (11.33) in the book).
% The transformation (T_n)^transpose c_n is already computed (and is available in
% in the variable y). So it remains to solve an upper triangular linear system
% to determine z_n.
for k = n:-1:1
sum = y(k);
for j = k+1:n
sum = sum - R(k,j)*z(j);
end
z(k) = sum/R(k,k);
end
% Now computation of the GMRES iterate x = x_n = Q z_n starts.
x = Q(:,1:n)*z(1:n);
res(n) = norm(A*x - b);
err(n) = norm(x - xstar);
if (go_ahead==1) n = n+1; end % End of n-th iteration step
end
% End of iterations
% Start output of data to file main.dat.
fid = fopen('main.dat','w');
fprintf(fid,' n || Ax_n-b ||_2 || x_n-x_* ||_2 \n\n');
for k = 1:n
fprintf(fid,' %2g %12.9f %12.9f \n',k,res(k),err(k));
end
fprintf(fid,'------------------------------------------------\n\n');
fclose(fid);
% (... End output of data ...)
%..........................................................
%
% Some additional comments:
% 1) The program runs under MATLAB 6 as well as under Octave.
% 2) Thanks are going to G. Fuss (TU Berlin). His comments
% encouraged me to revise an earlier version of this code.