LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ strevc()

subroutine strevc ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
real, dimension( ldt, * )  T,
integer  LDT,
real, dimension( ldvl, * )  VL,
integer  LDVL,
real, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
real, dimension( * )  WORK,
integer  INFO 
)

STREVC

Download STREVC + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STREVC computes some or all of the right and/or left eigenvectors of
 a real upper quasi-triangular matrix T.
 Matrices of this type are produced by the Schur factorization of
 a real general matrix:  A = Q*T*Q**T, as computed by SHSEQR.

 The right eigenvector x and the left eigenvector y of T corresponding
 to an eigenvalue w are defined by:

    T*x = w*x,     (y**H)*T = w*(y**H)

 where y**H denotes the conjugate transpose of y.
 The eigenvalues are not input to this routine, but are read directly
 from the diagonal blocks of T.

 This routine returns the matrices X and/or Y of right and left
 eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
 input matrix.  If Q is the orthogonal factor that reduces a matrix
 A to Schur form T, then Q*X and Q*Y are the matrices of right and
 left eigenvectors of A.
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'R':  compute right eigenvectors only;
          = 'L':  compute left eigenvectors only;
          = 'B':  compute both right and left eigenvectors.
[in]HOWMNY
          HOWMNY is CHARACTER*1
          = 'A':  compute all right and/or left eigenvectors;
          = 'B':  compute all right and/or left eigenvectors,
                  backtransformed by the matrices in VR and/or VL;
          = 'S':  compute selected right and/or left eigenvectors,
                  as indicated by the logical array SELECT.
[in,out]SELECT
          SELECT is LOGICAL array, dimension (N)
          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
          computed.
          If w(j) is a real eigenvalue, the corresponding real
          eigenvector is computed if SELECT(j) is .TRUE..
          If w(j) and w(j+1) are the real and imaginary parts of a
          complex eigenvalue, the corresponding complex eigenvector is
          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
          .FALSE..
          Not referenced if HOWMNY = 'A' or 'B'.
[in]N
          N is INTEGER
          The order of the matrix T. N >= 0.
[in]T
          T is REAL array, dimension (LDT,N)
          The upper quasi-triangular matrix T in Schur canonical form.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in,out]VL
          VL is REAL array, dimension (LDVL,MM)
          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
          contain an N-by-N matrix Q (usually the orthogonal matrix Q
          of Schur vectors returned by SHSEQR).
          On exit, if SIDE = 'L' or 'B', VL contains:
          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
          if HOWMNY = 'B', the matrix Q*Y;
          if HOWMNY = 'S', the left eigenvectors of T specified by
                           SELECT, stored consecutively in the columns
                           of VL, in the same order as their
                           eigenvalues.
          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part, and the second the imaginary part.
          Not referenced if SIDE = 'R'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.  LDVL >= 1, and if
          SIDE = 'L' or 'B', LDVL >= N.
[in,out]VR
          VR is REAL array, dimension (LDVR,MM)
          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
          contain an N-by-N matrix Q (usually the orthogonal matrix Q
          of Schur vectors returned by SHSEQR).
          On exit, if SIDE = 'R' or 'B', VR contains:
          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
          if HOWMNY = 'B', the matrix Q*X;
          if HOWMNY = 'S', the right eigenvectors of T specified by
                           SELECT, stored consecutively in the columns
                           of VR, in the same order as their
                           eigenvalues.
          A complex eigenvector corresponding to a complex eigenvalue
          is stored in two consecutive columns, the first holding the
          real part and the second the imaginary part.
          Not referenced if SIDE = 'L'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1, and if
          SIDE = 'R' or 'B', LDVR >= N.
[in]MM
          MM is INTEGER
          The number of columns in the arrays VL and/or VR. MM >= M.
[out]M
          M is INTEGER
          The number of columns in the arrays VL and/or VR actually
          used to store the eigenvectors.
          If HOWMNY = 'A' or 'B', M is set to N.
          Each selected real eigenvector occupies one column and each
          selected complex eigenvector occupies two columns.
