function [x, rel_res, iter, flag] = gmres_simple(A, x, b, M, restrt, max_it, tol)
%
%    A simplified program for GMRES, based on:
%    http://www.netlib.org/templates/matlab/gmres.m
%
%    Simplifed to make GMRES idea clearer.
%    Modified to do the least squares solve using Matlabs \ operator
%    at each iteration, instead of the Givens rotations.
%     
%
%    Faisal Saied
%    Computing Research Institute (CRI)
%    and Rosen Center for Adv. Computing (RCAC)
%    Purdue University
%
%    Oct 21, 2005
%

iter = 0;   flag = 0;                             % initialization

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

r = M \ (b-A*x);                                  % Matvec, preconditioning 
rel_res = norm(r) / bnrm2;
if (rel_res < tol) return, end

[n,n] = size(A);  m = restrt;                     % initialize workspace
V(1:n,1:m+1) = zeros(n,m+1);
H(1:m+1,1:m) = zeros(m+1,m);

e1 = zeros(n,1);   e1(1) = 1.0;

for iter = 1:max_it                               % begin iteration
   V(:,1) = r / norm(r);                          % v_1 for restart
   s = norm(r)*e1;                                % rhs for L.Sq. problem
   for j = 1:m                           % construct orthonormal
                                         % basis using modified Gram-Schmidt
      w = M \ (A*V(:,j));                % Matvec, preconditioning

      for i = 1:j 
        H(i,j)= w'*V(:,i);
        w = w - H(i,j)*V(:,i);
      end
      H(j+1,j) = norm(w);
      V(:,j+1) = w / H(j+1,j);

      Htilde = H(1:j+1,1:j); sj = s(1:j+1);
      y = Htilde \ sj;                             % Least squares solves 
                                                   % using Matlabs \ operator
      small_rel_res = norm(sj-Htilde*y)/bnrm2;
      fprintf(' iter = %d, j = %d, small rel res = %d \n', iter, j, small_rel_res);
      if (small_rel_res <= tol), break, end;
   end

   mj = length(y);
   x = x + V(:,1:mj)*y;                            % update approximation
   r = M \ (b-A*x);                                % compute res
   rel_res = norm(r) / bnrm2;                        % check convergence
   fprintf('rel res = %d \n', rel_res);
   if (rel_res <= tol), break, end;
end

if (rel_res > tol) flag = 1; end;                   % converged

% END of gmres