function [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M1,M2,x0,varargin)
%GMRES  Generalized Minimum Residual Method.
%   X = GMRES(A,B) and X = GMRES(A,B,RESTART) attempt to solve the system
%   of linear equations A*X = B for X.  The coefficient matrix A must be
%   square and the right hand side (column) vector B must have length N,
%   where A is N-by-N.  GMRES will start iterating from an initial guess
%   which by default is an all zero vector of length N.  GMRES will
%   restart itself at iteration RESTART using that iterate as its new
%   initial guess.  The default value of RESTART is N, in which case the
%   method never actually restarts.  Iterates are produced until the
%   method either converges, fails, or has computed the maximum number of
%   iterations.  Convergence is achieved when an iterate X has relative
%   residual NORM(B-A*X)/NORM(B) less than or equal to the tolerance of
%   the method.  The default tolerance is 1e-6.  The default maximum
%   number of iterations is the minimum of N/RESTART and 10.  No
%   preconditioning is used.
%
%   If GMRES converged, then a message to that effect is displayed.  
%   If GMRES failed to converge after the maximum number of iterations
%   or halted for any reason, then a warning message is printed
%   displaying the relative residual NORM(B-A*X)/NORM(B) and the
%   iteration number at which the method stopped or failed.
%
%   GMRES(A,B,RESTART,TOL) specifies the tolerance of the method.
%
%   GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of iterations.
%
%   GMRES(A,B,RESTART,TOL,MAXIT,M) and GMRES(A,B,RESTART,TOL,MAXIT,M1,M2)
%   use (left) preconditioner M or M=M1*M2 and effectively solve the
%   system inv(M)*A*X = inv(M)*B for X.  If any preconditioner is given as
%   the empty matrix ([]) then that preconditioner is considered to be the
%   identity matrix.  Since systems of the form M*Y = R are solved for Y
%   using \ (backslash), it is wise to factor a non-triangular
%   preconditioner M into its LU factors i.e. replace
%   GMRES(A,B,RESTART,TOL,MAXIT,M) with [M1,M2] = lu(M);
%   GMRES(A,B,RESTART,TOL,MAXIT,M1,M2).
%
%   GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial
%   guess.  If X0 is given as the empty matrix ([]) then the default all
%   zero vector is used.
%
%   Example:
%      load west0479
%      A = west0479;
%      b = sum(A,2);        % the true solution is the vector of all ones
%      x = gmres(A,b,5);
%   does not converge to a tolerance of 1e-6 within 20 iterations.
%   Try an incomplete factorization of A as a preconditioner
%      [L,U] = luinc(A,1e-6);
%      x1 = gmres(A,b,5,1e-15,10,L,U);
%   Now GMRES(5) will converge to within the requested tolerance of 1e-15
%   with a suitable preconditioner.
%
%   For more information and examples on the other outputs and on passing
%   functions as arguments, try: more on, type gmres, more off
%
%   See also BICG, BICGSTAB, CGS, PCG, QMR, SLASH, LUINC.

