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