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