1 SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
2 $ LDVR, MM, M, WORK, 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 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
16 $ WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DTREVC computes some or all of the right and/or left eigenvectors of
23 * a real upper quasi-triangular matrix T.
24 * Matrices of this type are produced by the Schur factorization of
25 * a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
26 *
27 * The right eigenvector x and the left eigenvector y of T corresponding
28 * to an eigenvalue w are defined by:
29 *
30 * T*x = w*x, (y**T)*T = w*(y**T)
31 *
32 * where y**T denotes the transpose of y.
33 * The eigenvalues are not input to this routine, but are read directly
34 * from the diagonal blocks 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 orthogonal factor that reduces a matrix
39 * A to Schur form T, then Q*X and Q*Y are the matrices of right and
40 * left eigenvectors of A.
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,
53 * backtransformed by the matrices in VR and/or VL;
54 * = 'S': compute selected right and/or left eigenvectors,
55 * as indicated by the logical array SELECT.
56 *
57 * SELECT (input/output) LOGICAL array, dimension (N)
58 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be
59 * computed.
60 * If w(j) is a real eigenvalue, the corresponding real
61 * eigenvector is computed if SELECT(j) is .TRUE..
62 * If w(j) and w(j+1) are the real and imaginary parts of a
63 * complex eigenvalue, the corresponding complex eigenvector is
64 * computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
65 * on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
66 * .FALSE..
67 * Not referenced if HOWMNY = 'A' or 'B'.
68 *
69 * N (input) INTEGER
70 * The order of the matrix T. N >= 0.
71 *
72 * T (input) DOUBLE PRECISION array, dimension (LDT,N)
73 * The upper quasi-triangular matrix T in Schur canonical form.
74 *
75 * LDT (input) INTEGER
76 * The leading dimension of the array T. LDT >= max(1,N).
77 *
78 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
79 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
80 * contain an N-by-N matrix Q (usually the orthogonal matrix Q
81 * of Schur vectors returned by DHSEQR).
82 * On exit, if SIDE = 'L' or 'B', VL contains:
83 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
84 * if HOWMNY = 'B', the matrix Q*Y;
85 * if HOWMNY = 'S', the left eigenvectors of T specified by
86 * SELECT, stored consecutively in the columns
87 * of VL, in the same order as their
88 * eigenvalues.
89 * A complex eigenvector corresponding to a complex eigenvalue
90 * is stored in two consecutive columns, the first holding the
91 * real part, and the second the imaginary part.
92 * Not referenced if SIDE = 'R'.
93 *
94 * LDVL (input) INTEGER
95 * The leading dimension of the array VL. LDVL >= 1, and if
96 * SIDE = 'L' or 'B', LDVL >= N.
97 *
98 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
99 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
100 * contain an N-by-N matrix Q (usually the orthogonal matrix Q
101 * of Schur vectors returned by DHSEQR).
102 * On exit, if SIDE = 'R' or 'B', VR contains:
103 * if HOWMNY = 'A', the matrix X of right eigenvectors of T;
104 * if HOWMNY = 'B', the matrix Q*X;
105 * if HOWMNY = 'S', the right eigenvectors of T specified by
106 * SELECT, stored consecutively in the columns
107 * of VR, in the same order as their
108 * eigenvalues.
109 * A complex eigenvector corresponding to a complex eigenvalue
110 * is stored in two consecutive columns, the first holding the
111 * real part and the second the imaginary part.
112 * Not referenced if SIDE = 'L'.
113 *
114 * LDVR (input) INTEGER
115 * The leading dimension of the array VR. LDVR >= 1, and if
116 * SIDE = 'R' or 'B', LDVR >= N.
117 *
118 * MM (input) INTEGER
119 * The number of columns in the arrays VL and/or VR. MM >= M.
120 *
121 * M (output) INTEGER
122 * The number of columns in the arrays VL and/or VR actually
123 * used to store the eigenvectors.
124 * If HOWMNY = 'A' or 'B', M is set to N.
125 * Each selected real eigenvector occupies one column and each
126 * selected complex eigenvector occupies two columns.
127 *
128 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
129 *
130 * INFO (output) INTEGER
131 * = 0: successful exit
132 * < 0: if INFO = -i, the i-th argument had an illegal value
133 *
134 * Further Details
135 * ===============
136 *
137 * The algorithm used in this program is basically backward (forward)
138 * substitution, with scaling to make the the code robust against
139 * possible overflow.
140 *
141 * Each eigenvector is normalized so that the element of largest
142 * magnitude has magnitude 1; here the magnitude of a complex number
143 * (x,y) is taken to be |x| + |y|.
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ZERO, ONE
149 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
153 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
154 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
155 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
156 $ XNORM
157 * ..
158 * .. External Functions ..
159 LOGICAL LSAME
160 INTEGER IDAMAX
161 DOUBLE PRECISION DDOT, DLAMCH
162 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
166 * ..
167 * .. Intrinsic Functions ..
168 INTRINSIC ABS, MAX, SQRT
169 * ..
170 * .. Local Arrays ..
171 DOUBLE PRECISION X( 2, 2 )
172 * ..
173 * .. Executable Statements ..
174 *
175 * Decode and test the input parameters
176 *
177 BOTHV = LSAME( SIDE, 'B' )
178 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
179 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
180 *
181 ALLV = LSAME( HOWMNY, 'A' )
182 OVER = LSAME( HOWMNY, 'B' )
183 SOMEV = LSAME( HOWMNY, 'S' )
184 *
185 INFO = 0
186 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
187 INFO = -1
188 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
189 INFO = -2
190 ELSE IF( N.LT.0 ) THEN
191 INFO = -4
192 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
193 INFO = -6
194 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
195 INFO = -8
196 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
197 INFO = -10
198 ELSE
199 *
200 * Set M to the number of columns required to store the selected
201 * eigenvectors, standardize the array SELECT if necessary, and
202 * test MM.
203 *
204 IF( SOMEV ) THEN
205 M = 0
206 PAIR = .FALSE.
207 DO 10 J = 1, N
208 IF( PAIR ) THEN
209 PAIR = .FALSE.
210 SELECT( J ) = .FALSE.
211 ELSE
212 IF( J.LT.N ) THEN
213 IF( T( J+1, J ).EQ.ZERO ) THEN
214 IF( SELECT( J ) )
215 $ M = M + 1
216 ELSE
217 PAIR = .TRUE.
218 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
219 SELECT( J ) = .TRUE.
220 M = M + 2
221 END IF
222 END IF
223 ELSE
224 IF( SELECT( N ) )
225 $ M = M + 1
226 END IF
227 END IF
228 10 CONTINUE
229 ELSE
230 M = N
231 END IF
232 *
233 IF( MM.LT.M ) THEN
234 INFO = -11
235 END IF
236 END IF
237 IF( INFO.NE.0 ) THEN
238 CALL XERBLA( 'DTREVC', -INFO )
239 RETURN
240 END IF
241 *
242 * Quick return if possible.
243 *
244 IF( N.EQ.0 )
245 $ RETURN
246 *
247 * Set the constants to control overflow.
248 *
249 UNFL = DLAMCH( 'Safe minimum' )
250 OVFL = ONE / UNFL
251 CALL DLABAD( UNFL, OVFL )
252 ULP = DLAMCH( 'Precision' )
253 SMLNUM = UNFL*( N / ULP )
254 BIGNUM = ( ONE-ULP ) / SMLNUM
255 *
256 * Compute 1-norm of each column of strictly upper triangular
257 * part of T to control overflow in triangular solver.
258 *
259 WORK( 1 ) = ZERO
260 DO 30 J = 2, N
261 WORK( J ) = ZERO
262 DO 20 I = 1, J - 1
263 WORK( J ) = WORK( J ) + ABS( T( I, J ) )
264 20 CONTINUE
265 30 CONTINUE
266 *
267 * Index IP is used to specify the real or complex eigenvalue:
268 * IP = 0, real eigenvalue,
269 * 1, first of conjugate complex pair: (wr,wi)
270 * -1, second of conjugate complex pair: (wr,wi)
271 *
272 N2 = 2*N
273 *
274 IF( RIGHTV ) THEN
275 *
276 * Compute right eigenvectors.
277 *
278 IP = 0
279 IS = M
280 DO 140 KI = N, 1, -1
281 *
282 IF( IP.EQ.1 )
283 $ GO TO 130
284 IF( KI.EQ.1 )
285 $ GO TO 40
286 IF( T( KI, KI-1 ).EQ.ZERO )
287 $ GO TO 40
288 IP = -1
289 *
290 40 CONTINUE
291 IF( SOMEV ) THEN
292 IF( IP.EQ.0 ) THEN
293 IF( .NOT.SELECT( KI ) )
294 $ GO TO 130
295 ELSE
296 IF( .NOT.SELECT( KI-1 ) )
297 $ GO TO 130
298 END IF
299 END IF
300 *
301 * Compute the KI-th eigenvalue (WR,WI).
302 *
303 WR = T( KI, KI )
304 WI = ZERO
305 IF( IP.NE.0 )
306 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
307 $ SQRT( ABS( T( KI-1, KI ) ) )
308 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
309 *
310 IF( IP.EQ.0 ) THEN
311 *
312 * Real right eigenvector
313 *
314 WORK( KI+N ) = ONE
315 *
316 * Form right-hand side
317 *
318 DO 50 K = 1, KI - 1
319 WORK( K+N ) = -T( K, KI )
320 50 CONTINUE
321 *
322 * Solve the upper quasi-triangular system:
323 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
324 *
325 JNXT = KI - 1
326 DO 60 J = KI - 1, 1, -1
327 IF( J.GT.JNXT )
328 $ GO TO 60
329 J1 = J
330 J2 = J
331 JNXT = J - 1
332 IF( J.GT.1 ) THEN
333 IF( T( J, J-1 ).NE.ZERO ) THEN
334 J1 = J - 1
335 JNXT = J - 2
336 END IF
337 END IF
338 *
339 IF( J1.EQ.J2 ) THEN
340 *
341 * 1-by-1 diagonal block
342 *
343 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
344 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
345 $ ZERO, X, 2, SCALE, XNORM, IERR )
346 *
347 * Scale X(1,1) to avoid overflow when updating
348 * the right-hand side.
349 *
350 IF( XNORM.GT.ONE ) THEN
351 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
352 X( 1, 1 ) = X( 1, 1 ) / XNORM
353 SCALE = SCALE / XNORM
354 END IF
355 END IF
356 *
357 * Scale if necessary
358 *
359 IF( SCALE.NE.ONE )
360 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
361 WORK( J+N ) = X( 1, 1 )
362 *
363 * Update right-hand side
364 *
365 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
366 $ WORK( 1+N ), 1 )
367 *
368 ELSE
369 *
370 * 2-by-2 diagonal block
371 *
372 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
373 $ T( J-1, J-1 ), LDT, ONE, ONE,
374 $ WORK( J-1+N ), N, WR, ZERO, X, 2,
375 $ SCALE, XNORM, IERR )
376 *
377 * Scale X(1,1) and X(2,1) to avoid overflow when
378 * updating the right-hand side.
379 *
380 IF( XNORM.GT.ONE ) THEN
381 BETA = MAX( WORK( J-1 ), WORK( J ) )
382 IF( BETA.GT.BIGNUM / XNORM ) THEN
383 X( 1, 1 ) = X( 1, 1 ) / XNORM
384 X( 2, 1 ) = X( 2, 1 ) / XNORM
385 SCALE = SCALE / XNORM
386 END IF
387 END IF
388 *
389 * Scale if necessary
390 *
391 IF( SCALE.NE.ONE )
392 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
393 WORK( J-1+N ) = X( 1, 1 )
394 WORK( J+N ) = X( 2, 1 )
395 *
396 * Update right-hand side
397 *
398 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
399 $ WORK( 1+N ), 1 )
400 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
401 $ WORK( 1+N ), 1 )
402 END IF
403 60 CONTINUE
404 *
405 * Copy the vector x or Q*x to VR and normalize.
406 *
407 IF( .NOT.OVER ) THEN
408 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
409 *
410 II = IDAMAX( KI, VR( 1, IS ), 1 )
411 REMAX = ONE / ABS( VR( II, IS ) )
412 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
413 *
414 DO 70 K = KI + 1, N
415 VR( K, IS ) = ZERO
416 70 CONTINUE
417 ELSE
418 IF( KI.GT.1 )
419 $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
420 $ WORK( 1+N ), 1, WORK( KI+N ),
421 $ VR( 1, KI ), 1 )
422 *
423 II = IDAMAX( N, VR( 1, KI ), 1 )
424 REMAX = ONE / ABS( VR( II, KI ) )
425 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
426 END IF
427 *
428 ELSE
429 *
430 * Complex right eigenvector.
431 *
432 * Initial solve
433 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
434 * [ (T(KI,KI-1) T(KI,KI) ) ]
435 *
436 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
437 WORK( KI-1+N ) = ONE
438 WORK( KI+N2 ) = WI / T( KI-1, KI )
439 ELSE
440 WORK( KI-1+N ) = -WI / T( KI, KI-1 )
441 WORK( KI+N2 ) = ONE
442 END IF
443 WORK( KI+N ) = ZERO
444 WORK( KI-1+N2 ) = ZERO
445 *
446 * Form right-hand side
447 *
448 DO 80 K = 1, KI - 2
449 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
450 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
451 80 CONTINUE
452 *
453 * Solve upper quasi-triangular system:
454 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
455 *
456 JNXT = KI - 2
457 DO 90 J = KI - 2, 1, -1
458 IF( J.GT.JNXT )
459 $ GO TO 90
460 J1 = J
461 J2 = J
462 JNXT = J - 1
463 IF( J.GT.1 ) THEN
464 IF( T( J, J-1 ).NE.ZERO ) THEN
465 J1 = J - 1
466 JNXT = J - 2
467 END IF
468 END IF
469 *
470 IF( J1.EQ.J2 ) THEN
471 *
472 * 1-by-1 diagonal block
473 *
474 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
475 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
476 $ X, 2, SCALE, XNORM, IERR )
477 *
478 * Scale X(1,1) and X(1,2) to avoid overflow when
479 * updating the right-hand side.
480 *
481 IF( XNORM.GT.ONE ) THEN
482 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
483 X( 1, 1 ) = X( 1, 1 ) / XNORM
484 X( 1, 2 ) = X( 1, 2 ) / XNORM
485 SCALE = SCALE / XNORM
486 END IF
487 END IF
488 *
489 * Scale if necessary
490 *
491 IF( SCALE.NE.ONE ) THEN
492 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
493 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
494 END IF
495 WORK( J+N ) = X( 1, 1 )
496 WORK( J+N2 ) = X( 1, 2 )
497 *
498 * Update the right-hand side
499 *
500 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
501 $ WORK( 1+N ), 1 )
502 CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
503 $ WORK( 1+N2 ), 1 )
504 *
505 ELSE
506 *
507 * 2-by-2 diagonal block
508 *
509 CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
510 $ T( J-1, J-1 ), LDT, ONE, ONE,
511 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
512 $ XNORM, IERR )
513 *
514 * Scale X to avoid overflow when updating
515 * the right-hand side.
516 *
517 IF( XNORM.GT.ONE ) THEN
518 BETA = MAX( WORK( J-1 ), WORK( J ) )
519 IF( BETA.GT.BIGNUM / XNORM ) THEN
520 REC = ONE / XNORM
521 X( 1, 1 ) = X( 1, 1 )*REC
522 X( 1, 2 ) = X( 1, 2 )*REC
523 X( 2, 1 ) = X( 2, 1 )*REC
524 X( 2, 2 ) = X( 2, 2 )*REC
525 SCALE = SCALE*REC
526 END IF
527 END IF
528 *
529 * Scale if necessary
530 *
531 IF( SCALE.NE.ONE ) THEN
532 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
533 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
534 END IF
535 WORK( J-1+N ) = X( 1, 1 )
536 WORK( J+N ) = X( 2, 1 )
537 WORK( J-1+N2 ) = X( 1, 2 )
538 WORK( J+N2 ) = X( 2, 2 )
539 *
540 * Update the right-hand side
541 *
542 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
543 $ WORK( 1+N ), 1 )
544 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
545 $ WORK( 1+N ), 1 )
546 CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
547 $ WORK( 1+N2 ), 1 )
548 CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
549 $ WORK( 1+N2 ), 1 )
550 END IF
551 90 CONTINUE
552 *
553 * Copy the vector x or Q*x to VR and normalize.
554 *
555 IF( .NOT.OVER ) THEN
556 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
557 CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
558 *
559 EMAX = ZERO
560 DO 100 K = 1, KI
561 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
562 $ ABS( VR( K, IS ) ) )
563 100 CONTINUE
564 *
565 REMAX = ONE / EMAX
566 CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
567 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
568 *
569 DO 110 K = KI + 1, N
570 VR( K, IS-1 ) = ZERO
571 VR( K, IS ) = ZERO
572 110 CONTINUE
573 *
574 ELSE
575 *
576 IF( KI.GT.2 ) THEN
577 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
578 $ WORK( 1+N ), 1, WORK( KI-1+N ),
579 $ VR( 1, KI-1 ), 1 )
580 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
581 $ WORK( 1+N2 ), 1, WORK( KI+N2 ),
582 $ VR( 1, KI ), 1 )
583 ELSE
584 CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
585 CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
586 END IF
587 *
588 EMAX = ZERO
589 DO 120 K = 1, N
590 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
591 $ ABS( VR( K, KI ) ) )
592 120 CONTINUE
593 REMAX = ONE / EMAX
594 CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
595 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
596 END IF
597 END IF
598 *
599 IS = IS - 1
600 IF( IP.NE.0 )
601 $ IS = IS - 1
602 130 CONTINUE
603 IF( IP.EQ.1 )
604 $ IP = 0
605 IF( IP.EQ.-1 )
606 $ IP = 1
607 140 CONTINUE
608 END IF
609 *
610 IF( LEFTV ) THEN
611 *
612 * Compute left eigenvectors.
613 *
614 IP = 0
615 IS = 1
616 DO 260 KI = 1, N
617 *
618 IF( IP.EQ.-1 )
619 $ GO TO 250
620 IF( KI.EQ.N )
621 $ GO TO 150
622 IF( T( KI+1, KI ).EQ.ZERO )
623 $ GO TO 150
624 IP = 1
625 *
626 150 CONTINUE
627 IF( SOMEV ) THEN
628 IF( .NOT.SELECT( KI ) )
629 $ GO TO 250
630 END IF
631 *
632 * Compute the KI-th eigenvalue (WR,WI).
633 *
634 WR = T( KI, KI )
635 WI = ZERO
636 IF( IP.NE.0 )
637 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
638 $ SQRT( ABS( T( KI+1, KI ) ) )
639 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
640 *
641 IF( IP.EQ.0 ) THEN
642 *
643 * Real left eigenvector.
644 *
645 WORK( KI+N ) = ONE
646 *
647 * Form right-hand side
648 *
649 DO 160 K = KI + 1, N
650 WORK( K+N ) = -T( KI, K )
651 160 CONTINUE
652 *
653 * Solve the quasi-triangular system:
654 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
655 *
656 VMAX = ONE
657 VCRIT = BIGNUM
658 *
659 JNXT = KI + 1
660 DO 170 J = KI + 1, N
661 IF( J.LT.JNXT )
662 $ GO TO 170
663 J1 = J
664 J2 = J
665 JNXT = J + 1
666 IF( J.LT.N ) THEN
667 IF( T( J+1, J ).NE.ZERO ) THEN
668 J2 = J + 1
669 JNXT = J + 2
670 END IF
671 END IF
672 *
673 IF( J1.EQ.J2 ) THEN
674 *
675 * 1-by-1 diagonal block
676 *
677 * Scale if necessary to avoid overflow when forming
678 * the right-hand side.
679 *
680 IF( WORK( J ).GT.VCRIT ) THEN
681 REC = ONE / VMAX
682 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
683 VMAX = ONE
684 VCRIT = BIGNUM
685 END IF
686 *
687 WORK( J+N ) = WORK( J+N ) -
688 $ DDOT( J-KI-1, T( KI+1, J ), 1,
689 $ WORK( KI+1+N ), 1 )
690 *
691 * Solve (T(J,J)-WR)**T*X = WORK
692 *
693 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
694 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
695 $ ZERO, X, 2, SCALE, XNORM, IERR )
696 *
697 * Scale if necessary
698 *
699 IF( SCALE.NE.ONE )
700 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
701 WORK( J+N ) = X( 1, 1 )
702 VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
703 VCRIT = BIGNUM / VMAX
704 *
705 ELSE
706 *
707 * 2-by-2 diagonal block
708 *
709 * Scale if necessary to avoid overflow when forming
710 * the right-hand side.
711 *
712 BETA = MAX( WORK( J ), WORK( J+1 ) )
713 IF( BETA.GT.VCRIT ) THEN
714 REC = ONE / VMAX
715 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
716 VMAX = ONE
717 VCRIT = BIGNUM
718 END IF
719 *
720 WORK( J+N ) = WORK( J+N ) -
721 $ DDOT( J-KI-1, T( KI+1, J ), 1,
722 $ WORK( KI+1+N ), 1 )
723 *
724 WORK( J+1+N ) = WORK( J+1+N ) -
725 $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
726 $ WORK( KI+1+N ), 1 )
727 *
728 * Solve
729 * [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
730 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
731 *
732 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
733 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
734 $ ZERO, X, 2, SCALE, XNORM, IERR )
735 *
736 * Scale if necessary
737 *
738 IF( SCALE.NE.ONE )
739 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
740 WORK( J+N ) = X( 1, 1 )
741 WORK( J+1+N ) = X( 2, 1 )
742 *
743 VMAX = MAX( ABS( WORK( J+N ) ),
744 $ ABS( WORK( J+1+N ) ), VMAX )
745 VCRIT = BIGNUM / VMAX
746 *
747 END IF
748 170 CONTINUE
749 *
750 * Copy the vector x or Q*x to VL and normalize.
751 *
752 IF( .NOT.OVER ) THEN
753 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
754 *
755 II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
756 REMAX = ONE / ABS( VL( II, IS ) )
757 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
758 *
759 DO 180 K = 1, KI - 1
760 VL( K, IS ) = ZERO
761 180 CONTINUE
762 *
763 ELSE
764 *
765 IF( KI.LT.N )
766 $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
767 $ WORK( KI+1+N ), 1, WORK( KI+N ),
768 $ VL( 1, KI ), 1 )
769 *
770 II = IDAMAX( N, VL( 1, KI ), 1 )
771 REMAX = ONE / ABS( VL( II, KI ) )
772 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
773 *
774 END IF
775 *
776 ELSE
777 *
778 * Complex left eigenvector.
779 *
780 * Initial solve:
781 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
782 * ((T(KI+1,KI) T(KI+1,KI+1)) )
783 *
784 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
785 WORK( KI+N ) = WI / T( KI, KI+1 )
786 WORK( KI+1+N2 ) = ONE
787 ELSE
788 WORK( KI+N ) = ONE
789 WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
790 END IF
791 WORK( KI+1+N ) = ZERO
792 WORK( KI+N2 ) = ZERO
793 *
794 * Form right-hand side
795 *
796 DO 190 K = KI + 2, N
797 WORK( K+N ) = -WORK( KI+N )*T( KI, K )
798 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
799 190 CONTINUE
800 *
801 * Solve complex quasi-triangular system:
802 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
803 *
804 VMAX = ONE
805 VCRIT = BIGNUM
806 *
807 JNXT = KI + 2
808 DO 200 J = KI + 2, N
809 IF( J.LT.JNXT )
810 $ GO TO 200
811 J1 = J
812 J2 = J
813 JNXT = J + 1
814 IF( J.LT.N ) THEN
815 IF( T( J+1, J ).NE.ZERO ) THEN
816 J2 = J + 1
817 JNXT = J + 2
818 END IF
819 END IF
820 *
821 IF( J1.EQ.J2 ) THEN
822 *
823 * 1-by-1 diagonal block
824 *
825 * Scale if necessary to avoid overflow when
826 * forming the right-hand side elements.
827 *
828 IF( WORK( J ).GT.VCRIT ) THEN
829 REC = ONE / VMAX
830 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
831 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
832 VMAX = ONE
833 VCRIT = BIGNUM
834 END IF
835 *
836 WORK( J+N ) = WORK( J+N ) -
837 $ DDOT( J-KI-2, T( KI+2, J ), 1,
838 $ WORK( KI+2+N ), 1 )
839 WORK( J+N2 ) = WORK( J+N2 ) -
840 $ DDOT( J-KI-2, T( KI+2, J ), 1,
841 $ WORK( KI+2+N2 ), 1 )
842 *
843 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
844 *
845 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
846 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
847 $ -WI, X, 2, SCALE, XNORM, IERR )
848 *
849 * Scale if necessary
850 *
851 IF( SCALE.NE.ONE ) THEN
852 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
853 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
854 END IF
855 WORK( J+N ) = X( 1, 1 )
856 WORK( J+N2 ) = X( 1, 2 )
857 VMAX = MAX( ABS( WORK( J+N ) ),
858 $ ABS( WORK( J+N2 ) ), VMAX )
859 VCRIT = BIGNUM / VMAX
860 *
861 ELSE
862 *
863 * 2-by-2 diagonal block
864 *
865 * Scale if necessary to avoid overflow when forming
866 * the right-hand side elements.
867 *
868 BETA = MAX( WORK( J ), WORK( J+1 ) )
869 IF( BETA.GT.VCRIT ) THEN
870 REC = ONE / VMAX
871 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
872 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
873 VMAX = ONE
874 VCRIT = BIGNUM
875 END IF
876 *
877 WORK( J+N ) = WORK( J+N ) -
878 $ DDOT( J-KI-2, T( KI+2, J ), 1,
879 $ WORK( KI+2+N ), 1 )
880 *
881 WORK( J+N2 ) = WORK( J+N2 ) -
882 $ DDOT( J-KI-2, T( KI+2, J ), 1,
883 $ WORK( KI+2+N2 ), 1 )
884 *
885 WORK( J+1+N ) = WORK( J+1+N ) -
886 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
887 $ WORK( KI+2+N ), 1 )
888 *
889 WORK( J+1+N2 ) = WORK( J+1+N2 ) -
890 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
891 $ WORK( KI+2+N2 ), 1 )
892 *
893 * Solve 2-by-2 complex linear equation
894 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
895 * ([T(j+1,j) T(j+1,j+1)] )
896 *
897 CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
898 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
899 $ -WI, X, 2, SCALE, XNORM, IERR )
900 *
901 * Scale if necessary
902 *
903 IF( SCALE.NE.ONE ) THEN
904 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
905 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
906 END IF
907 WORK( J+N ) = X( 1, 1 )
908 WORK( J+N2 ) = X( 1, 2 )
909 WORK( J+1+N ) = X( 2, 1 )
910 WORK( J+1+N2 ) = X( 2, 2 )
911 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
912 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
913 VCRIT = BIGNUM / VMAX
914 *
915 END IF
916 200 CONTINUE
917 *
918 * Copy the vector x or Q*x to VL and normalize.
919 *
920 IF( .NOT.OVER ) THEN
921 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
922 CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
923 $ 1 )
924 *
925 EMAX = ZERO
926 DO 220 K = KI, N
927 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
928 $ ABS( VL( K, IS+1 ) ) )
929 220 CONTINUE
930 REMAX = ONE / EMAX
931 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
932 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
933 *
934 DO 230 K = 1, KI - 1
935 VL( K, IS ) = ZERO
936 VL( K, IS+1 ) = ZERO
937 230 CONTINUE
938 ELSE
939 IF( KI.LT.N-1 ) THEN
940 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
941 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
942 $ VL( 1, KI ), 1 )
943 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
944 $ LDVL, WORK( KI+2+N2 ), 1,
945 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
946 ELSE
947 CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
948 CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
949 END IF
950 *
951 EMAX = ZERO
952 DO 240 K = 1, N
953 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
954 $ ABS( VL( K, KI+1 ) ) )
955 240 CONTINUE
956 REMAX = ONE / EMAX
957 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
958 CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
959 *
960 END IF
961 *
962 END IF
963 *
964 IS = IS + 1
965 IF( IP.NE.0 )
966 $ IS = IS + 1
967 250 CONTINUE
968 IF( IP.EQ.-1 )
969 $ IP = 0
970 IF( IP.EQ.1 )
971 $ IP = -1
972 *
973 260 CONTINUE
974 *
975 END IF
976 *
977 RETURN
978 *
979 * End of DTREVC
980 *
981 END
2 $ LDVR, MM, M, WORK, 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 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
16 $ WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DTREVC computes some or all of the right and/or left eigenvectors of
23 * a real upper quasi-triangular matrix T.
24 * Matrices of this type are produced by the Schur factorization of
25 * a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
26 *
27 * The right eigenvector x and the left eigenvector y of T corresponding
28 * to an eigenvalue w are defined by:
29 *
30 * T*x = w*x, (y**T)*T = w*(y**T)
31 *
32 * where y**T denotes the transpose of y.
33 * The eigenvalues are not input to this routine, but are read directly
34 * from the diagonal blocks 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 orthogonal factor that reduces a matrix
39 * A to Schur form T, then Q*X and Q*Y are the matrices of right and
40 * left eigenvectors of A.
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,
53 * backtransformed by the matrices in VR and/or VL;
54 * = 'S': compute selected right and/or left eigenvectors,
55 * as indicated by the logical array SELECT.
56 *
57 * SELECT (input/output) LOGICAL array, dimension (N)
58 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be
59 * computed.
60 * If w(j) is a real eigenvalue, the corresponding real
61 * eigenvector is computed if SELECT(j) is .TRUE..
62 * If w(j) and w(j+1) are the real and imaginary parts of a
63 * complex eigenvalue, the corresponding complex eigenvector is
64 * computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
65 * on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
66 * .FALSE..
67 * Not referenced if HOWMNY = 'A' or 'B'.
68 *
69 * N (input) INTEGER
70 * The order of the matrix T. N >= 0.
71 *
72 * T (input) DOUBLE PRECISION array, dimension (LDT,N)
73 * The upper quasi-triangular matrix T in Schur canonical form.
74 *
75 * LDT (input) INTEGER
76 * The leading dimension of the array T. LDT >= max(1,N).
77 *
78 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
79 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
80 * contain an N-by-N matrix Q (usually the orthogonal matrix Q
81 * of Schur vectors returned by DHSEQR).
82 * On exit, if SIDE = 'L' or 'B', VL contains:
83 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
84 * if HOWMNY = 'B', the matrix Q*Y;
85 * if HOWMNY = 'S', the left eigenvectors of T specified by
86 * SELECT, stored consecutively in the columns
87 * of VL, in the same order as their
88 * eigenvalues.
89 * A complex eigenvector corresponding to a complex eigenvalue
90 * is stored in two consecutive columns, the first holding the
91 * real part, and the second the imaginary part.
92 * Not referenced if SIDE = 'R'.
93 *
94 * LDVL (input) INTEGER
95 * The leading dimension of the array VL. LDVL >= 1, and if
96 * SIDE = 'L' or 'B', LDVL >= N.
97 *
98 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
99 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
100 * contain an N-by-N matrix Q (usually the orthogonal matrix Q
101 * of Schur vectors returned by DHSEQR).
102 * On exit, if SIDE = 'R' or 'B', VR contains:
103 * if HOWMNY = 'A', the matrix X of right eigenvectors of T;
104 * if HOWMNY = 'B', the matrix Q*X;
105 * if HOWMNY = 'S', the right eigenvectors of T specified by
106 * SELECT, stored consecutively in the columns
107 * of VR, in the same order as their
108 * eigenvalues.
109 * A complex eigenvector corresponding to a complex eigenvalue
110 * is stored in two consecutive columns, the first holding the
111 * real part and the second the imaginary part.
112 * Not referenced if SIDE = 'L'.
113 *
114 * LDVR (input) INTEGER
115 * The leading dimension of the array VR. LDVR >= 1, and if
116 * SIDE = 'R' or 'B', LDVR >= N.
117 *
118 * MM (input) INTEGER
119 * The number of columns in the arrays VL and/or VR. MM >= M.
120 *
121 * M (output) INTEGER
122 * The number of columns in the arrays VL and/or VR actually
123 * used to store the eigenvectors.
124 * If HOWMNY = 'A' or 'B', M is set to N.
125 * Each selected real eigenvector occupies one column and each
126 * selected complex eigenvector occupies two columns.
127 *
128 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
129 *
130 * INFO (output) INTEGER
131 * = 0: successful exit
132 * < 0: if INFO = -i, the i-th argument had an illegal value
133 *
134 * Further Details
135 * ===============
136 *
137 * The algorithm used in this program is basically backward (forward)
138 * substitution, with scaling to make the the code robust against
139 * possible overflow.
140 *
141 * Each eigenvector is normalized so that the element of largest
142 * magnitude has magnitude 1; here the magnitude of a complex number
143 * (x,y) is taken to be |x| + |y|.
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148 DOUBLE PRECISION ZERO, ONE
149 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
150 * ..
151 * .. Local Scalars ..
152 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
153 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
154 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
155 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
156 $ XNORM
157 * ..
158 * .. External Functions ..
159 LOGICAL LSAME
160 INTEGER IDAMAX
161 DOUBLE PRECISION DDOT, DLAMCH
162 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
166 * ..
167 * .. Intrinsic Functions ..
168 INTRINSIC ABS, MAX, SQRT
169 * ..
170 * .. Local Arrays ..
171 DOUBLE PRECISION X( 2, 2 )
172 * ..
173 * .. Executable Statements ..
174 *
175 * Decode and test the input parameters
176 *
177 BOTHV = LSAME( SIDE, 'B' )
178 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
179 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
180 *
181 ALLV = LSAME( HOWMNY, 'A' )
182 OVER = LSAME( HOWMNY, 'B' )
183 SOMEV = LSAME( HOWMNY, 'S' )
184 *
185 INFO = 0
186 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
187 INFO = -1
188 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
189 INFO = -2
190 ELSE IF( N.LT.0 ) THEN
191 INFO = -4
192 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
193 INFO = -6
194 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
195 INFO = -8
196 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
197 INFO = -10
198 ELSE
199 *
200 * Set M to the number of columns required to store the selected
201 * eigenvectors, standardize the array SELECT if necessary, and
202 * test MM.
203 *
204 IF( SOMEV ) THEN
205 M = 0
206 PAIR = .FALSE.
207 DO 10 J = 1, N
208 IF( PAIR ) THEN
209 PAIR = .FALSE.
210 SELECT( J ) = .FALSE.
211 ELSE
212 IF( J.LT.N ) THEN
213 IF( T( J+1, J ).EQ.ZERO ) THEN
214 IF( SELECT( J ) )
215 $ M = M + 1
216 ELSE
217 PAIR = .TRUE.
218 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
219 SELECT( J ) = .TRUE.
220 M = M + 2
221 END IF
222 END IF
223 ELSE
224 IF( SELECT( N ) )
225 $ M = M + 1
226 END IF
227 END IF
228 10 CONTINUE
229 ELSE
230 M = N
231 END IF
232 *
233 IF( MM.LT.M ) THEN
234 INFO = -11
235 END IF
236 END IF
237 IF( INFO.NE.0 ) THEN
238 CALL XERBLA( 'DTREVC', -INFO )
239 RETURN
240 END IF
241 *
242 * Quick return if possible.
243 *
244 IF( N.EQ.0 )
245 $ RETURN
246 *
247 * Set the constants to control overflow.
248 *
249 UNFL = DLAMCH( 'Safe minimum' )
250 OVFL = ONE / UNFL
251 CALL DLABAD( UNFL, OVFL )
252 ULP = DLAMCH( 'Precision' )
253 SMLNUM = UNFL*( N / ULP )
254 BIGNUM = ( ONE-ULP ) / SMLNUM
255 *
256 * Compute 1-norm of each column of strictly upper triangular
257 * part of T to control overflow in triangular solver.
258 *
259 WORK( 1 ) = ZERO
260 DO 30 J = 2, N
261 WORK( J ) = ZERO
262 DO 20 I = 1, J - 1
263 WORK( J ) = WORK( J ) + ABS( T( I, J ) )
264 20 CONTINUE
265 30 CONTINUE
266 *
267 * Index IP is used to specify the real or complex eigenvalue:
268 * IP = 0, real eigenvalue,
269 * 1, first of conjugate complex pair: (wr,wi)
270 * -1, second of conjugate complex pair: (wr,wi)
271 *
272 N2 = 2*N
273 *
274 IF( RIGHTV ) THEN
275 *
276 * Compute right eigenvectors.
277 *
278 IP = 0
279 IS = M
280 DO 140 KI = N, 1, -1
281 *
282 IF( IP.EQ.1 )
283 $ GO TO 130
284 IF( KI.EQ.1 )
285 $ GO TO 40
286 IF( T( KI, KI-1 ).EQ.ZERO )
287 $ GO TO 40
288 IP = -1
289 *
290 40 CONTINUE
291 IF( SOMEV ) THEN
292 IF( IP.EQ.0 ) THEN
293 IF( .NOT.SELECT( KI ) )
294 $ GO TO 130
295 ELSE
296 IF( .NOT.SELECT( KI-1 ) )
297 $ GO TO 130
298 END IF
299 END IF
300 *
301 * Compute the KI-th eigenvalue (WR,WI).
302 *
303 WR = T( KI, KI )
304 WI = ZERO
305 IF( IP.NE.0 )
306 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
307 $ SQRT( ABS( T( KI-1, KI ) ) )
308 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
309 *
310 IF( IP.EQ.0 ) THEN
311 *
312 * Real right eigenvector
313 *
314 WORK( KI+N ) = ONE
315 *
316 * Form right-hand side
317 *
318 DO 50 K = 1, KI - 1
319 WORK( K+N ) = -T( K, KI )
320 50 CONTINUE
321 *
322 * Solve the upper quasi-triangular system:
323 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
324 *
325 JNXT = KI - 1
326 DO 60 J = KI - 1, 1, -1
327 IF( J.GT.JNXT )
328 $ GO TO 60
329 J1 = J
330 J2 = J
331 JNXT = J - 1
332 IF( J.GT.1 ) THEN
333 IF( T( J, J-1 ).NE.ZERO ) THEN
334 J1 = J - 1
335 JNXT = J - 2
336 END IF
337 END IF
338 *
339 IF( J1.EQ.J2 ) THEN
340 *
341 * 1-by-1 diagonal block
342 *
343 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
344 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
345 $ ZERO, X, 2, SCALE, XNORM, IERR )
346 *
347 * Scale X(1,1) to avoid overflow when updating
348 * the right-hand side.
349 *
350 IF( XNORM.GT.ONE ) THEN
351 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
352 X( 1, 1 ) = X( 1, 1 ) / XNORM
353 SCALE = SCALE / XNORM
354 END IF
355 END IF
356 *
357 * Scale if necessary
358 *
359 IF( SCALE.NE.ONE )
360 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
361 WORK( J+N ) = X( 1, 1 )
362 *
363 * Update right-hand side
364 *
365 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
366 $ WORK( 1+N ), 1 )
367 *
368 ELSE
369 *
370 * 2-by-2 diagonal block
371 *
372 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
373 $ T( J-1, J-1 ), LDT, ONE, ONE,
374 $ WORK( J-1+N ), N, WR, ZERO, X, 2,
375 $ SCALE, XNORM, IERR )
376 *
377 * Scale X(1,1) and X(2,1) to avoid overflow when
378 * updating the right-hand side.
379 *
380 IF( XNORM.GT.ONE ) THEN
381 BETA = MAX( WORK( J-1 ), WORK( J ) )
382 IF( BETA.GT.BIGNUM / XNORM ) THEN
383 X( 1, 1 ) = X( 1, 1 ) / XNORM
384 X( 2, 1 ) = X( 2, 1 ) / XNORM
385 SCALE = SCALE / XNORM
386 END IF
387 END IF
388 *
389 * Scale if necessary
390 *
391 IF( SCALE.NE.ONE )
392 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
393 WORK( J-1+N ) = X( 1, 1 )
394 WORK( J+N ) = X( 2, 1 )
395 *
396 * Update right-hand side
397 *
398 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
399 $ WORK( 1+N ), 1 )
400 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
401 $ WORK( 1+N ), 1 )
402 END IF
403 60 CONTINUE
404 *
405 * Copy the vector x or Q*x to VR and normalize.
406 *
407 IF( .NOT.OVER ) THEN
408 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
409 *
410 II = IDAMAX( KI, VR( 1, IS ), 1 )
411 REMAX = ONE / ABS( VR( II, IS ) )
412 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
413 *
414 DO 70 K = KI + 1, N
415 VR( K, IS ) = ZERO
416 70 CONTINUE
417 ELSE
418 IF( KI.GT.1 )
419 $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
420 $ WORK( 1+N ), 1, WORK( KI+N ),
421 $ VR( 1, KI ), 1 )
422 *
423 II = IDAMAX( N, VR( 1, KI ), 1 )
424 REMAX = ONE / ABS( VR( II, KI ) )
425 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
426 END IF
427 *
428 ELSE
429 *
430 * Complex right eigenvector.
431 *
432 * Initial solve
433 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
434 * [ (T(KI,KI-1) T(KI,KI) ) ]
435 *
436 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
437 WORK( KI-1+N ) = ONE
438 WORK( KI+N2 ) = WI / T( KI-1, KI )
439 ELSE
440 WORK( KI-1+N ) = -WI / T( KI, KI-1 )
441 WORK( KI+N2 ) = ONE
442 END IF
443 WORK( KI+N ) = ZERO
444 WORK( KI-1+N2 ) = ZERO
445 *
446 * Form right-hand side
447 *
448 DO 80 K = 1, KI - 2
449 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
450 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
451 80 CONTINUE
452 *
453 * Solve upper quasi-triangular system:
454 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
455 *
456 JNXT = KI - 2
457 DO 90 J = KI - 2, 1, -1
458 IF( J.GT.JNXT )
459 $ GO TO 90
460 J1 = J
461 J2 = J
462 JNXT = J - 1
463 IF( J.GT.1 ) THEN
464 IF( T( J, J-1 ).NE.ZERO ) THEN
465 J1 = J - 1
466 JNXT = J - 2
467 END IF
468 END IF
469 *
470 IF( J1.EQ.J2 ) THEN
471 *
472 * 1-by-1 diagonal block
473 *
474 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
475 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
476 $ X, 2, SCALE, XNORM, IERR )
477 *
478 * Scale X(1,1) and X(1,2) to avoid overflow when
479 * updating the right-hand side.
480 *
481 IF( XNORM.GT.ONE ) THEN
482 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
483 X( 1, 1 ) = X( 1, 1 ) / XNORM
484 X( 1, 2 ) = X( 1, 2 ) / XNORM
485 SCALE = SCALE / XNORM
486 END IF
487 END IF
488 *
489 * Scale if necessary
490 *
491 IF( SCALE.NE.ONE ) THEN
492 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
493 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
494 END IF
495 WORK( J+N ) = X( 1, 1 )
496 WORK( J+N2 ) = X( 1, 2 )
497 *
498 * Update the right-hand side
499 *
500 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
501 $ WORK( 1+N ), 1 )
502 CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
503 $ WORK( 1+N2 ), 1 )
504 *
505 ELSE
506 *
507 * 2-by-2 diagonal block
508 *
509 CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
510 $ T( J-1, J-1 ), LDT, ONE, ONE,
511 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
512 $ XNORM, IERR )
513 *
514 * Scale X to avoid overflow when updating
515 * the right-hand side.
516 *
517 IF( XNORM.GT.ONE ) THEN
518 BETA = MAX( WORK( J-1 ), WORK( J ) )
519 IF( BETA.GT.BIGNUM / XNORM ) THEN
520 REC = ONE / XNORM
521 X( 1, 1 ) = X( 1, 1 )*REC
522 X( 1, 2 ) = X( 1, 2 )*REC
523 X( 2, 1 ) = X( 2, 1 )*REC
524 X( 2, 2 ) = X( 2, 2 )*REC
525 SCALE = SCALE*REC
526 END IF
527 END IF
528 *
529 * Scale if necessary
530 *
531 IF( SCALE.NE.ONE ) THEN
532 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
533 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
534 END IF
535 WORK( J-1+N ) = X( 1, 1 )
536 WORK( J+N ) = X( 2, 1 )
537 WORK( J-1+N2 ) = X( 1, 2 )
538 WORK( J+N2 ) = X( 2, 2 )
539 *
540 * Update the right-hand side
541 *
542 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
543 $ WORK( 1+N ), 1 )
544 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
545 $ WORK( 1+N ), 1 )
546 CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
547 $ WORK( 1+N2 ), 1 )
548 CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
549 $ WORK( 1+N2 ), 1 )
550 END IF
551 90 CONTINUE
552 *
553 * Copy the vector x or Q*x to VR and normalize.
554 *
555 IF( .NOT.OVER ) THEN
556 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
557 CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
558 *
559 EMAX = ZERO
560 DO 100 K = 1, KI
561 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
562 $ ABS( VR( K, IS ) ) )
563 100 CONTINUE
564 *
565 REMAX = ONE / EMAX
566 CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
567 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
568 *
569 DO 110 K = KI + 1, N
570 VR( K, IS-1 ) = ZERO
571 VR( K, IS ) = ZERO
572 110 CONTINUE
573 *
574 ELSE
575 *
576 IF( KI.GT.2 ) THEN
577 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
578 $ WORK( 1+N ), 1, WORK( KI-1+N ),
579 $ VR( 1, KI-1 ), 1 )
580 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
581 $ WORK( 1+N2 ), 1, WORK( KI+N2 ),
582 $ VR( 1, KI ), 1 )
583 ELSE
584 CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
585 CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
586 END IF
587 *
588 EMAX = ZERO
589 DO 120 K = 1, N
590 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
591 $ ABS( VR( K, KI ) ) )
592 120 CONTINUE
593 REMAX = ONE / EMAX
594 CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
595 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
596 END IF
597 END IF
598 *
599 IS = IS - 1
600 IF( IP.NE.0 )
601 $ IS = IS - 1
602 130 CONTINUE
603 IF( IP.EQ.1 )
604 $ IP = 0
605 IF( IP.EQ.-1 )
606 $ IP = 1
607 140 CONTINUE
608 END IF
609 *
610 IF( LEFTV ) THEN
611 *
612 * Compute left eigenvectors.
613 *
614 IP = 0
615 IS = 1
616 DO 260 KI = 1, N
617 *
618 IF( IP.EQ.-1 )
619 $ GO TO 250
620 IF( KI.EQ.N )
621 $ GO TO 150
622 IF( T( KI+1, KI ).EQ.ZERO )
623 $ GO TO 150
624 IP = 1
625 *
626 150 CONTINUE
627 IF( SOMEV ) THEN
628 IF( .NOT.SELECT( KI ) )
629 $ GO TO 250
630 END IF
631 *
632 * Compute the KI-th eigenvalue (WR,WI).
633 *
634 WR = T( KI, KI )
635 WI = ZERO
636 IF( IP.NE.0 )
637 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
638 $ SQRT( ABS( T( KI+1, KI ) ) )
639 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
640 *
641 IF( IP.EQ.0 ) THEN
642 *
643 * Real left eigenvector.
644 *
645 WORK( KI+N ) = ONE
646 *
647 * Form right-hand side
648 *
649 DO 160 K = KI + 1, N
650 WORK( K+N ) = -T( KI, K )
651 160 CONTINUE
652 *
653 * Solve the quasi-triangular system:
654 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
655 *
656 VMAX = ONE
657 VCRIT = BIGNUM
658 *
659 JNXT = KI + 1
660 DO 170 J = KI + 1, N
661 IF( J.LT.JNXT )
662 $ GO TO 170
663 J1 = J
664 J2 = J
665 JNXT = J + 1
666 IF( J.LT.N ) THEN
667 IF( T( J+1, J ).NE.ZERO ) THEN
668 J2 = J + 1
669 JNXT = J + 2
670 END IF
671 END IF
672 *
673 IF( J1.EQ.J2 ) THEN
674 *
675 * 1-by-1 diagonal block
676 *
677 * Scale if necessary to avoid overflow when forming
678 * the right-hand side.
679 *
680 IF( WORK( J ).GT.VCRIT ) THEN
681 REC = ONE / VMAX
682 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
683 VMAX = ONE
684 VCRIT = BIGNUM
685 END IF
686 *
687 WORK( J+N ) = WORK( J+N ) -
688 $ DDOT( J-KI-1, T( KI+1, J ), 1,
689 $ WORK( KI+1+N ), 1 )
690 *
691 * Solve (T(J,J)-WR)**T*X = WORK
692 *
693 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
694 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
695 $ ZERO, X, 2, SCALE, XNORM, IERR )
696 *
697 * Scale if necessary
698 *
699 IF( SCALE.NE.ONE )
700 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
701 WORK( J+N ) = X( 1, 1 )
702 VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
703 VCRIT = BIGNUM / VMAX
704 *
705 ELSE
706 *
707 * 2-by-2 diagonal block
708 *
709 * Scale if necessary to avoid overflow when forming
710 * the right-hand side.
711 *
712 BETA = MAX( WORK( J ), WORK( J+1 ) )
713 IF( BETA.GT.VCRIT ) THEN
714 REC = ONE / VMAX
715 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
716 VMAX = ONE
717 VCRIT = BIGNUM
718 END IF
719 *
720 WORK( J+N ) = WORK( J+N ) -
721 $ DDOT( J-KI-1, T( KI+1, J ), 1,
722 $ WORK( KI+1+N ), 1 )
723 *
724 WORK( J+1+N ) = WORK( J+1+N ) -
725 $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
726 $ WORK( KI+1+N ), 1 )
727 *
728 * Solve
729 * [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
730 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
731 *
732 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
733 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
734 $ ZERO, X, 2, SCALE, XNORM, IERR )
735 *
736 * Scale if necessary
737 *
738 IF( SCALE.NE.ONE )
739 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
740 WORK( J+N ) = X( 1, 1 )
741 WORK( J+1+N ) = X( 2, 1 )
742 *
743 VMAX = MAX( ABS( WORK( J+N ) ),
744 $ ABS( WORK( J+1+N ) ), VMAX )
745 VCRIT = BIGNUM / VMAX
746 *
747 END IF
748 170 CONTINUE
749 *
750 * Copy the vector x or Q*x to VL and normalize.
751 *
752 IF( .NOT.OVER ) THEN
753 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
754 *
755 II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
756 REMAX = ONE / ABS( VL( II, IS ) )
757 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
758 *
759 DO 180 K = 1, KI - 1
760 VL( K, IS ) = ZERO
761 180 CONTINUE
762 *
763 ELSE
764 *
765 IF( KI.LT.N )
766 $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
767 $ WORK( KI+1+N ), 1, WORK( KI+N ),
768 $ VL( 1, KI ), 1 )
769 *
770 II = IDAMAX( N, VL( 1, KI ), 1 )
771 REMAX = ONE / ABS( VL( II, KI ) )
772 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
773 *
774 END IF
775 *
776 ELSE
777 *
778 * Complex left eigenvector.
779 *
780 * Initial solve:
781 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
782 * ((T(KI+1,KI) T(KI+1,KI+1)) )
783 *
784 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
785 WORK( KI+N ) = WI / T( KI, KI+1 )
786 WORK( KI+1+N2 ) = ONE
787 ELSE
788 WORK( KI+N ) = ONE
789 WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
790 END IF
791 WORK( KI+1+N ) = ZERO
792 WORK( KI+N2 ) = ZERO
793 *
794 * Form right-hand side
795 *
796 DO 190 K = KI + 2, N
797 WORK( K+N ) = -WORK( KI+N )*T( KI, K )
798 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
799 190 CONTINUE
800 *
801 * Solve complex quasi-triangular system:
802 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
803 *
804 VMAX = ONE
805 VCRIT = BIGNUM
806 *
807 JNXT = KI + 2
808 DO 200 J = KI + 2, N
809 IF( J.LT.JNXT )
810 $ GO TO 200
811 J1 = J
812 J2 = J
813 JNXT = J + 1
814 IF( J.LT.N ) THEN
815 IF( T( J+1, J ).NE.ZERO ) THEN
816 J2 = J + 1
817 JNXT = J + 2
818 END IF
819 END IF
820 *
821 IF( J1.EQ.J2 ) THEN
822 *
823 * 1-by-1 diagonal block
824 *
825 * Scale if necessary to avoid overflow when
826 * forming the right-hand side elements.
827 *
828 IF( WORK( J ).GT.VCRIT ) THEN
829 REC = ONE / VMAX
830 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
831 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
832 VMAX = ONE
833 VCRIT = BIGNUM
834 END IF
835 *
836 WORK( J+N ) = WORK( J+N ) -
837 $ DDOT( J-KI-2, T( KI+2, J ), 1,
838 $ WORK( KI+2+N ), 1 )
839 WORK( J+N2 ) = WORK( J+N2 ) -
840 $ DDOT( J-KI-2, T( KI+2, J ), 1,
841 $ WORK( KI+2+N2 ), 1 )
842 *
843 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
844 *
845 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
846 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
847 $ -WI, X, 2, SCALE, XNORM, IERR )
848 *
849 * Scale if necessary
850 *
851 IF( SCALE.NE.ONE ) THEN
852 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
853 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
854 END IF
855 WORK( J+N ) = X( 1, 1 )
856 WORK( J+N2 ) = X( 1, 2 )
857 VMAX = MAX( ABS( WORK( J+N ) ),
858 $ ABS( WORK( J+N2 ) ), VMAX )
859 VCRIT = BIGNUM / VMAX
860 *
861 ELSE
862 *
863 * 2-by-2 diagonal block
864 *
865 * Scale if necessary to avoid overflow when forming
866 * the right-hand side elements.
867 *
868 BETA = MAX( WORK( J ), WORK( J+1 ) )
869 IF( BETA.GT.VCRIT ) THEN
870 REC = ONE / VMAX
871 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
872 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
873 VMAX = ONE
874 VCRIT = BIGNUM
875 END IF
876 *
877 WORK( J+N ) = WORK( J+N ) -
878 $ DDOT( J-KI-2, T( KI+2, J ), 1,
879 $ WORK( KI+2+N ), 1 )
880 *
881 WORK( J+N2 ) = WORK( J+N2 ) -
882 $ DDOT( J-KI-2, T( KI+2, J ), 1,
883 $ WORK( KI+2+N2 ), 1 )
884 *
885 WORK( J+1+N ) = WORK( J+1+N ) -
886 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
887 $ WORK( KI+2+N ), 1 )
888 *
889 WORK( J+1+N2 ) = WORK( J+1+N2 ) -
890 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
891 $ WORK( KI+2+N2 ), 1 )
892 *
893 * Solve 2-by-2 complex linear equation
894 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
895 * ([T(j+1,j) T(j+1,j+1)] )
896 *
897 CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
898 $ LDT, ONE, ONE, WORK( J+N ), N, WR,
899 $ -WI, X, 2, SCALE, XNORM, IERR )
900 *
901 * Scale if necessary
902 *
903 IF( SCALE.NE.ONE ) THEN
904 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
905 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
906 END IF
907 WORK( J+N ) = X( 1, 1 )
908 WORK( J+N2 ) = X( 1, 2 )
909 WORK( J+1+N ) = X( 2, 1 )
910 WORK( J+1+N2 ) = X( 2, 2 )
911 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
912 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
913 VCRIT = BIGNUM / VMAX
914 *
915 END IF
916 200 CONTINUE
917 *
918 * Copy the vector x or Q*x to VL and normalize.
919 *
920 IF( .NOT.OVER ) THEN
921 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
922 CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
923 $ 1 )
924 *
925 EMAX = ZERO
926 DO 220 K = KI, N
927 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
928 $ ABS( VL( K, IS+1 ) ) )
929 220 CONTINUE
930 REMAX = ONE / EMAX
931 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
932 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
933 *
934 DO 230 K = 1, KI - 1
935 VL( K, IS ) = ZERO
936 VL( K, IS+1 ) = ZERO
937 230 CONTINUE
938 ELSE
939 IF( KI.LT.N-1 ) THEN
940 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
941 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
942 $ VL( 1, KI ), 1 )
943 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
944 $ LDVL, WORK( KI+2+N2 ), 1,
945 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
946 ELSE
947 CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
948 CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
949 END IF
950 *
951 EMAX = ZERO
952 DO 240 K = 1, N
953 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
954 $ ABS( VL( K, KI+1 ) ) )
955 240 CONTINUE
956 REMAX = ONE / EMAX
957 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
958 CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
959 *
960 END IF
961 *
962 END IF
963 *
964 IS = IS + 1
965 IF( IP.NE.0 )
966 $ IS = IS + 1
967 250 CONTINUE
968 IF( IP.EQ.-1 )
969 $ IP = 0
970 IF( IP.EQ.1 )
971 $ IP = -1
972 *
973 260 CONTINUE
974 *
975 END IF
976 *
977 RETURN
978 *
979 * End of DTREVC
980 *
981 END