[nm,n,a,b,matz,z]=qzhes(nm,n,a,b,matz,z);

function [nm,n,a,b,matz,z]=qzhes(nm,n,a,b,matz,z);
%***BEGIN PROLOGUE  QZHES
%***PURPOSE  The first step of the QZ algorithm for solving generalized
%            matrix eigenproblems.  Accepts a pair of real general
%            matrices and reduces one of them to upper Hessenberg
%            and the other to upper triangular form using orthogonal
%            transformations. Usually followed by QZIT, QZVAL, QZVEC.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4C1B3
%***TYPE      SINGLE PRECISION (QZHES-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine is the first step of the QZ algorithm
%     for solving generalized matrix eigenvalue problems,
%     SIAM J. NUMER. ANAL. 10, 241-256(1973) by MOLER and STEWART.
%
%     This subroutine accepts a pair of REAL GENERAL matrices and
%     reduces one of them to upper Hessenberg form and the other
%     to upper triangular form using orthogonal transformations.
%     It is usually followed by  QZIT,  QZVAL  and, possibly,  QZVEC.
%
%     On Input
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameters, A, B, and Z, as declared in the calling
%          program dimension statement.  NM is an INTEGER variable.
%
%        N is the order of the matrices A and B.  N is an INTEGER
%          variable.  N must be less than or equal to NM.
%
%        A contains a real general matrix.  A is a two-dimensional
%          REAL array, dimensioned A(NM,N).
%
%        B contains a real general matrix.  B is a two-dimensional
%          REAL array, dimensioned B(NM,N).
%
%        MATZ should be set to true if the right hand transformations
%          are to be accumulated for later use in computing
%          eigenvectors, and to false otherwise.  MATZ is a LOGICAL
%          variable.
%
%     On Output
%
%        A has been reduced to upper Hessenberg form.  The elements
%          below the first subdiagonal have been set to zero.
%
%        B has been reduced to upper triangular form.  The elements
%          below the main diagonal have been set to zero.
%
%        Z contains the product of the right hand transformations if
%          MATZ has been set to true  Otherwise, Z is not referenced.
%          Z is a two-dimensional REAL array, dimensioned Z(NM,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
%   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  QZHES
%
persistent i j k l l1 lb nk1 nm1 nm2 r rho s t u1 u2 v1 v2 ; 

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(lb), lb=0; end;
if isempty(l1), l1=0; end;
if isempty(nk1), nk1=0; end;
if isempty(nm1), nm1=0; end;
if isempty(nm2), nm2=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([nm])).*prod([nm])-numel(b))],nm,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
if isempty(u1), u1=0; end;
if isempty(u2), u2=0; end;
if isempty(v1), v1=0; end;
if isempty(v2), v2=0; end;
if isempty(rho), rho=0; end;
%
%     .......... INITIALIZE Z ..........
%***FIRST EXECUTABLE STATEMENT  QZHES
if( matz )
%
for i = 1 : n;
%
for j = 1 : n;
z(i,j) = 0.0e0;
end; j = fix(n+1);
%
z(i,i) = 1.0e0;
end; i = fix(n+1);
end;
%     .......... REDUCE B TO UPPER TRIANGULAR FORM ..........
if( n>1 )
nm1 = fix(n - 1);
%
for l = 1 : nm1;
l1 = fix(l + 1);
s = 0.0e0;
%
for i = l1 : n;
s = s + abs(b(i,l));
end; i = fix(n+1);
%
if( s~=0.0e0 )
s = s + abs(b(l,l));
r = 0.0e0;
%
for i = l : n;
b(i,l) = b(i,l)./s;
r = r + b(i,l).^2;
end; i = fix(n+1);
%
r = (abs(sqrt(r)).*sign(b(l,l)));
b(l,l) = b(l,l) + r;
rho = r.*b(l,l);
%
for j = l1 : n;
t = 0.0e0;
%
for i = l : n;
t = t + b(i,l).*b(i,j);
end; i = fix(n+1);
%
t = -t./rho;
%
for i = l : n;
b(i,j) = b(i,j) + t.*b(i,l);
end; i = fix(n+1);
%
end; j = fix(n+1);
%
for j = 1 : n;
t = 0.0e0;
%
for i = l : n;
t = t + b(i,l).*a(i,j);
end; i = fix(n+1);
%
t = -t./rho;
%
for i = l : n;
a(i,j) = a(i,j) + t.*b(i,l);
end; i = fix(n+1);
%
end; j = fix(n+1);
%
b(l,l) = -s.*r;
%
for i = l1 : n;
b(i,l) = 0.0e0;
end; i = fix(n+1);
end;
%
end; l = fix(nm1+1);
%     .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE
%                KEEPING B TRIANGULAR ..........
if( n~=2 )
nm2 = fix(n - 2);
%
for k = 1 : nm2;
nk1 = fix(nm1 - k);
%     .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- ..........
for lb = 1 : nk1;
l = fix(n - lb);
l1 = fix(l + 1);
%     .......... ZERO A(L+1,K) ..........
s = abs(a(l,k)) + abs(a(l1,k));
if( s~=0.0e0 )
u1 = a(l,k)./s;
u2 = a(l1,k)./s;
r = (abs(sqrt(u1.*u1+u2.*u2)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
u2 = v2./v1;
%
for j = k : n;
t = a(l,j) + u2.*a(l1,j);
a(l,j) = a(l,j) + t.*v1;
a(l1,j) = a(l1,j) + t.*v2;
end; j = fix(n+1);
%
a(l1,k) = 0.0e0;
%
for j = l : n;
t = b(l,j) + u2.*b(l1,j);
b(l,j) = b(l,j) + t.*v1;
b(l1,j) = b(l1,j) + t.*v2;
end; j = fix(n+1);
%     .......... ZERO B(L+1,L) ..........
s = abs(b(l1,l1)) + abs(b(l1,l));
if( s~=0.0e0 )
u1 = b(l1,l1)./s;
u2 = b(l1,l)./s;
r = (abs(sqrt(u1.*u1+u2.*u2)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
u2 = v2./v1;
%
for i = 1 : l1;
t = b(i,l1) + u2.*b(i,l);
b(i,l1) = b(i,l1) + t.*v1;
b(i,l) = b(i,l) + t.*v2;
end; i = fix(l1+1);
%
b(l1,l) = 0.0e0;
%
for i = 1 : n;
t = a(i,l1) + u2.*a(i,l);
a(i,l1) = a(i,l1) + t.*v1;
a(i,l) = a(i,l) + t.*v2;
end; i = fix(n+1);
%
if( matz )
%
for i = 1 : n;
t = z(i,l1) + u2.*z(i,l);
z(i,l1) = z(i,l1) + t.*v1;
z(i,l) = z(i,l) + t.*v2;
end; i = fix(n+1);
end;
end;
end;
%
end; lb = fix(nk1+1);
%
end; k = fix(nm2+1);
end;
end;
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end
%DECK QZIT