function varargout = eigsb(A, nev, m, tol, maxiter)
%eigsb executes the Arnoldi process with implicit restarts to
%compute a specified number of the largest eigenvalues of the sparse
%nonsymmetric matrix A.
%
%   [V D] = eigsb(A,v0,nev,m,tol,maxiter) returns a matrix V of
%   eigenvectors and diagonal matrix D of eigenvalues for the sparse
%   nonsymmetric A. m is the size of the Krylov subspace to be
%   used and should be as small as possible. nev is the number of
%   eigenvalues desired and should also be small. The error tolerance
%   tol defaults to 1.0e-6, maxiter defaults to 100.
%
%   E = eigsb(A, nev, m, tol, maxiter) - assigns E a column
%   vector containing the eigenvalues. The default values of tol and
%   maxiter are the same as the previous calling sequence.

% initialization
n = size(A, 1);

if nargin < 3
	error('eigssymb requires at least 3 arguments.');
end

% take care of defaults
if nargin == 3
	tol = 1.0e-6;
	maxiter = 100;
elseif nargin == 4
	maxiter = 100;
end

% we use deflation as k increases
k = 1;
I = eye(n);

V = zeros(n, m);
H = zeros(m, m);
% start with random f
f = randn(n,1);
f = f/norm(f);
%peform an m-step Arnoldi factorization
[V, H, f] = arnoldif(A,V,H,f,k,m);

% we only allow maxiter iterations
iter = 0;                                                               

while true
   iter = iter + 1;

% find eigenvalues and eigenvectors for H.
	[UH, DH] = eig(H);  
	sigma = diag(DH);

% reorder the eigenvalues and the corresponding eigenvectors.
	[~,ind] = sort(abs(sigma),1,'descend');
   sigma = sigma(ind);

	VHtmp = UH;
	for j =1:nev
      UH(:,j) = VHtmp(:, ind(j));
   end
   
   Q = eye(m);
   % use sigma(nev+1), ..., sigma(m) as the shifts for
   % the QR iteration
   j = m;
   while j >= nev+1
      lambda = sigma(j);
      alpha = imag(lambda);
      if abs(alpha) > 0 
         % lambda is complex. use a implicit double shift
         beta = real(lambda);
         [Qj,~] = qr((H - beta*eye(m))^2 + alpha^2*eye(m));
         j = j - 2;
         H = Qj'*H*Qj;
      else
         % lambda is real. use an implicit single shift.

         % NOTE: If you don't want to use chase modified for this
         % iteration, the same thing can be accomplished with the
         % statements
      
         %       [Qj,~] = qr(H - lambda*eye(m));
         %       j = j - 1;
         %       H = Qj'*H*Qj;
         [Qj,H] = chase2(H,lambda);
         j = j-1;
     end
      Q = Q*Qj;
   end
       
   % compute the residual norm for the kth eigenpair
   u = UH(:,k);
   resid_norm = norm(f)*abs(u(m));
   % lock v_k if the tolerance is obtained
   if(resid_norm < tol)
      if(k < nev)
         k = k+1;
      else
         break;
      end
   end
    
   % build an m-step factorization from the nev step one
   betak = H(nev+1,nev);
   sigmak = Q(m,nev);
   fk = V(:,nev+1)*betak + f*sigmak;
   V(:,1:nev) = V(:,1:m)*Q(:,1:nev);
   [V, H, f] = arnoldif(A,V(:,1:nev),H(1:nev,1:nev),fk,nev+1,m);
   
   if(iter >= maxiter)
      error('Cannot compute the requested eigenvalues in the specfied number of iterations');
   end
end

if nargout == 2
   V(:,1:nev) = V(:,1:m)*UH(:,1:nev);
   varargout{1} =  V(:,1:nev);
   varargout{2} =  diag(sigma(1:nev));
else
   varargout{1} = sigma(1:nev);
end
