1DGGSVD(1) LAPACK driver routine (version 3.1) DGGSVD(1)
2
6 DGGSVD - the generalized singular value decomposition (GSVD) of an M-
7 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 )
29 where U, V and Q are orthogonal matrices, and Z' is the transpose of Z.
31 Let K+L = the effective numerical rank of the matrix (A',B')', then R
32 is a K+L-by-K+L nonsingular upper triangular matrix, D1 and D2 are M-
33 by-(K+L) and P-by-(K+L) "diagonal" matrices and of the following struc‐
34 tures, respectively:
35 If M-K-L >= 0,
37 K L
39 D1 = K ( I 0 )
40 L ( 0 C )
41 M-K-L ( 0 0 )
42 K L
44 D2 = L ( 0 S )
45 P-L ( 0 0 )
46 N-K-L K L
48 ( 0 R ) = K ( 0 R11 R12 )
49 L ( 0 0 R22 )
50 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 orthogonal transforma‐
86 tion 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 orthonormal 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, 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': Orthogonal matrix U is computed;
108 = 'N': U is not computed.
109 JOBV (input) CHARACTER*1
111 = 'V': Orthogonal matrix V is computed;
112 = 'N': V is not computed.
113 JOBQ (input) CHARACTER*1
115 = 'Q': Orthogonal 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 the Purpose section. K + L =
130 effective numerical rank of (A',B')'.
131 A (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDB,N)
140 On entry, the P-by-N matrix B. On exit, B contains the trian‐
141 gular 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) DOUBLE PRECISION array, dimension (LDU,M)
157 If JOBU = 'U', U contains the M-by-M orthogonal matrix U. If
158 JOBU = '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) DOUBLE PRECISION array, dimension (LDV,P)
165 If JOBV = 'V', V contains the P-by-P orthogonal matrix V. If
166 JOBV = '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) DOUBLE PRECISION array, dimension (LDQ,N)
173 If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. If
174 JOBQ = '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) DOUBLE PRECISION array,
181 dimension (max(3*N,M,P)+N)
182 IWORK (workspace/output) INTEGER array, dimension (N)
184 On exit, IWORK stores the sorting information. More precisely,
185 the following loop will sort ALPHA for I = K+1, min(M,K+L) swap
186 ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >=
187 ALPHA(2) >= ... >= ALPHA(N).
188 INFO (output) INTEGER
190 = 0: successful exit
191 < 0: if INFO = -i, the i-th argument had an illegal value.
192 > 0: if INFO = 1, the Jacobi-type procedure failed to con‐
193 verge. For further details, see subroutine DTGSJA.
194
196 TOLA DOUBLE PRECISION
197 TOLB DOUBLE PRECISION TOLA and TOLB are the thresholds to
198 determine the effective rank of (A',B')'. Generally, they are
199 set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB =
200 MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect
201 the size of backward errors of the decomposition.
202 Further Details ===============
204 2-96 Based on modifications by Ming Gu and Huan Ren, Computer
206 Science Division, University of California at Berkeley, USA
207 LAPACK driver routine (version 3.N1o)vember 2006 DGGSVD(1)