function [xv,lmb,iresult] = sptarn(A,B,lb,ub, varargin)
%SPTARN   Solve generalized sparse eigenvalue problem.
%
%       [XV,LMB,IRESULT] = SPTARN(A,B,LB,UB,SPD,TOLCONV,JMAX,MAXMUL)
%       Finds eigenvalues of pencil (A-LMB*B)*XV=0 in interval LB < LMB <= UB.
%
%       Input data:
%       A, B sparse matrices
%       LB lower bound if -Inf all eigenvalues to the left of UB sought
%       UB upper bound if +Inf all eigenvalues to the right of LB sought
%       One of LB, UB must be finite. Narrower interval gives faster result.
%       In complex case real parts of LMB are compared to LB and UB.
%
%       Output results:
%       XV eigenvectors, ordered so that NORM(A*XV-B*XV*DIAG(LMB)) small
%       LMB eigenvalues in interval, sorted
%       IRESULT flag ABS(IRESULT) number of eigenvalues found
%       IRESULT =>0 successful termination, all eigenvalues in interval found
%       IRESULT < 0 not yet successful, there may be more eigenvalues. Try
%       again with a smaller interval!
%
%       Parameters that may be set, but have defaults:
%       SPD   true if pencil known to be symmetric pos definite
%       default false
%
%       TOLCONV   expected relative accuracy
%       default 100*eps machine precision
%
%       JMAX  maximum number of basis vectors. You need JMAX*N working space
%       so a small value may be justified on a small computer, otherwise
%       let it be the default value JMAX=100. Normally the algorithm stops
%       earlier when enough eigenvalues have converged.
%
%       MAXMUL number of Arnoldi runs tried. Must at least be as large
%       as maximum multiplicity of any eigenvalue. If a small value of
%       JMAX is given, many Arnoldi runs are necessary.
%       Default value is MAXMUL=N, which will be needed when all the
%       eigenvalues of the unit matrix are sought!
%
%       Algorithm description:
%       Uses Arnoldi with spectral transformation. The shift is chosen
%       at UB, LB or at a random point in the interval (LB,UB) when both
%       bounds are finite.
%       Number of steps J in the Arnoldi run depends on how many
%       eigenvalues there are in the interval, but stops at J=MIN(JMAX,N).
%       After a
%       stop, restart to find more Schur vectors in orthogonal complement to
%       all those already found. When no more eigenvalues are found in
%       LB < LMB <= UB, algorithm stops. For small values of JMAX, several
%       restarts may be needed before a certain eigenvalue has converged,
%       algorithm works for JMAX > #eigenvalues in interval + 1, but then
%       many restarts are needed. For large values of JMAX, which is the
%       preferred choice, MUL+1 runs are needed, where MUL is the maximum
%       multiplicity of an eigenvalue in the interval.
%         Works on nonsymmetric as well as symmetric pencils, but then accuracy
%       is approximately TOL*(Henrici departure from normality). The
%       parameter SPD is used only to choose between SYMAMD and COLAMD when
%       factorizing, the former being marginally better for symmetric matrices
%       close to the lower end of the spectrum.
%        case of trouble:
%       If convergence is too slow, try (in this order of priority)
%       1) a smaller interval LB,UB
%       2) a larger JMAX
%       3) a larger MAXMUL.
%         If factorization fails, try again with LB and UB finite.
%       Then shift is chosen at random and hopefully not at an eigenvalue.
%       If it fails again, check whether pencil is singular.
%         If it goes on forever, there may be too many eigenvalues in the
%       strip. Try with a small value MAXMUL=2 and see which eigenvalues you
%       get! Those you get are some of the evs, but a negative iresult tells
%       you that you have not gotten them all.
%         If memory overflows, try smaller JMAX.
%         Algorithm is designed for eigenvalues close to the real axis. If you
%       want those close to the imaginary axis, try A=i*A!
%         When SPD is true, shift is at LB so that advantage is taken of the
%       faster factorization for symmetric positive definite matrices.
%       No harm done but slower execution if LB above lowest eigenvalue.

%       Copyright 1994-2017 The MathWorks, Inc.

if nargin > 4
    [varargin{:}] = convertStringsToChars(varargin{:});
end

cpustart=cputime;
n=size(A,1);

parser = inputParser;
parser.addRequired('A', @isnumeric);
parser.addRequired('B', @isnumeric);
parser.addRequired('lb', @isnumeric);
parser.addRequired('ub', @isnumeric);
parser.addOptional('spd', 0, @isnumeric);
parser.addOptional('tolconv', 100*eps, @isnumeric);
parser.addOptional('jmax', min(100,n), @isnumeric);
parser.addOptional('maxmul', n, @isnumeric);
parser.addParamValue('Stats', 'on', @pdeIsValid.onOrOff);
parser.parse(A,B,lb,ub, varargin{:});
spd = parser.Results.spd;
tolconv = parser.Results.tolconv;
jmax = parser.Results.jmax;
maxmul = parser.Results.maxmul;
stats = strcmpi(parser.Results.Stats, 'on');


% Choose shift
if lb==-inf && ub==inf
  % No info on eigenvalue location, pick shift in the dark
  shift=rand(1);
elseif lb==-inf
  shift=ub;
elseif ub==inf
  shift=lb;
elseif spd
  shift=lb;
else
  ranv=rand(2,1);
  shift=[lb ub]*ranv/sum(ranv);
end

% Never more than n basis vectors!
jmax=min(jmax,n);

