1CGGSVD(1) LAPACK driver routine (version 3.2) CGGSVD(1)
2
6 CGGSVD - 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 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 )
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 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 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) REAL array, dimension (N)
128 BETA (output) REAL array, dimension (N) On exit, ALPHA and
129 BETA contain the generalized singular value pairs of A and B;
130 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 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 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 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 array, dimension (max(3*N,M,P)+N)
162 RWORK (workspace) REAL 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 CTGSJA.
176
178 TOLA REAL
179 TOLB REAL TOLA and TOLB are the thresholds to determine the
180 effective rank of (A',B')'. Generally, they are set to TOLA =
181 MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS. The
182 size of TOLA and TOLB may affect the size of backward errors of
183 the decomposition. Further Details =============== 2-96 Based
184 on modifications by Ming Gu and Huan Ren, Computer Science
185 Division, University of California at Berkeley, USA
186 LAPACK driver routine (version 3.N2o)vember 2008 CGGSVD(1)