  function [v,vals,res] = arnR(A,v,m,order,tol,maxits)
% function [v,vals,res] = arnR(A,v,m,order,tol,maxits)
% restarted  Arnoldi process 
% m = dimension of Krylov subspace
% order = -1 --> leftmost eigenvalue
%       = +1 --> rightmost eigenvalue 
% tol   = tolerance for stopping criterion
% maxits= max number of (inner) Arnoldi steps;
%
% OUT:
%  v = eigenvector
%  vals = eigenvalues of last hessenberg matrix
%  res  = residual norm history
%
%%--------------------
 n = size(A,1);
 res = 1.0; its = 0; tol1 = 0.0;
%%-------------------- restart loop 
while (its*m < maxits  & res > tol1) 
  its = its+1;
  ro = norm(v,2)  ; 
%%-------------------- if ro == 0 ...
  t = 1.0/ ro   ;
  V(:,1) = v*t; 
  v  = A*V(:,1) ;
  val  = v'*V(:,1) ;
  ro = norm(v - val*V(:,1)); 
  if (its == 1) 
    tol1 = tol*res;
  end
  res(its)  = ro;
%%-------------------- This is the Anoldi loop 
for i=1:m  
   i1 = i + 1  ; 
   V(:,i1) = A*V(:,i)  ;
%%-------------------- modified GS  ;
 for j=1:i 
   t = V(:,j)'*V(:,i1)  ;
   H(j,i) = t  ;
   V(:,i1) = V(:,i1) - t*V(:,j)  ; 
  end  ;
%
  t = norm(V(:,i1),2)  ;
  H(i1,i) = t  ;
  if (t  ~= 0.0)  
   t = 1.0 / t  ;
   V(:,i1) = V(:,i1)*t  ;  
  end  
%%-------------------- Arnoldi Objects 
  Vm = V(:,1:i);
  Hm = H(1:i,1:i);
%%-------------------- Get Ritz values
  [X,vals] = eig(Hm);
  vals = diag(vals);
%%-------------------- get index of wanted eigenvalue
  [t, k] = max(order*real(vals));
%%-------------------- restart with this vector:
  v = V(:,1:i)*X(:,k) ;
 end
end
