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