%   [X,FLAG] = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) returns a solution
%   X and a FLAG which describes the convergence of GMRES.  If FLAG is
%    0 then GMRES converged to the desired tolerance TOL within MAXIT
%      iterations without failing for any reason.
%    1 then GMRES iterated MAXIT times but did not converge.
%    2 then a system of equations of the form M*Y = R was ill-conditioned and
%      did not return a useable result when solved by \ (backslash).
%    3 then GMRES stagnated (two consecutive iterates were the same).
%   Whenever FLAG is not 0, the solution X returned is that with
%   minimal norm residual computed over all the iterations.  No
%   messages are displayed if the FLAG output is specified.
%
%   [X,FLAG,RELRES] = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) also returns
%   the relative residual NORM(B-A*X)/NORM(B).  If FLAG is 0, RELRES <= TOL.
%
%   [X,FLAG,RELRES,ITER] = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) also
%   returns both the outer and inner iteration numbers at which X was
%   computed: 0 <= ITER(1) <= MAXIT and 0 <= ITER(2) <= RESTART.
%
%   [X,FLAG,RELRES,ITER,RESVEC] = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) also
%   returns a vector of the residual norms at each inner iteration, starting
%   from RESVEC(1) = NORM(B-A*X0).  If FLAG is 0 and ITER = [I J], then RESVEC
%   is of length (I-1)*RESTART+J+1 and RESVEC(END) <= TOL*NORM(B).
%
%   When A is not explicitly available as a matrix, express A as an operator
%   that accepts vector input X and returns the matrix-vector product A*X.
%   Similarly the preconditioner operators should return the result of
%   solving M\X.
%   These operators may be functions, string expressions or inline objects.
%
%   Extra inputs to the functions or inline objects representing the matrix
%   and the preconditioner operators are passed as extra arguments to pcg.
%
%   Example:
%      A = gallery('wilk',21);  b = sum(A,2);
%      tol = 1e-12;  maxit = 15;
%      M1 = diag([10:-1:1 1 1:10]); M2 = []; x0 = [];
%      x = gmres(A,b,10,tol,maxit,M1,M2,x0);
%   Alternatively, use this matrix-vector product function
%      function y = afun(x,n)
%      y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'] .*x + [x(2:n); 0];
%   and this preconditioner backsolve function
%      function y = mfun(r,n)
%      y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%   as inputs to gmres
%      x1 = gmres('afun',b,10,tol,maxit,'mfun',M2,x0,21);
%   Note that both afun and mfun had to accept gmres's extra input n=21.

%   Penny Anderson, 1996.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.13 $ $Date: 1998/10/21 20:51:36 $

if (nargin < 2)
   error('Not enough input arguments.');
end

% Determine whether A is a matrix, a string expression,
% the name of a function or an inline object.
[atype,afun,afcnstr] = iterchk(A);
if isequal(atype,'matrix')
   % Check matrix and right hand side vector inputs have appropriate sizes
   [m,n] = size(A);
   if (m ~= n)
      error('Matrix must be square.');
   end
   if ~isequal(size(b),[m,1])
      es = sprintf(['Right hand side must be a column vector of' ...
            ' length %d to match the coefficient matrix.'],m);
      error(es);
   end
else
   m = size(b,1);
   n = m;
   if (size(b,2) ~= 1)
      error('Right hand side must be a column vector.');
   end
end

% Assign default values to unspecified parameters
if nargin < 3
   restart = n;
end
if nargin < 4
   tol = 1e-6;
end
if nargin < 5
   maxit = min(ceil(n/restart),10);
end

% Check for all zero right hand side vector => all zero solution
n2b = norm(b);                      % Norm of rhs vector, b
if (n2b == 0)                       % if    rhs vector is all zeros
   x = zeros(n,1);                  % then  solution is all zeros
   flag = 0;                        % a valid solution has been obtained
   relres = 0;                      % the relative residual is actually 0/0
   iter = [0 0];                    % no iterations need be performed
   resvec = 0;                      % resvec(1) = norm(b-A*x) = norm(0)      
   if (nargout < 2)
      itermsg('gmres',tol,maxit,0,flag,iter,NaN);
   end
   return
end

if ((nargin >= 6) & ~isempty(M1))
   existM1 = 1;
   [m1type,m1fun,m1fcnstr] = iterchk(M1);
   if isequal(m1type,'matrix')
      if ~isequal(size(M1),[m,m])
         es = sprintf(['Preconditioner must be a square matrix' ...
               ' of size %d to match the problem size.'],m);
         error(es);
      end      
   end   
else
   existM1 = 0;
   m1type = 'matrix';
end

if ((nargin >= 7) & ~isempty(M2))
   existM2 = 1;
   [m2type,m2fun,m2fcnstr] = iterchk(M2);
   if isequal(m2type,'matrix')
      if ~isequal(size(M2),[m,m])
         es = sprintf(['Preconditioner must be a square matrix' ...
               ' of size %d to match the problem size.'],m);
         error(es);
      end
   end
else
   existM2 = 0;
   m2type = 'matrix';
end

