ZTGSJA(l)		LAPACK routine (version	1.1)		    ZTGSJA(l)

NAME
  ZTGSJA - compute the generalized singular value decomposition	(GSVD) of two
  complex upper	triangular (or trapezoidal) matrices A and B

SYNOPSIS

  SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N,	K, L, A, LDA, B, LDB, TOLA,
		     TOLB, ALPHA, BETA,	U, LDU,	V, LDV,	Q, LDQ,	WORK, NCYCLE,
		     INFO )

      CHARACTER	     JOBQ, JOBU, JOBV

      INTEGER	     INFO, K, L, LDA, LDB, LDQ,	LDU, LDV, M, N,	NCYCLE,	P

      DOUBLE	     PRECISION TOLA, TOLB

      DOUBLE	     PRECISION ALPHA( *	), BETA( * )

      COMPLEX*16     A(	LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ), V(
		     LDV, * ), WORK( * )

PURPOSE
  ZTGSJA computes the generalized singular value decomposition (GSVD) of two
  complex upper	triangular (or trapezoidal) matrices A and B.

  On entry, it is assumed that matrices	A and B	have the following forms,
  which	may be obtained	by the preprocessing subroutine	ZGGSVP for two gen-
  eral M-by-N matrix A and P-by-N matrix B:

  If M-K-L >= 0

     A = ( 0	A12  A13 ) K	,  B = ( 0     0   B13 ) L
	 ( 0	 0   A23 ) L	       ( 0     0    0  ) P-L
	 ( 0	 0    0	 ) M-K-L	N-K-L  K    L
	  N-K-L	 K    L

  if M-K-L < 0

      A	= ( 0	 A12  A13 ) K  ,   B = ( 0     0   B13 ) L
	  ( 0	  0   A23 ) M-K	       ( 0     0    0  ) P-L
	   N-K-L  K    L		N-K-L  K    L

  where	K-by-K matrix A12 and L-by-L matrix B13	are nonsingular	upper tri-
  angular. A23 is L-by-L upper triangular if M-K-L > 0,	otherwise A23 is L-
  by-(M-K) upper trapezoidal.

  On exit,

	 U'*A*Q	= D1*( 0 R ),	 V'*B*Q	= D2*( 0 R ),

  where	U, V and Q are unitary matrices, Z' denotes the	conjugate transpose
  of Z,	R is a nonsingular upper triangular matrix, and	D1 and D2 are ``diag-
  onal'' matrices, which are of	the following structures:

  If M-K-L >= 0,

     U'*A*Q = D1*( 0 R )

	    = K	    ( I	 0 ) * (  0   R11  R12 ) K
	      L	    ( 0	 C )   (  0    0   R22 ) L
	      M-K-L ( 0	 0 )	N-K-L  K    L
		      K	 L

     V'*B*Q = D2*( 0 R )

	    = L	    ( 0	 S ) * (  0   R11  R12 ) K
	      P-L   ( 0	 0 )   (  0    0   R22 ) L
		      K	 L	N-K-L  K    L
  where

    C =	diag( ALPHA(K+1), ... ,	ALPHA(K+L) ),
    S =	diag( BETA(K+1),  ... ,	BETA(K+L) ),
    C**2 + S**2	= I.

    The	nonsingular triangular matrix R	= ( R11	R12 ) is stored
					  (  0	R22 )
    in A(1:K+L,N-K-L+1:N) on exit.

  If M-K-L < 0,

  U'*A*Q = D1*(	0 R )

	 = K   ( I  0	 0   ) * ( 0	R11  R12  R13  ) K
	   M-K ( 0  C	 0   )	 ( 0	 0   R22  R23  ) M-K
		 K M-K K+L-M	 ( 0	 0    0	  R33  ) K+L-M
				  N-K-L	 K   M-K  K+L-M

  V'*B*Q = D2*(	0 R )

	 = M-K	 ( 0  S	   0   ) * ( 0	  R11  R12  R13	 ) K
	   K+L-M ( 0  0	   I   )   ( 0	   0   R22  R23	 ) M-K
	   P-L	 ( 0  0	   0   )   ( 0	   0	0   R33	 ) K+L-M
		   K M-K K+L-M	    N-K-L  K   M-K  K+L-M

  where
  C = diag( ALPHA(K+1),	... , ALPHA(M) ),
  S = diag( BETA(K+1),	... , BETA(M) ),
  C**2 + S**2 =	I.

  R = (	R11 R12	R13 ) is a nonsingular upper triangular	matrix,	the
      (	 0  R22	R23 )
      (	 0   0	R33 )
  first	M rows of R are	stored in A(1:M, N-K-L+1:N) and	R33 is stored in
  B(M-K+1:L,N+M-K-L+1:N) on exit.

  The computations of the unitary transformation matrices U, V and Q are
  optional and may also	be applied to the input	unitary	matrices U, V and Q.

