function X = gmresid(A,b,n)
% usage: X = gmresid(A,b,n)
% description: given mxm A and mx1 b != 0, this function
% calculates n steps of GMRES, provided termination
% does not occur before n, and returns mxn X containing
% the first n approximates to solution of A*x = b.
% We do not count the 0th, which is x^(0)=b.

% local variables:
% ii,jj: index variables

m = length(b);
Q = zeros(m,n+1);
H = zeros(n+1,n);
X = zeros(m,n);
Aeps = norm(A,inf)*eps;
normb = norm(b);
e1 = [normb;zeros(n,1)];
v = b/normb;
Q(:,1) = v;
for jj = 1:n
  v = A*Q(:,jj);
  for ii = 1:jj
    H(ii,jj) = Q(:,ii)'*v;
    v = v - H(ii,jj)*Q(:,ii);
  end
  xtmp = norm(v);
  if (xtmp <= Aeps)
    jj = n;
  else
    H(jj+1,jj) = xtmp;
    Q(:,jj+1) = v/xtmp;
  end
  y = H(1:jj+1,1:jj)\e1(1:jj+1);
  X(:,jj) = Q(:,1:jj)*y;
end