function [x, error, iter] = gmres_m(x, b, atv, errtol, nstart, maxit, varargin)
% GMRES(m) linear equation solver
% based on a Matlab script by C. T. Kelley, July 10, 1994
% function [x, error, iter] = gmres_m(x, b, atv, errtol, nstart, maxit, varargin)
% input parameters:
% x      initial iterate
% b      right hand side
% atv    character name of a matrix-vector product routine returning x+(b-Ax)
%        when x is input. The format for atv is "function x = atv(x,b,p1,p2,...)"
%        where p1 and p2 are optional parameters
% errtol relative residual reduction factor
% nstart restarts the method every nstart iterations
% maxit  maximum number of iterations
% varargin optional parameters (p1,p2,...) for atv
% output parameters:
% x      solution of the linear system
% error  vector of residual norms for the history of the iteration
% iter   number of iterations
% (c) 2007 Alain Hebert, Ecole Polytechnique de Montreal
  n=length(b);
  errtol=errtol*norm(b);
  rho=Inf; error=[ ];
%
% global GMRES(m) iteration
  iter=0;
  while((rho > errtol) && (iter < maxit))
    r=feval(atv,x,b,varargin{:})-x ;
    rho=norm(r); error=[error;rho];
%
%   test for termination on entry
    if(rho < errtol), return, end
%
    h=zeros(nstart); v=zeros(n,nstart+1);
    c=zeros(nstart+1,1); s=zeros(nstart+1,1);
    g=rho*eye(nstart+1,1);
    v(:,1)=r/rho;
%
%   GMRES(1) iteration
    k=0;
    while((rho > errtol) && (k < nstart) && (iter < maxit))
      k=k+1;
      iter=iter+1;
      gar=feval(atv,v(:,k),zeros(n,1),varargin{:}) ;
      v(:,k+1)=v(:,k)-gar(:) ;
%
%     Modified Gram-Schmidt
      for j=1:k
        h(j,k)=v(:,j)'*v(:,k+1);
        v(:,k+1)=v(:,k+1)-h(j,k)*v(:,j);
      end
      h(k+1,k)=norm(v(:,k+1));
%
%     Reorthogonalize
      for j=1:k
        hr=v(:,j)'*v(:,k+1);
        h(j,k)=h(j,k)+hr;
        v(:,k+1)=v(:,k+1)-hr*v(:,j);
      end
      h(k+1,k)=norm(v(:,k+1));
%
%     watch out for happy breakdown 
      if(h(k+1,k) ~= 0)
        v(:,k+1)=v(:,k+1)/h(k+1,k);
      end
%
%     Form and store the information for the new Givens rotation
      if k > 1
        for i=1:k-1
          w1=c(i)*h(i,k)-s(i)*h(i+1,k); w2=s(i)*h(i,k)+c(i)*h(i+1,k);
          h(i:i+1,k)=[w1,w2];
        end
      end
      nu=norm(h(k:k+1,k));
      if nu~=0
        c(k)=h(k,k)/nu;
        s(k)=-h(k+1,k)/nu;
        h(k,k)=c(k)*h(k,k)-s(k)*h(k+1,k);
        h(k+1,k)=0;
        w1=c(k)*g(k)-s(k)*g(k+1); w2=s(k)*g(k)+c(k)*g(k+1);
        g(k:k+1)=[w1,w2];
      end
%
%     Update the residual norm
      rho=abs(g(k+1)); error=[error;rho];
    end
%
%   At this point either k > nstart or rho < errtol.
%   It's time to compute x and cycle.
    x=x+v(1:n,1:k)*(h(1:k,1:k)\g(1:k));
  end