1DTGSJA(1) LAPACK routine (version 3.1) DTGSJA(1)
2
6 DTGSJA - the generalized singular value decomposition (GSVD) of two
7 real upper triangular (or trapezoidal) matrices A and B
8
10 SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
11 TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
12 WORK, NCYCLE, INFO )
13 CHARACTER JOBQ, JOBU, JOBV
15 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, NCYCLE, P
17 DOUBLE PRECISION TOLA, TOLB
19 DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ),
21 BETA( * ), Q( LDQ, * ), U( LDU, * ), V( LDV, * ),
22 WORK( * )
23
25 DTGSJA computes the generalized singular value decomposition (GSVD) of
26 two real upper triangular (or trapezoidal) matrices A and B.
27 On entry, it is assumed that matrices A and B have the following forms,
29 which may be obtained by the preprocessing subroutine DGGSVP from a
30 general M-by-N matrix A and P-by-N matrix B:
31 N-K-L K L
33 A = K ( 0 A12 A13 ) if M-K-L >= 0;
34 L ( 0 0 A23 )
35 M-K-L ( 0 0 0 )
36 N-K-L K L
38 A = K ( 0 A12 A13 ) if M-K-L < 0;
39 M-K ( 0 0 A23 )
40 N-K-L K L
42 B = L ( 0 0 B13 )
43 P-L ( 0 0 0 )
44 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper
46 triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23
47 is (M-K)-by-L upper trapezoidal.
48 On exit,
50 U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
52 where U, V and Q are orthogonal matrices, Z' denotes the transpose of
54 Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
55 ``diagonal'' matrices, which are of the following structures:
56 If M-K-L >= 0,
58 K L
60 D1 = K ( I 0 )
61 L ( 0 C )
62 M-K-L ( 0 0 )
63 K L
65 D2 = L ( 0 S )
66 P-L ( 0 0 )
67 N-K-L K L
69 ( 0 R ) = K ( 0 R11 R12 ) K
70 L ( 0 0 R22 ) L
71 where
73 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
75 S = diag( BETA(K+1), ... , BETA(K+L) ),
76 C**2 + S**2 = I.
77 R is stored in A(1:K+L,N-K-L+1:N) on exit.
79 If M-K-L < 0,
81 K M-K K+L-M
83 D1 = K ( I 0 0 )
84 M-K ( 0 C 0 )
85 K M-K K+L-M
87 D2 = M-K ( 0 S 0 )
88 K+L-M ( 0 0 I )
89 P-L ( 0 0 0 )
90 N-K-L K M-K K+L-M
92 M-K ( 0 0 R22 R23 )
94 K+L-M ( 0 0 0 R33 )
95 where
97 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
98 S = diag( BETA(K+1), ... , BETA(M) ),
99 C**2 + S**2 = I.
100 R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
102 ( 0 R22 R23 )
103 in B(M-K+1:L,N+M-K-L+1:N) on exit.
104 The computation of the orthogonal transformation matrices U, V or Q is
106 optional. These matrices may either be formed explicitly, or they may
107 be postmultiplied into input matrices U1, V1, or Q1.
108
111 JOBU (input) CHARACTER*1
112 = 'U': U must contain an orthogonal matrix U1 on entry, and
113 the product U1*U is returned; = 'I': U is initialized to the
114 unit matrix, and the orthogonal matrix U is returned; = 'N': U
115 is not computed.
116 JOBV (input) CHARACTER*1
118 = 'V': V must contain an orthogonal matrix V1 on entry, and
119 the product V1*V is returned; = 'I': V is initialized to the
120 unit matrix, and the orthogonal matrix V is returned; = 'N': V
121 is not computed.
122 JOBQ (input) CHARACTER*1
124 = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
125 the product Q1*Q is returned; = 'I': Q is initialized to the
126 unit matrix, and the orthogonal matrix Q is returned; = 'N': Q
127 is not computed.
128 M (input) INTEGER
130 The number of rows of the matrix A. M >= 0.
131 P (input) INTEGER
133 The number of rows of the matrix B. P >= 0.
134 N (input) INTEGER
136 The number of columns of the matrices A and B. N >= 0.
137 K (input) INTEGER
139 L (input) INTEGER K and L specify the subblocks in the
140 input matrices A and B:
141 A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N) of A
142 and B, whose GSVD is going to be computed by DTGSJA. See Fur‐
143 ther details.
144 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
146 On entry, the M-by-N matrix A. On exit, A(N-K+1:N,1:MIN(K+L,M)
147 ) contains the triangular matrix R or part of R. See Purpose
148 for details.
149 LDA (input) INTEGER
151 The leading dimension of the array A. LDA >= max(1,M).
152 B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
154 On entry, the P-by-N matrix B. On exit, if necessary, B(M-
155 K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for
156 details.
157 LDB (input) INTEGER
159 The leading dimension of the array B. LDB >= max(1,P).
160 TOLA (input) DOUBLE PRECISION
162 TOLB (input) DOUBLE PRECISION TOLA and TOLB are the conver‐
163 gence criteria for the Jacobi- Kogbetliantz iteration proce‐
164 dure. Generally, they are the same as used in the preprocessing
165 step, say TOLA = max(M,N)*norm(A)*MAZHEPS, TOLB =
166 max(P,N)*norm(B)*MAZHEPS.
167 ALPHA (output) DOUBLE PRECISION array, dimension (N)
169 BETA (output) DOUBLE PRECISION array, dimension (N) On exit,
170 ALPHA and BETA contain the generalized singular value pairs of
171 A and B; ALPHA(1:K) = 1,
172 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = diag(C),
173 BETA(K+1:K+L) = diag(S), or if M-K-L < 0, ALPHA(K+1:M)= C,
174 ALPHA(M+1:K+L)= 0
175 BETA(K+1:M) = S, BETA(M+1:K+L) = 1. Furthermore, if K+L < N,
176 ALPHA(K+L+1:N) = 0 and
177 BETA(K+L+1:N) = 0.
178 U (input/output) DOUBLE PRECISION array, dimension (LDU,M)
180 On entry, if JOBU = 'U', U must contain a matrix U1 (usually
181 the orthogonal matrix returned by DGGSVP). On exit, if JOBU =
182 'I', U contains the orthogonal matrix U; if JOBU = 'U', U con‐
183 tains the product U1*U. If JOBU = 'N', U is not referenced.
184 LDU (input) INTEGER
186 The leading dimension of the array U. LDU >= max(1,M) if JOBU =
187 'U'; LDU >= 1 otherwise.
188 V (input/output) DOUBLE PRECISION array, dimension (LDV,P)
190 On entry, if JOBV = 'V', V must contain a matrix V1 (usually
191 the orthogonal matrix returned by DGGSVP). On exit, if JOBV =
192 'I', V contains the orthogonal matrix V; if JOBV = 'V', V con‐
193 tains the product V1*V. If JOBV = 'N', V is not referenced.
194 LDV (input) INTEGER
196 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
197 'V'; LDV >= 1 otherwise.
198 Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
200 On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
201 the orthogonal matrix returned by DGGSVP). On exit, if JOBQ =
202 'I', Q contains the orthogonal matrix Q; if JOBQ = 'Q', Q con‐
203 tains the product Q1*Q. If JOBQ = 'N', Q is not referenced.
204 LDQ (input) INTEGER
206 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
207 'Q'; LDQ >= 1 otherwise.
208 WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
210 NCYCLE (output) INTEGER
212 The number of cycles required for convergence.
213 INFO (output) INTEGER
215 = 0: successful exit
216 < 0: if INFO = -i, the i-th argument had an illegal value.
217 = 1: the procedure does not converge after MAXIT cycles.
218
220 MAXIT INTEGER
221 MAXIT specifies the total loops that the iterative procedure
222 may take. If after MAXIT cycles, the routine fails to converge,
223 we return INFO = 1.
224 Further Details ===============
226 DTGSJA essentially uses a variant of Kogbetliantz algorithm to
228 reduce min(L,M-K)-by-L triangular (or trapezoidal) matrix A23
229 and L-by-L matrix B13 to the form:
230 U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
232 where U1, V1 and Q1 are orthogonal matrix, and Z' is the trans‐
234 pose of Z. C1 and S1 are diagonal matrices satisfying
235 C1**2 + S1**2 = I,
237 and R1 is an L-by-L nonsingular upper triangular matrix.
239 LAPACK routine (version 3.1) November 2006 DTGSJA(1)