1ZGGSVD(1) LAPACK driver routine (version 3.1) ZGGSVD(1)
2
6 ZGGSVD - the generalized singular value decomposition (GSVD) of an M-
7 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 )
30 where U, V and Q are unitary matrices, and Z' means the conjugate
32 transpose of Z. Let K+L = the effective numerical rank of the matrix
33 (A',B')', then R is a (K+L)-by-(K+L) nonsingular upper triangular
34 matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and
35 of the following structures, respectively:
36 If M-K-L >= 0,
38 K L
40 D1 = K ( I 0 )
41 L ( 0 C )
42 M-K-L ( 0 0 )
43 K L
45 D2 = L ( 0 S )
46 P-L ( 0 0 )
47 N-K-L K L
49 ( 0 R ) = K ( 0 R11 R12 )
50 L ( 0 0 R22 )
51 where
52 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
54 S = diag( BETA(K+1), ... , BETA(K+L) ),
55 C**2 + S**2 = I.
56 R is stored in A(1:K+L,N-K-L+1:N) on exit.
58 If M-K-L < 0,
60 K M-K K+L-M
62 D1 = K ( I 0 0 )
63 M-K ( 0 C 0 )
64 K M-K K+L-M
66 D2 = M-K ( 0 S 0 )
67 K+L-M ( 0 0 I )
68 P-L ( 0 0 0 )
69 N-K-L K M-K K+L-M
71 ( 0 R ) = K ( 0 R11 R12 R13 )
72 M-K ( 0 0 R22 R23 )
73 K+L-M ( 0 0 0 R33 )
74 where
76 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
78 S = diag( BETA(K+1), ... , BETA(M) ),
79 C**2 + S**2 = I.
80 (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
82 ( 0 R22 R23 )
83 in B(M-K+1:L,N+M-K-L+1:N) on exit.
84 The routine computes C, S, R, and optionally the unitary
86 transformation matrices U, V and Q.
87 In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A
89 and B implicitly gives the SVD of A*inv(B):
90 A*inv(B) = U*(D1*inv(D2))*V'.
91 If ( A',B')' has orthnormal columns, then the GSVD of A and B is also
92 equal to the CS decomposition of A and B. Furthermore, the GSVD can be
93 used to derive the solution of the eigenvalue problem:
94 A'*A x = lambda* B'*B x.
95 In some literature, the GSVD of A and B is presented in the form
96 U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
97 where U and V are orthogonal and X is nonsingular, and D1 and D2 are
98 ``diagonal''. The former GSVD form can be converted to the latter form
99 by taking the nonsingular matrix X as
100 X = Q*( I 0 )
102 ( 0 inv(R) )
103
106 JOBU (input) CHARACTER*1
107 = 'U': Unitary matrix U is computed;
108 = 'N': U is not computed.
109 JOBV (input) CHARACTER*1
111 = 'V': Unitary matrix V is computed;
112 = 'N': V is not computed.
113 JOBQ (input) CHARACTER*1
115 = 'Q': Unitary matrix Q is computed;
116 = 'N': Q is not computed.
117 M (input) INTEGER
119 The number of rows of the matrix A. M >= 0.
120 N (input) INTEGER
122 The number of columns of the matrices A and B. N >= 0.
123 P (input) INTEGER
125 The number of rows of the matrix B. P >= 0.
126 K (output) INTEGER
128 L (output) INTEGER On exit, K and L specify the dimension
129 of the subblocks described in Purpose. K + L = effective
130 numerical rank of (A',B')'.
131 A (input/output) COMPLEX*16 array, dimension (LDA,N)
133 On entry, the M-by-N matrix A. On exit, A contains the trian‐
134 gular matrix R, or part of R. See Purpose for details.
135 LDA (input) INTEGER
137 The leading dimension of the array A. LDA >= max(1,M).
138 B (input/output) COMPLEX*16 array, dimension (LDB,N)
140 On entry, the P-by-N matrix B. On exit, B contains part of the
141 triangular matrix R if M-K-L < 0. See Purpose for details.
142 LDB (input) INTEGER
144 The leading dimension of the array B. LDB >= max(1,P).
145 ALPHA (output) DOUBLE PRECISION array, dimension (N)
147 BETA (output) DOUBLE PRECISION array, dimension (N) On exit,
148 ALPHA and BETA contain the generalized singular value pairs of
149 A and B; ALPHA(1:K) = 1,
150 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C,
151 BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)= C,
152 ALPHA(M+1:K+L)= 0
153 BETA(K+1:M) = S, BETA(M+1:K+L) = 1 and ALPHA(K+L+1:N) = 0
154 BETA(K+L+1:N) = 0
155 U (output) COMPLEX*16 array, dimension (LDU,M)
157 If JOBU = 'U', U contains the M-by-M unitary matrix U. If JOBU
158 = 'N', U is not referenced.
159 LDU (input) INTEGER
161 The leading dimension of the array U. LDU >= max(1,M) if JOBU =
162 'U'; LDU >= 1 otherwise.
163 V (output) COMPLEX*16 array, dimension (LDV,P)
165 If JOBV = 'V', V contains the P-by-P unitary matrix V. If JOBV
166 = 'N', V is not referenced.
167 LDV (input) INTEGER
169 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
170 'V'; LDV >= 1 otherwise.
171 Q (output) COMPLEX*16 array, dimension (LDQ,N)
173 If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. If JOBQ
174 = 'N', Q is not referenced.
175 LDQ (input) INTEGER
177 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
178 'Q'; LDQ >= 1 otherwise.
179 WORK (workspace) COMPLEX*16 array, dimension (max(3*N,M,P)+N)
181 RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
183 IWORK (workspace/output) INTEGER array, dimension (N)
185 On exit, IWORK stores the sorting information. More precisely,
186 the following loop will sort ALPHA for I = K+1, min(M,K+L) swap
187 ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >=
188 ALPHA(2) >= ... >= ALPHA(N).
189 INFO (output) INTEGER
191 = 0: successful exit.
192 < 0: if INFO = -i, the i-th argument had an illegal value.
193 > 0: if INFO = 1, the Jacobi-type procedure failed to con‐
194 verge. For further details, see subroutine ZTGSJA.
195
197 TOLA DOUBLE PRECISION
198 TOLB DOUBLE PRECISION TOLA and TOLB are the thresholds to
199 determine the effective rank of (A',B')'. Generally, they are
200 set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB =
201 MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect
202 the size of backward errors of the decomposition.
203 Further Details ===============
205 2-96 Based on modifications by Ming Gu and Huan Ren, Computer
207 Science Division, University of California at Berkeley, USA
208 LAPACK driver routine (version 3.N1o)vember 2006 ZGGSVD(1)