% ...................   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.