/*--------------------------------------------------------------------------
**  ZtgSen
**
**  Purpose:   Calculate a sorted generalized complex Schur decomposition
**             of A*x=lambda*Bx, where A and B are complex matrices.
**             It generates Q,Z,S, and P in
**                                         A = Q*S*ZH
**                                         B = Q*P*ZH
**             where ZH is the transpose of the complex conjugate of Z.
**             S and P are upper triangular complex matrices,
**             and Q and Z are unitary (i.e. QH = inv(Q) and ZH=inv(Z)).
**             The decomposition is sorted so that abs(eigenvalues) less than
**             a given number comes first.
**
**  Usage:     dlibrary -a PsZtgSen;     (only once)
**             { S,P,Q,Z } = ZtgSen( A,B,cutoff,pii_sii,GtIt )
**
**  Input:     A            NxN complex matrix
**             B              ""
**             Cutoff       scalar, cutoff point. See below
**             pii_sii      scalar, governs the logical condition. See below
**             GtIt             ""
**
**                      pii_sii = 1                  pii_sii < 0 
**            -----|---------------------------------------------------------------- 
**        GtIt = 1 |  abs(pii/sii)>Cutoff first      abs(sii/pii)>Cutoff first
**                 |
**        GtIt < 0 |  abs(pii/sii)<Cutoff first      abs(sii/pii)<Cutoff first
**                 |  
**
**        
**
**  Output:    S            NxN complex matrix (see above)
**             P             ""
**             Q             ""
**             Z             ""
**
**
**  Note:  This procedure calls on a dynamic library (dll file) created
**         by compiling and linking my FORTRAN code PsZtgSen.f with the
**         FORTRAN library Lapack.lib created from LAPACK 2.1 
**         (a pre-release of Lapack 3.0). The calculations are made in
**         complex*16 (complex(8) in FORTRAN 90 languange). In the FORTRAN 
**         code, very small denominators are replaced by 1E-100 to avoid 
**         overflow.
**         
**
**
**
**  Paul Soderlind (nepso@hhs.se), 5 August 1997
----------------------------------------------------------------------------*/

proc (4) = ZtgSen( A,B,Cutoff,pii_sii,GtIt );

	local Ar,Br,N,Ai,Bi,Qr,Qi,Zr,Zi,S,P,Q,Z;


  Ar = real(A');   /*real part of A: easier to pass doubles than complex*16*/
  Br = real(B');   /*to FORTRAN*/
  Ai = imag(A');   /*imaginary part of A*/
  Bi = imag(B');

  N = rows(ar);

  Qr = zeros(N,N);   /*empty matrices to be filled by FORTRAN*/
  Qi = zeros(N,N);
  Zr = zeros(N,N);
  Zi = zeros(N,N);

  dllcall PSZTGSEN(Ar,Ai,Br,Bi,Qr,Qi,Zr,Zi,N,Cutoff,pii_sii,GtIt);

  S = complex( Ar',Ai' ); /*transpose since FORTRAN has column-major order*/
  P = complex( Br',Bi' ); /*converting to complex matrices again*/
  Q = complex( Qr',Qi' );
  Z = complex( Zr',Zi' );

	retp( S,P,Q,Z );

endp;
/*------------------------------------------------------------------------*/