% Factorize and check if problem it not too singular
ntries=0;
maxtries=2;
while ntries<=maxtries

  C=A-shift*B;
  if spd
    [R,p,pr]=chol(C, 'vector');
    if p==0
      L=R';
      U=R;
      pc = pr;
    else
      %      fprintf('Problem not symmetric positive definite\n');
      [L, U, pr, pc] = lu(C, 'vector');
    end
  else
    [L, U, pr, pc] = lu(C, 'vector');
  end

  % Check for singularity
  v=U\(L\rand(n,1));
  fact=norm(v,inf)*norm(C,inf);

  if fact>0.01/eps
    % Try another shift
    if lb==-inf && ub==inf
      % No info on eigenvalue location, pick shift in the dark
      shift=rand(1,1);
    elseif lb==-inf
      if ub==0
        % No info here either
        shift=-rand(1);
      else
        shift=ub-abs(ub)*rand(1);
      end
    elseif ub==inf
      if lb==0
        % No info here either
        shift=rand(1);
      else
        shift=lb+abs(lb)*rand(1);
      end
    else
      ranv=rand(2,1);
      shift=[lb ub]*ranv/sum(ranv);
    end
  else
    clear C;
    B=B(pr,pc);
    break;
  end

  ntries=ntries+1;
end

if ntries>maxtries
  error(message('pde:sptarn:EigenvalueShiftNotFound'))
end

% tol is tolerance for convergence
% tolspur tol to determine spurious vectors
% tolstab decides when eigenvalue appr is stabilized
tol=tolconv;
tolspur=100*tolconv;
tolstab=1e-5;
% check for convergence first after jf steps, then every js steps
% The larger n the cheaper is the test compared to an extra step
js=ceil(4000/n);
jf=10;
%
normhj=0;
%
% nt is total number of converged values, both inside and outside interval
% US is n*nt matrix of Schur vectors
% T is nt*nt upper triangular Schur form
nt=0;
US=[];
T=[];
%
% Starting vector
vstart=randn(n,1);
terminatealg=0;
% Finally, it is time to start
fprintf('\n')
for mul=1:maxmul
  % Starting vector orthogonalized against already converged
  v=[US pdegrmsc(vstart,US)];
  h=[T;zeros(1,nt)];
  convrun=0;
  % First testing point
  jt=min(nt+jf,jmax);
  for j=nt+1:jmax
    r=U\(L\(B*v(:,j)));
    [vj1,hj]=pdegrmsc(r,v);
    h=[h hj(1:j)];
    normhj=max(normhj,norm(hj));
    % Test if linear dependence, signals decrease of band width or
    % convergence
    if hj(j+1) < tol*normhj
      convrun=1;
      jt=j;
    else
      v=[v vj1];
    end
    h=[h;zeros(1,j-1) hj(j+1)];

    % One vector added, is it time to test?
    if j==jt
      iv=nt+1:jt;
      [us,ts]=schur(h(iv,iv));
      [uc,tc]=rsf2csf(us,ts);
      lms=shift+1 ./diag(tc);
      sij=abs((h(j+1,j)/normhj)*us(j-nt,:))';

      % Test if it is time to stop
      convoutside=(lms<=lb|lms>ub) & sij<=tolstab;
      nonconvinside= lb<lms & lms<=ub & sij>tol & abs(us(1,:))'>tolspur;
      convrun= convrun | (any(convoutside) & ~any(nonconvinside));
% Report progress if wanted
   if any(nonconvinside)
       incinside=find(nonconvinside);
       iintr=incinside(1);
   else
       iintr=1;
   end
   if(stats)
     fprintf('              Basis=%3.0f,  Time=%7.2f,  New conv eig=%3.0f\n',j,cputime-cpustart,sum(sij<=tol));
     %   fprintf(1,'%4.0f %7.2f %3.0f ',j, cputime-cpustart,sum(sij<=tol));
     %   fprintf(1,'%12.5e +i%12.5e %10.3e \n',real(lms(iintr)),imag(lms(iintr)), sij(iintr));
   end % report part
      jt=min(jt+js,jmax);
    end
    if convrun
      break;
    end
  end % for j=nt+1:jmax,

  % nc is number of converged eigenvalues above the first nonconverged along h
  nc=sum(cumprod(double(sij<=tol)));
  if(stats)
    % Report progress of sweep
    fprintf('End of sweep: Basis=%3.0f,  Time=%7.2f,  New conv eig=%3.0f\n',j,cputime-cpustart,nc);
  end %report progress
  if nc>0
    T=[T h(1:nt,iv)*us(:,1:nc); zeros(nc,nt) ts(1:nc,1:nc)];
    US=[US v(:,iv)*us(:,1:nc)];
    nt=nt+nc;
    if nt>=jmax
      terminatealg=1;
    end
  end
  newinside=lb<lms & lms<=ub;
  if nc>0 && any(newinside & sij<=tol)
    vstart=randn(n,1);
  elseif any(newinside)
    vstart=v(:,iv)*us(:,find(newinside, 1 ));
  else
    terminatealg=1;
  end
  if terminatealg
    break;
  end
end % for mul=1:maxmul

% Triangular matrix is complete. Get eigensolution!
%[s,dt]=eig(T,'nobalance');
[s,dt]=eig(T.*(abs(T)>tol*normhj),'nobalance');
ev=shift+1./diag(dt);

% Deliver eigenvalues in interval as result
inside=find(lb<ev & ev<=ub);
iresult=length(inside);
[lmb,ievsrt]=sort(ev(inside));
xv=US*s(:,inside(ievsrt));
if size(xv,2)>0
  xv(pc,:)=xv;
end

if ~terminatealg
  iresult=-iresult;
end

