function [x,itno,lrrnrm]  = gmres514( A, b, x, m, M, tol, itmax) 
%
% Implements standard Gram-Schmidt GMRES(m) to produce x such that
% ||b-Ax|| <= tol*||initial residual||.  
% 
% The inputs are 
% Required:
%   A       = nxn coefficient matrix. This can also be a function handle (@A) for 
%                producing products of the matrix A with vectors; the form is w = A(v).
%   b       = nx1 right-hand side vector. 
% Optional:
%   x       = initial iterate; default 0.
%   M       = left preconditioner, i.e., the system solved is M^{-1}Ax = M^{-1}b. 
%                This can also be a function handle (@M) for producing solutions
%                 of Mw = v for given v; the form is w = M(v). The default is M = []
%                 (no preconditioning). 
%   tol     = stopping tolerance; default 1.e-6.
%   m       = maximum number of iterations before restarting; default min(n,20).
%   itmax   = maximum allowable number of iterations; default min(n,2000). 
%
% The returns are 
%  x = final iterate. 
%   itno    = vector of iteration numbers. 
%   lrrnrm  = vector of log relative residual norms. 
%

%
% Check the inputs and initialize. 
%

if (nargin<2 | nargin>7) 
  error('gmres514: The number of arguments must be between 2 and 7.');
end
if ( ~isreal(A) & ~isa(A,'function_handle') );
  error('gmres514: A must be either a matrix or a function handle.');
end
if (size(b,2) ~= 1)
  error('gmres514: b must be a column vector.');
end
n = size(b,1);
if (nargin<3 | isempty(x)),  x = zeros(n,1); end
if ( (n~= size(x,1)|1~=size(x,2)) |  (isreal(A) & (n~= size(A,1)|n~=size(A,2)) ) )
  error('gmres514: A, x, and b must have consistent dimensions.')
end
if (nargin<4 | isempty(m)), m = min(n,20); end
if nargin < 5,   M = []; end
if ( ~isreal(M) & ~isa(M,'function_handle') )
  error('gmres514: M must be either a matrix or a function handle.');
end
if (isreal(M) & ~isempty(M) & (n~= size(M,1)|n~=size(M,2)) )
  error('gmres514: M must have the same dimensions as A.');
end 
if (nargin<6 | isempty(tol)), tol = 1.e-6; end
if (nargin<7 | isempty(itmax)), itmax = min(n,2000); end
  
if isreal(A)
  r = b - A*x;
else
  r = b - A(x);
end
if ~isempty(M)
  if isreal(M) 
    r = M\r;
  else
    r = M(r);
  end
end

rho = norm(r);
atol = tol*rho; 
if atol<eps
  warning('gmres514: absolute stopping tolerance < eps; resetting to 100*eps.'); 
end

V = zeros(length(b),m+1);
R = zeros(m,m);
c = zeros(m,1);
s = c;
w = zeros(m+1,1);
ittot = 0;

itno = [0];
lrrnrm = [0];
rho0 = rho;

%
%  Top of outer loop.
%
while rho > atol 
  V(:,1) = r/rho;
  w(1) = rho;

  %
  % Top of inner loop. 
  %
  for k = 1:m
    ittot = ittot + 1;
    if ittot > itmax, 
      warning('gmres514: Maximum number of iterations reached without convergence.'); 
      break; 
    end 
    if isreal(A)
      V(:,k+1) = A*V(:,k);
    else
      V(:,k+1) = A(V(:,k));
    end
    if ~isempty(M)
      if isreal(M) 
        V(:,k+1) = M\V(:,k+1);
      else
        V(:,k+1) = M(V(:,k+1));
      end
    end
    for i = 1:k  
      R(i,k) = V(:,k+1)'*V(:,i);
      V(:,k+1) = V(:,k+1) - R(i,k)*V(:,i);
    end
    for i = 1:k-1
      temp = R(i,k);
      R(i,k) = c(i)*temp + s(i)*R(i+1,k);
      R(i+1,k) = -s(i)*temp + c(i)*R(i+1,k);
    end
    tempnorm = norm(V(:,k+1)); 
    temp = sqrt(R(k,k)^2 + tempnorm^2);
    c(k) = R(k,k)/temp;
    s(k) = tempnorm/temp;
    R(k,k) = temp; 
    w(k+1) = -s(k)*w(k);
    w(k) = c(k)*w(k);
    rho = abs(w(k+1));
    itno = [itno, ittot];
    lrrnrm = [lrrnrm, log10(rho0\rho)];
    if rho <= atol, break, end
    V(:,k+1) = V(:,k+1)/tempnorm; 
  end
  %
  % Bottom of the inner loop.
  %

  %
  % Form the solution.
  %
%    for i = k:-1:1 
%      w(i) = w(i)/R(i,i);
%      if i>1
%        w(1:i-1) = w(1:i-1) - w(i)*R(1:i-1,i);
%      end
%    end
  w(1:k) = R(1:k,1:k)\w(1:k);
  x = x + V(:,1:k)*w(1:k);
  if ittot > itmax, break, end 
  if isreal(A)
    r = b - A*x;
  else
    r = b - A(x);
  end
  if ~isempty(M)
    if isreal(M) 
      r = M\r;
    else
      r = M(r);
    end
  end
  rho = norm(r); 
end
%
% Bottom of the outer loop.
%