if ((nargin >= 8) & ~isempty(x0))
   if ~isequal(size(x0),[n,1])
      es = sprintf(['Initial guess must be a column vector of' ...
            ' length %d to match the problem size.'],n);
      error(es);
   end
else
   x0 = zeros(n,1);
end
x = x0;

if ((nargin > 8) & isequal(atype,'matrix') & ...
      isequal(m1type,'matrix') & isequal(m2type,'matrix'))
   error('Too many input arguments.');
end

% Set up for the method
flag = 1;
xmin = x;                          % Iterate which has minimal residual so far
imin = 0;                          % "Outer" iteration at which xmin was computed
jmin = 0;                          % "Inner" iteration at which xmin was computed
tolb = tol * n2b;                  % Relative tolerance
if isequal(atype,'matrix')
   r = b - A * x;                  % Zero-th residual
else
   r = b - iterapp(afun,atype,afcnstr,x,varargin{:});
end
normr = norm(r);                   % Norm of residual

if (normr <= tolb)                 % Initial guess is a good enough solution
   flag = 0;
   relres = normr / n2b;
   iter = [0 0];
   resvec = normr;
   if (nargout < 2)
      itermsg('gmres',tol,maxit,[0 0],flag,iter,relres);
   end
   return
end

resvec = zeros(restart*maxit+1,1); % Preallocate vector for norm of residuals
resvec(1) = normr;                 % resvec(1) = norm(b-A*x0)
normrmin = normr;                  % Norm of residual from xmin
rho = 1;
stag = 0;                          % stagnation of the method

% loop over maxit "outer" iterations (unless convergence or failure)

