[n,d,e2,ierr]=tevls(n,d,e2,ierr);

function [n,d,e2,ierr]=tevls(n,d,e2,ierr);
persistent b c dhold f g h i igo ii j l l1 m mml nhalf ntop p r s ; 

global cblkt_4; if isempty(cblkt_4), cblkt_4=0; end;
if isempty(dhold), dhold=0; end;
global cblkt_7; if isempty(cblkt_7), cblkt_7=0; end;
global cblkt_2; if isempty(cblkt_2), cblkt_2=0; end;
global cblkt_6; if isempty(cblkt_6), cblkt_6=0; end;
if isempty(nhalf), nhalf=0; end;
global cblkt_5; if isempty(cblkt_5), cblkt_5=0; end;
global cblkt_1; if isempty(cblkt_1), cblkt_1=0; end;
if isempty(ntop), ntop=0; end;
if isempty(igo), igo=0; end;
%***BEGIN PROLOGUE  TEVLS
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BLKTRI
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (TEVLS-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     This subroutine finds the eigenvalues of a symmetric tridiagonal
%     matrix by the rational QL method.
%
%     On Input-
%
%        N is the order of the matrix,
%
%        D contains the diagonal elements of the input matrix,
%
%        E2 contains the subdiagonal elements of the input matrix
%           in its last N-1 positions.  E2(1) is arbitrary.
%
%      On Output-
%
%        D 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,
%
%        E2 has been destroyed,
%
%        IERR is set to
%          ZERO       for normal return,
%          J          if the J-th eigenvalue has not been
%                     determined after 30 iterations.
%
%***SEE ALSO  BLKTRI
%***REFERENCES  C. H. Reinsch, Eigenvalues of a real, symmetric, tri-
%                 diagonal matrix, Algorithm 464, Communications of the
%                 ACM 16, 11 (November 1973), pp. 689.
%***ROUTINES CALLED  (NONE)
%***COMMON BLOCKS    CBLKT
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%   920528  DESCRIPTION revised and REFERENCES section added.  (WRB)
%***end PROLOGUE  TEVLS
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
if isempty(ii), ii=0; end;
if isempty(l1), l1=0; end;
if isempty(mml), mml=0; end;
d_shape=size(d);d=reshape(d,1,[]);
e2_shape=size(e2);e2=reshape(e2,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(h), h=0; end;
if isempty(p), p=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
global cblkt_3; if isempty(cblkt_3), cblkt_3=0; end;
%
% common :: ;
%% common /cblkt / npp , k , machep , cnv , nm , ncmplx , ik;
%% common /cblkt / cblkt_1 , cblkt_2 , cblkt_3 , cblkt_4 , cblkt_5 , cblkt_6 , cblkt_7;
%***FIRST EXECUTABLE STATEMENT  TEVLS
ierr = 0;
if( n~=1 )
%
for i = 2 : n;
e2(i-1) = e2(i).*e2(i);
end; i = fix(n+1);
%
f = 0.0;
b = 0.0;
e2(n) = 0.0;
%
for l = 1 : n;
j = 0;
h = cblkt_3.*(abs(d(l))+sqrt(e2(l)));
if( b<=h )
b = h;
c = b.*b;
end;
%
%     ********** LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT **********
%
for m = l : n;
if( e2(m)<=c )
break;
end;
%
%     ********** E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
%                THROUGH THE BOTTOM OF THE LOOP **********
%
end;
%
if( m~=l )
while (1);
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;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
return;
else;
j = fix(j + 1);
%
%     ********** FORM SHIFT **********
%
l1 = fix(l + 1);
s = sqrt(e2(l));
g = d(l);
p =(d(l1)-g)./(2.0.*s);
r = sqrt(p.*p+1.0);
d(l) = s./(p+(abs(r).*sign(p)));
h = g - d(l);
%
for i = l1 : n;
d(i) = d(i) - h;
end; i = fix(n+1);
%
f = f + h;
%
%     ********** RATIONAL QL TRANSFORMATION **********
%
g = d(m);
if( g==0.0 )
g = b;
end;
h = g;
s = 0.0;
mml = fix(m - l);
%
%     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
%
for ii = 1 : mml;
i = fix(m - ii);
p = g.*h;
r = p + e2(i);
e2(i+1) = s.*r;
s = e2(i)./r;
d(i+1) = h + s.*(h+d(i));
g = d(i) - e2(i)./g;
if( g==0.0 )
g = b;
end;
h = g.*p./r;
end; ii = fix(mml+1);
%
e2(l) = s.*g;
d(l) = h;
%
%     ********** GUARD AGAINST UNDERFLOWED H **********
%
if( h~=0.0 )
if( abs(e2(l))>abs(c./h) )
e2(l) = h.*e2(l);
if( e2(l)~=0.0 )
continue;
end;
end;
end;
end;
break;
end;
end;
p = d(l) + f;
%
%     ********** ORDER EIGENVALUES **********
%
if( l~=1 )
%
%     ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
%
igo=0;
for ii = 2 : l;
i = fix(l + 2 - ii);
if( p>=d(i-1) )
igo=1;
break;
end;
d(i) = d(i-1);
end;
end;
%
if(igo==0)
i = 1;
end;
d(i) = p;
end;
%
if( abs(d(n))<abs(d(1)) )
nhalf = fix(fix(n./2));
for i = 1 : nhalf;
ntop = fix(n - i);
dhold = d(i);
d(i) = d(ntop+1);
d(ntop+1) = dhold;
end; i = fix(nhalf+1);
end;
end;
%
%     ********** LAST CARD OF TQLRAT **********
%
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e2_shape=zeros(e2_shape);e2_shape(:)=e2(1:numel(e2_shape));e2=e2_shape;
end %subroutine tevls