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