[out]WORK
          WORK is REAL array, dimension (3*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  The algorithm used in this program is basically backward (forward)
  substitution, with scaling to make the the code robust against
  possible overflow.

  Each eigenvector is normalized so that the element of largest
  magnitude has magnitude 1; here the magnitude of a complex number
  (x,y) is taken to be |x| + |y|.

Definition at line 220 of file strevc.f.

222 *
223 * -- LAPACK computational routine --
224 * -- LAPACK is a software package provided by Univ. of Tennessee, --
225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226 *
227 * .. Scalar Arguments ..
228  CHARACTER HOWMNY, SIDE
229  INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
230 * ..
231 * .. Array Arguments ..
232  LOGICAL SELECT( * )
233  REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
234  $ WORK( * )
235 * ..
236 *
237 * =====================================================================
238 *
239 * .. Parameters ..
240  REAL ZERO, ONE
241  parameter( zero = 0.0e+0, one = 1.0e+0 )
242 * ..
243 * .. Local Scalars ..
244  LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
245  INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
246  REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
247  $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
248  $ XNORM
249 * ..
250 * .. External Functions ..
251  LOGICAL LSAME
252  INTEGER ISAMAX
253  REAL SDOT, SLAMCH
254  EXTERNAL lsame, isamax, sdot, slamch
255 * ..
256 * .. External Subroutines ..
257  EXTERNAL saxpy, scopy, sgemv, slabad, slaln2, sscal,
258  $ xerbla
259 * ..
260 * .. Intrinsic Functions ..
261  INTRINSIC abs, max, sqrt
262 * ..
263 * .. Local Arrays ..
264  REAL X( 2, 2 )
265 * ..
266 * .. Executable Statements ..
267 *
268 * Decode and test the input parameters
269 *
270  bothv = lsame( side, 'B' )
271  rightv = lsame( side, 'R' ) .OR. bothv
272  leftv = lsame( side, 'L' ) .OR. bothv
273 *
274  allv = lsame( howmny, 'A' )
275  over = lsame( howmny, 'B' )
276  somev = lsame( howmny, 'S' )
277 *
278  info = 0
279  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
280  info = -1
281  ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
282  info = -2
283  ELSE IF( n.LT.0 ) THEN
284  info = -4
285  ELSE IF( ldt.LT.max( 1, n ) ) THEN
286  info = -6
287  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
288  info = -8
289  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
290  info = -10
291  ELSE
292 *
293 * Set M to the number of columns required to store the selected
294 * eigenvectors, standardize the array SELECT if necessary, and
295 * test MM.
296 *
297  IF( somev ) THEN
298  m = 0
299  pair = .false.
300  DO 10 j = 1, n
301  IF( pair ) THEN
302  pair = .false.
303  SELECT( j ) = .false.
304  ELSE
305  IF( j.LT.n ) THEN
306  IF( t( j+1, j ).EQ.zero ) THEN
307  IF( SELECT( j ) )
308  $ m = m + 1
309  ELSE
310  pair = .true.
311  IF( SELECT( j ) .OR. SELECT( j+1 ) ) THEN
312  SELECT( j ) = .true.
313  m = m + 2
314  END IF
315  END IF
316  ELSE
317  IF( SELECT( n ) )
318  $ m = m + 1
319  END IF
320  END IF
321  10 CONTINUE
322  ELSE
323  m = n
324  END IF
325 *
326  IF( mm.LT.m ) THEN
327  info = -11
328  END IF
329  END IF
330  IF( info.NE.0 ) THEN
331  CALL xerbla( 'STREVC', -info )
332  RETURN
333  END IF
334 *
335 * Quick return if possible.
336 *
337  IF( n.EQ.0 )
338  $ RETURN
339 *
340 * Set the constants to control overflow.
341 *
342  unfl = slamch( 'Safe minimum' )
343  ovfl = one / unfl
344  CALL slabad( unfl, ovfl )
345  ulp = slamch( 'Precision' )
346  smlnum = unfl*( n / ulp )
347  bignum = ( one-ulp ) / smlnum
348 *
349 * Compute 1-norm of each column of strictly upper triangular
350 * part of T to control overflow in triangular solver.
351 *
352  work( 1 ) = zero
353  DO 30 j = 2, n
354  work( j ) = zero
355  DO 20 i = 1, j - 1
356  work( j ) = work( j ) + abs( t( i, j ) )
357  20 CONTINUE
358  30 CONTINUE
359 *
360 * Index IP is used to specify the real or complex eigenvalue:
361 * IP = 0, real eigenvalue,
362 * 1, first of conjugate complex pair: (wr,wi)
363 * -1, second of conjugate complex pair: (wr,wi)
364 *
365  n2 = 2*n
366 *
367  IF( rightv ) THEN
368 *
369 * Compute right eigenvectors.
370 *
371  ip = 0
372  is = m
373  DO 140 ki = n, 1, -1
374 *
375  IF( ip.EQ.1 )
376  $ GO TO 130
377  IF( ki.EQ.1 )
378  $ GO TO 40
379  IF( t( ki, ki-1 ).EQ.zero )
380  $ GO TO 40
381  ip = -1
382 *
383  40 CONTINUE
384  IF( somev ) THEN
385  IF( ip.EQ.0 ) THEN
386  IF( .NOT.SELECT( ki ) )
387  $ GO TO 130
388  ELSE
389  IF( .NOT.SELECT( ki-1 ) )
390  $ GO TO 130
391  END IF
392  END IF
393 *
394 * Compute the KI-th eigenvalue (WR,WI).
395 *
396  wr = t( ki, ki )
397  wi = zero
398  IF( ip.NE.0 )
399  $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
400  $ sqrt( abs( t( ki-1, ki ) ) )
401  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
402 *
403  IF( ip.EQ.0 ) THEN
404 *
405 * Real right eigenvector
406 *
407  work( ki+n ) = one
408 *
409 * Form right-hand side
410 *
411  DO 50 k = 1, ki - 1
412  work( k+n ) = -t( k, ki )
413  50 CONTINUE
414 *
415 * Solve the upper quasi-triangular system:
416 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
417 *
418  jnxt = ki - 1
419  DO 60 j = ki - 1, 1, -1
420  IF( j.GT.jnxt )
421  $ GO TO 60
422  j1 = j
423  j2 = j
424  jnxt = j - 1
425  IF( j.GT.1 ) THEN
426  IF( t( j, j-1 ).NE.zero ) THEN
427  j1 = j - 1
428  jnxt = j - 2
429  END IF
430  END IF
431 *
432  IF( j1.EQ.j2 ) THEN
433 *
434 * 1-by-1 diagonal block
435 *
436  CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
437  $ ldt, one, one, work( j+n ), n, wr,
438  $ zero, x, 2, scale, xnorm, ierr )
439 *
440 * Scale X(1,1) to avoid overflow when updating
441 * the right-hand side.
442 *
443  IF( xnorm.GT.one ) THEN
444  IF( work( j ).GT.bignum / xnorm ) THEN
445  x( 1, 1 ) = x( 1, 1 ) / xnorm
446  scale = scale / xnorm
447  END IF
448  END IF
449 *
450 * Scale if necessary
451 *
452  IF( scale.NE.one )
453  $ CALL sscal( ki, scale, work( 1+n ), 1 )
454  work( j+n ) = x( 1, 1 )
455 *
456 * Update right-hand side
457 *
458  CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
459  $ work( 1+n ), 1 )
460 *
461  ELSE
462 *
463 * 2-by-2 diagonal block
464 *
465  CALL slaln2( .false., 2, 1, smin, one,
466  $ t( j-1, j-1 ), ldt, one, one,
467  $ work( j-1+n ), n, wr, zero, x, 2,
468  $ scale, xnorm, ierr )
469 *
470 * Scale X(1,1) and X(2,1) to avoid overflow when
471 * updating the right-hand side.
472 *
473  IF( xnorm.GT.one ) THEN
474  beta = max( work( j-1 ), work( j ) )
475  IF( beta.GT.bignum / xnorm ) THEN
476  x( 1, 1 ) = x( 1, 1 ) / xnorm
477  x( 2, 1 ) = x( 2, 1 ) / xnorm
478  scale = scale / xnorm
479  END IF
480  END IF
481 *
482 * Scale if necessary
483 *
484  IF( scale.NE.one )
485  $ CALL sscal( ki, scale, work( 1+n ), 1 )
486  work( j-1+n ) = x( 1, 1 )
487  work( j+n ) = x( 2, 1 )
488 *
489 * Update right-hand side
490 *
491  CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
492  $ work( 1+n ), 1 )
493  CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
494  $ work( 1+n ), 1 )
495  END IF
496  60 CONTINUE
497 *
498 * Copy the vector x or Q*x to VR and normalize.
499 *
500  IF( .NOT.over ) THEN
501  CALL scopy( ki, work( 1+n ), 1, vr( 1, is ), 1 )
502 *
503  ii = isamax( ki, vr( 1, is ), 1 )
504  remax = one / abs( vr( ii, is ) )
505  CALL sscal( ki, remax, vr( 1, is ), 1 )
506 *
507  DO 70 k = ki + 1, n
508  vr( k, is ) = zero
509  70 CONTINUE
510  ELSE
511  IF( ki.GT.1 )
512  $ CALL sgemv( 'N', n, ki-1, one, vr, ldvr,
513  $ work( 1+n ), 1, work( ki+n ),
514  $ vr( 1, ki ), 1 )
515 *
516  ii = isamax( n, vr( 1, ki ), 1 )
517  remax = one / abs( vr( ii, ki ) )
518  CALL sscal( n, remax, vr( 1, ki ), 1 )
519  END IF
520 *
521  ELSE
522 *
523 * Complex right eigenvector.
524 *
525 * Initial solve
526 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
527 * [ (T(KI,KI-1) T(KI,KI) ) ]
528 *
529  IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) ) THEN
530  work( ki-1+n ) = one
531  work( ki+n2 ) = wi / t( ki-1, ki )
532  ELSE
533  work( ki-1+n ) = -wi / t( ki, ki-1 )
534  work( ki+n2 ) = one
535  END IF
536  work( ki+n ) = zero
537  work( ki-1+n2 ) = zero
538 *
539 * Form right-hand side
540 *
541  DO 80 k = 1, ki - 2
542  work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
543  work( k+n2 ) = -work( ki+n2 )*t( k, ki )
544  80 CONTINUE
545 *
546 * Solve upper quasi-triangular system:
547 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
548 *
549  jnxt = ki - 2
550  DO 90 j = ki - 2, 1, -1
551  IF( j.GT.jnxt )
552  $ GO TO 90
553  j1 = j
554  j2 = j
555  jnxt = j - 1
556  IF( j.GT.1 ) THEN
557  IF( t( j, j-1 ).NE.zero ) THEN
558  j1 = j - 1
559  jnxt = j - 2
560  END IF
561  END IF
562 *
563  IF( j1.EQ.j2 ) THEN
564 *
565 * 1-by-1 diagonal block
566 *
567  CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
568  $ ldt, one, one, work( j+n ), n, wr, wi,
569  $ x, 2, scale, xnorm, ierr )
570 *
571 * Scale X(1,1) and X(1,2) to avoid overflow when
572 * updating the right-hand side.
573 *
574  IF( xnorm.GT.one ) THEN
575  IF( work( j ).GT.bignum / xnorm ) THEN
576  x( 1, 1 ) = x( 1, 1 ) / xnorm
577  x( 1, 2 ) = x( 1, 2 ) / xnorm
578  scale = scale / xnorm
579  END IF
580  END IF
581 *
582 * Scale if necessary
583 *
584  IF( scale.NE.one ) THEN
585  CALL sscal( ki, scale, work( 1+n ), 1 )
586  CALL sscal( ki, scale, work( 1+n2 ), 1 )
587  END IF
588  work( j+n ) = x( 1, 1 )
589  work( j+n2 ) = x( 1, 2 )
590 *
591 * Update the right-hand side
592 *
593  CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
594  $ work( 1+n ), 1 )
595  CALL saxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
596  $ work( 1+n2 ), 1 )
597 *
598  ELSE
599 *
600 * 2-by-2 diagonal block
601 *
602  CALL slaln2( .false., 2, 2, smin, one,
603  $ t( j-1, j-1 ), ldt, one, one,
604  $ work( j-1+n ), n, wr, wi, x, 2, scale,
605  $ xnorm, ierr )
606 *
607 * Scale X to avoid overflow when updating
608 * the right-hand side.
609 *
610  IF( xnorm.GT.one ) THEN
611  beta = max( work( j-1 ), work( j ) )
612  IF( beta.GT.bignum / xnorm ) THEN
613  rec = one / xnorm
614  x( 1, 1 ) = x( 1, 1 )*rec
615  x( 1, 2 ) = x( 1, 2 )*rec
616  x( 2, 1 ) = x( 2, 1 )*rec
617  x( 2, 2 ) = x( 2, 2 )*rec
618  scale = scale*rec
619  END IF
620  END IF
621 *
622 * Scale if necessary
623 *
624  IF( scale.NE.one ) THEN
625  CALL sscal( ki, scale, work( 1+n ), 1 )
626  CALL sscal( ki, scale, work( 1+n2 ), 1 )
627  END IF
628  work( j-1+n ) = x( 1, 1 )
629  work( j+n ) = x( 2, 1 )
630  work( j-1+n2 ) = x( 1, 2 )
631  work( j+n2 ) = x( 2, 2 )
632 *
633 * Update the right-hand side
634 *
635  CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
636  $ work( 1+n ), 1 )
637  CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
638  $ work( 1+n ), 1 )
639  CALL saxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
640  $ work( 1+n2 ), 1 )
641  CALL saxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
642  $ work( 1+n2 ), 1 )
643  END IF
644  90 CONTINUE
645 *
646 * Copy the vector x or Q*x to VR and normalize.
647 *
648  IF( .NOT.over ) THEN
649  CALL scopy( ki, work( 1+n ), 1, vr( 1, is-1 ), 1 )
650  CALL scopy( ki, work( 1+n2 ), 1, vr( 1, is ), 1 )
651 *
652  emax = zero
653  DO 100 k = 1, ki
654  emax = max( emax, abs( vr( k, is-1 ) )+
655  $ abs( vr( k, is ) ) )
656  100 CONTINUE
657 *
658  remax = one / emax
659  CALL sscal( ki, remax, vr( 1, is-1 ), 1 )
660  CALL sscal( ki, remax, vr( 1, is ), 1 )
661 *
662  DO 110 k = ki + 1, n
663  vr( k, is-1 ) = zero
664  vr( k, is ) = zero
665  110 CONTINUE
666 *
667  ELSE
668 *
669  IF( ki.GT.2 ) THEN
670  CALL sgemv( 'N', n, ki-2, one, vr, ldvr,
671  $ work( 1+n ), 1, work( ki-1+n ),
672  $ vr( 1, ki-1 ), 1 )
673  CALL sgemv( 'N', n, ki-2, one, vr, ldvr,
674  $ work( 1+n2 ), 1, work( ki+n2 ),
675  $ vr( 1, ki ), 1 )
676  ELSE
677  CALL sscal( n, work( ki-1+n ), vr( 1, ki-1 ), 1 )
678  CALL sscal( n, work( ki+n2 ), vr( 1, ki ), 1 )
679  END IF
680 *
681  emax = zero
682  DO 120 k = 1, n
683  emax = max( emax, abs( vr( k, ki-1 ) )+
684  $ abs( vr( k, ki ) ) )
685  120 CONTINUE
686  remax = one / emax
687  CALL sscal( n, remax, vr( 1, ki-1 ), 1 )
688  CALL sscal( n, remax, vr( 1, ki ), 1 )
689  END IF
690  END IF
691 *
692  is = is - 1
693  IF( ip.NE.0 )
694  $ is = is - 1
695  130 CONTINUE
696  IF( ip.EQ.1 )
697  $ ip = 0
698  IF( ip.EQ.-1 )
699  $ ip = 1
700  140 CONTINUE
701  END IF
702 *
703  IF( leftv ) THEN
704 *
705 * Compute left eigenvectors.
706 *
707  ip = 0
708  is = 1
709  DO 260 ki = 1, n
710 *
711  IF( ip.EQ.-1 )
712  $ GO TO 250
713  IF( ki.EQ.n )
714  $ GO TO 150
715  IF( t( ki+1, ki ).EQ.zero )
716  $ GO TO 150
717  ip = 1
718 *
719  150 CONTINUE
720  IF( somev ) THEN
721  IF( .NOT.SELECT( ki ) )
722  $ GO TO 250
723  END IF
724 *
725 * Compute the KI-th eigenvalue (WR,WI).
726 *
727  wr = t( ki, ki )
728  wi = zero
729  IF( ip.NE.0 )
730  $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
731  $ sqrt( abs( t( ki+1, ki ) ) )
732  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
733 *
734  IF( ip.EQ.0 ) THEN
735 *
736 * Real left eigenvector.
737 *
738  work( ki+n ) = one
739 *
740 * Form right-hand side
741 *
742  DO 160 k = ki + 1, n
743  work( k+n ) = -t( ki, k )
744  160 CONTINUE
745 *
746 * Solve the quasi-triangular system:
747 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
748 *
749  vmax = one
750  vcrit = bignum
751 *
752  jnxt = ki + 1
753  DO 170 j = ki + 1, n
754  IF( j.LT.jnxt )
755  $ GO TO 170
756  j1 = j
757  j2 = j
758  jnxt = j + 1
759  IF( j.LT.n ) THEN
760  IF( t( j+1, j ).NE.zero ) THEN
761  j2 = j + 1
762  jnxt = j + 2
763  END IF
764  END IF
765 *
766  IF( j1.EQ.j2 ) THEN
767 *
768 * 1-by-1 diagonal block
769 *
770 * Scale if necessary to avoid overflow when forming
771 * the right-hand side.
772 *
773  IF( work( j ).GT.vcrit ) THEN
774  rec = one / vmax
775  CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
776  vmax = one
777  vcrit = bignum
778  END IF
779 *
780  work( j+n ) = work( j+n ) -
781  $ sdot( j-ki-1, t( ki+1, j ), 1,
782  $ work( ki+1+n ), 1 )
783 *
784 * Solve (T(J,J)-WR)**T*X = WORK
785 *
786  CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
787  $ ldt, one, one, work( j+n ), n, wr,
788  $ zero, x, 2, scale, xnorm, ierr )
789 *
790 * Scale if necessary
791 *
792  IF( scale.NE.one )
793  $ CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
794  work( j+n ) = x( 1, 1 )
795  vmax = max( abs( work( j+n ) ), vmax )
796  vcrit = bignum / vmax
797 *
798  ELSE
799 *
800 * 2-by-2 diagonal block
801 *
802 * Scale if necessary to avoid overflow when forming
803 * the right-hand side.
804 *
805  beta = max( work( j ), work( j+1 ) )
806  IF( beta.GT.vcrit ) THEN
807  rec = one / vmax
808  CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
809  vmax = one
810  vcrit = bignum
811  END IF
812 *
813  work( j+n ) = work( j+n ) -
814  $ sdot( j-ki-1, t( ki+1, j ), 1,
815  $ work( ki+1+n ), 1 )
816 *
817  work( j+1+n ) = work( j+1+n ) -
818  $ sdot( j-ki-1, t( ki+1, j+1 ), 1,
819  $ work( ki+1+n ), 1 )
820 *
821 * Solve
822 * [T(J,J)-WR T(J,J+1) ]**T* X = SCALE*( WORK1 )
823 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
824 *
825  CALL slaln2( .true., 2, 1, smin, one, t( j, j ),
826  $ ldt, one, one, work( j+n ), n, wr,
827  $ zero, x, 2, scale, xnorm, ierr )
828 *
829 * Scale if necessary
830 *
831  IF( scale.NE.one )
832  $ CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
833  work( j+n ) = x( 1, 1 )
834  work( j+1+n ) = x( 2, 1 )
835 *
836  vmax = max( abs( work( j+n ) ),
837  $ abs( work( j+1+n ) ), vmax )
838  vcrit = bignum / vmax
839 *
840  END IF
841  170 CONTINUE
842 *
843 * Copy the vector x or Q*x to VL and normalize.
844 *
845  IF( .NOT.over ) THEN
846  CALL scopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
847 *
848  ii = isamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
849  remax = one / abs( vl( ii, is ) )
850  CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
851 *
852  DO 180 k = 1, ki - 1
853  vl( k, is ) = zero
854  180 CONTINUE
855 *
856  ELSE
857 *
858  IF( ki.LT.n )
859  $ CALL sgemv( 'N', n, n-ki, one, vl( 1, ki+1 ), ldvl,
860  $ work( ki+1+n ), 1, work( ki+n ),
861  $ vl( 1, ki ), 1 )
862 *
863  ii = isamax( n, vl( 1, ki ), 1 )
864  remax = one / abs( vl( ii, ki ) )
865  CALL sscal( n, remax, vl( 1, ki ), 1 )
866 *
867  END IF
868 *
869  ELSE
870 *
871 * Complex left eigenvector.
872 *
873 * Initial solve:
874 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
875 * ((T(KI+1,KI) T(KI+1,KI+1)) )
876 *
877  IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) ) THEN
878  work( ki+n ) = wi / t( ki, ki+1 )
879  work( ki+1+n2 ) = one
880  ELSE
881  work( ki+n ) = one
882  work( ki+1+n2 ) = -wi / t( ki+1, ki )
883  END IF
884  work( ki+1+n ) = zero
885  work( ki+n2 ) = zero
886 *
887 * Form right-hand side
888 *
889  DO 190 k = ki + 2, n
890  work( k+n ) = -work( ki+n )*t( ki, k )
891  work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
892  190 CONTINUE
893 *
894 * Solve complex quasi-triangular system:
895 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
896 *
897  vmax = one
898  vcrit = bignum
899 *
900  jnxt = ki + 2
901  DO 200 j = ki + 2, n
902  IF( j.LT.jnxt )
903  $ GO TO 200
904  j1 = j
905  j2 = j
906  jnxt = j + 1
907  IF( j.LT.n ) THEN
908  IF( t( j+1, j ).NE.zero ) THEN
909  j2 = j + 1
910  jnxt = j + 2
911  END IF
912  END IF
913 *
914  IF( j1.EQ.j2 ) THEN
915 *
916 * 1-by-1 diagonal block
917 *
918 * Scale if necessary to avoid overflow when
919 * forming the right-hand side elements.
920 *
921  IF( work( j ).GT.vcrit ) THEN
922  rec = one / vmax
923  CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
924  CALL sscal( n-ki+1, rec, work( ki+n2 ), 1 )
925  vmax = one
926  vcrit = bignum
927  END IF
928 *
929  work( j+n ) = work( j+n ) -
930  $ sdot( j-ki-2, t( ki+2, j ), 1,
931  $ work( ki+2+n ), 1 )
932  work( j+n2 ) = work( j+n2 ) -
933  $ sdot( j-ki-2, t( ki+2, j ), 1,
934  $ work( ki+2+n2 ), 1 )
935 *
936 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
937 *
938  CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
939  $ ldt, one, one, work( j+n ), n, wr,
940  $ -wi, x, 2, scale, xnorm, ierr )
941 *
942 * Scale if necessary
943 *
944  IF( scale.NE.one ) THEN
945  CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
946  CALL sscal( n-ki+1, scale, work( ki+n2 ), 1 )
947  END IF
948  work( j+n ) = x( 1, 1 )
949  work( j+n2 ) = x( 1, 2 )
950  vmax = max( abs( work( j+n ) ),
951  $ abs( work( j+n2 ) ), vmax )
952  vcrit = bignum / vmax
953 *
954  ELSE
955 *
956 * 2-by-2 diagonal block
957 *
958 * Scale if necessary to avoid overflow when forming
959 * the right-hand side elements.
960 *
961  beta = max( work( j ), work( j+1 ) )
962  IF( beta.GT.vcrit ) THEN
963  rec = one / vmax
964  CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
965  CALL sscal( n-ki+1, rec, work( ki+n2 ), 1 )
966  vmax = one
967  vcrit = bignum
968  END IF
969 *
970  work( j+n ) = work( j+n ) -
971  $ sdot( j-ki-2, t( ki+2, j ), 1,
972  $ work( ki+2+n ), 1 )
973 *
974  work( j+n2 ) = work( j+n2 ) -
975  $ sdot( j-ki-2, t( ki+2, j ), 1,
976  $ work( ki+2+n2 ), 1 )
977 *
978  work( j+1+n ) = work( j+1+n ) -
979  $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
980  $ work( ki+2+n ), 1 )
981 *
982  work( j+1+n2 ) = work( j+1+n2 ) -
983  $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
984  $ work( ki+2+n2 ), 1 )
985 *
986 * Solve 2-by-2 complex linear equation
987 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
988 * ([T(j+1,j) T(j+1,j+1)] )
989 *
990  CALL slaln2( .true., 2, 2, smin, one, t( j, j ),
991  $ ldt, one, one, work( j+n ), n, wr,
992  $ -wi, x, 2, scale, xnorm, ierr )
993 *
994 * Scale if necessary
995 *
996  IF( scale.NE.one ) THEN
997  CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
998  CALL sscal( n-ki+1, scale, work( ki+n2 ), 1 )
999  END IF
1000  work( j+n ) = x( 1, 1 )
1001  work( j+n2 ) = x( 1, 2 )
1002  work( j+1+n ) = x( 2, 1 )
1003  work( j+1+n2 ) = x( 2, 2 )
1004  vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1005  $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ), vmax )
1006  vcrit = bignum / vmax
1007 *
1008  END IF
1009  200 CONTINUE
1010 *
1011 * Copy the vector x or Q*x to VL and normalize.
1012 *
1013  IF( .NOT.over ) THEN
1014  CALL scopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
1015  CALL scopy( n-ki+1, work( ki+n2 ), 1, vl( ki, is+1 ),
1016  $ 1 )
1017 *
1018  emax = zero
1019  DO 220 k = ki, n
1020  emax = max( emax, abs( vl( k, is ) )+
1021  $ abs( vl( k, is+1 ) ) )
1022  220 CONTINUE
1023  remax = one / emax
1024  CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1025  CALL sscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1026 *
1027  DO 230 k = 1, ki - 1
1028  vl( k, is ) = zero
1029  vl( k, is+1 ) = zero
1030  230 CONTINUE
1031  ELSE
1032  IF( ki.LT.n-1 ) THEN
1033  CALL sgemv( 'N', n, n-ki-1, one, vl( 1, ki+2 ),
1034  $ ldvl, work( ki+2+n ), 1, work( ki+n ),
1035  $ vl( 1, ki ), 1 )
1036  CALL sgemv( 'N', n, n-ki-1, one, vl( 1, ki+2 ),
1037  $ ldvl, work( ki+2+n2 ), 1,
1038  $ work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1039  ELSE
1040  CALL sscal( n, work( ki+n ), vl( 1, ki ), 1 )
1041  CALL sscal( n, work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1042  END IF
1043 *
1044  emax = zero
1045  DO 240 k = 1, n
1046  emax = max( emax, abs( vl( k, ki ) )+
1047  $ abs( vl( k, ki+1 ) ) )
1048  240 CONTINUE
1049  remax = one / emax
1050  CALL sscal( n, remax, vl( 1, ki ), 1 )
1051  CALL sscal( n, remax, vl( 1, ki+1 ), 1 )
1052 *
1053  END IF
1054 *
1055  END IF
1056 *
1057  is = is + 1
1058  IF( ip.NE.0 )
1059  $ is = is + 1
1060  250 CONTINUE
1061  IF( ip.EQ.-1 )
1062  $ ip = 0
1063  IF( ip.EQ.1 )
1064  $ ip = -1
1065 *
1066  260 CONTINUE
1067 *
1068  END IF
1069 *
1070  RETURN
1071 *
1072 * End of STREVC
1073 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition: slaln2.f:218
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: