[nm,n,low,igh,ar,ai,intv]=comhes(nm,n,low,igh,ar,ai,intv);

function [nm,n,low,igh,ar,ai,intv]=comhes(nm,n,low,igh,ar,ai,intv);
%***BEGIN PROLOGUE  COMHES
%***PURPOSE  Reduce a complex general matrix to complex upper Hessenberg
%            form using stabilized elementary similarity
%            transformations.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C1B2
%***TYPE      COMPLEX (ELMHES-S, COMHES-C)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is a translation of the ALGOL procedure COMHES,
%     NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson.
%     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
%
%     Given a COMPLEX GENERAL matrix, this subroutine
%     reduces a submatrix situated in rows and columns
%     LOW through IGH to upper Hessenberg form by
%     stabilized elementary similarity transformations.
%
%     On INPUT
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameters, AR and AI, as declared in the calling
%          program dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrix A=(AR,AI).  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.
%
%        AR and AI contain the real and imaginary parts, respectively,
%          of the complex input matrix.  AR and AI are two-dimensional
%          REAL arrays, dimensioned AR(NM,N) and AI(NM,N).
%
%     On OUTPUT
%
%        AR and AI contain the real and imaginary parts, respectively,
%          of the upper Hessenberg matrix.  The multipliers which
%          were used in the reduction are stored in the remaining
%          triangles under the Hessenberg matrix.
%
%        INT contains information on the rows and columns
%          interchanged in the reduction.  Only elements LOW through
%          IGH are used.  INT is a one-dimensional INTEGER array,
%          dimensioned INT(IGH).
%
%     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
%***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  COMHES
%
persistent i j kp1 la m mm1 mp1 xi xr yi yr ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(m), m=0; end;
if isempty(la), la=0; end;
if isempty(kp1), kp1=0; end;
if isempty(mm1), mm1=0; end;
if isempty(mp1), mp1=0; end;
ar_shape=size(ar);ar=reshape([ar(:).',zeros(1,ceil(numel(ar)./prod([nm])).*prod([nm])-numel(ar))],nm,[]);
ai_shape=size(ai);ai=reshape([ai(:).',zeros(1,ceil(numel(ai)./prod([nm])).*prod([nm])-numel(ai))],nm,[]);
if isempty(xr), xr=0; end;
if isempty(xi), xi=0; end;
if isempty(yr), yr=0; end;
if isempty(yi), yi=0; end;
intv_shape=size(intv);intv=reshape(intv,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  COMHES
la = fix(igh - 1);
kp1 = fix(low + 1);
if( la>=kp1 )
%
for m = kp1 : la;
mm1 = fix(m - 1);
xr = 0.0e0;
xi = 0.0e0;
i = fix(m);
%
for j = m : igh;
if( abs(ar(j,mm1))+abs(ai(j,mm1))>abs(xr)+abs(xi) )
xr = ar(j,mm1);
xi = ai(j,mm1);
i = fix(j);
end;
end; j = fix(igh+1);
%
intv(m) = fix(i);
if( i~=m )
%     .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
for j = mm1 : n;
yr = ar(i,j);
ar(i,j) = ar(m,j);
ar(m,j) = yr;
yi = ai(i,j);
ai(i,j) = ai(m,j);
ai(m,j) = yi;
end; j = fix(n+1);
%
for j = 1 : igh;
yr = ar(j,i);
ar(j,i) = ar(j,m);
ar(j,m) = yr;
yi = ai(j,i);
ai(j,i) = ai(j,m);
ai(j,m) = yi;
end; j = fix(igh+1);
end;
%     .......... end INTERCHANGE ..........
if( xr~=0.0e0 || xi~=0.0e0 )
mp1 = fix(m + 1);
%
for i = mp1 : igh;
yr = ar(i,mm1);
yi = ai(i,mm1);
if( yr~=0.0e0 || yi~=0.0e0 )
yr_orig=yr; yi_orig=yi;    [yr,yi,xr,xi,dumvar5,dumvar6]=cdiv(yr,yi,xr,xi,yr,yi);    yi(dumvar6~=yi_orig)=dumvar6(dumvar6~=yi_orig); yr(dumvar5~=yr_orig)=dumvar5(dumvar5~=yr_orig);
ar(i,mm1) = yr;
ai(i,mm1) = yi;
%
for j = m : n;
ar(i,j) = ar(i,j) - yr.*ar(m,j) + yi.*ai(m,j);
ai(i,j) = ai(i,j) - yr.*ai(m,j) - yi.*ar(m,j);
end; j = fix(n+1);
%
for j = 1 : igh;
ar(j,m) = ar(j,m) + yr.*ar(j,i) - yi.*ai(j,i);
ai(j,m) = ai(j,m) + yr.*ai(j,i) + yi.*ar(j,i);
end; j = fix(igh+1);
end;
%
end; i = fix(igh+1);
end;
%
end; m = fix(la+1);
end;
%
ar_shape=zeros(ar_shape);ar_shape(:)=ar(1:numel(ar_shape));ar=ar_shape;
ai_shape=zeros(ai_shape);ai_shape(:)=ai(1:numel(ai_shape));ai=ai_shape;
intv_shape=zeros(intv_shape);intv_shape(:)=intv(1:numel(intv_shape));intv=intv_shape;
end
%DECK COMLR2