sggsvd(3)

NAME

SGGSVD - compute the generalized singular value decomposi

tion (GSVD) of an M-by-N real matrix A and P-by-N real matrix B

SYNOPSIS

SUBROUTINE  SGGSVD(  JOBU,  JOBV,  JOBQ, M, N, P, K, L, A,
LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q,  LDQ,  WORK,  IWORK,
INFO )
    CHARACTER      JOBQ, JOBU, JOBV
    INTEGER        INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M,
N, P
    INTEGER        IWORK( * )
    REAL           A( LDA, * ), ALPHA( * ), B( LDB,  *  ),
BETA( * ), Q( LDQ, * ), U( LDU, * ), V( LDV, * ), WORK( * )

PURPOSE

SGGSVD computes the generalized singular value decomposi

tion (GSVD) of an M-by-N real matrix A and P-by-N real matrix B:

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

where U, V and Q are orthogonal matrices, and Z' is the

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 trian

gular matrix, 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,

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 orthogo

nal 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 orthonormal 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, 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': Orthogonal matrix U is computed;
= 'N': U is not computed.

JOBV (input) CHARACTER*1

= 'V': Orthogonal matrix V is computed;
= 'N': V is not computed.

JOBQ (input) CHARACTER*1

= 'Q': Orthogonal 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 the Purpose section.

K + L = effective numerical rank of (A',B')'.

A (input/output) REAL 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) REAL array, dimension (LDB,N)

On entry, the P-by-N matrix B. On exit, B con

tains the triangular matrix R if M-K-L < 0. See Purpose for de

tails.

LDB (input) INTEGER

The leading dimension of the array B. LDA >=

max(1,P).

ALPHA (output) REAL array, dimension (N)

BETA (output) REAL array, dimension (N) On ex

it, ALPHA and BETA contain the generalized singular value 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, AL

PHA(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) REAL array, dimension (LDU,M)

If JOBU = 'U', U contains the M-by-M orthogonal

matrix 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) REAL array, dimension (LDV,P)

If JOBV = 'V', V contains the P-by-P orthogonal

matrix 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) REAL array, dimension (LDQ,N)

If JOBQ = 'Q', Q contains the N-by-N orthogonal

matrix 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) REAL array,

dimension (max(3*N,M,P)+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 STGSJA.

PARAMETERS

TOLA REAL

TOLB REAL TOLA and TOLB are the thresholds to

determine the effective rank of (A',B')'. Generally, they are set

to TOLA = MAX(M,N)*norm(A)*MACHEPS, TOLB =

MAX(P,N)*norm(B)*MACHEPS. 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

SGGSVD(3)