1CGGSVD(1) LAPACK driver routine (version 3.1) CGGSVD(1)
2
6 CGGSVD - 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 CGGSVD( 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 REAL ALPHA( * ), BETA( * ), RWORK( * )
21 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ),
23 V( LDV, * ), WORK( * )
24
26 CGGSVD 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 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 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) REAL array, dimension (N)
147 BETA (output) REAL array, dimension (N) On exit, ALPHA and
148 BETA contain the generalized singular value pairs of A and B;
149 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 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 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 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 array, dimension (max(3*N,M,P)+N)
181 RWORK (workspace) REAL 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 CTGSJA.
195
197 TOLA REAL
198 TOLB REAL TOLA and TOLB are the thresholds to determine the
199 effective rank of (A',B')'. Generally, they are set to TOLA =
200 MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS. The
201 size of TOLA and TOLB may affect the size of backward errors of
202 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 CGGSVD(1)