function [u,uk,res]=GMRES(A,f,u0,tol,m)
% GMRES Generalized Minimal Residual
%   [u,uk,res]=GMRES(A,f,u0,tol,m) solves Au=f using the GMRES 
%   method starting at the initial guess u0 up to a tolerance tol using 
%   at most m iterations. It computes u in the affine Krylov space
%   u0+(r0,Ar0...,A^(k-1)r0) by minimizing the norm ||f-A*u||_2 where k is
%   the smallest integer such that ||f-A*u||_2<tol.
%   GMRES returns in the matrix uk the iterates, in u the solution
%   computed, and in res the history of the norm of the residuals.

if nargin<5, m=100; end            % default values
if nargin<4, tol=1e-6; end
k=1;                               % Initialize the iteration index
uk(:,1)=u0;                        
r0=f-A*u0;                         % Compute the initial residual
res=norm(r0);                      % Initialize the vercor of residuals
rhs=res;                           % Initialize the right-hand side
Q(:,1)=r0/res;                     % First column of Q
while res(end)/norm(f)>tol && k<=m % GMRES iterations
  Q(:,k+1)=A*Q(:,k);               % PART 1: Arnoldi: (A Q_k = Q_{k+1} H_k)
  for j=1:k                        % modified Gram-Schmidt on v=A*Q(:,k)
    H(j,k)=Q(:,k+1)'*Q(:,j);              
    Q(:,k+1)=Q(:,k+1)-H(j,k)*Q(:,j);
  end
  H(k+1,k)=norm(Q(:,k+1));         % Lower-diagonal element
  Q(:,k+1)=Q(:,k+1)/H(k+1,k);      % Compute the vector Q(:,k+1)  
  R(k+1,k)=0;                      % PART 2: QR by Givens rotations
  R(:,k)=H(:,k);                   % New column of R
  for j=1:k-1                      % Former Givens rotations on last column
    Rk=c(j)*R(j,k)+s(j)*R(j+1,k);
    R(j+1,k)=-s(j)*R(j,k)+c(j)*R(j+1,k);
    R(j,k)=Rk;
  end                              % Next apply the new Givens rotation  
  c(k)=R(k,k)/sqrt(R(k,k)^2+R(k+1,k)^2);
  s(k)=R(k+1,k)/sqrt(R(k,k)^2+R(k+1,k)^2);
  R(k,k)=c(k)*R(k,k)+s(k)*R(k+1,k);
  R(k+1,k)=0;                      
  rhs(k+1)=-s(k)*rhs(k);           % PART 3: Update rhs and residuals
  rhs(k)=c(k)*rhs(k);              % New Givens rotation on rhs
  res(k+1)=abs(rhs(k+1));          % Compute the new norm of the residual
  k=k+1;                           % Update iteration index
  if nargout>1
    y=R\rhs';                      % Compute the coefficients y_j
    uk(:,k)=u0+Q(:,1:end-1)*y;     % Compute u as sum_j y_j q_j 
  end
end
y=R\rhs';                          % Compute the coefficients y_j
u=u0+Q(:,1:end-1)*y;               % Compute u as u0+sum_j y_j q_j 
