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