ARGUMENTS

  JOBU	  (input) CHARACTER*1
	  = 'U':  U is overwritten on the input	unitary	matrix U;
	  = 'I':  U is initialized to the identity matrix;
	  = 'N':  U is not computed.

  JOBV	  (input) CHARACTER*1
	  = 'V':  V is overwritten on the input	unitary	matrix V;
	  = 'I':  V is initialized to the identity matrix;
	  = 'N':  V is not computed.

  JOBQ	  (input) CHARACTER*1
	  = 'Q':  Q is overwritten on the input	unitary	matrix Q;
	  = 'I':  Q is initialized to the identity matrix;
	  = 'N':  Q is not computed.

  M	  (input) INTEGER
	  The number of	rows of	the matrix A.  M >= 0.

  P	  (input) INTEGER
	  The number of	rows of	the matrix B.  P >= 0.

  N	  (input) INTEGER
	  The number of	columns	of the matrices	A and B.  N >= 0.

  K	  (input) INTEGER
	  L	  (input) INTEGER K and	L specify the subblocks	in the input
	  matrices A and B:
	  A23 =	A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) of A and
	  B, whose GSVD	is going to be computed	by ZTGSJA.  See	Further
	  details.

  A	  (input/output) COMPLEX*16 array, dimension (LDA,N)
	  On entry, the	M-by-N matrix A.  On exit, A(N-K+1:N,1:MIN(K+L,M) )
	  contains the triangular matrix R or part of R.  See Purpose for
	  details.

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

  B	  (input/output) COMPLEX*16 array, dimension (LDB,N)
	  On entry, the	P-by-N matrix B.  On exit, if necessary, B(M-
	  K+1:L,N+M-K-L+1:N) contains a	part of	R.  See	Purpose	for details.

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

  TOLA	  (input) DOUBLE PRECISION
	  TOLB	  (input) DOUBLE PRECISION TOLA	and TOLB are the convergence
	  criteria for the Jacobi- Kogbetliantz	iteration procedure. Gen-
	  erally, they are the same as used in the preprocessing step, say
	  TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB	= MAX(P,N)*norm(B)*MAZHEPS.

  ALPHA	  (output) DOUBLE PRECISION array, dimension (N)
	  BETA	  (output) DOUBLE PRECISION array, dimension (N) On exit,
	  ALPHA	and BETA contain the generalized singular value	pairs of A
	  and B; If M-K-L >= 0,	ALPHA(1:K) = ONE,  ALPHA(K+1:K+L) = diag(C),
	  BETA(1:K)  = ZERO, BETA(K+1:K+L)  = diag(S), and if M-K-L < 0,
	  ALPHA(1:K)= ONE,  ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= ZERO
	  BETA(1:K) = ZERO, BETA(K+1:M)	= S, BETA(M+1:K+L) = ONE.  Further-
	  more,	if K+L < N, ALPHA(K+L+1:N) = ZERO
	  BETA(K+L+1:N)	 = ZERO.

  U	  (input/output) COMPLEX*16 array, dimension (LDU,M)
	  On entry, if JOBU = 'U', U contains the unitary matrix U, On exit,
	  if JOBU = 'U', U  is overwritten on the input	unitary	matrix U. If
	  JOBU = 'I', U	is first set to	the identity matrix.  If JOBU =	'N',
	  U is not referenced.

  LDU	  (input) INTEGER
	  The leading dimension	of the array U.	LDU >= max(1,M).

  V	  (input/output) COMPLEX*16 array, dimension (LDV,P)
	  On entry, if JOBV = 'V', V contains the unitary matrix V.  On	exit,
	  if JOBV = 'V', V  is overwritten on the input	unitary	matrix V. If
	  JOBV = 'I', U	is first set to	the identity matrix.  If JOBV =	'N',
	  V is not referenced.

  LDV	  (input) INTEGER
	  The leading dimension	of the array V.	LDV >= max(1,P).

  Q	  (input/output) COMPLEX*16 array, dimension (LDQ,N)
	  On entry, if JOBQ = 'Q', Q contains the unitary matrix Q.  On	exit,
	  if JOBQ = 'Q', Q  is overwritten on the input	unitary	matrix Q. If
	  JOBQ = 'I', Q	is first set to	the identity matrix.  If JOBQ =	'N',
	  Q is not referenced.

  LDQ	  (input) INTEGER
	  The leading dimension	of the array Q.	LDQ >= MAX(1,N).

  WORK	  (workspace) COMPLEX*16 array,	dimension (2*N)

  NCYCLE  (output) INTEGER
	  The number of	cycles required	for convergence.

  INFO	  (output) INTEGER
	  = 0:	successful exit
	  < 0:	if INFO	= -i, the i-th argument	had an illegal value.
	  = 1:	the procedure does not converge	after MAXIT cycles.

PARAMETERS

  MAXIT	  INTEGER
	  MAXIT	specifies the total loops that the iterative procedure may
	  take.	If after MAXIT cycles, the routine fails to converge, we
	  return INFO =	1.

	  Further Details ===============

	  ZTGSJA essentially uses a variant of Kogbetliantz algorithm to
	  reduce min(L,M-K)-by-L triangular (or	trapezoidal) matrix A23	and
	  L-by-L matrix	B13 to the form:

	  U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,

	  where	U1, V1 and Q1 are unitary matrix, and Z' is the	conjugate
	  transpose of Z.  C1 and S1 are diagonal matrices satisfying

	  C1**2	+ S1**2	= I,

	  and R1 is an L-by-L nonsingular upper	triangular matrix.