[nm,n,low,igh,h,wr,wi,ierr]=hqr(nm,n,low,igh,h,wr,wi,ierr);

function [nm,n,low,igh,h,wr,wi,ierr]=hqr(nm,n,low,igh,h,wr,wi,ierr);
%***BEGIN PROLOGUE  HQR
%***PURPOSE  Compute the eigenvalues of a real upper Hessenberg matrix
%            using the QR method.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C2B
%***TYPE      SINGLE PRECISION (HQR-S, COMQR-C)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure HQR,
%     NUM. MATH. 14, 219-231(1970) by Martin, Peters, and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
%
%     This subroutine finds the eigenvalues of a REAL
%     UPPER Hessenberg matrix by the QR method.
%
%     On INPUT
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameter, H, as declared in the calling program
%          dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrix H.  N is an INTEGER variable.
%          N must be less than or equal to NM.
%
%        LOW and IGH are two INTEGER variables determined by the
%          balancing subroutine  BALANC.  If  BALANC  has not been
%          used, set LOW=1 and IGH equal to the order of the matrix, N.
%
%        H contains the upper Hessenberg matrix.  Information about
%          the transformations used in the reduction to Hessenberg
%          form by  ELMHES  or  ORTHES, if performed, is stored
%          in the remaining triangle under the Hessenberg matrix.
%          H is a two-dimensional REAL array, dimensioned H(NM,N).
%
%     On OUTPUT
%
%        H has been destroyed.  Therefore, it must be saved before
%          calling  HQR  if subsequent calculation and back
%          transformation of eigenvectors is to be performed.
%
%        WR and WI contain the real and imaginary parts, respectively,
%          of the eigenvalues.  The eigenvalues are unordered except
%          that complex conjugate pairs of values appear consecutively
%          with the eigenvalue having the positive imaginary part first.
%          If an error exit is made, the eigenvalues should be correct
%          for indices IERR+1, IERR+2, ..., N.  WR and WI are one-
%          dimensional REAL arrays, dimensioned WR(N) and WI(N).
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          J          if the J-th eigenvalue has not been
%                     determined after a total of 30*N iterations.
%                     The eigenvalues should be correct for indices
%                     IERR+1, IERR+2, ..., N.
%
%     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  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   760101  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   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  HQR
%
persistent en enm2 gt i itn its j k l ll m mm mp2 na norm notlas p q r s s1 s2 t w x y zz ; 

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(en), en=0; end;
if isempty(ll), ll=0; end;
if isempty(mm), mm=0; end;
if isempty(na), na=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(mp2), mp2=0; end;
if isempty(enm2), enm2=0; end;
if isempty(gt), gt=0; end;
h_shape=size(h);h=reshape([h(:).',zeros(1,ceil(numel(h)./prod([nm])).*prod([nm])-numel(h))],nm,[]);
wr_shape=size(wr);wr=reshape(wr,1,[]);
wi_shape=size(wi);wi=reshape(wi,1,[]);
if isempty(p), p=0; end;
if isempty(q), q=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
if isempty(w), w=0; end;
if isempty(x), x=0; end;
if isempty(y), y=0; end;
if isempty(zz), zz=0; end;
if isempty(norm), norm=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(notlas), notlas=false; end;
%
%***FIRST EXECUTABLE STATEMENT  HQR
gt=0;
ierr = 0;
norm = 0.0e0;
k = 1;
%     .......... STORE ROOTS ISOLATED BY BALANC
%                AND COMPUTE MATRIX NORM ..........
for i = 1 : n;
%
for j = k : n;
norm = norm + abs(h(i,j));
end; j = fix(n+1);
%
k = fix(i);
if( i<low || i>igh )
wr(i) = h(i,i);
wi(i) = 0.0e0;
end;
end; i = fix(n+1);
%
en = fix(igh);
t = 0.0e0;
itn = fix(30.*n);
%     .......... SEARCH FOR NEXT EIGENVALUES ..........
while( en>=low );
gt=0;
its = 0;
na = fix(en - 1);
enm2 = fix(na - 1);
%     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
%                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
while( true );
for ll = low : en;
l = fix(en + low - ll);
if( l==low )
break;
end;
s = abs(h(l-1,l-1)) + abs(h(l,l));
if( s==0.0e0 )
s = norm;
end;
s2 = s + abs(h(l,l-1));
if( s2==s )
break;
end;
end;
%     .......... FORM SHIFT ..........
x = h(en,en);
if( l==en )
%     .......... ONE ROOT FOUND ..........
wr(en) = x + t;
wi(en) = 0.0e0;
en = fix(na);
gt=1;
break;
else;
y = h(na,na);
w = h(en,na).*h(na,en);
if( l==na )
break;
end;
if( itn==0 )
%     .......... SET ERROR -- NO CONVERGENCE TO AN
%                EIGENVALUE AFTER 30*N ITERATIONS ..........
ierr = fix(en);
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
wr_shape=zeros(wr_shape);wr_shape(:)=wr(1:numel(wr_shape));wr=wr_shape;
wi_shape=zeros(wi_shape);wi_shape(:)=wi(1:numel(wi_shape));wi=wi_shape;
return;
else;
if( its==10 || its==20 )
%     .......... FORM EXCEPTIONAL SHIFT ..........
t = t + x;
%
for i = low : en;
h(i,i) = h(i,i) - x;
end; i = fix(en+1);
%
s = abs(h(en,na)) + abs(h(na,enm2));
x = 0.75e0.*s;
y = x;
w = -0.4375e0.*s.*s;
end;
its = fix(its + 1);
itn = fix(itn - 1);
%     .......... LOOK FOR TWO CONSECUTIVE SMALL
%                SUB-DIAGONAL ELEMENTS.
%                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
for mm = l : enm2;
m = fix(enm2 + l - mm);
zz = h(m,m);
r = x - zz;
s = y - zz;
p =(r.*s-w)./h(m+1,m) + h(m,m+1);
q = h(m+1,m+1) - zz - r - s;
r = h(m+2,m+1);
s = abs(p) + abs(q) + abs(r);
p = p./s;
q = q./s;
r = r./s;
if( m==l )
break;
end;
s1 = abs(p).*(abs(h(m-1,m-1))+abs(zz)+abs(h(m+1,m+1)));
s2 = s1 + abs(h(m,m-1)).*(abs(q)+abs(r));
if( s2==s1 )
break;
end;
end;
%
mp2 = fix(m + 2);
%
for i = mp2 : en;
h(i,i-2) = 0.0e0;
if( i~=mp2 )
h(i,i-3) = 0.0e0;
end;
end; i = fix(en+1);
%     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
%                COLUMNS M TO EN ..........
for k = m : na;
notlas = k~=na;
if( k~=m )
p = h(k,k-1);
q = h(k+1,k-1);
r = 0.0e0;
if( notlas )
r = h(k+2,k-1);
end;
x = abs(p) + abs(q) + abs(r);
if( x==0.0e0 )
continue;
end;
p = p./x;
q = q./x;
r = r./x;
end;
s = (abs(sqrt(p.*p+q.*q+r.*r)).*sign(p));
if( k==m )
if( l~=m )
h(k,k-1) = -h(k,k-1);
end;
else;
h(k,k-1) = -s.*x;
end;
p = p + s;
x = p./s;
y = q./s;
zz = r./s;
q = q./p;
r = r./p;
%     .......... ROW MODIFICATION ..........
for j = k : en;
p = h(k,j) + q.*h(k+1,j);
if( notlas )
p = p + r.*h(k+2,j);
h(k+2,j) = h(k+2,j) - p.*zz;
end;
h(k+1,j) = h(k+1,j) - p.*y;
h(k,j) = h(k,j) - p.*x;
end; j = fix(en+1);
%
j = fix(min(en,k+3));
%     .......... COLUMN MODIFICATION ..........
for i = l : j;
%
p = x.*h(i,k) + y.*h(i,k+1);
if( notlas )
p = p + zz.*h(i,k+2);
h(i,k+2) = h(i,k+2) - p.*r;
end;
h(i,k+1) = h(i,k+1) - p.*q;
h(i,k) = h(i,k) - p;
%
end; i = fix(j+1);
end; k = fix(na+1);
end;
end;
end;
%     .......... TWO ROOTS FOUND ..........
if(gt==0)
p =(y-x)./2.0e0;
q = p.*p + w;
zz = sqrt(abs(q));
x = x + t;
if( q<0.0e0 )
%     .......... COMPLEX PAIR ..........
wr(na) = x + p;
wr(en) = x + p;
wi(na) = zz;
wi(en) = -zz;
else;
%     .......... REAL PAIR ..........
zz = p + (abs(zz).*sign(p));
wr(na) = x + zz;
wr(en) = wr(na);
if( zz~=0.0e0 )
wr(en) = x - w./zz;
end;
wi(na) = 0.0e0;
wi(en) = 0.0e0;
end;
en = fix(enm2);
end;
end;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
wr_shape=zeros(wr_shape);wr_shape(:)=wr(1:numel(wr_shape));wr=wr_shape;
wi_shape=zeros(wi_shape);wi_shape(:)=wi(1:numel(wi_shape));wi=wi_shape;
return;
end %subroutine hqr
%DECK HSTART