 function [sol,res]=gmres0(A,rhs,x0,tol,im,maxits) 
%function [sol,res]=gmres0(A,rhs,x0,tol,im,maxits) 
%-------------------------------------------------
% restarted gmres with Krylov subspace of dim = im 
% A      = matrix
% im     = krylov subspace dimension
% rhs    = right-hand side
% x0     = initial guess
% tol    = tolerance for stopping criterion
% maxits = max number of matrix - vector products 
%          allowed (= max tot. # of gmres steps)
%-------------------------------------------------
      tolmac = eps; 
      n = size(A,1);
      vv = zeros(n,im+1);
      its = 1 ;
      ierr = 0 ;
      sol = x0; 
%%
  while (its < maxits) 
      vv(:,1) = rhs - A*sol ;
      %      
      ro = norm(vv(1:n,1),2) ; 
      res(its) = ro;
      if (its  == 1) 
	tol1=tol*ro ;
      end 
      if (ro < tol1 | its > maxits)  
	return
      end
      t = 1.0/ ro ;
      vv(1:n,1) = vv(1:n,1) * t ;
%%---initialize 1-st term of rhs of hess.system..
      rs(1) = ro ;
      i = 0 ;
 while (i < im & (ro > tol1)  &  its < maxits)
      i=i+1 ;
      its = its + 1 ;
      i1 = i + 1 ;
      vv(:,i1) = A*vv(1:n,i);
%---------- modified GS ;
      for j=1:i
         t = vv(1:n,j)'*vv(1:n,i1) ;
         hh(j,i) = t ;
	 vv(1:n,i1) = 	 vv(1:n,i1) - t*vv(1:n,j) ;
      end  
      t = norm(vv(1:n,i1),2) ;
      hh(i1,i) = t ;
      if (t  ~= 0.0)  
        t = 1.0 / t ;
        vv(1:n,i1) = vv(1:n,i1)*t ;
      end  
%%
    if (i ~= 1) 
%---------- apply previous transf on i-th column h ;
       for k=2:i 
         k1 = k-1 ;
         t = hh(k1,i) ;
         hh(k1,i) = c(k1)*t + s(k1)*hh(k,i) ;
         hh(k,i) = -s(k1)*t + c(k1)*hh(k,i) ;
       end  
     end  %% IF 
%%
     gam = sqrt(hh(i,i)^2 + hh(i1,i)^2) ;
     gam = max(tolmac,gam);
%---------- get plane rotation + update rhs of LS Pb
     c(i) = hh(i,i)/gam ;
     s(i) = hh(i1,i)/gam ;
     rs(i1) = -s(i)*rs(i) ;
     rs(i) =  c(i)*rs(i) ;
%---------- get res. norm. and test for convergence
     hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i) ;
     ro = abs(rs(i1)) ;
     res(its) = ro ; 
   end 
%---------- solve upper triangular system
   rs(i) = rs(i)/hh(i,i) ;
   for  k=i-1:-1:1 ;
      t=rs(k) ;
      for j=k+1:i ;
         t = t-hh(k,j)*rs(j) ;
      end  
      rs(k) = t/hh(k,k) ;
   end  
%---------- form linear combination to get solution 
   for j=1:i 
      sol = sol +rs(j)* vv(1:n,j) ;
   end 
 end
