1DTGSJA(1) LAPACK routine (version 3.2) DTGSJA(1)
2
6 DTGSJA - computes the generalized singular value decomposition (GSVD)
7 of two 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. On entry,
27 it is assumed that matrices A and B have the following forms, which may
28 be obtained by the preprocessing subroutine DGGSVP from a general M-by-
29 N matrix A and P-by-N matrix B:
30 N-K-L K L
31 A = K ( 0 A12 A13 ) if M-K-L >= 0;
32 L ( 0 0 A23 )
33 M-K-L ( 0 0 0 )
34 N-K-L K L
35 A = K ( 0 A12 A13 ) if M-K-L < 0;
36 M-K ( 0 0 A23 )
37 N-K-L K L
38 B = L ( 0 0 B13 )
39 P-L ( 0 0 0 )
40 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper
41 triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23
42 is (M-K)-by-L upper trapezoidal.
43 On exit,
44 U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
45 where U, V and Q are orthogonal matrices, Z' denotes the transpose of
46 Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
47 ``diagonal'' matrices, which are of the following structures: If M-K-L
48 >= 0,
49 K L
50 D1 = K ( I 0 )
51 L ( 0 C )
52 M-K-L ( 0 0 )
53 K L
54 D2 = L ( 0 S )
55 P-L ( 0 0 )
56 N-K-L K L
57 ( 0 R ) = K ( 0 R11 R12 ) K
58 L ( 0 0 R22 ) L
59 where
60 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
61 S = diag( BETA(K+1), ... , BETA(K+L) ),
62 C**2 + S**2 = I.
63 R is stored in A(1:K+L,N-K-L+1:N) on exit.
64 If M-K-L < 0,
65 K M-K K+L-M
66 D1 = K ( I 0 0 )
67 M-K ( 0 C 0 )
68 K M-K K+L-M
69 D2 = M-K ( 0 S 0 )
70 K+L-M ( 0 0 I )
71 P-L ( 0 0 0 )
72 N-K-L K M-K K+L-M
73 M-K ( 0 0 R22 R23 )
74 K+L-M ( 0 0 0 R33 )
75 where
76 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
77 S = diag( BETA(K+1), ... , BETA(M) ),
78 C**2 + S**2 = I.
79 R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
80 ( 0 R22 R23 )
81 in B(M-K+1:L,N+M-K-L+1:N) on exit.
82 The computation of the orthogonal transformation matrices U, V or Q is
83 optional. These matrices may either be formed explicitly, or they may
84 be postmultiplied into input matrices U1, V1, or Q1.
85
87 JOBU (input) CHARACTER*1
88 = 'U': U must contain an orthogonal matrix U1 on entry, and
89 the product U1*U is returned; = 'I': U is initialized to the
90 unit matrix, and the orthogonal matrix U is returned; = 'N': U
91 is not computed.
92 JOBV (input) CHARACTER*1
94 = 'V': V must contain an orthogonal matrix V1 on entry, and
95 the product V1*V is returned; = 'I': V is initialized to the
96 unit matrix, and the orthogonal matrix V is returned; = 'N': V
97 is not computed.
98 JOBQ (input) CHARACTER*1
100 = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
101 the product Q1*Q is returned; = 'I': Q is initialized to the
102 unit matrix, and the orthogonal matrix Q is returned; = 'N': Q
103 is not computed.
104 M (input) INTEGER
106 The number of rows of the matrix A. M >= 0.
107 P (input) INTEGER
109 The number of rows of the matrix B. P >= 0.
110 N (input) INTEGER
112 The number of columns of the matrices A and B. N >= 0.
113 K (input) INTEGER
115 L (input) INTEGER K and L specify the subblocks in the
116 input matrices A and B:
117 A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N) of A
118 and B, whose GSVD is going to be computed by DTGSJA. See Fur‐
119 ther Details. A (input/output) DOUBLE PRECISION array,
120 dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A(N-
121 K+1:N,1:MIN(K+L,M) ) contains the triangular matrix R or part
122 of R. See Purpose for details.
123 LDA (input) INTEGER
125 The leading dimension of the array A. LDA >= max(1,M).
126 B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
128 On entry, the P-by-N matrix B. On exit, if necessary, B(M-
129 K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for
130 details.
131 LDB (input) INTEGER
133 The leading dimension of the array B. LDB >= max(1,P).
134 TOLA (input) DOUBLE PRECISION
136 TOLB (input) DOUBLE PRECISION TOLA and TOLB are the conver‐
137 gence criteria for the Jacobi- Kogbetliantz iteration proce‐
138 dure. Generally, they are the same as used in the preprocessing
139 step, say TOLA = max(M,N)*norm(A)*MAZHEPS, TOLB =
140 max(P,N)*norm(B)*MAZHEPS.
141 ALPHA (output) DOUBLE PRECISION array, dimension (N)
143 BETA (output) DOUBLE PRECISION array, dimension (N) On exit,
144 ALPHA and BETA contain the generalized singular value pairs of
145 A and B; ALPHA(1:K) = 1,
146 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = diag(C),
147 BETA(K+1:K+L) = diag(S), or if M-K-L < 0, ALPHA(K+1:M)= C,
148 ALPHA(M+1:K+L)= 0
149 BETA(K+1:M) = S, BETA(M+1:K+L) = 1. Furthermore, if K+L < N,
150 ALPHA(K+L+1:N) = 0 and
151 BETA(K+L+1:N) = 0.
152 U (input/output) DOUBLE PRECISION array, dimension (LDU,M)
154 On entry, if JOBU = 'U', U must contain a matrix U1 (usually
155 the orthogonal matrix returned by DGGSVP). On exit, if JOBU =
156 'I', U contains the orthogonal matrix U; if JOBU = 'U', U con‐
157 tains the product U1*U. If 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 (input/output) DOUBLE PRECISION array, dimension (LDV,P)
164 On entry, if JOBV = 'V', V must contain a matrix V1 (usually
165 the orthogonal matrix returned by DGGSVP). On exit, if JOBV =
166 'I', V contains the orthogonal matrix V; if JOBV = 'V', V con‐
167 tains the product V1*V. If JOBV = 'N', V is not referenced.
168 LDV (input) INTEGER
170 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
171 'V'; LDV >= 1 otherwise.
172 Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
174 On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
175 the orthogonal matrix returned by DGGSVP). On exit, if JOBQ =
176 'I', Q contains the orthogonal matrix Q; if JOBQ = 'Q', Q con‐
177 tains the product Q1*Q. If JOBQ = 'N', Q is not referenced.
178 LDQ (input) INTEGER
180 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
181 'Q'; LDQ >= 1 otherwise.
182 WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
184 NCYCLE (output) INTEGER
186 The number of cycles required for convergence.
187 INFO (output) INTEGER
189 = 0: successful exit
190 < 0: if INFO = -i, the i-th argument had an illegal value.
191 = 1: the procedure does not converge after MAXIT cycles.
192
194 MAXIT INTEGER
195 MAXIT specifies the total loops that the iterative procedure
196 may take. If after MAXIT cycles, the routine fails to converge,
197 we return INFO = 1. Further Details =============== DTGSJA
198 essentially uses a variant of Kogbetliantz algorithm to reduce
199 min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-
200 by-L matrix B13 to the form: U1'*A13*Q1 = C1*R1; V1'*B13*Q1 =
201 S1*R1, where U1, V1 and Q1 are orthogonal matrix, and Z' is the
202 transpose of Z. C1 and S1 are diagonal matrices satisfying
203 C1**2 + S1**2 = I, and R1 is an L-by-L nonsingular upper trian‐
204 gular matrix.
205 LAPACK routine (version 3.2) November 2008 DTGSJA(1)