%% gmres_k.m --- Illustration of GMRES(k) -- Chapter 3
%%
%%       Solve: A*x = b for x
%%
%% Usage:     x=gmres_k(A,b,k,tol,maxit,x0)
%%    or      x=gmres_k(A,b,k)
%%    or      x=gmres_k(A,b)
%% Sample data:  rand('seed',3);  n=16,
%                A=rand(n)+2*eye(n); A=diag(diag(A))\A; x=(n:-1:1)'; b=A*x;
%             x=gmres_k(A,b,7,1.0e-2,3);

function x=gmres_k(A,b,k,tol,maxit,x0)
if nargin<2, help gmres_k; return; end
   [m n]=size(A);   if m~=n, error('Check dimension of A!'); end
if nargin<3,  k=min(5,m); end
if nargin<4,  tol=1.0e-4; end %=========================%
if nargin<5,     maxit=m; end
if nargin<6, x0=zeros(m,1); end
x=x0;
           r = b - A*x;  norm_0=norm(r); iter=0; ok=0;

while iter<maxit

% Step (1)
          x0 = x;  r = b - A*x;  r2=norm(r); Q=[];R=[];c=[];s=[];
% Step (2)
          q1 = r/r2; Q=q1;   rhs=zeros(m,1);rhs(1)=r2;
   for i=1:k
% Step (3)
          w = A*Q(:,i);
% Step (4)
   for L=1:i
     R(L,i) = w'*Q(:,L);  w=w-R(L,i)*Q(:,L);
   end 
% Step (5)
   R(i+1,i) = norm(w);   Q(:,i+1)=w/R(i+1,i);
% Step (6) 
   for L=1:(i-1)
     t=c(L)*R(L,i) + s(L)*R(L+1,i);
     R(L+1,i)=s(L)*R(L,i) - c(L)*R(L+1,i);
     R(L,i)=t;
   end
% Step (7)
   if abs(R(i+1,i))<1.0e-16, c(i)=1; s(i)=0;
   elseif abs(R(i+1,i))>abs(R(i,i))
      t=R(i,i)/R(i+1,i); s(i)=1/sqrt(1+t^2);  c(i)=s(i)*t; else
      t=R(i+1,i)/R(i,i); c(i)=1/sqrt(1+t^2);  s(i)=c(i)*t; 
   end
% Step (8)
   t=c(i)*rhs(i);  rhs(i+1)=s(i)*rhs(i);  rhs(i)=t;
   R(i,i)=c(i)*R(i,i) + s(i)*R(i+1,i);
   R(i+1,i) = 0; 
% Step (9)
   y = R(1:i,1:i) \ rhs(1:i,1);
% Step (10)
   x = x0 + Q(:,1:i)*y;
% Step (11)
   r=b-A*x;   r2=norm(r)/norm_0;
   fprintf('%s iter=%3d, res=%e\n',mfilename,iter*k+i,r2)
   if r2 <= tol, ok=1; break; disp('gmres_k done...'); end
% Step (12)
end % i
   if ok==1, break; end
   iter=iter+1;
end %while
   if iter>=maxit, iter=iter-1; end
   fprintf('%s: Total GMRES(%d) Steps = %d to achieve %e [N=%d]\n',...
            mfilename, k, iter*k+i, norm(r)/norm_0, m);
