 function varargout = eigsbexplicit(A, x0, nev, m, tol, maxiter, reorthog)
%EIGSBEXPLICIT executes the Arnoldi process with explicit restarts to
%compute a specified number of the largest eigenvalues of the sparse
%nonsymmetric matrix A with whose nev largest eigenvalues in magnitude
%are real.
%
%   [V D] = eigsbexplicit(A,x0,nev,m,tol,maxiter) returns an orthogonal
%   matrix V of eigenvectors and diagonal matrix D of eigenvalues
%   for the sparse nonsymmetric A. x0 is the starting vector and
%   is normally random. 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, and default
%   value of reorthog is 1 (to reorthogonalize).
%
%   E = eigsb(A, x0, nev, m, tol, maxiter) - assigns E a column
%   vector containing the eigenvalues. The default values of tol,
%   maxiter, and reorthog are the same as the previous calling sequence.

% initialization
n = size(A, 1);

if nargin < 4
	error('eigssymb requires at least 4 arguments.');
end
if nargin == 4
	tol = 1.0e-6;
	maxiter = 100;
	reorthog = 1;
elseif nargin == 5
	maxiter = 100;
	reorthog = 1;
elseif nargin == 6
	reorthog = 1;
end

k = 1;
V = zeros(n, m+1);
H = zeros(m+1, m);
V(:, 1) = x0/norm(x0);

iter = 0;
number_computed = 0;

% restart Arnoldi or begin a new factorization
while true
	iter = iter + 1;
% execute m-k+1 steps of Arnoldi
	for j = k:m
		w = A*V(:,j);
		for i = 1:j
			H(i,j) = V(:,i)' * w;
			w = w - H(i,j)*V(:,i);
		end
	if reorthog == 1
       for p=1:j,
          tmp = V(:,p)'*w;
          w = w - tmp*V(:,p);
          H(p,j)=H(p,j)+tmp;
		 end
	end
		H(j+1,j) = norm(w);
		if H(j+1,j) == 0
			% very good! break out and find eigenvalues/eigenvectors
			break;
		end
		V(:,j+1) = w/H(j+1,j);
	end
	
	% find eigenvalues and eigenvectors for H(1:m,1:m)
	[UH, DH] = eig(H(1:m, 1:m));  
	d = diag(DH);

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

	VHtmp = UH;
	for j =1:nev
		UH(:,j) = VHtmp(:, ind(j));
		DH(j,j) = d(j);
	end

	% compute the Ritz vector.
	V(:,k) = V(:,1:m)*UH(:,k);
	   
%	orthogonalize V(:,k) with respect to the
%	previously estimated Ritz vectors.
	if(k > 1)
		V(:,k) = V(:,k) - V(:,1:(k-1))*(V(:,1:(k-1))'*V(:,k));
	end

  % normalize V(:,k)
	V(:,k) = V(:,k)/norm(V(:,k));

	% compute the residual norm for the kth eigenpair
	u = UH(:,k);
   resid_norm = H(m+1,m)*abs(u(m));

    % lock v_k if the tolerance is obtained
	if(resid_norm < tol)
		H(1:k,1:k) = V(:,1:k)'*A*V(:,1:k);
		number_computed = number_computed+ 1;
		if k < nev
			k = k+1;
      else
			break;
      end      
   end
   
   if(iter >= maxiter)
      break;
   end
end

if number_computed == 0
	error('No eigenvalues found.');
end

% setup the return values
if nargout == 2
	varargout{1} =  V(:,1:number_computed);
   varargout{2} =  diag(diag(H(1:number_computed,1:number_computed)));
else
    varargout{1} = diag(H(1:number_computed,1:number_computed));
end