LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dtrevc ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
double precision, dimension( * )  WORK,
integer  INFO 
)

DTREVC

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

Purpose:
 DTREVC 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 DHSEQR.
 
 The right eigenvector x and the left eigenvector y of T corresponding
 to an eigenvalue w are defined by:
 
    T*x = w*x,     (y**T)*T = w*(y**T)
 
 where y**T denotes the 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DHSEQR).
          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 DOUBLE PRECISION 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 DHSEQR).
          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 DOUBLE PRECISION 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.
Date
November 2011
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 224 of file dtrevc.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: