1CTGSJA(1) LAPACK routine (version 3.1) CTGSJA(1)
2
6 CTGSJA - the generalized singular value decomposition (GSVD) of two
7 complex upper triangular (or trapezoidal) matrices A and B
8
10 SUBROUTINE CTGSJA( 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 ALPHA( * ), BETA( * )
21 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ),
23 V( LDV, * ), WORK( * )
24
26 CTGSJA computes the generalized singular value decomposition (GSVD) of
27 two complex upper triangular (or trapezoidal) matrices A and B.
28 On entry, it is assumed that matrices A and B have the following forms,
30 which may be obtained by the preprocessing subroutine CGGSVP from a
31 general M-by-N matrix A and P-by-N matrix B:
32 N-K-L K L
34 A = K ( 0 A12 A13 ) if M-K-L >= 0;
35 L ( 0 0 A23 )
36 M-K-L ( 0 0 0 )
37 N-K-L K L
39 A = K ( 0 A12 A13 ) if M-K-L < 0;
40 M-K ( 0 0 A23 )
41 N-K-L K L
43 B = L ( 0 0 B13 )
44 P-L ( 0 0 0 )
45 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper
47 triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23
48 is (M-K)-by-L upper trapezoidal.
49 On exit,
51 U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
53 where U, V and Q are unitary matrices, Z' denotes the conjugate trans‐
55 pose of Z, R is a nonsingular upper triangular matrix, and D1 and D2
56 are ``diagonal'' matrices, which are of the following structures:
57 If M-K-L >= 0,
59 K L
61 D1 = K ( I 0 )
62 L ( 0 C )
63 M-K-L ( 0 0 )
64 K L
66 D2 = L ( 0 S )
67 P-L ( 0 0 )
68 N-K-L K L
70 ( 0 R ) = K ( 0 R11 R12 ) K
71 L ( 0 0 R22 ) L
72 where
74 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
76 S = diag( BETA(K+1), ... , BETA(K+L) ),
77 C**2 + S**2 = I.
78 R is stored in A(1:K+L,N-K-L+1:N) on exit.
80 If M-K-L < 0,
82 K M-K K+L-M
84 D1 = K ( I 0 0 )
85 M-K ( 0 C 0 )
86 K M-K K+L-M
88 D2 = M-K ( 0 S 0 )
89 K+L-M ( 0 0 I )
90 P-L ( 0 0 0 )
91 N-K-L K M-K K+L-M
93 M-K ( 0 0 R22 R23 )
95 K+L-M ( 0 0 0 R33 )
96 where
98 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
99 S = diag( BETA(K+1), ... , BETA(M) ),
100 C**2 + S**2 = I.
101 R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
103 ( 0 R22 R23 )
104 in B(M-K+1:L,N+M-K-L+1:N) on exit.
105 The computation of the unitary transformation matrices U, V or Q is
107 optional. These matrices may either be formed explicitly, or they may
108 be postmultiplied into input matrices U1, V1, or Q1.
109
112 JOBU (input) CHARACTER*1
113 = 'U': U must contain a unitary matrix U1 on entry, and the
114 product U1*U is returned; = 'I': U is initialized to the unit
115 matrix, and the unitary matrix U is returned; = 'N': U is not
116 computed.
117 JOBV (input) CHARACTER*1
119 = 'V': V must contain a unitary matrix V1 on entry, and the
120 product V1*V is returned; = 'I': V is initialized to the unit
121 matrix, and the unitary matrix V is returned; = 'N': V is not
122 computed.
123 JOBQ (input) CHARACTER*1
125 = 'Q': Q must contain a unitary matrix Q1 on entry, and the
126 product Q1*Q is returned; = 'I': Q is initialized to the unit
127 matrix, and the unitary matrix Q is returned; = 'N': Q is not
128 computed.
129 M (input) INTEGER
131 The number of rows of the matrix A. M >= 0.
132 P (input) INTEGER
134 The number of rows of the matrix B. P >= 0.
135 N (input) INTEGER
137 The number of columns of the matrices A and B. N >= 0.
138 K (input) INTEGER
140 L (input) INTEGER K and L specify the subblocks in the
141 input matrices A and B:
142 A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) of A
143 and B, whose GSVD is going to be computed by CTGSJA. See Fur‐
144 ther details.
145 A (input/output) COMPLEX array, dimension (LDA,N)
147 On entry, the M-by-N matrix A. On exit, A(N-K+1:N,1:MIN(K+L,M)
148 ) contains the triangular matrix R or part of R. See Purpose
149 for details.
150 LDA (input) INTEGER
152 The leading dimension of the array A. LDA >= max(1,M).
153 B (input/output) COMPLEX array, dimension (LDB,N)
155 On entry, the P-by-N matrix B. On exit, if necessary, B(M-
156 K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for
157 details.
158 LDB (input) INTEGER
160 The leading dimension of the array B. LDB >= max(1,P).
161 TOLA (input) REAL
163 TOLB (input) REAL TOLA and TOLB are the convergence criteria
164 for the Jacobi- Kogbetliantz iteration procedure. Generally,
165 they are the same as used in the preprocessing step, say TOLA =
166 MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS.
167 ALPHA (output) REAL array, dimension (N)
169 BETA (output) REAL array, dimension (N) On exit, ALPHA and
170 BETA contain the generalized singular value pairs of A and B;
171 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
177 BETA(K+L+1:N) = 0.
178 U (input/output) COMPLEX array, dimension (LDU,M)
180 On entry, if JOBU = 'U', U must contain a matrix U1 (usually
181 the unitary matrix returned by CGGSVP). On exit, if JOBU =
182 'I', U contains the unitary matrix U; if JOBU = 'U', U contains
183 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) COMPLEX array, dimension (LDV,P)
190 On entry, if JOBV = 'V', V must contain a matrix V1 (usually
191 the unitary matrix returned by CGGSVP). On exit, if JOBV =
192 'I', V contains the unitary matrix V; if JOBV = 'V', V contains
193 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) COMPLEX array, dimension (LDQ,N)
200 On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
201 the unitary matrix returned by CGGSVP). On exit, if JOBQ =
202 'I', Q contains the unitary matrix Q; if JOBQ = 'Q', Q contains
203 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) COMPLEX 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 CTGSJA 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 unitary matrix, and Z' is the conjugate
234 transpose 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 CTGSJA(1)