for i = 1 : maxit
   V = zeros(n,restart+1);          % Arnoldi vectors
   h = zeros(restart+1,1);          % upper Hessenberg st A*V = V*H ...
   QT = zeros(restart+1,restart+1); % orthogonal factor st QT*H = R
   R = zeros(restart,restart);      % upper triangular factor st H = Q*R
   f = zeros(restart,1);            % y = R\f => x = x0 + V*y
   W = zeros(n,restart);            % W = V*inv(R)
   j = 0;                           % "inner" iteration counter
   
   if existM1
      if isequal(m1type,'matrix')
         vh1 = M1 \ r;
      else
         vh1 = iterapp(m1fun,m1type,m1fcnstr,r,varargin{:});
      end
      if isinf(norm(vh1,inf))
         flag = 2;
         break
      end
   else
      vh1 = r;
   end
   if existM2
      if isequal(m2type,'matrix')
         vh = M2 \ vh1;
      else
         vh = iterapp(m2fun,m2type,m2fcnstr,vh1,varargin{:});
      end
      if isinf(norm(vh,inf))
         flag = 2;
         break
      end
   else
      vh = vh1;
   end
   h(1) = norm(vh);
   V(:,1) = vh / h(1);
   QT(1,1) = 1;
   phibar = h(1);
   
   for j = 1 : restart
      if isequal(atype,'matrix')
         u2 = A * V(:,j);
      else
         u2 = iterapp(afun,atype,afcnstr,V(:,j),varargin{:});         
      end
      if existM1
         if isequal(m1type,'matrix')
            u1 = M1 \ u2;
         else
            u1 = iterapp(m1fun,m1type,m1fcnstr,u2,varargin{:});
         end
         if isinf(norm(u1,inf))
            flag = 2;
            break
         end
      else
         u1 = u2;
      end
      if existM2
         if isequal(m2type,'matrix')
            u = M2 \ u1;
         else
            u = iterapp(m2fun,m2type,m2fcnstr,u1,varargin{:});
         end
         if isinf(norm(u,inf))
            flag = 2;
            break
         end
      else
         u = u1;
      end
      for k = 1 : j
         h(k) = V(:,k)' * u;
         u = u - h(k) * V(:,k);
      end
      h(j+1) = norm(u);
      V(:,j+1) = u / h(j+1);
      R(1:j,j) = QT(1:j,1:j) * h(1:j);
      rt = R(j,j);
      
      % find cos(theta) and sin(theta) of Givens rotation
      if h(j+1) == 0
         c = 1.0;                      % theta = 0
         s = 0.0;
      elseif abs(h(j+1)) > abs(rt)
         temp = rt / h(j+1);
         s = 1.0 / sqrt(1.0 + abs(temp)^2); % pi/4 < theta < 3pi/4
         c = - temp * s;
      else
         temp = h(j+1) / rt;
         c = 1.0 / sqrt(1.0 + abs(temp)^2); % -pi/4 <= theta < 0 < theta <= pi/4
         s = - temp * c;
      end
      R(j,j) = c' * rt - s' * h(j+1);
      %   zero = s * rt + c * h(j+1);
      
      q = QT(j,1:j);
      QT(j,1:j) = c' * q;
      QT(j+1,1:j) = s * q;
      QT(j,j+1) = -s';
      QT(j+1,j+1) = c;
      f(j) = c' * phibar;
      phibar = s * phibar;
      
      if j < restart
         if f(j) == 0                 % stagnation of the method
            stag = 1;
         end
         W(:,j) = (V(:,j) - W(:,1:j-1) * R(1:j-1,j))/ R(j,j);
         % Check for stagnation of the method
         if stag == 0
            stagtest = zeros(n,1);
            ind = (x ~= 0);
            if ~(i==1 & j==1)
               stagtest(ind) = W(ind,j) ./ x(ind);
               stagtest(~ind & W(:,j) ~= 0) = Inf;
               if abs(f(j))*norm(stagtest,inf) < eps
                  stag = 1;
               end
            end
         end
         x = x + f(j) * W(:,j);       % form the new inner iterate
      else % j == restart
         vrf = V(:,1:j)*(R(1:j,1:j)\f(1:j));
         % Check for stagnation of the method
         if stag == 0
            stagtest = zeros(n,1);
            ind = (x0 ~= 0);
            stagtest(ind) = vrf(ind) ./ x0(ind);
            stagtest(~ind & vrf ~= 0) = Inf;
            if norm(stagtest,inf) < eps
               stag = 1;
            end
         end
         x = x0 + vrf;                % form the new outer iterate
      end
      
      if isequal(atype,'matrix')
         normr = norm(b - A * x);
      else
         normr = norm(b - iterapp(afun,atype,afcnstr,x,varargin{:}));   
      end
      resvec((i-1)*restart+j+1) = normr;
      
      if normr <= tolb               % check for convergence
         if j < restart
            y = R(1:j,1:j) \ f(1:j);
            x = x0 + V(:,1:j) * y;     % more accurate computation of xj
            if isequal(atype,'matrix')
               r = b - A * x;
            else
               r = b - iterapp(afun,atype,afcnstr,x,varargin{:});
            end
            normr = norm(r);
            resvec((i-1)*restart+j+1) = normr;
         end
         if normr <= tolb             % check using more accurate xj
            flag = 0;
            iter = [i j];
            break
         end
      end
      
      if stag == 1
         flag = 3;
         break
      end
      
      if normr < normrmin            % update minimal norm quantities
         normrmin = normr;
         xmin = x;
         imin = i;
         jmin = j;
      end
   end                              % for j = 1 : restart
   
   if flag == 1
      x0 = x;                        % save for the next outer iteration
      if isequal(atype,'matrix')
         r = b - A * x0;
      else
         r = b - iterapp(afun,atype,afcnstr,x0,varargin{:});
      end
   else
      break
   end
   
end                                % for i = 1 : maxit

% returned solution is that with minimum residual
if flag == 0
   relres = normr / n2b;
else
   x = xmin;
   iter = [imin jmin];
   relres = normrmin / n2b;
end

% truncate the zeros from resvec
if flag <= 1 | flag == 3
   resvec = resvec(1:(i-1)*restart+j+1);
else
   if j == 0
      resvec = resvec(1:(i-1)*restart+1);
   else
      resvec = resvec(1:(i-1)*restart+j);
   end
end

% only display a message if the output flag is not used
if nargout < 2
   itermsg(sprintf(['gmres(%d)'],restart),tol,maxit,[i j],flag,iter,relres);
end
