ZHGEQZ(l)		LAPACK routine (version	1.1)		    ZHGEQZ(l)

NAME
  ZHGEQZ - implement a single-shift version of the QZ method for finding the
  generalized eigenvalues w(i)=ALPHA(i)/BETA(i)	of the equation	  det( A -
  w(i) B ) = 0	If JOB='S', then the pair (A,B)	is simultaneously reduced to
  Schur	form (i.e., A and B are	both upper triangular) using one unitary
  tranformation	(usually called	Q) on the left and another (usually called Z)
  on the right

SYNOPSIS

  SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ,	N, ILO,	IHI, A,	LDA, B,	LDB, ALPHA,
		     BETA, Q, LDQ, Z, LDZ, WORK, LWORK,	RWORK, INFO )

      CHARACTER	     COMPQ, COMPZ, JOB

      INTEGER	     IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK,	N

      DOUBLE	     PRECISION RWORK( *	)

      COMPLEX*16     A(	LDA, * ), ALPHA( * ), B( LDB, *	), BETA( * ), Q( LDQ,
		     * ), WORK(	* ), Z(	LDZ, * )

PURPOSE
  ZHGEQZ implements a single-shift version of the QZ method for	finding	the
  generalized eigenvalues w(i)=ALPHA(i)/BETA(i)	of the equation	A are then
  ALPHA(1),...,ALPHA(N), and of	B are BETA(1),...,BETA(N).

  If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the unitary transforma-
  tions	used to	reduce (A,B) are accumulated into the arrays Q and Z s.t.:

       Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
       Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*

  Ref: C.B. Moler & G.W. Stewart, "An Algorithm	for Generalized	Matrix
       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
       pp. 241--256.

ARGUMENTS

  JOB	  (input) CHARACTER*1
	  = 'E': compute only ALPHA and	BETA.  A and B will not	necessarily
	  be put into generalized Schur	form.  = 'S': put A and	B into gen-
	  eralized Schur form, as well as computing ALPHA and BETA.

  COMPQ	  (input) CHARACTER*1
	  = 'N': do not	modify Q.
	  = 'V': multiply the array Q on the right by the conjugate transpose
	  of the unitary tranformation that is applied to the left side	of A
	  and B	to reduce them to Schur	form.  = 'I': like COMPQ='V', except
	  that Q will be initialized to	the identity first.

  COMPZ	  (input) CHARACTER*1
	  = 'N': do not	modify Z.
	  = 'V': multiply the array Z on the right by the unitary tranforma-
	  tion that is applied to the right side of A and B to reduce them to
	  Schur	form.  = 'I': like COMPZ='V', except that Z will be initial-
	  ized to the identity first.

  N	  (input) INTEGER
	  The number of	rows and columns in the	matrices A, B, Q, and Z.  N
	  must be at least 0.

  ILO	  (input) INTEGER
	  Columns 1 through ILO-1 are assumed to be in triangular form
	  already, and will not	be modified.  ILO must be at least 1.

  IHI	  (input) INTEGER
	  Rows IHI+1 through N are assumed to be in triangular form already,
	  and will not be touched.  IHI	may not	be greater than	N.

  A	  (input/output) COMPLEX*16 array, dimension (LDA, N)
	  On entry, the	N x N upper Hessenberg matrix A.  Entries below	the
	  subdiagonal must be zero.  If	JOB='S', then on exit A	and B will
	  have been simultaneously reduced to upper triangular form.  If
	  JOB='E', then	on exit	A will have been destroyed.

  LDA	  (input) INTEGER
	  The leading dimension	of the array A.	 LDA >=	max( 1,	N ).

  B	  (input/output) COMPLEX*16 array, dimension (LDB, N)
	  On entry, the	N x N upper triangular matrix B.  Entries below	the
	  diagonal must	be zero.  If JOB='S', then on exit A and B will	have
	  been simultaneously reduced to upper triangular form.	 If JOB='E',
	  then on exit B will have been	destroyed.

  LDB	  (input) INTEGER
	  The leading dimension	of the array B.	 LDB >=	max( 1,	N ).

  ALPHA	  (output) COMPLEX*16 array, dimension (N)
	  The diagonal elements	of A when the pair (A,B) has been reduced to
	  Schur	form.  ALPHA(i)/BETA(i)	i=1,...,N are the generalized eigen-
	  values.

  BETA	  (output) COMPLEX*16 array, dimension (N)
	  The diagonal elements	of B when the pair (A,B) has been reduced to
	  Schur	form.  ALPHA(i)/BETA(i)	i=1,...,N are the generalized eigen-
	  values.  A and B are normalized so that BETA(1),...,BETA(N) are
	  non-negative real numbers.

  Q	  (input/output) COMPLEX*16 array, dimension (LDQ, N)
	  If COMPQ='N',	then Q will not	be referenced.	If COMPQ='V' or	'I',
	  then the conjugate transpose of the unitary transformations which
	  are applied to A and B on the	left will be applied to	the array Q
	  on the right.

  LDQ	  (input) INTEGER
	  The leading dimension	of the array Q.	 LDQ must be at	least 1.  If
	  COMPQ='V' or 'I', then LDQ must also be at least N.

  Z	  (input/output) COMPLEX*16 array, dimension (LDZ, N)
	  If COMPZ='N',	then Z will not	be referenced.	If COMPZ='V' or	'I',
	  then the unitary transformations which are applied to	A and B	on
	  the right will be applied to the array Z on the right.

  LDZ	  (input) INTEGER
	  The leading dimension	of the array Z.	 LDZ must be at	least 1.  If
	  COMPZ='V' or 'I', then LDZ must also be at least N.

  WORK	  (workspace) COMPLEX*16 array,	dimension (LWORK)
	  On exit, if INFO is not negative, WORK(1) will be set	to the
	  optimal size of the array WORK.

  LWORK	  (input) INTEGER
	  The number of	elements in WORK.  It must be at least 1.

  RWORK	  (workspace) DOUBLE PRECISION array, dimension	(N)

  INFO	  (output) INTEGER
	  < 0: if INFO = -i, the i-th argument had an illegal value
	  = 0: successful exit.
	  = 1,...,N: the QZ iteration did not converge.	 (A,B) is not in
	  Schur	form, but ALPHA(i) and BETA(i),	i=INFO+1,...,N should be
	  correct.  = N+1,...,2*N: the shift calculation failed.  (A,B)	is
	  not in Schur form, but ALPHA(i) and BETA(i), i=INFO-N+1,...,N
	  should be correct.  >	2*N:	 various "impossible" errors.

FURTHER	DETAILS
  We assume that complex ABS works as long as its value	is less	than over-
  flow.