CTGSJA(l)		LAPACK routine (version	1.1)		    CTGSJA(l)

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

SYNOPSIS

  SUBROUTINE CTGSJA( 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

      REAL	     TOLA, TOLB

      REAL	     ALPHA( * ), BETA( * )

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

PURPOSE
  CTGSJA 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	CGGSVP 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 CTGSJA.  See	Further
	  details.

  A	  (input/output) COMPLEX 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 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) REAL
	  TOLB	  (input) REAL TOLA and	TOLB are the convergence criteria for
	  the Jacobi- Kogbetliantz iteration procedure.	Generally, they	are
	  the same as used in the preprocessing	step, say TOLA =
	  MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS.

  ALPHA	  (output) REAL	array, dimension (N)
	  BETA	  (output) REAL	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 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 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 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 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 ===============

	  CTGSJA 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.