1ZGGSVD(1) LAPACK driver routine (version 3.2) ZGGSVD(1)
2
6 ZGGSVD - computes the generalized singular value decomposition (GSVD)
7 of an M-by-N complex matrix A and P-by-N complex matrix B
8
10 SUBROUTINE ZGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB,
11 ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK,
12 IWORK, INFO )
13 CHARACTER JOBQ, JOBU, JOBV
15 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
17 INTEGER IWORK( * )
19 DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * )
21 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ),
23 V( LDV, * ), WORK( * )
24
26 ZGGSVD computes the generalized singular value decomposition (GSVD) of
27 an M-by-N complex matrix A and P-by-N complex matrix B:
28 U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
29 where U, V and Q are unitary matrices, and Z' means the conjugate
30 transpose of Z. Let K+L = the effective numerical rank of the matrix
31 (A',B')', then R is a (K+L)-by-(K+L) nonsingular upper triangular
32 matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and
33 of the following structures, respectively:
34 If M-K-L >= 0,
35 K L
36 D1 = K ( I 0 )
37 L ( 0 C )
38 M-K-L ( 0 0 )
39 K L
40 D2 = L ( 0 S )
41 P-L ( 0 0 )
42 N-K-L K L
43 ( 0 R ) = K ( 0 R11 R12 )
44 L ( 0 0 R22 )
45 where
46 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
47 S = diag( BETA(K+1), ... , BETA(K+L) ),
48 C**2 + S**2 = I.
49 R is stored in A(1:K+L,N-K-L+1:N) on exit.
50 If M-K-L < 0,
51 K M-K K+L-M
52 D1 = K ( I 0 0 )
53 M-K ( 0 C 0 )
54 K M-K K+L-M
55 D2 = M-K ( 0 S 0 )
56 K+L-M ( 0 0 I )
57 P-L ( 0 0 0 )
58 N-K-L K M-K K+L-M
59 ( 0 R ) = K ( 0 R11 R12 R13 )
60 M-K ( 0 0 R22 R23 )
61 K+L-M ( 0 0 0 R33 )
62 where
63 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
64 S = diag( BETA(K+1), ... , BETA(M) ),
65 C**2 + S**2 = I.
66 (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
67 ( 0 R22 R23 )
68 in B(M-K+1:L,N+M-K-L+1:N) on exit.
69 The routine computes C, S, R, and optionally the unitary
70 transformation matrices U, V and Q.
71 In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A
72 and B implicitly gives the SVD of A*inv(B):
73 A*inv(B) = U*(D1*inv(D2))*V'.
74 If ( A',B')' has orthnormal columns, then the GSVD of A and B is also
75 equal to the CS decomposition of A and B. Furthermore, the GSVD can be
76 used to derive the solution of the eigenvalue problem:
77 A'*A x = lambda* B'*B x.
78 In some literature, the GSVD of A and B is presented in the form
79 U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
80 where U and V are orthogonal and X is nonsingular, and D1 and D2 are
81 ``diagonal''. The former GSVD form can be converted to the latter form
82 by taking the nonsingular matrix X as
83 X = Q*( I 0 )
84 ( 0 inv(R) )
85
87 JOBU (input) CHARACTER*1
88 = 'U': Unitary matrix U is computed;
89 = 'N': U is not computed.
90 JOBV (input) CHARACTER*1
92 = 'V': Unitary matrix V is computed;
93 = 'N': V is not computed.
94 JOBQ (input) CHARACTER*1
96 = 'Q': Unitary matrix Q is computed;
97 = 'N': Q is not computed.
98 M (input) INTEGER
100 The number of rows of the matrix A. M >= 0.
101 N (input) INTEGER
103 The number of columns of the matrices A and B. N >= 0.
104 P (input) INTEGER
106 The number of rows of the matrix B. P >= 0.
107 K (output) INTEGER
109 L (output) INTEGER On exit, K and L specify the dimension
110 of the subblocks described in Purpose. K + L = effective
111 numerical rank of (A',B')'.
112 A (input/output) COMPLEX*16 array, dimension (LDA,N)
114 On entry, the M-by-N matrix A. On exit, A contains the trian‐
115 gular matrix R, or part of R. See Purpose for details.
116 LDA (input) INTEGER
118 The leading dimension of the array A. LDA >= max(1,M).
119 B (input/output) COMPLEX*16 array, dimension (LDB,N)
121 On entry, the P-by-N matrix B. On exit, B contains part of the
122 triangular matrix R if M-K-L < 0. See Purpose for details.
123 LDB (input) INTEGER
125 The leading dimension of the array B. LDB >= max(1,P).
126 ALPHA (output) DOUBLE PRECISION array, dimension (N)
128 BETA (output) DOUBLE PRECISION array, dimension (N) On exit,
129 ALPHA and BETA contain the generalized singular value pairs of
130 A and B; ALPHA(1:K) = 1,
131 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C,
132 BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)= C,
133 ALPHA(M+1:K+L)= 0
134 BETA(K+1:M) = S, BETA(M+1:K+L) = 1 and ALPHA(K+L+1:N) = 0
135 BETA(K+L+1:N) = 0
136 U (output) COMPLEX*16 array, dimension (LDU,M)
138 If JOBU = 'U', U contains the M-by-M unitary matrix U. If JOBU
139 = 'N', U is not referenced.
140 LDU (input) INTEGER
142 The leading dimension of the array U. LDU >= max(1,M) if JOBU =
143 'U'; LDU >= 1 otherwise.
144 V (output) COMPLEX*16 array, dimension (LDV,P)
146 If JOBV = 'V', V contains the P-by-P unitary matrix V. If JOBV
147 = 'N', V is not referenced.
148 LDV (input) INTEGER
150 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
151 'V'; LDV >= 1 otherwise.
152 Q (output) COMPLEX*16 array, dimension (LDQ,N)
154 If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. If JOBQ
155 = 'N', Q is not referenced.
156 LDQ (input) INTEGER
158 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
159 'Q'; LDQ >= 1 otherwise.
160 WORK (workspace) COMPLEX*16 array, dimension (max(3*N,M,P)+N)
162 RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
164 IWORK (workspace/output) INTEGER array, dimension (N)
166 On exit, IWORK stores the sorting information. More precisely,
167 the following loop will sort ALPHA for I = K+1, min(M,K+L) swap
168 ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >=
169 ALPHA(2) >= ... >= ALPHA(N).
170 INFO (output) INTEGER
172 = 0: successful exit.
173 < 0: if INFO = -i, the i-th argument had an illegal value.
174 > 0: if INFO = 1, the Jacobi-type procedure failed to con‐
175 verge. For further details, see subroutine ZTGSJA.
176
178 TOLA DOUBLE PRECISION
179 TOLB DOUBLE PRECISION TOLA and TOLB are the thresholds to
180 determine the effective rank of (A',B')'. Generally, they are
181 set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB =
182 MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect
183 the size of backward errors of the decomposition. Further
184 Details =============== 2-96 Based on modifications by Ming Gu
185 and Huan Ren, Computer Science Division, University of Califor‐
186 nia at Berkeley, USA
187 LAPACK driver routine (version 3.N2o)vember 2008 ZGGSVD(1)