1ZTGSJA(1) LAPACK routine (version 3.2) ZTGSJA(1)
2
6 ZTGSJA - computes the generalized singular value decomposition (GSVD)
7 of two complex upper triangular (or trapezoidal) matrices A and B
8
10 SUBROUTINE ZTGSJA( 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 ALPHA( * ), BETA( * )
21 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ),
23 V( LDV, * ), WORK( * )
24
26 ZTGSJA 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 ZGGSVP 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 ZTGSJA. See Fur‐
120 ther Details. A (input/output) COMPLEX*16 array, dimen‐
121 sion (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*16 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) DOUBLE PRECISION
137 TOLB (input) DOUBLE PRECISION TOLA and TOLB are the conver‐
138 gence criteria for the Jacobi- Kogbetliantz iteration proce‐
139 dure. Generally, they are the same as used in the preprocessing
140 step, say TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB =
141 MAX(P,N)*norm(B)*MAZHEPS.
142 ALPHA (output) DOUBLE PRECISION array, dimension (N)
144 BETA (output) DOUBLE PRECISION array, dimension (N) On exit,
145 ALPHA and BETA contain the generalized singular value pairs of
146 A and B; ALPHA(1:K) = 1,
147 BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = diag(C),
148 BETA(K+1:K+L) = diag(S), or if M-K-L < 0, ALPHA(K+1:M)= C,
149 ALPHA(M+1:K+L)= 0
150 BETA(K+1:M) = S, BETA(M+1:K+L) = 1. Furthermore, if K+L < N,
151 ALPHA(K+L+1:N) = 0
152 BETA(K+L+1:N) = 0.
153 U (input/output) COMPLEX*16 array, dimension (LDU,M)
155 On entry, if JOBU = 'U', U must contain a matrix U1 (usually
156 the unitary matrix returned by ZGGSVP). On exit, if JOBU =
157 'I', U contains the unitary matrix U; if JOBU = 'U', U contains
158 the product U1*U. If JOBU = 'N', U is not referenced.
159 LDU (input) INTEGER
161 The leading dimension of the array U. LDU >= max(1,M) if JOBU =
162 'U'; LDU >= 1 otherwise.
163 V (input/output) COMPLEX*16 array, dimension (LDV,P)
165 On entry, if JOBV = 'V', V must contain a matrix V1 (usually
166 the unitary matrix returned by ZGGSVP). On exit, if JOBV =
167 'I', V contains the unitary matrix V; if JOBV = 'V', V contains
168 the product V1*V. If JOBV = 'N', V is not referenced.
169 LDV (input) INTEGER
171 The leading dimension of the array V. LDV >= max(1,P) if JOBV =
172 'V'; LDV >= 1 otherwise.
173 Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
175 On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
176 the unitary matrix returned by ZGGSVP). On exit, if JOBQ =
177 'I', Q contains the unitary matrix Q; if JOBQ = 'Q', Q contains
178 the product Q1*Q. If JOBQ = 'N', Q is not referenced.
179 LDQ (input) INTEGER
181 The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
182 'Q'; LDQ >= 1 otherwise.
183 WORK (workspace) COMPLEX*16 array, dimension (2*N)
185 NCYCLE (output) INTEGER
187 The number of cycles required for convergence.
188 INFO (output) INTEGER
190 = 0: successful exit
191 < 0: if INFO = -i, the i-th argument had an illegal value.
192 = 1: the procedure does not converge after MAXIT cycles.
193
195 MAXIT INTEGER
196 MAXIT specifies the total loops that the iterative procedure
197 may take. If after MAXIT cycles, the routine fails to converge,
198 we return INFO = 1. Further Details =============== ZTGSJA
199 essentially uses a variant of Kogbetliantz algorithm to reduce
200 min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-
201 by-L matrix B13 to the form: U1'*A13*Q1 = C1*R1; V1'*B13*Q1 =
202 S1*R1, where U1, V1 and Q1 are unitary matrix, and Z' is the
203 conjugate transpose of Z. C1 and S1 are diagonal matrices sat‐
204 isfying C1**2 + S1**2 = I, and R1 is an L-by-L nonsingular
205 upper triangular matrix.
206 LAPACK routine (version 3.2) November 2008 ZTGSJA(1)