DTGSJA(l)		LAPACK routine (version	1.1)		    DTGSJA(l)

NAME
  DTGSJA - compute the generalized singular value decomposition	(GSVD) of two
  real upper ``triangular (or trapezoidal)'' matrices A	and B

SYNOPSIS

  SUBROUTINE DTGSJA( 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 A( LDA, * ), ALPHA( * ),	B( LDB,	* ), BETA( *
		     ),	Q( LDQ,	* ), U(	LDU, * ), V( LDV, * ), WORK( * )

PURPOSE
  DTGSJA computes the generalized singular value decomposition (GSVD) of two
  real 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	DGGSVP 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 orthogonal matrices, Z' denotes the transpose of	Z, R
  is a nonsingular upper triangular matrix, and	D1 and D2 are ``diagonal''
  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 orthogonal transformation matrices U,	V and Q	are
  optional and may also	be applied to the input	orthogonal matrices U, V and
  Q.

ARGUMENTS

  JOBU	  (input) CHARACTER*1
	  = 'U':  U is overwritten on the input	orthogonal 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	orthogonal 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	orthogonal 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 DTGSJA.  See	Further
	  details.

  A	  (input/output) DOUBLE	PRECISION 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) DOUBLE	PRECISION 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) DOUBLE	PRECISION array, dimension (LDU,M)
	  On entry, if JOBU = 'U', U contains the orthogonal matrix U, On
	  exit,	if JOBU	= 'U', U  is overwritten on the	input orthogonal
	  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) DOUBLE	PRECISION array, dimension (LDV,P)
	  On entry, if JOBV = 'V', V contains the orthogonal matrix V.	On
	  exit,	if JOBV	= 'V', V  is overwritten on the	input orthogonal
	  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) DOUBLE	PRECISION array, dimension (LDQ,N)
	  On entry, if JOBQ = 'Q', Q contains the orthogonal matrix Q, On
	  exit,	if JOBQ	= 'Q', Q  is overwritten on the	input orthogonal
	  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) DOUBLE PRECISION 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 ===============

	  DTGSJA 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 orthogonal matrix, and Z' is the 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.