STGSJA(l)		LAPACK routine (version	1.1)		    STGSJA(l)

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

SYNOPSIS

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

PURPOSE
  STGSJA 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	SGGSVP 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 STGSJA.  See	Further
	  details.

  A	  (input/output) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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) REAL 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 ===============

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