ZGGSVD(l)	     LAPACK driver routine (version 1.1)	    ZGGSVD(l)

NAME
  ZGGSVD - compute the generalized singular value decomposition	(GSVD) of the
  M-by-N complex matrix	A and P-by-N complex matrix B

SYNOPSIS

  SUBROUTINE ZGGSVD( JOBU, JOBV, JOBQ, M, N, P,	K, L, A, LDA, B, LDB, ALPHA,
		     BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,	IWORK, INFO )

      CHARACTER	     JOBQ, JOBU, JOBV

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

      INTEGER	     IWORK( * )

      DOUBLE	     PRECISION ALPHA( *	), BETA( * ), RWORK( * )

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

PURPOSE
  ZGGSVD computes the generalized singular value decomposition (GSVD) of the
  M-by-N complex matrix	A and P-by-N complex matrix B:

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

  where	U, V and Q are unitary matrices, R is an upper triangular matrix, and
  Z' means the conjugate transpose of Z.  Let K+L = the	numerical effective
  rank of the matrix (A',B')', then D1 and D2 are M-by-(K+L) and P-by-(K+L)
  "diagonal" matrices and of the following structures, respectively:

  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,
	(  0  R22 R23 )
	(  0   0  R33 )
    (R11 R12 R13 ) is stored in	A(1:M, N-K-L+1:N), and R33 is stored
    ( 0	 R22 R23 )
    in B(M-K+1:L,N+M-K-L+1:N) on exit.

  The routine computes C, S, R,	and optionally the unitary
  transformation matrices U, V and Q.

  In particular, if B is an N-by-N nonsingular matrix, then the	GSVD of	A and
  B implicitly gives the SVD of	the matrix A*inv(B):
		       A*inv(B)	= U*(D1*inv(D2))*V'.
  If ( A',B')' has orthnormal columns, then the	GSVD of	A and B	is also	equal
  to the CS decomposition of A and B. Furthermore, the GSVD can	be used	to
  derive the solution of the eigenvalue	problem:
		       A'*A x =	lambda*	B'*B x.
  In some literature, the GSVD of A and	B is presented in the form
		   U'*A*X = ( 0	D1 ),	V'*B*X = ( 0 D2	)	   (2) where
  U and	V are orthogonal and X is nonsingular, and D1 and D2 are ``diago-
  nal''.  It is	easy to	see that the GSVD form (1) can be converted to the
  form (2) by taking the nonsingular matrix X as

			X = Q*(	 I   0	  )
			      (	 0 inv(R) )

ARGUMENTS

  JOBU	  (input) CHARACTER*1
	  = 'U':  Unitary matrix U is computed;
	  = 'N':  U is not computed.

  JOBV	  (input) CHARACTER*1
	  = 'V':  Unitary matrix V is computed;
	  = 'N':  V is not computed.

  JOBQ	  (input) CHARACTER*1
	  = 'Q':  Unitary matrix Q is computed;
	  = 'N':  Q is not computed.

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

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

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

  K	  (output) INTEGER
	  L	  (output) INTEGER On exit, K and L specify the	dimension of
	  the subblocks	described in Purpose.  K + L = effective numerical
	  rank of (A',B')'.

  A	  (input/output) COMPLEX*16 array, dimension (LDA,N)
	  On entry, the	M-by-N matrix A.  On exit, A 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, B contains part of the
	  triangular matrix R if M-K-L < 0.  See Purpose for details.

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

  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) = C,
	  BETA(1:K)  = ZERO, BETA(K+1:K+L)  = S; or 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.  and
	  ALPHA(K+L+1:N) = ZERO
	  BETA(K+L+1:N)	 = ZERO

  U	  (output) COMPLEX*16 array, dimension (LDU,M)
	  If JOBU = 'U', U contains the	M-by-M unitary matrix U.  If JOBU =
	  'N', U is not	referenced.

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

  V	  (output) COMPLEX*16 array, dimension (LDV,P)
	  If JOBU = 'V', V contains the	P-by-P unitary matrix V.  If JOBV =
	  'N', V is not	referenced.

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

  Q	  (output) COMPLEX*16 array, dimension (LDQ,N)
	  If JOBU = 'Q', Q contains the	N-by-N unitary matrix Q.  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 (MAX(3*N,M,P)+N)

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

  IWORK	  (workspace) INTEGER array, dimension (N)

  INFO	  (output)INTEGER
	  = 0:	successful exit.
	  < 0:	if INFO	= -i, the i-th argument	had an illegal value.
	  > 0:	if INFO	= 1, the Jacobi-type procedure failed to converge.
	  For further details, see subroutine ZTGSJA.

PARAMETERS

  TOLA	  DOUBLE PRECISION
	  TOLB	  DOUBLE PRECISION TOLA	and TOLB are the thresholds to deter-
	  mine the effective rank of (A',B')'. Generally, they are set to
	  TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB	= MAX(P,N)*norm(B)*MAZHEPS.
	  The size of TOLA and TOLB may	affect the size	of backward errors of
	  the decomposition.


Back to the listing of simple driver routines