function [x,r] = arngmres(A,b,M)
% GMRES   GMRES for a linear system (demo only).
% Input:
%   A       square matrix (n by n)
%   b       right-hand side
%   M       number of iterations
% Output: 
%   x       approximate solution
%   r       history of norms of the residuals

n = length(A);
Q = zeros(n,M);  H = zeros(M,M-1);

% Initial "solution" is zero.
r(1) = norm(b);

for m = 1:M

  % Use Arnoldi iteration to get the mth column of Q.
  if m==1
    Q(:,1) = b/norm(b);
  else
    v = A*Q(:,m-1);
    for j = 1:m-1    
      H(j,m-1) = Q(:,j)'*v;
      v = v - H(j,m-1)*Q(:,j);
    end
    H(m,m-1) = norm(v);
    Q(:,m) = v/H(m,m-1);
  end
  
  % Solve the GMRES problem with the current Q.
  z = (A*Q(:,1:m)) \ b;
  x = Q(:,1:m)*z;
  r(m+1) = norm( A*x - b );
  
end