1CTGSJA(1) LAPACK routine (version 3.2) CTGSJA(1)
2
6 CTGSJA - computes the generalized singular value decomposition (GSVD)
7 of two 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. On
28 entry, it is assumed that matrices A and B have the following forms,
29 which may be obtained by the preprocessing subroutine CGGSVP from a
30 general M-by-N matrix A and P-by-N matrix B:
31 N-K-L K L
32 A = K ( 0 A12 A13 ) if M-K-L >= 0;
33 L ( 0 0 A23 )
34 M-K-L ( 0 0 0 )
35 N-K-L K L
36 A = K ( 0 A12 A13 ) if M-K-L < 0;
37 M-K ( 0 0 A23 )
38 N-K-L K L
39 B = L ( 0 0 B13 )
40 P-L ( 0 0 0 )
41 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper
42 triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23
43 is (M-K)-by-L upper trapezoidal.
44 On exit,
45 U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ),
46 where U, V and Q are unitary matrices, Z' denotes the conjugate trans‐
47 pose of Z, R is a nonsingular upper triangular matrix, and D1 and D2
48 are ``diagonal'' matrices, which are of the following structures:
49 If M-K-L >= 0,
50 K L
51 D1 = K ( I 0 )
52 L ( 0 C )
53 M-K-L ( 0 0 )
54 K L
55 D2 = L ( 0 S )
56 P-L ( 0 0 )
57 N-K-L K L
58 ( 0 R ) = K ( 0 R11 R12 ) K
59 L ( 0 0 R22 ) L
60 where
61 C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
62 S = diag( BETA(K+1), ... , BETA(K+L) ),
63 C**2 + S**2 = I.
64 R is stored in A(1:K+L,N-K-L+1:N) on exit.
65 If M-K-L < 0,
66 K M-K K+L-M
67 D1 = K ( I 0 0 )
68 M-K ( 0 C 0 )
69 K M-K K+L-M
70 D2 = M-K ( 0 S 0 )
71 K+L-M ( 0 0 I )
72 P-L ( 0 0 0 )
73 N-K-L K M-K K+L-M
74 M-K ( 0 0 R22 R23 )
75 K+L-M ( 0 0 0 R33 )
76 where
77 C = diag( ALPHA(K+1), ... , ALPHA(M) ),
78 S = diag( BETA(K+1), ... , BETA(M) ),
79 C**2 + S**2 = I.
80 R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
81 ( 0 R22 R23 )
82 in B(M-K+1:L,N+M-K-L+1:N) on exit.
83 The computation of the unitary transformation matrices U, V or Q is
84 optional. These matrices may either be formed explicitly, or they may
85 be postmultiplied into input matrices U1, V1, or Q1.
86
88 JOBU (input) CHARACTER*1
89 = 'U': U must contain a unitary matrix U1 on entry, and the
90 product U1*U is returned; = 'I': U is initialized to the unit
91 matrix, and the unitary matrix U is returned; = 'N': U is not
92 computed.
93 JOBV (input) CHARACTER*1
95 = 'V': V must contain a unitary matrix V1 on entry, and the
96 product V1*V is returned; = 'I': V is initialized to the unit
97 matrix, and the unitary matrix V is returned; = 'N': V is not
98 computed.
99 JOBQ (input) CHARACTER*1
101 = 'Q': Q must contain a unitary matrix Q1 on entry, and the
102 product Q1*Q is returned; = 'I': Q is initialized to the unit
103 matrix, and the unitary matrix Q is returned; = 'N': Q is not
104 computed.
105 M (input) INTEGER
107 The number of rows of the matrix A. M >= 0.
108 P (input) INTEGER
110 The number of rows of the matrix B. P >= 0.
111 N (input) INTEGER
113 The number of columns of the matrices A and B. N >= 0.
114 K (input) INTEGER
116 L (input) INTEGER K and L specify the subblocks in the
117 input matrices A and B:
118 A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) of A
119 and B, whose GSVD is going to be computed by CTGSJA. See Fur‐
120 ther Details. A (input/output) COMPLEX array, dimension
121 (LDA,N) On entry, the M-by-N matrix A. On exit, A(N-
122 K+1:N,1:MIN(K+L,M) ) contains the triangular matrix R or part
123 of R. See Purpose for details.
124 LDA (input) INTEGER
126 The leading dimension of the array A. LDA >= max(1,M).
127 B (input/output) COMPLEX array, dimension (LDB,N)
129 On entry, the P-by-N matrix B. On exit, if necessary, B(M-
130 K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for
131 details.
132 LDB (input) INTEGER
134 The leading dimension of the array B. LDB >= max(1,P).
135 TOLA (input) REAL
137 TOLB (input) REAL TOLA and TOLB are the convergence criteria
138 for the Jacobi- Kogbetliantz iteration procedure. Generally,
139 they are the same as used in the preprocessing step, say TOLA =
140 MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS.
141 ALPHA (output) REAL array, dimension (N)
143 BETA (output) REAL array, dimension (N) On exit, ALPHA and
144 BETA contain the generalized singular value pairs of A and B;
145 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
151 BETA(K+L+1:N) = 0.
152 U (input/output) COMPLEX array, dimension (LDU,M)
154 On entry, if JOBU = 'U', U must contain a matrix U1 (usually
155 the unitary matrix returned by CGGSVP). On exit, if JOBU =
156 'I', U contains the unitary matrix U; if JOBU = 'U', U contains
157 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) COMPLEX array, dimension (LDV,P)
164 On entry, if JOBV = 'V', V must contain a matrix V1 (usually
165 the unitary matrix returned by CGGSVP). On exit, if JOBV =
166 'I', V contains the unitary matrix V; if JOBV = 'V', V contains
167 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) COMPLEX array, dimension (LDQ,N)
174 On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
175 the unitary matrix returned by CGGSVP). On exit, if JOBQ =
176 'I', Q contains the unitary matrix Q; if JOBQ = 'Q', Q contains
177 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) COMPLEX 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 =============== CTGSJA
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 unitary matrix, and Z' is the
202 conjugate transpose of Z. C1 and S1 are diagonal matrices sat‐
203 isfying C1**2 + S1**2 = I, and R1 is an L-by-L nonsingular
204 upper triangular matrix.
205 LAPACK routine (version 3.2) November 2008 CTGSJA(1)