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