|
2329
|
1 SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, |
|
|
2 $ LDVR, MM, M, WORK, RWORK, INFO ) |
|
|
3 * |
|
7034
|
4 * -- LAPACK routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 CHARACTER HOWMNY, SIDE |
|
|
10 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N |
|
|
11 * .. |
|
|
12 * .. Array Arguments .. |
|
|
13 LOGICAL SELECT( * ) |
|
|
14 DOUBLE PRECISION RWORK( * ) |
|
|
15 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), |
|
|
16 $ WORK( * ) |
|
|
17 * .. |
|
|
18 * |
|
|
19 * Purpose |
|
|
20 * ======= |
|
|
21 * |
|
|
22 * ZTREVC computes some or all of the right and/or left eigenvectors of |
|
|
23 * a complex upper triangular matrix T. |
|
7034
|
24 * Matrices of this type are produced by the Schur factorization of |
|
|
25 * a complex general matrix: A = Q*T*Q**H, as computed by ZHSEQR. |
|
|
26 * |
|
2329
|
27 * The right eigenvector x and the left eigenvector y of T corresponding |
|
|
28 * to an eigenvalue w are defined by: |
|
7034
|
29 * |
|
|
30 * T*x = w*x, (y**H)*T = w*(y**H) |
|
|
31 * |
|
|
32 * where y**H denotes the conjugate transpose of the vector y. |
|
|
33 * The eigenvalues are not input to this routine, but are read directly |
|
|
34 * from the diagonal of T. |
|
|
35 * |
|
|
36 * This routine returns the matrices X and/or Y of right and left |
|
|
37 * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an |
|
|
38 * input matrix. If Q is the unitary factor that reduces a matrix A to |
|
|
39 * Schur form T, then Q*X and Q*Y are the matrices of right and left |
|
|
40 * eigenvectors of A. |
|
2329
|
41 * |
|
|
42 * Arguments |
|
|
43 * ========= |
|
|
44 * |
|
|
45 * SIDE (input) CHARACTER*1 |
|
|
46 * = 'R': compute right eigenvectors only; |
|
|
47 * = 'L': compute left eigenvectors only; |
|
|
48 * = 'B': compute both right and left eigenvectors. |
|
|
49 * |
|
|
50 * HOWMNY (input) CHARACTER*1 |
|
|
51 * = 'A': compute all right and/or left eigenvectors; |
|
|
52 * = 'B': compute all right and/or left eigenvectors, |
|
7034
|
53 * backtransformed using the matrices supplied in |
|
|
54 * VR and/or VL; |
|
2329
|
55 * = 'S': compute selected right and/or left eigenvectors, |
|
7034
|
56 * as indicated by the logical array SELECT. |
|
2329
|
57 * |
|
|
58 * SELECT (input) LOGICAL array, dimension (N) |
|
|
59 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be |
|
|
60 * computed. |
|
7034
|
61 * The eigenvector corresponding to the j-th eigenvalue is |
|
|
62 * computed if SELECT(j) = .TRUE.. |
|
|
63 * Not referenced if HOWMNY = 'A' or 'B'. |
|
2329
|
64 * |
|
|
65 * N (input) INTEGER |
|
|
66 * The order of the matrix T. N >= 0. |
|
|
67 * |
|
|
68 * T (input/output) COMPLEX*16 array, dimension (LDT,N) |
|
|
69 * The upper triangular matrix T. T is modified, but restored |
|
|
70 * on exit. |
|
|
71 * |
|
|
72 * LDT (input) INTEGER |
|
|
73 * The leading dimension of the array T. LDT >= max(1,N). |
|
|
74 * |
|
|
75 * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) |
|
|
76 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must |
|
|
77 * contain an N-by-N matrix Q (usually the unitary matrix Q of |
|
|
78 * Schur vectors returned by ZHSEQR). |
|
|
79 * On exit, if SIDE = 'L' or 'B', VL contains: |
|
|
80 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; |
|
|
81 * if HOWMNY = 'B', the matrix Q*Y; |
|
|
82 * if HOWMNY = 'S', the left eigenvectors of T specified by |
|
|
83 * SELECT, stored consecutively in the columns |
|
|
84 * of VL, in the same order as their |
|
|
85 * eigenvalues. |
|
7034
|
86 * Not referenced if SIDE = 'R'. |
|
2329
|
87 * |
|
|
88 * LDVL (input) INTEGER |
|
7034
|
89 * The leading dimension of the array VL. LDVL >= 1, and if |
|
|
90 * SIDE = 'L' or 'B', LDVL >= N. |
|
2329
|
91 * |
|
|
92 * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) |
|
|
93 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must |
|
|
94 * contain an N-by-N matrix Q (usually the unitary matrix Q of |
|
|
95 * Schur vectors returned by ZHSEQR). |
|
|
96 * On exit, if SIDE = 'R' or 'B', VR contains: |
|
|
97 * if HOWMNY = 'A', the matrix X of right eigenvectors of T; |
|
|
98 * if HOWMNY = 'B', the matrix Q*X; |
|
|
99 * if HOWMNY = 'S', the right eigenvectors of T specified by |
|
|
100 * SELECT, stored consecutively in the columns |
|
|
101 * of VR, in the same order as their |
|
|
102 * eigenvalues. |
|
7034
|
103 * Not referenced if SIDE = 'L'. |
|
2329
|
104 * |
|
|
105 * LDVR (input) INTEGER |
|
7034
|
106 * The leading dimension of the array VR. LDVR >= 1, and if |
|
|
107 * SIDE = 'R' or 'B'; LDVR >= N. |
|
2329
|
108 * |
|
|
109 * MM (input) INTEGER |
|
|
110 * The number of columns in the arrays VL and/or VR. MM >= M. |
|
|
111 * |
|
|
112 * M (output) INTEGER |
|
|
113 * The number of columns in the arrays VL and/or VR actually |
|
|
114 * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M |
|
|
115 * is set to N. Each selected eigenvector occupies one |
|
|
116 * column. |
|
|
117 * |
|
|
118 * WORK (workspace) COMPLEX*16 array, dimension (2*N) |
|
|
119 * |
|
|
120 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) |
|
|
121 * |
|
|
122 * INFO (output) INTEGER |
|
|
123 * = 0: successful exit |
|
|
124 * < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
125 * |
|
|
126 * Further Details |
|
|
127 * =============== |
|
|
128 * |
|
|
129 * The algorithm used in this program is basically backward (forward) |
|
|
130 * substitution, with scaling to make the the code robust against |
|
|
131 * possible overflow. |
|
|
132 * |
|
|
133 * Each eigenvector is normalized so that the element of largest |
|
|
134 * magnitude has magnitude 1; here the magnitude of a complex number |
|
|
135 * (x,y) is taken to be |x| + |y|. |
|
|
136 * |
|
|
137 * ===================================================================== |
|
|
138 * |
|
|
139 * .. Parameters .. |
|
|
140 DOUBLE PRECISION ZERO, ONE |
|
|
141 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
|
|
142 COMPLEX*16 CMZERO, CMONE |
|
|
143 PARAMETER ( CMZERO = ( 0.0D+0, 0.0D+0 ), |
|
|
144 $ CMONE = ( 1.0D+0, 0.0D+0 ) ) |
|
|
145 * .. |
|
|
146 * .. Local Scalars .. |
|
|
147 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV |
|
|
148 INTEGER I, II, IS, J, K, KI |
|
|
149 DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL |
|
|
150 COMPLEX*16 CDUM |
|
|
151 * .. |
|
|
152 * .. External Functions .. |
|
|
153 LOGICAL LSAME |
|
|
154 INTEGER IZAMAX |
|
|
155 DOUBLE PRECISION DLAMCH, DZASUM |
|
|
156 EXTERNAL LSAME, IZAMAX, DLAMCH, DZASUM |
|
|
157 * .. |
|
|
158 * .. External Subroutines .. |
|
3333
|
159 EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS |
|
2329
|
160 * .. |
|
|
161 * .. Intrinsic Functions .. |
|
|
162 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX |
|
|
163 * .. |
|
|
164 * .. Statement Functions .. |
|
|
165 DOUBLE PRECISION CABS1 |
|
|
166 * .. |
|
|
167 * .. Statement Function definitions .. |
|
|
168 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) |
|
|
169 * .. |
|
|
170 * .. Executable Statements .. |
|
|
171 * |
|
|
172 * Decode and test the input parameters |
|
|
173 * |
|
|
174 BOTHV = LSAME( SIDE, 'B' ) |
|
|
175 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV |
|
|
176 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV |
|
|
177 * |
|
|
178 ALLV = LSAME( HOWMNY, 'A' ) |
|
3333
|
179 OVER = LSAME( HOWMNY, 'B' ) |
|
2329
|
180 SOMEV = LSAME( HOWMNY, 'S' ) |
|
|
181 * |
|
|
182 * Set M to the number of columns required to store the selected |
|
|
183 * eigenvectors. |
|
|
184 * |
|
|
185 IF( SOMEV ) THEN |
|
|
186 M = 0 |
|
|
187 DO 10 J = 1, N |
|
|
188 IF( SELECT( J ) ) |
|
|
189 $ M = M + 1 |
|
|
190 10 CONTINUE |
|
|
191 ELSE |
|
|
192 M = N |
|
|
193 END IF |
|
|
194 * |
|
|
195 INFO = 0 |
|
|
196 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN |
|
|
197 INFO = -1 |
|
|
198 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN |
|
|
199 INFO = -2 |
|
|
200 ELSE IF( N.LT.0 ) THEN |
|
|
201 INFO = -4 |
|
|
202 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN |
|
|
203 INFO = -6 |
|
|
204 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN |
|
|
205 INFO = -8 |
|
|
206 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN |
|
|
207 INFO = -10 |
|
|
208 ELSE IF( MM.LT.M ) THEN |
|
|
209 INFO = -11 |
|
|
210 END IF |
|
|
211 IF( INFO.NE.0 ) THEN |
|
|
212 CALL XERBLA( 'ZTREVC', -INFO ) |
|
|
213 RETURN |
|
|
214 END IF |
|
|
215 * |
|
|
216 * Quick return if possible. |
|
|
217 * |
|
|
218 IF( N.EQ.0 ) |
|
|
219 $ RETURN |
|
|
220 * |
|
|
221 * Set the constants to control overflow. |
|
|
222 * |
|
|
223 UNFL = DLAMCH( 'Safe minimum' ) |
|
|
224 OVFL = ONE / UNFL |
|
|
225 CALL DLABAD( UNFL, OVFL ) |
|
|
226 ULP = DLAMCH( 'Precision' ) |
|
|
227 SMLNUM = UNFL*( N / ULP ) |
|
|
228 * |
|
|
229 * Store the diagonal elements of T in working array WORK. |
|
|
230 * |
|
|
231 DO 20 I = 1, N |
|
|
232 WORK( I+N ) = T( I, I ) |
|
|
233 20 CONTINUE |
|
|
234 * |
|
|
235 * Compute 1-norm of each column of strictly upper triangular |
|
|
236 * part of T to control overflow in triangular solver. |
|
|
237 * |
|
|
238 RWORK( 1 ) = ZERO |
|
|
239 DO 30 J = 2, N |
|
|
240 RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 ) |
|
|
241 30 CONTINUE |
|
|
242 * |
|
|
243 IF( RIGHTV ) THEN |
|
|
244 * |
|
|
245 * Compute right eigenvectors. |
|
|
246 * |
|
|
247 IS = M |
|
|
248 DO 80 KI = N, 1, -1 |
|
|
249 * |
|
|
250 IF( SOMEV ) THEN |
|
|
251 IF( .NOT.SELECT( KI ) ) |
|
|
252 $ GO TO 80 |
|
|
253 END IF |
|
|
254 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) |
|
|
255 * |
|
|
256 WORK( 1 ) = CMONE |
|
|
257 * |
|
|
258 * Form right-hand side. |
|
|
259 * |
|
|
260 DO 40 K = 1, KI - 1 |
|
|
261 WORK( K ) = -T( K, KI ) |
|
|
262 40 CONTINUE |
|
|
263 * |
|
|
264 * Solve the triangular system: |
|
|
265 * (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. |
|
|
266 * |
|
|
267 DO 50 K = 1, KI - 1 |
|
|
268 T( K, K ) = T( K, K ) - T( KI, KI ) |
|
|
269 IF( CABS1( T( K, K ) ).LT.SMIN ) |
|
|
270 $ T( K, K ) = SMIN |
|
|
271 50 CONTINUE |
|
|
272 * |
|
|
273 IF( KI.GT.1 ) THEN |
|
|
274 CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', |
|
|
275 $ KI-1, T, LDT, WORK( 1 ), SCALE, RWORK, |
|
|
276 $ INFO ) |
|
|
277 WORK( KI ) = SCALE |
|
|
278 END IF |
|
|
279 * |
|
|
280 * Copy the vector x or Q*x to VR and normalize. |
|
|
281 * |
|
|
282 IF( .NOT.OVER ) THEN |
|
|
283 CALL ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 ) |
|
|
284 * |
|
|
285 II = IZAMAX( KI, VR( 1, IS ), 1 ) |
|
|
286 REMAX = ONE / CABS1( VR( II, IS ) ) |
|
|
287 CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 ) |
|
|
288 * |
|
|
289 DO 60 K = KI + 1, N |
|
|
290 VR( K, IS ) = CMZERO |
|
|
291 60 CONTINUE |
|
|
292 ELSE |
|
|
293 IF( KI.GT.1 ) |
|
|
294 $ CALL ZGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ), |
|
|
295 $ 1, DCMPLX( SCALE ), VR( 1, KI ), 1 ) |
|
|
296 * |
|
|
297 II = IZAMAX( N, VR( 1, KI ), 1 ) |
|
|
298 REMAX = ONE / CABS1( VR( II, KI ) ) |
|
|
299 CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 ) |
|
|
300 END IF |
|
|
301 * |
|
|
302 * Set back the original diagonal elements of T. |
|
|
303 * |
|
|
304 DO 70 K = 1, KI - 1 |
|
|
305 T( K, K ) = WORK( K+N ) |
|
|
306 70 CONTINUE |
|
|
307 * |
|
|
308 IS = IS - 1 |
|
|
309 80 CONTINUE |
|
|
310 END IF |
|
|
311 * |
|
|
312 IF( LEFTV ) THEN |
|
|
313 * |
|
|
314 * Compute left eigenvectors. |
|
|
315 * |
|
|
316 IS = 1 |
|
|
317 DO 130 KI = 1, N |
|
|
318 * |
|
|
319 IF( SOMEV ) THEN |
|
|
320 IF( .NOT.SELECT( KI ) ) |
|
|
321 $ GO TO 130 |
|
|
322 END IF |
|
|
323 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) |
|
|
324 * |
|
|
325 WORK( N ) = CMONE |
|
|
326 * |
|
|
327 * Form right-hand side. |
|
|
328 * |
|
|
329 DO 90 K = KI + 1, N |
|
|
330 WORK( K ) = -DCONJG( T( KI, K ) ) |
|
|
331 90 CONTINUE |
|
|
332 * |
|
|
333 * Solve the triangular system: |
|
|
334 * (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. |
|
|
335 * |
|
|
336 DO 100 K = KI + 1, N |
|
|
337 T( K, K ) = T( K, K ) - T( KI, KI ) |
|
|
338 IF( CABS1( T( K, K ) ).LT.SMIN ) |
|
|
339 $ T( K, K ) = SMIN |
|
|
340 100 CONTINUE |
|
|
341 * |
|
|
342 IF( KI.LT.N ) THEN |
|
|
343 CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', |
|
|
344 $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, |
|
|
345 $ WORK( KI+1 ), SCALE, RWORK, INFO ) |
|
|
346 WORK( KI ) = SCALE |
|
|
347 END IF |
|
|
348 * |
|
|
349 * Copy the vector x or Q*x to VL and normalize. |
|
|
350 * |
|
|
351 IF( .NOT.OVER ) THEN |
|
|
352 CALL ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 ) |
|
|
353 * |
|
|
354 II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 |
|
|
355 REMAX = ONE / CABS1( VL( II, IS ) ) |
|
|
356 CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) |
|
|
357 * |
|
|
358 DO 110 K = 1, KI - 1 |
|
|
359 VL( K, IS ) = CMZERO |
|
|
360 110 CONTINUE |
|
|
361 ELSE |
|
|
362 IF( KI.LT.N ) |
|
|
363 $ CALL ZGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL, |
|
|
364 $ WORK( KI+1 ), 1, DCMPLX( SCALE ), |
|
|
365 $ VL( 1, KI ), 1 ) |
|
|
366 * |
|
|
367 II = IZAMAX( N, VL( 1, KI ), 1 ) |
|
|
368 REMAX = ONE / CABS1( VL( II, KI ) ) |
|
|
369 CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 ) |
|
|
370 END IF |
|
|
371 * |
|
|
372 * Set back the original diagonal elements of T. |
|
|
373 * |
|
|
374 DO 120 K = KI + 1, N |
|
|
375 T( K, K ) = WORK( K+N ) |
|
|
376 120 CONTINUE |
|
|
377 * |
|
|
378 IS = IS + 1 |
|
|
379 130 CONTINUE |
|
|
380 END IF |
|
|
381 * |
|
|
382 RETURN |
|
|
383 * |
|
|
384 * End of ZTREVC |
|
|
385 * |
|
|
386 END |
