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: *     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:       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
015:      $                   WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  STREVC 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 SHSEQR.
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) REAL 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) REAL 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 SHSEQR).
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) REAL 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 SHSEQR).
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) REAL 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:       REAL               ZERO, ONE
148:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+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:       REAL               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            ISAMAX
160:       REAL               SDOT, SLAMCH
161:       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
162: *     ..
163: *     .. External Subroutines ..
164:       EXTERNAL           SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL,
165:      $                   XERBLA
166: *     ..
167: *     .. Intrinsic Functions ..
168:       INTRINSIC          ABS, MAX, SQRT
169: *     ..
170: *     .. Local Arrays ..
171:       REAL               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( 'STREVC', -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 = SLAMCH( 'Safe minimum' )
250:       OVFL = ONE / UNFL
251:       CALL SLABAD( UNFL, OVFL )
252:       ULP = SLAMCH( '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 SLALN2( .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 SSCAL( KI, SCALE, WORK( 1+N ), 1 )
361:                      WORK( J+N ) = X( 1, 1 )
362: *
363: *                    Update right-hand side
364: *
365:                      CALL SAXPY( 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 SLALN2( .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 SSCAL( 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 SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
399:      $                           WORK( 1+N ), 1 )
400:                      CALL SAXPY( 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 SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
409: *
410:                   II = ISAMAX( KI, VR( 1, IS ), 1 )
411:                   REMAX = ONE / ABS( VR( II, IS ) )
412:                   CALL SSCAL( 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 SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
420:      $                           WORK( 1+N ), 1, WORK( KI+N ),
421:      $                           VR( 1, KI ), 1 )
422: *
423:                   II = ISAMAX( N, VR( 1, KI ), 1 )
424:                   REMAX = ONE / ABS( VR( II, KI ) )
425:                   CALL SSCAL( 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 SLALN2( .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 SSCAL( KI, SCALE, WORK( 1+N ), 1 )
493:                         CALL SSCAL( 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 SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
501:      $                           WORK( 1+N ), 1 )
502:                      CALL SAXPY( 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 SLALN2( .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 SSCAL( KI, SCALE, WORK( 1+N ), 1 )
533:                         CALL SSCAL( 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 SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
543:      $                           WORK( 1+N ), 1 )
544:                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
545:      $                           WORK( 1+N ), 1 )
546:                      CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
547:      $                           WORK( 1+N2 ), 1 )
548:                      CALL SAXPY( 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 SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
557:                   CALL SCOPY( 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 SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
567:                   CALL SSCAL( 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 SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
578:      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
579:      $                           VR( 1, KI-1 ), 1 )
580:                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
581:      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
582:      $                           VR( 1, KI ), 1 )
583:                   ELSE
584:                      CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
585:                      CALL SSCAL( 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 SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
595:                   CALL SSCAL( 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)'*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 SSCAL( 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:      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
689:      $                             WORK( KI+1+N ), 1 )
690: *
691: *                    Solve (T(J,J)-WR)'*X = WORK
692: *
693:                      CALL SLALN2( .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 SSCAL( 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 SSCAL( 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:      $                             SDOT( 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:      $                               SDOT( 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)     ]'* X = SCALE*( WORK1 )
730: *                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
731: *
732:                      CALL SLALN2( .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 SSCAL( 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 SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
754: *
755:                   II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
756:                   REMAX = ONE / ABS( VL( II, IS ) )
757:                   CALL SSCAL( 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 SGEMV( '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 = ISAMAX( N, VL( 1, KI ), 1 )
771:                   REMAX = ONE / ABS( VL( II, KI ) )
772:                   CALL SSCAL( 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) )' - (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 SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
831:                         CALL SSCAL( 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:      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
838:      $                             WORK( KI+2+N ), 1 )
839:                      WORK( J+N2 ) = WORK( J+N2 ) -
840:      $                              SDOT( 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 SLALN2( .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 SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
853:                         CALL SSCAL( 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 SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
872:                         CALL SSCAL( 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:      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
879:      $                             WORK( KI+2+N ), 1 )
880: *
881:                      WORK( J+N2 ) = WORK( J+N2 ) -
882:      $                              SDOT( 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:      $                               SDOT( 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:      $                                SDOT( 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)  ]'-(wr-i*wi)*I)*X = SCALE*B
895: *                      ([T(j+1,j) T(j+1,j+1)]             )
896: *
897:                      CALL SLALN2( .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 SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
905:                         CALL SSCAL( 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 SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
922:                   CALL SCOPY( 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 SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
932:                   CALL SSCAL( 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 SGEMV( '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 SGEMV( '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 SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
948:                      CALL SSCAL( 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 SSCAL( N, REMAX, VL( 1, KI ), 1 )
958:                   CALL SSCAL( 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 STREVC
980: *
981:       END
982: