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+1);  H = zeros(m+1,m);
r(1) = norm(b);
Q(:,1) = b/r(1);

for k = 1:m

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

end