[nm,n,ar,ai,wr,wi,matz,zr,zi,fv1,fv2,fv3,ierr]=cg(nm,n,ar,ai,wr,wi,matz,zr,zi,fv1,fv2,fv3,ierr);

function [nm,n,ar,ai,wr,wi,matz,zr,zi,fv1,fv2,fv3,ierr]=cg(nm,n,ar,ai,wr,wi,matz,zr,zi,fv1,fv2,fv3,ierr);
%***BEGIN PROLOGUE  CG
%***PURPOSE  Compute the eigenvalues and, optionally, the eigenvectors
%            of a complex general matrix.
%***LIBRARY   SLATEC (EISPACK)
%***CATEGORY  D4A4
%***TYPE      COMPLEX (RG-S, CG-C)
%***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)
%     of a COMPLEX GENERAL matrix.
%
%     On INPUT
%
%        NM must be set to the row dimension of the two-dimensional
%          array parameters, AR, AI, ZR and ZI, 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.
%
%        AR and AI contain the real and imaginary parts, respectively,
%          of the complex general matrix.  AR and AI are two-dimensional
%          REAL arrays, dimensioned AR(NM,N) and AI(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
%
%        WR and WI contain the real and imaginary parts, respectively,
%          of the eigenvalues.  WR and WI are one-dimensional REAL
%          arrays, dimensioned WR(N) and WI(N).
%
%        ZR and ZI contain the real and imaginary parts, respectively,
%          of the eigenvectors if MATZ is not zero.  ZR and ZI are
%          two-dimensional REAL arrays, dimensioned ZR(NM,N) and
%          ZI(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 iterations.
%                     The eigenvalues should be correct for indices
%                     IERR+1, IERR+2, ..., N, but no eigenvectors are
%                     computed.
%
%        FV1, FV2, and FV3 are one-dimensional REAL arrays used for
%          temporary storage, dimensioned FV1(N), FV2(N), and FV3(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  CBABK2, CBAL, COMQR, COMQR2, CORTH
%***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  CG
%
persistent is1 is2 ; 

if isempty(is1), is1=0; end;
if isempty(is2), is2=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,[]);
wr_shape=size(wr);wr=reshape(wr,1,[]);
wi_shape=size(wi);wi=reshape(wi,1,[]);
zr_shape=size(zr);zr=reshape([zr(:).',zeros(1,ceil(numel(zr)./prod([nm])).*prod([nm])-numel(zr))],nm,[]);
zi_shape=size(zi);zi=reshape([zi(:).',zeros(1,ceil(numel(zi)./prod([nm])).*prod([nm])-numel(zi))],nm,[]);
fv1_shape=size(fv1);fv1=reshape(fv1,1,[]);
fv2_shape=size(fv2);fv2=reshape(fv2,1,[]);
fv3_shape=size(fv3);fv3=reshape(fv3,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  CG
if( n<=nm )
%
[nm,n,ar,ai,is1,is2,fv1]=cbal(nm,n,ar,ai,is1,is2,fv1);
[nm,n,is1,is2,ar,ai,fv2,fv3]=corth(nm,n,is1,is2,ar,ai,fv2,fv3);
if( matz~=0 )
%     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
[nm,n,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr]=comqr2(nm,n,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr);
if( ierr==0 )
n_orig=n;    [nm,n,is1,is2,fv1,dumvar6,zr,zi]=cbabk2(nm,n,is1,is2,fv1,n,zr,zi);    n(dumvar6~=n_orig)=dumvar6(dumvar6~=n_orig);
end;
else;
%     .......... FIND EIGENVALUES ONLY ..........
[nm,n,is1,is2,ar,ai,wr,wi,ierr]=comqr(nm,n,is1,is2,ar,ai,wr,wi,ierr);
end;
else;
ierr = fix(10.*n);
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;
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;
zr_shape=zeros(zr_shape);zr_shape(:)=zr(1:numel(zr_shape));zr=zr_shape;
zi_shape=zeros(zi_shape);zi_shape(:)=zi(1:numel(zi_shape));zi=zi_shape;
fv1_shape=zeros(fv1_shape);fv1_shape(:)=fv1(1:numel(fv1_shape));fv1=fv1_shape;
fv2_shape=zeros(fv2_shape);fv2_shape(:)=fv2(1:numel(fv2_shape));fv2=fv2_shape;
fv3_shape=zeros(fv3_shape);fv3_shape(:)=fv3(1:numel(fv3_shape));fv3=fv3_shape;
end
%DECK CGTSL