1SGGSVD(1) LAPACK driver routine (version 3.1) SGGSVD(1)
2
6 SGGSVD - 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 SGGSVD( 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 REAL A( LDA, * ), ALPHA( * ), B( LDB, * ), BETA( * ), Q(
21 LDQ, * ), U( LDU, * ), V( LDV, * ), WORK( * )
22
24 SGGSVD computes the generalized singular value decomposition (GSVD) of
25 an M-by-N real matrix A and P-by-N real matrix B:
26 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.
30 Let K+L = the effective numerical rank of the matrix (A',B')', then R
31 is a K+L-by-K+L nonsingular upper triangular matrix, D1 and D2 are M-
32 by-(K+L) and P-by-(K+L) "diagonal" matrices and of the following struc‐
33 tures, respectively:
34 If M-K-L >= 0,
36 K L
38 D1 = K ( I 0 )
39 L ( 0 C )
40 M-K-L ( 0 0 )
41 K L
43 D2 = L ( 0 S )
44 P-L ( 0 0 )
45 N-K-L K L
47 ( 0 R ) = K ( 0 R11 R12 )
48 L ( 0 0 R22 )
49 where
51 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
53 S = diag( BETA(K+1), ... , BETA(K+L) ),
54 C**2 + S**2 = I.
55 R is stored in A(1:K+L,N-K-L+1:N) on exit.
57 If M-K-L < 0,
59 K M-K K+L-M
61 D1 = K ( I 0 0 )
62 M-K ( 0 C 0 )
63 K M-K K+L-M
65 D2 = M-K ( 0 S 0 )
66 K+L-M ( 0 0 I )
67 P-L ( 0 0 0 )
68 N-K-L K M-K K+L-M
70 ( 0 R ) = K ( 0 R11 R12 R13 )
71 M-K ( 0 0 R22 R23 )
72 K+L-M ( 0 0 0 R33 )
73 where
75 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
77 S = diag( BETA(K+1), ... , BETA(M) ),
78 C**2 + S**2 = I.
79 (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
81 ( 0 R22 R23 )
82 in B(M-K+1:L,N+M-K-L+1:N) on exit.
83 The routine computes C, S, R, and optionally the orthogonal transforma‐
85 tion matrices U, V and Q.
86 In particular, if B is an N-by-N nonsingular matrix, then the GSVD of A
88 and B implicitly gives the SVD of A*inv(B):
89 A*inv(B) = U*(D1*inv(D2))*V'.
90 If ( A',B')' has orthonormal columns, then the GSVD of A and B is also
91 equal to the CS decomposition of A and B. Furthermore, the GSVD can be
92 used to derive the solution of the eigenvalue problem:
93 A'*A x = lambda* B'*B x.
94 In some literature, the GSVD of A and B is presented in the form
95 U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
96 where U and V are orthogonal and X is nonsingular, D1 and D2 are
97 ``diagonal''. The former GSVD form can be converted to the latter form
98 by taking the nonsingular matrix X as
99 X = Q*( I 0 )
101 ( 0 inv(R) ).
102
105 JOBU (input) CHARACTER*1
106 = 'U': Orthogonal matrix U is computed;
107 = 'N': U is not computed.
108 JOBV (input) CHARACTER*1
110 = 'V': Orthogonal matrix V is computed;
111 = 'N': V is not computed.
112 JOBQ (input) CHARACTER*1
114 = 'Q': Orthogonal matrix Q is computed;
115 = 'N': Q is not computed.
116 M (input) INTEGER
118 The number of rows of the matrix A. M >= 0.
119 N (input) INTEGER
121 The number of columns of the matrices A and B. N >= 0.
122 P (input) INTEGER
124 The number of rows of the matrix B. P >= 0.
125 K (output) INTEGER
127 L (output) INTEGER On exit, K and L specify the dimension
128 of the subblocks described in the Purpose section. K + L =
129 effective numerical rank of (A',B')'.
130 A (input/output) REAL array, dimension (LDA,N)
132 On entry, the M-by-N matrix A. On exit, A contains the trian‐
133 gular matrix R, or part of R. See Purpose for details.
134 LDA (input) INTEGER
136 The leading dimension of the array A. LDA >= max(1,M).
137 B (input/output) REAL array, dimension (LDB,N)
139 On entry, the P-by-N matrix B. On exit, B contains the trian‐
140 gular matrix R if M-K-L < 0. See Purpose for details.
141 LDB (input) INTEGER
143 The leading dimension of the array B. LDB >= max(1,P).
144 ALPHA (output) REAL array, dimension (N)
146 BETA (output) REAL array, dimension (N) On exit, ALPHA and
147 BETA contain the generalized singular value pairs of A and B;
148 ALPHA(1:K) = 1,
149 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C,
150 BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)=C,
151 ALPHA(M+1:K+L)=0
152 BETA(K+1:M) =S, BETA(M+1:K+L) =1 and ALPHA(K+L+1:N) = 0
153 BETA(K+L+1:N) = 0
154 U (output) REAL array, dimension (LDU,M)
156 If JOBU = 'U', U contains the M-by-M orthogonal matrix U. If
157 JOBU = 'N', U is not referenced.
158 LDU (input) INTEGER
160 The leading dimension of the array U. LDU >= max(1,M) if JOBU =
161 'U'; LDU >= 1 otherwise.
162 V (output) REAL array, dimension (LDV,P)
164 If JOBV = 'V', V contains the P-by-P orthogonal matrix V. If
165 JOBV = 'N', V is not referenced.
166 LDV (input) INTEGER
168 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
169 'V'; LDV >= 1 otherwise.
170 Q (output) REAL array, dimension (LDQ,N)
172 If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. If
173 JOBQ = 'N', Q is not referenced.
174 LDQ (input) INTEGER
176 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
177 'Q'; LDQ >= 1 otherwise.
178 WORK (workspace) REAL array,
180 dimension (max(3*N,M,P)+N)
181 IWORK (workspace/output) INTEGER array, dimension (N)
183 On exit, IWORK stores the sorting information. More precisely,
184 the following loop will sort ALPHA for I = K+1, min(M,K+L) swap
185 ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >=
186 ALPHA(2) >= ... >= ALPHA(N).
187 INFO (output) INTEGER
189 = 0: successful exit
190 < 0: if INFO = -i, the i-th argument had an illegal value.
191 > 0: if INFO = 1, the Jacobi-type procedure failed to con‐
192 verge. For further details, see subroutine STGSJA.
193
195 TOLA REAL
196 TOLB REAL TOLA and TOLB are the thresholds to determine the
197 effective rank of (A',B')'. Generally, they are set to TOLA =
198 MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS. The
199 size of TOLA and TOLB may affect the size of backward errors of
200 the decomposition.
201 Further Details ===============
203 2-96 Based on modifications by Ming Gu and Huan Ren, Computer
205 Science Division, University of California at Berkeley, USA
206 LAPACK driver routine (version 3.N1o)vember 2006 SGGSVD(1)