  function [W,R,res] = arnRD(A,m,v,nev,order,tol,maxits)
% function [W,R,res] = arnR(A,v,m,nev,order,tol,maxits)
% restarted  Arnoldi process with explicit deflation.
% m      = dimension of Krylov subspace
% order  = 1='descend' or -1= 'ascend' -- compute
%   eigenvalues of largest (resp. largest) real parts
% tol    = tolerance for stopping criterion
% maxits = max number of (oiuter) Arnoldi iterations
%
%  OUT:
%  W = sets of schur vectors 
%  R = related Schur form 
%  res  = residual norm history [outer iterations]
%%--------------------
 W = [];
%%R = zeros(m,m);
 R = []; 
 n = size(A,1);
 res = 1.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;
 its = 0;
 while (its < maxits) 
   its = its+1;
%%-------------------- anoldi loop 
  for i=1:m  
   i1 = i + 1  ; 
   v  = A*V(:,i)  ;
%%-------------------- Deflation step 
   if  (size(W,2) > 0) 
     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:m);
  Hm = H(1:m,1:m);
%%--------------------  get Ritz values and vectors
  [X,vals] = eig(Hm);
  vals = diag(vals);
  [t, k] = max(order*real(vals));
  v = V(:,1:m)*X(:,k) ;
  %%
  Av = A*v;
%%-------------------- nevc = # converged eigenvalues
  if  (nevc >1) 
    Av = Av - W*(R*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;
    end
    v = v / norm(v);
    W(:,nevc) = v;
    R(1:nevc,nevc) = W'*(A*v);
  end
  if (nevc >= nev);
    disp(' ----------converged ')
    [nevc, nev]
    return
  end
  %%-------------------- else restart
  V(:,1) = v;
end
