[n,d,e,e2,w,ind,ierr,rv1]=imtqlv(n,d,e,e2,w,ind,ierr,rv1);

function [n,d,e,e2,w,ind,ierr,rv1]=imtqlv(n,d,e,e2,w,ind,ierr,rv1);
%***BEGIN PROLOGUE  IMTQLV
%***PURPOSE  Compute the eigenvalues of a symmetric tridiagonal matrix
%            using the implicit QL method.  Eigenvectors may be computed
%            later.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4A5, D4C2A
%***TYPE      SINGLE PRECISION (IMTQLV-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a variant of  IMTQL1  which is a translation of
%     ALGOL procedure IMTQL1, NUM. MATH. 12, 377-383(1968) by Martin and
%     Wilkinson, as modified in NUM. MATH. 15, 450(1970) by Dubrulle.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
%
%     This subroutine finds the eigenvalues of a SYMMETRIC TRIDIAGONAL
%     matrix by the implicit QL method and associates with them
%     their corresponding submatrix indices.
%
%     On INPUT
%
%        N is the order of the matrix.  N is an INTEGER variable.
%
%        D contains the diagonal elements of the symmetric tridiagonal
%          matrix.  D is a one-dimensional REAL array, dimensioned D(N).
%
%        E contains the subdiagonal elements of the symmetric
%          tridiagonal matrix in its last N-1 positions.  E(1) is
%          arbitrary.  E is a one-dimensional REAL array, dimensioned
%          E(N).
%
%        E2 contains the squares of the corresponding elements of E in
%          its last N-1 positions.  E2(1) is arbitrary.  E2 is a one-
%          dimensional REAL array, dimensioned E2(N).
%
%     On OUTPUT
%
%        D and E are unaltered.
%
%        Elements of E2, corresponding to elements of E regarded as
%          negligible, have been replaced by zero causing the matrix to
%          split into a direct sum of submatrices.  E2(1) is also set
%          to zero.
%
%        W contains the eigenvalues in ascending order.  If an error
%          exit is made, the eigenvalues are correct and ordered for
%          indices 1, 2, ..., IERR-1, but may not be the smallest
%          eigenvalues.  W is a one-dimensional REAL array, dimensioned
%          W(N).
%
%        IND contains the submatrix indices associated with the
%          corresponding eigenvalues in W -- 1 for eigenvalues belonging
%          to the first submatrix from the top, 2 for those belonging to
%          the second submatrix, etc.  IND is a one-dimensional REAL
%          array, dimensioned IND(N).
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          J          if the J-th eigenvalue has not been
%                     determined after 30 iterations.
%                     The eigenvalues should be correct for indices
%                     1, 2, ..., IERR-1.  These eigenvalues are
%                     ordered, but are not necessarily the smallest.
%
%        RV1 is a one-dimensional REAL array used for temporary storage,
%          dimensioned RV1(N).
%
%     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
%
%     Questions and comments should be directed to B. S. Garbow,
%     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
%     ------------------------------------------------------------------
%
%***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
%                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
%                 system Routines - EISPACK Guide, Springer-Verlag,
%                 1976.
%***ROUTINES CALLED  PYTHAG
%***REVISION HISTORY  (YYMMDD)
%   760101  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  IMTQLV
%
persistent b c f g gt i ii j k l m mml p r s s1 s2 tag ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
if isempty(ii), ii=0; end;
if isempty(mml), mml=0; end;
if isempty(tag), tag=0; end;
if isempty(gt), gt=0; end;
d_shape=size(d);d=reshape(d,1,[]);
e_shape=size(e);e=reshape(e,1,[]);
e2_shape=size(e2);e2=reshape(e2,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
rv1_shape=size(rv1);rv1=reshape(rv1,1,[]);
if isempty(b), b=0; end;
if isempty(c), c=0; end;
if isempty(f), f=0; end;
if isempty(g), g=0; end;
if isempty(p), p=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
ind_shape=size(ind);ind=reshape(ind,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  IMTQLV
ierr = 0;
k = 0;
tag = 0;
%
for i = 1 : n;
w(i) = d(i);
if( i~=1 )
rv1(i-1) = e(i);
end;
end; i = fix(n+1);
%
e2(1) = 0.0e0;
rv1(n) = 0.0e0;
%
gt=0;
for l = 1 : n;
j = 0;
%     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
while( true );
for m = l : n;
if( m==n )
break;
end;
s1 = abs(w(m)) + abs(w(m+1));
s2 = s1 + abs(rv1(m));
if( s2==s1 )
break;
end;
%     .......... GUARD AGAINST UNDERFLOWED ELEMENT OF E2 ..........
if( e2(m+1)==0.0e0 )
k = fix(m);
tag = fix(tag + 1);
gt=1;
end;
end;
%
if( m>k && gt==0)
if( m~=n )
e2(m+1) = 0.0e0;
end;
k = fix(m);
tag = fix(tag + 1);
end;
p = w(l);
if( m==l )
break;
end;
%
if( j==30 )
%     .......... SET ERROR -- NO CONVERGENCE TO AN
%                EIGENVALUE AFTER 30 ITERATIONS ..........
ierr = fix(l);
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
return;
else;
j = fix(j + 1);
%     .......... FORM SHIFT ..........
g =(w(l+1)-p)./(2.0e0.*rv1(l));
[r ,g]=pythag(g,1.0e0);
g = w(m) - p + rv1(l)./(g+(abs(r).*sign(g)));
s = 1.0e0;
c = 1.0e0;
p = 0.0e0;
mml = fix(m - l);
%     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
for ii = 1 : mml;
i = fix(m - ii);
f = s.*rv1(i);
b = c.*rv1(i);
if( abs(f)<abs(g) )
s = f./g;
r = sqrt(s.*s+1.0e0);
rv1(i+1) = g.*r;
c = 1.0e0./r;
s = s.*c;
else;
c = g./f;
r = sqrt(c.*c+1.0e0);
rv1(i+1) = f.*r;
s = 1.0e0./r;
c = c.*s;
end;
g = w(i+1) - p;
r =(w(i)-g).*s + 2.0e0.*c.*b;
p = s.*r;
w(i+1) = g + p;
g = c.*r - b;
end; ii = fix(mml+1);
%
w(l) = w(l) - p;
rv1(l) = g;
rv1(m) = 0.0e0;
end;
end;
%     .......... ORDER EIGENVALUES ..........
gt=0;
if( l~=1 )
%     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
for ii = 2 : l;
i = fix(l + 2 - ii);
if( p>=w(i-1) )
gt=1;
break;
end;
w(i) = w(i-1);
ind(i) = fix(ind(i-1));
end;
end;
%
if(gt==0)
i = 1;
end;
w(i) = p;
ind(i) = fix(tag);
end;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
ind_shape=zeros(ind_shape);ind_shape(:)=ind(1:numel(ind_shape));ind=ind_shape;
end %subroutine imtqlv
%DECK INDXA