  function [W,R,res] = arnRD(A,v,m,nev,order,tol,maxits)
% function [W,R,res] = arnR(A,v,m,nev,order,tol,maxits)
% restarted  Arnoldi process with deflation.
% m = dimension of Krylov subspace
% order = 'ascend' or 'descend' -- computes eigenvalues of
% smallest (resp. largest)  real parts
% tol   = tolerance for stopping criterion
% maxits= max number of (inner) Arnoldi steps;
%
%  OUT:
%  W = sets of schur vectors 
%  R = related Schur form 
%  res  = residual norm history
%%--------------------
W = [];
%%R = zeros(m,m);
R = []; 
 n = size(A,1);
 res = 1.0; its = 0; tol1 = 0.0;
%%-------------------- restart loop 
 nevc = 0;
 ro = norm(v)  ; 
%%-------------------- if ro == 0 ...
 if (ro == 0), error('zero vector for v'), end
 V(:,1) = v/ro;
while (its*m < maxits) 
   its = its+1;
  %% restart 
  %%-------------------- anoldi loop 
  for i=1:m  
   i1 = i + 1  ; 
   v  = A*V(:,i)  ;
%%-------------------- Deflation step 
  if  (nevc >1) 
     v = v - W*(R*W'*v);
   end
   V(:,i1) = v;
   %%-------------------- 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),break,end
   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 and vectors
  [X,vals] = eig(Hm);
  vals = diag(vals);
  [t, k] = max(order*real(vals));
  v = V(:,1:i)*X(:,k) ;
  %%
  Av = A*v;
%%-------------------- nevc = number of converged eigenvalues --   
  if  (nevc >1) 
    %%%%% THE R IS NOT NEEDED
    %%   Av = Av - W*(R*W'*Av);
       Av = Av - W*(W'*Av);
  end
%%--------------------
  lam = v'*Av;
  ro = norm(Av - lam*v);
  if (its == 1) 
    tol1 = tol*res;
  end
  res(its)  = ro;
  if (ro < tol1) 
    nevc = nevc+1;
    if (nevc >1) 
      t = W'*v;
      v = v - W*t;
      v = v / norm(v);
    end
    nevc
    W(:,nevc) = v;
    R(1:nevc,nevc) = W'*(A*v);
  end
  if (nevc >= nev);
    return
  end
  %%-------------------- else restart
  V(:,1) = v;
end
