[nm,n,d,e,z,ierr]=tql2(nm,n,d,e,z,ierr);

function [nm,n,d,e,z,ierr]=tql2(nm,n,d,e,z,ierr);
%***BEGIN PROLOGUE  TQL2
%***PURPOSE  Compute the eigenvalues and eigenvectors of symmetric
%            tridiagonal matrix.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4A5, D4C2A
%***TYPE      SINGLE PRECISION (TQL2-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure TQL2,
%     NUM. MATH. 11, 293-306(1968) by Bowdler, Martin, Reinsch, and
%     Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
%
%     This subroutine finds the eigenvalues and eigenvectors
%     of a SYMMETRIC TRIDIAGONAL matrix by the QL method.
%     The eigenvectors of a FULL SYMMETRIC matrix can also
%     be found if  TRED2  has been used to reduce this
%     full matrix to tridiagonal form.
%
%     On Input
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameter, Z, as declared in the calling program
%          dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrix.  N is an INTEGER variable.
%          N must be less than or equal to NM.
%
%        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).
%
%        Z contains the transformation matrix produced in the
%          reduction by  TRED2, if performed.  If the eigenvectors
%          of the tridiagonal matrix are desired, Z must contain
%          the identity matrix.  Z is a two-dimensional REAL array,
%          dimensioned Z(NM,N).
%
%      On Output
%
%        D contains the eigenvalues in ascending order.  If an
%          error exit is made, the eigenvalues are correct but
%          unordered for indices 1, 2, ..., IERR-1.
%
%        E has been destroyed.
%
%        Z contains orthonormal eigenvectors of the symmetric
%          tridiagonal (or full) matrix.  If an error exit is made,
%          Z contains the eigenvectors associated with the stored
%          eigenvalues.
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          J          if the J-th eigenvalue has not been
%                     determined after 30 iterations.
%
%     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  TQL2
%
persistent b c c2 c3 dl1 el1 f g h i ii j k l l1 l2 m mml p r s s2 ; 

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(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(mml), mml=0; end;
d_shape=size(d);d=reshape(d,1,[]);
e_shape=size(e);e=reshape(e,1,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(b), b=0; end;
if isempty(c), c=0; end;
if isempty(c2), c2=0; end;
if isempty(c3), c3=0; end;
if isempty(dl1), dl1=0; end;
if isempty(el1), el1=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;
if isempty(s2), s2=0; end;
%
%***FIRST EXECUTABLE STATEMENT  TQL2
ierr = 0;
if( n~=1 )
%
for i = 2 : n;
e(i-1) = e(i);
end; i = fix(n+1);
%
f = 0.0e0;
b = 0.0e0;
e(n) = 0.0e0;
%
for l = 1 : n;
j = 0;
h = abs(d(l)) + abs(e(l));
if( b<h )
b = h;
end;
%     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
for m = l : n;
if( b+abs(e(m))==b )
break;
end;
%     .......... E(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;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
j = fix(j + 1);
%     .......... FORM SHIFT ..........
l1 = fix(l + 1);
l2 = fix(l1 + 1);
g = d(l);
p =(d(l1)-g)./(2.0e0.*e(l));
[r ,p]=pythag(p,1.0e0);
d(l) = e(l)./(p+(abs(r).*sign(p)));
d(l1) = e(l).*(p+(abs(r).*sign(p)));
dl1 = d(l1);
h = g - d(l);
if( l2<=n )
%
for i = l2 : n;
d(i) = d(i) - h;
end; i = fix(n+1);
end;
%
f = f + h;
%     .......... QL TRANSFORMATION ..........
p = d(m);
c = 1.0e0;
c2 = c;
el1 = e(l1);
s = 0.0e0;
mml = fix(m - l);
%     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
for ii = 1 : mml;
c3 = c2;
c2 = c;
s2 = s;
i = fix(m - ii);
g = c.*e(i);
h = c.*p;
if( abs(p)<abs(e(i)) )
c = p./e(i);
r = sqrt(c.*c+1.0e0);
e(i+1) = s.*e(i).*r;
s = 1.0e0./r;
c = c.*s;
else;
c = e(i)./p;
r = sqrt(c.*c+1.0e0);
e(i+1) = s.*p.*r;
s = c./r;
c = 1.0e0./r;
end;
p = c.*d(i) - s.*g;
d(i+1) = h + s.*(c.*g+s.*d(i));
%     .......... FORM VECTOR ..........
for k = 1 : n;
h = z(k,i+1);
z(k,i+1) = s.*z(k,i) + c.*h;
z(k,i) = c.*z(k,i) - s.*h;
end; k = fix(n+1);
%
end; ii = fix(mml+1);
%
p = -s.*s2.*c3.*el1.*e(l)./dl1;
e(l) = s.*p;
d(l) = c.*p;
if( b+abs(e(l))>b )
continue;
end;
end;
break;
end;
end;
d(l) = d(l) + f;
end;
%     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
for ii = 2 : n;
i = fix(ii - 1);
k = fix(i);
p = d(i);
%
for j = ii : n;
if( d(j)<p )
k = fix(j);
p = d(j);
end;
end; j = fix(n+1);
%
if( k~=i )
d(k) = d(i);
d(i) = p;
%
for j = 1 : n;
p = z(j,i);
z(j,i) = z(j,k);
z(j,k) = p;
end; j = fix(n+1);
end;
%
%
end; ii = fix(n+1);
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine tql2
%DECK TQLRAT