function [x, resnrm ] = gmres( A, x0, b, maxit, tol )

% [x, resnrm ] = gmres( A, x, b, maxit, tol )
%
% gmres.m solves the linear system Ax=b
% using the Generalized Minimal residual (GMRES) method.
%
% input   A        REAL nonsymmetric positive definite matrix
%         x0       REAL initial guess vector
%         b        REAL right hand side vector
%         maxit    INTEGER maximum number of iterations
%         tol      REAL error tolerance
%
% output  x        REAL solution vector
%         rsnrm    REAL residual norm

   iter = 0;                                         % initialization
   flag = 0;
   resnrm = zeros(maxit+1,1);

   bnrm2 = norm( b );
   if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end
   tol = tol*bnrm2;

   x = x0;
   r = ( b-A*x );
   resnrm(1) = norm( r );
   if ( resnrm(1) < tol ) return, end

   [n,n] = size(A);                                  % initialize workspace
   V(1:n,1:maxit+1) = zeros(n,maxit+1);
   H(1:maxit+1,1:maxit) = zeros(maxit+1,maxit);
   e1    = zeros(maxit+1,1);
   e1(1) = 1.0;

   s = norm( r )*e1;
   V(:,1) = r / norm( r );

   for i = 1:maxit,                                  % begin iteration

       w = A*V(:,i);                                 % basis using Gram-Schmidt
       for k = 1:i,
           H(k,i)= w'*V(:,k);
           w = w - H(k,i)*V(:,k);
       end
       H(i+1,i) = norm( w );
       V(:,i+1) = w / H(i+1,i);
       
       y = H(1:i+1,1:i)\s(1:i+1);                    % solve projected system 
       resnrm(i+1) = norm( H(1:i+1,1:i)*y - s(1:i+1) );
       if ( resnrm(i+1) <= tol ),                    % update approximation
          resnrm = resnrm(1:i+1);
	  break;
       end
   end
   x = x + V(:,1:i)*y;


% END of gmres.m
