zggsvd(3)

NAME

ZGGSVD - compute the generalized singular value decomposi

tion (GSVD) of an M-by-N complex matrix A and P-by-N complex ma

trix 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 decomposi

tion (GSVD) of an M-by-N complex matrix A and P-by-N complex ma

trix B:

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

where U, V and Q are unitary matrices, and Z' means the

conjugate transpose of Z. Let K+L = the effective numerical rank

of the matrix (A',B')', then R is a (K+L)-by-(K+L) nonsingular

upper triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L)

"diagonal" matrices and of the following structures, respective

ly:

If M-K-L >= 0,

K L

D1 = K ( I 0 )

L ( 0 C )

M-K-L ( 0 0 )

K L

D2 = L ( 0 S )

P-L ( 0 0 )

N-K-L K L

( 0 R ) = K ( 0 R11 R12 )

L ( 0 0 R22 )

where

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

R is stored in A(1:K+L,N-K-L+1:N) on exit.

If M-K-L < 0,

K M-K K+L-M

D1 = K ( I 0 0 )

M-K ( 0 C 0 )

K M-K K+L-M

D2 = M-K ( 0 S 0 )

K+L-M ( 0 0 I )

P-L ( 0 0 0 )

N-K-L K M-K K+L-M

( 0 R ) = K ( 0 R11 R12 R13 )

M-K ( 0 0 R22 R23 )

K+L-M ( 0 0 0 R33 )

where

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

(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 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. Further

more, the GSVD can be used to derive the solution of the eigen

value 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 )

where U and V are orthogonal and X is nonsingular, and D1

and D2 are ``diagonal''. The former GSVD form can be converted

to the latter form 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 = ef

fective 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 con

tains the triangular matrix R, or part of R. See Purpose for de

tails.

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 con

tains 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 val

ue pairs of A and B; ALPHA(1:K) = 1,
BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L)

= C,
BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)=

C, ALPHA(M+1:K+L)= 0
BETA(K+1:M) = S, BETA(M+1:K+L) = 1 and AL

PHA(K+L+1:N) = 0
BETA(K+L+1:N) = 0

U (output) COMPLEX*16 array, dimension (LDU,M)

If JOBU = 'U', U contains the M-by-M unitary ma

trix U. If JOBU = 'N', U is not referenced.

LDU (input) INTEGER

The leading dimension of the array U. LDU >=

max(1,M) if JOBU = 'U'; LDU >= 1 otherwise.

V (output) COMPLEX*16 array, dimension (LDV,P)

If JOBV = 'V', V contains the P-by-P unitary ma

trix V. If JOBV = 'N', V is not referenced.

LDV (input) INTEGER

The leading dimension of the array V. LDV >=

max(1,P) if JOBV = 'V'; LDV >= 1 otherwise.

Q (output) COMPLEX*16 array, dimension (LDQ,N)

If JOBQ = 'Q', Q contains the N-by-N unitary ma

trix Q. If JOBQ = 'N', Q is not referenced.

LDQ (input) INTEGER

The leading dimension of the array Q. LDQ >=

max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise.

WORK (workspace) COMPLEX*16 array, dimension

(max(3*N,M,P)+N)

RWORK (workspace) DOUBLE PRECISION array, dimension

(2*N)

IWORK (workspace/output) INTEGER array, dimension (N)

On exit, IWORK stores the sorting information.

More precisely, the following loop will sort ALPHA for I = K+1,

min(M,K+L) swap ALPHA(I) and ALPHA(IWORK(I)) endfor such that AL

PHA(1) >= ALPHA(2) >= ... >= ALPHA(N).

INFO (output)INTEGER

= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an ille

gal 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 determine the effective rank of (A',B')'. General

ly, 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.

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

2-96 Based on modifications by Ming Gu and Huan

Ren, Computer Science Division, University of California at

Berkeley, USA

LAPACK version 3.0 15 June 2000

ZGGSVD(3)