[nm,n,a,b,alfr,alfi,beta,matz,z,ierr]=rgg(nm,n,a,b,alfr,alfi,beta,matz,z,ierr);

function [nm,n,a,b,alfr,alfi,beta,matz,z,ierr]=rgg(nm,n,a,b,alfr,alfi,beta,matz,z,ierr);
%***BEGIN PROLOGUE  RGG
%***PURPOSE  Compute the eigenvalues and eigenvectors for a real
%            generalized eigenproblem.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4B2
%***TYPE      SINGLE PRECISION (RGG-S)
%***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR  Smith, B. T., et al.
%***DESCRIPTION
%
%     This subroutine calls the recommended sequence of
%     subroutines from the eigensystem subroutine package (EISPACK)
%     to find the eigenvalues and eigenvectors (if desired)
%     for the REAL GENERAL GENERALIZED eigenproblem  Ax = (LAMBDA)Bx.
%
%     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 is an INTEGER variable set equal to zero if only
%          eigenvalues are desired.  Otherwise, it is set to any
%          non-zero integer for both eigenvalues and eigenvectors.
%
%     On Output
%
%        A and B have been destroyed.
%
%        ALFR and ALFI contain the real and imaginary parts,
%          respectively, of the numerators of the eigenvalues.
%          ALFR and ALFI are one-dimensional REAL arrays,
%          dimensioned ALFR(N) and ALFI(N).
%
%        BETA contains the denominators of the eigenvalues,
%          which are thus given by the ratios  (ALFR+I*ALFI)/BETA.
%          Complex conjugate pairs of eigenvalues appear consecutively
%          with the eigenvalue having the positive imaginary part first.
%          BETA is a one-dimensional REAL array, dimensioned BETA(N).
%
%        Z contains the real and imaginary parts of the eigenvectors
%          if MATZ is not zero.  If the J-th eigenvalue is real, the
%          J-th column of  Z  contains its eigenvector.  If the J-th
%          eigenvalue is complex with positive imaginary part, the
%          J-th and (J+1)-th columns of  Z  contain the real and
%          imaginary parts of its eigenvector.  The conjugate of this
%          vector is the eigenvector for the conjugate eigenvalue.
%          Z is a two-dimensional REAL array, dimensioned Z(NM,N).
%
%        IERR is an INTEGER flag set to
%          Zero       for normal return,
%          10*N       if N is greater than NM,
%          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, but no eigenvectors are
%                     computed.
%
%     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  QZHES, QZIT, QZVAL, QZVEC
%***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  RGG
%
persistent tf ; 

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,[]);
alfr_shape=size(alfr);alfr=reshape(alfr,1,[]);
alfi_shape=size(alfi);alfi=reshape(alfi,1,[]);
beta_shape=size(beta);beta=reshape(beta,1,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(tf), tf=false; end;
%
%***FIRST EXECUTABLE STATEMENT  RGG
if( n>nm )
ierr = fix(10.*n);
%
elseif( matz~=0 ) ;
%     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
tf = true;
[nm,n,a,b,tf,z]=qzhes(nm,n,a,b,tf,z);
[nm,n,a,b,dumvar5,tf,z,ierr]=qzit(nm,n,a,b,0.0e0,tf,z,ierr);
[nm,n,a,b,alfr,alfi,beta,tf,z]=qzval(nm,n,a,b,alfr,alfi,beta,tf,z);
if( ierr==0 )
[nm,n,a,b,alfr,alfi,beta,z]=qzvec(nm,n,a,b,alfr,alfi,beta,z);
end;
else;
%     .......... FIND EIGENVALUES ONLY ..........
tf = false;
[nm,n,a,b,tf,z]=qzhes(nm,n,a,b,tf,z);
[nm,n,a,b,dumvar5,tf,z,ierr]=qzit(nm,n,a,b,0.0e0,tf,z,ierr);
[nm,n,a,b,alfr,alfi,beta,tf,z]=qzval(nm,n,a,b,alfr,alfi,beta,tf,z);
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;
alfr_shape=zeros(alfr_shape);alfr_shape(:)=alfr(1:numel(alfr_shape));alfr=alfr_shape;
alfi_shape=zeros(alfi_shape);alfi_shape(:)=alfi(1:numel(alfi_shape));alfi=alfi_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end
%DECK RJ