[nm,n,low,igh,hr,hi,wr,wi,ierr]=comlr(nm,n,low,igh,hr,hi,wr,wi,ierr);

function [nm,n,low,igh,hr,hi,wr,wi,ierr]=comlr(nm,n,low,igh,hr,hi,wr,wi,ierr);
%***BEGIN PROLOGUE  COMLR
%***PURPOSE  Compute the eigenvalues of a complex upper Hessenberg
%            matrix using the modified LR method.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C2B
%***TYPE      COMPLEX (COMLR-C)
%***KEYWORDS  EIGENVALUES, EISPACK, LR METHOD
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure COMLR,
%     NUM. MATH. 12, 369-376(1968) by Martin and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
%
%     This subroutine finds the eigenvalues of a COMPLEX
%     UPPER Hessenberg matrix by the modified LR method.
%
%     On INPUT
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameters, HR and HI, as declared in the calling
%          program dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrix H=(HR,HI).  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  CBAL.  If  CBAL  has not been used,
%          set LOW=1 and IGH equal to the order of the matrix, N.
%
%        HR and HI contain the real and imaginary parts, respectively,
%          of the complex upper Hessenberg matrix.  Their lower
%          triangles below the subdiagonal contain the multipliers
%          which were used in the reduction by  COMHES, if performed.
%          HR and HI are two-dimensional REAL arrays, dimensioned
%          HR(NM,N) and HI(NM,N).
%
%     On OUTPUT
%
%        The upper Hessenberg portions of HR and HI have been
%          destroyed.  Therefore, they must be saved before calling
%          COMLR  if subsequent calculation of eigenvectors is to
%          be performed.
%
%        WR and WI contain the real and imaginary parts, respectively,
%          of the eigenvalues of the upper Hessenberg matrix.  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.
%
%     Calls CSROOT for complex square root.
%     Calls CDIV for complex division.
%
%     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  CDIV, CSROOT
%***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  COMLR
%
persistent en enm1 i im1 itn its j l ll m mm mp1 s1 s2 si sr ti tr xi xr yi yr zzi zzr ; 

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(en), en=0; end;
if isempty(ll), ll=0; end;
if isempty(mm), mm=0; end;
if isempty(im1), im1=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(mp1), mp1=0; end;
if isempty(enm1), enm1=0; end;
hr_shape=size(hr);hr=reshape([hr(:).',zeros(1,ceil(numel(hr)./prod([nm])).*prod([nm])-numel(hr))],nm,[]);
hi_shape=size(hi);hi=reshape([hi(:).',zeros(1,ceil(numel(hi)./prod([nm])).*prod([nm])-numel(hi))],nm,[]);
wr_shape=size(wr);wr=reshape(wr,1,[]);
wi_shape=size(wi);wi=reshape(wi,1,[]);
if isempty(si), si=0; end;
if isempty(sr), sr=0; end;
if isempty(ti), ti=0; end;
if isempty(tr), tr=0; end;
if isempty(xi), xi=0; end;
if isempty(xr), xr=0; end;
if isempty(yi), yi=0; end;
if isempty(yr), yr=0; end;
if isempty(zzi), zzi=0; end;
if isempty(zzr), zzr=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
%
%***FIRST EXECUTABLE STATEMENT  COMLR
ierr = 0;
%     .......... STORE ROOTS ISOLATED BY CBAL ..........
for i = 1 : n;
if( i<low || i>igh )
wr(i) = hr(i,i);
wi(i) = hi(i,i);
end;
end; i = fix(n+1);
%
en = fix(igh);
tr = 0.0e0;
ti = 0.0e0;
itn = fix(30.*n);
%     .......... SEARCH FOR NEXT EIGENVALUE ..........
while( en>=low );
its = 0;
enm1 = fix(en - 1);
%     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
%                FOR L=EN STEP -1 UNTIL LOW E0 -- ..........
while( true );
for ll = low : en;
l = fix(en + low - ll);
if( l==low )
break;
end;
s1 = abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) + abs(hr(l,l))+ abs(hi(l,l));
s2 = s1 + abs(hr(l,l-1)) + abs(hi(l,l-1));
if( s2==s1 )
break;
end;
end;
%     .......... FORM SHIFT ..........
if( l==en )
break;
end;
if( itn==0 )
%     .......... SET ERROR -- NO CONVERGENCE TO AN
%                EIGENVALUE AFTER 30*N ITERATIONS ..........
ierr = fix(en);
hr_shape=zeros(hr_shape);hr_shape(:)=hr(1:numel(hr_shape));hr=hr_shape;
hi_shape=zeros(hi_shape);hi_shape(:)=hi(1:numel(hi_shape));hi=hi_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 ..........
sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2));
si = abs(hi(en,enm1)) + abs(hi(enm1,en-2));
else;
sr = hr(en,en);
si = hi(en,en);
xr = hr(enm1,en).*hr(en,enm1) - hi(enm1,en).*hi(en,enm1);
xi = hr(enm1,en).*hi(en,enm1) + hi(enm1,en).*hr(en,enm1);
if( xr~=0.0e0 || xi~=0.0e0 )
yr =(hr(enm1,enm1)-sr)./2.0e0;
yi =(hi(enm1,enm1)-si)./2.0e0;
[dumvar1,dumvar2,zzr,zzi]=csroot(yr.^2-yi.^2+xr,2.0e0.*yr.*yi+xi,zzr,zzi);
if( yr.*zzr+yi.*zzi<0.0e0 )
zzr = -zzr;
zzi = -zzi;
end;
xr_orig=xr; xi_orig=xi;    [xr,xi,dumvar3,dumvar4,dumvar5,dumvar6]=cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi);    xi(dumvar6~=xi_orig)=dumvar6(dumvar6~=xi_orig); xr(dumvar5~=xr_orig)=dumvar5(dumvar5~=xr_orig);
sr = sr - xr;
si = si - xi;
end;
end;
%
for i = low : en;
hr(i,i) = hr(i,i) - sr;
hi(i,i) = hi(i,i) - si;
end; i = fix(en+1);
%
tr = tr + sr;
ti = ti + si;
its = fix(its + 1);
itn = fix(itn - 1);
%     .......... LOOK FOR TWO CONSECUTIVE SMALL
%                SUB-DIAGONAL ELEMENTS ..........
xr = abs(hr(enm1,enm1)) + abs(hi(enm1,enm1));
yr = abs(hr(en,enm1)) + abs(hi(en,enm1));
zzr = abs(hr(en,en)) + abs(hi(en,en));
%     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- ..........
for mm = l : enm1;
m = fix(enm1 + l - mm);
if( m==l )
break;
end;
yi = yr;
yr = abs(hr(m,m-1)) + abs(hi(m,m-1));
xi = zzr;
zzr = xr;
xr = abs(hr(m-1,m-1)) + abs(hi(m-1,m-1));
s1 = zzr./yi.*(zzr+xr+xi);
s2 = s1 + yr;
if( s2==s1 )
break;
end;
end;
%     .......... TRIANGULAR DECOMPOSITION H=L*R ..........
mp1 = fix(m + 1);
%
for i = mp1 : en;
im1 = fix(i - 1);
xr = hr(im1,im1);
xi = hi(im1,im1);
yr = hr(i,im1);
yi = hi(i,im1);
if( abs(xr)+abs(xi)>=abs(yr)+abs(yi) )
[yr,yi,xr,xi,zzr,zzi]=cdiv(yr,yi,xr,xi,zzr,zzi);
wr(i) = -1.0e0;
else;
%     .......... INTERCHANGE ROWS OF HR AND HI ..........
for j = im1 : en;
zzr = hr(im1,j);
hr(im1,j) = hr(i,j);
hr(i,j) = zzr;
zzi = hi(im1,j);
hi(im1,j) = hi(i,j);
hi(i,j) = zzi;
end; j = fix(en+1);
%
[xr,xi,yr,yi,zzr,zzi]=cdiv(xr,xi,yr,yi,zzr,zzi);
wr(i) = 1.0e0;
end;
hr(i,im1) = zzr;
hi(i,im1) = zzi;
%
for j = i : en;
%
hr(i,j) = hr(i,j) - zzr.*hr(im1,j) + zzi.*hi(im1,j);
hi(i,j) = hi(i,j) - zzr.*hi(im1,j) - zzi.*hr(im1,j);
end; j = fix(en+1);
end; i = fix(en+1);
%     .......... COMPOSITION R*L=H ..........
for j = mp1 : en;
xr = hr(j,j-1);
xi = hi(j,j-1);
hr(j,j-1) = 0.0e0;
hi(j,j-1) = 0.0e0;
%     .......... INTERCHANGE COLUMNS OF HR AND HI,
%                IF NECESSARY ..........
if( wr(j)>0.0e0 )
%
for i = l : j;
zzr = hr(i,j-1);
hr(i,j-1) = hr(i,j);
hr(i,j) = zzr;
zzi = hi(i,j-1);
hi(i,j-1) = hi(i,j);
hi(i,j) = zzi;
end; i = fix(j+1);
end;
%
for i = l : j;
%
hr(i,j-1) = hr(i,j-1) + xr.*hr(i,j) - xi.*hi(i,j);
hi(i,j-1) = hi(i,j-1) + xr.*hi(i,j) + xi.*hr(i,j);
%
end; i = fix(j+1);
end; j = fix(en+1);
end;
end;
%     .......... A ROOT FOUND ..........
wr(en) = hr(en,en) + tr;
wi(en) = hi(en,en) + ti;
en = fix(enm1);
end;
hr_shape=zeros(hr_shape);hr_shape(:)=hr(1:numel(hr_shape));hr=hr_shape;
hi_shape=zeros(hi_shape);hi_shape(:)=hi(1:numel(hi_shape));hi=hi_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 comlr
%DECK COMPB