LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dtgevc ( character  SIDE,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
double precision, dimension( lds, * )  S,
integer  LDS,
double precision, dimension( ldp, * )  P,
integer  LDP,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
integer  MM,
integer  M,
double precision, dimension( * )  WORK,
integer  INFO 
)

DTGEVC

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

Purpose:
 DTGEVC computes some or all of the right and/or left eigenvectors of
 a pair of real matrices (S,P), where S is a quasi-triangular matrix
 and P is upper triangular.  Matrix pairs of this type are produced by
 the generalized Schur factorization of a matrix pair (A,B):

    A = Q*S*Z**T,  B = Q*P*Z**T

 as computed by DGGHRD + DHGEQZ.

 The right eigenvector x and the left eigenvector y of (S,P)
 corresponding to an eigenvalue w are defined by:
 
    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
 
 where y**H denotes the conjugate tranpose of y.
 The eigenvalues are not input to this routine, but are computed
 directly from the diagonal blocks of S and P.
 
 This routine returns the matrices X and/or Y of right and left
 eigenvectors of (S,P), or the products Z*X and/or Q*Y,
 where Z and Q are input matrices.
 If Q and Z are the orthogonal factors from the generalized Schur
 factorization of a matrix pair (A,B), then Z*X and Q*Y
 are the matrices of right and left eigenvectors of (A,B).
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,
                 specified by the logical array SELECT.
[in]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 matrices S and P.  N >= 0.
[in]S
          S is DOUBLE PRECISION array, dimension (LDS,N)
          The upper quasi-triangular matrix S from a generalized Schur
          factorization, as computed by DHGEQZ.
[in]LDS
          LDS is INTEGER
          The leading dimension of array S.  LDS >= max(1,N).
[in]P
          P is DOUBLE PRECISION array, dimension (LDP,N)
          The upper triangular matrix P from a generalized Schur
          factorization, as computed by DHGEQZ.
          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
          of S must be in positive diagonal form.
[in]LDP
          LDP is INTEGER
          The leading dimension of array P.  LDP >= 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 left Schur vectors returned by DHGEQZ).
          On exit, if SIDE = 'L' or 'B', VL contains:
          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
          if HOWMNY = 'B', the matrix Q*Y;
          if HOWMNY = 'S', the left eigenvectors of (S,P) 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 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 Z (usually the orthogonal matrix Z
          of right Schur vectors returned by DHGEQZ).

          On exit, if SIDE = 'R' or 'B', VR contains:
          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
          if HOWMNY = 'B' or 'b', the matrix Z*X;
          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
                      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 (6*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
                eigenvalue.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  Allocation of workspace:
  ---------- -- ---------

     WORK( j ) = 1-norm of j-th column of A, above the diagonal
     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
     WORK( 2*N+1:3*N ) = real part of eigenvector
     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector

  Rowwise vs. columnwise solution methods:
  ------- --  ---------- -------- -------

  Finding a generalized eigenvector consists basically of solving the
  singular triangular system

   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)

  Consider finding the i-th right eigenvector (assume all eigenvalues
  are real). The equation to be solved is:
       n                   i
  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
      k=j                 k=j

  where  C = (A - w B)  (The components v(i+1:n) are 0.)

  The "rowwise" method is:

  (1)  v(i) := 1
  for j = i-1,. . .,1:
                          i
      (2) compute  s = - sum C(j,k) v(k)   and
                        k=j+1

      (3) v(j) := s / C(j,j)

  Step 2 is sometimes called the "dot product" step, since it is an
  inner product between the j-th row and the portion of the eigenvector
  that has been computed so far.

  The "columnwise" method consists basically in doing the sums
  for all the rows in parallel.  As each v(j) is computed, the
  contribution of v(j) times the j-th column of C is added to the
  partial sums.  Since FORTRAN arrays are stored columnwise, this has
  the advantage that at each step, the elements of C that are accessed
  are adjacent to one another, whereas with the rowwise method, the
  elements accessed at a step are spaced LDS (and LDP) words apart.

  When finding left eigenvectors, the matrix in question is the
  transpose of the one in storage, so the rowwise method then
  actually accesses columns of A and B at each step, and so is the
  preferred method.

Definition at line 297 of file dtgevc.f.

297 *
298 * -- LAPACK computational routine (version 3.4.0) --
299 * -- LAPACK is a software package provided by Univ. of Tennessee, --
300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301 * November 2011
302 *
303 * .. Scalar Arguments ..
304  CHARACTER howmny, side
305  INTEGER info, ldp, lds, ldvl, ldvr, m, mm, n
306 * ..
307 * .. Array Arguments ..
308  LOGICAL select( * )
309  DOUBLE PRECISION p( ldp, * ), s( lds, * ), vl( ldvl, * ),
310  $ vr( ldvr, * ), work( * )
311 * ..
312 *
313 *
314 * =====================================================================
315 *
316 * .. Parameters ..
317  DOUBLE PRECISION zero, one, safety
318  parameter ( zero = 0.0d+0, one = 1.0d+0,
319  $ safety = 1.0d+2 )
320 * ..
321 * .. Local Scalars ..
322  LOGICAL compl, compr, il2by2, ilabad, ilall, ilback,
323  $ ilbbad, ilcomp, ilcplx, lsa, lsb
324  INTEGER i, ibeg, ieig, iend, ihwmny, iinfo, im, iside,
325  $ j, ja, jc, je, jr, jw, na, nw
326  DOUBLE PRECISION acoef, acoefa, anorm, ascale, bcoefa, bcoefi,
327  $ bcoefr, big, bignum, bnorm, bscale, cim2a,
328  $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
329  $ crealb, dmin, safmin, salfar, sbeta, scale,
330  $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
331  $ xscale
332 * ..
333 * .. Local Arrays ..
334  DOUBLE PRECISION bdiag( 2 ), sum( 2, 2 ), sums( 2, 2 ),
335  $ sump( 2, 2 )
336 * ..
337 * .. External Functions ..
338  LOGICAL lsame
339  DOUBLE PRECISION dlamch
340  EXTERNAL lsame, dlamch
341 * ..
342 * .. External Subroutines ..
343  EXTERNAL dgemv, dlabad, dlacpy, dlag2, dlaln2, xerbla
344 * ..
345 * .. Intrinsic Functions ..
346  INTRINSIC abs, max, min
347 * ..
348 * .. Executable Statements ..
349 *
350 * Decode and Test the input parameters
351 *
352  IF( lsame( howmny, 'A' ) ) THEN
353  ihwmny = 1
354  ilall = .true.
355  ilback = .false.
356  ELSE IF( lsame( howmny, 'S' ) ) THEN
357  ihwmny = 2
358  ilall = .false.
359  ilback = .false.
360  ELSE IF( lsame( howmny, 'B' ) ) THEN
361  ihwmny = 3
362  ilall = .true.
363  ilback = .true.
364  ELSE
365  ihwmny = -1
366  ilall = .true.
367  END IF
368 *
369  IF( lsame( side, 'R' ) ) THEN
370  iside = 1
371  compl = .false.
372  compr = .true.
373  ELSE IF( lsame( side, 'L' ) ) THEN
374  iside = 2
375  compl = .true.
376  compr = .false.
377  ELSE IF( lsame( side, 'B' ) ) THEN
378  iside = 3
379  compl = .true.
380  compr = .true.
381  ELSE
382  iside = -1
383  END IF
384 *
385  info = 0
386  IF( iside.LT.0 ) THEN
387  info = -1
388  ELSE IF( ihwmny.LT.0 ) THEN
389  info = -2
390  ELSE IF( n.LT.0 ) THEN
391  info = -4
392  ELSE IF( lds.LT.max( 1, n ) ) THEN
393  info = -6
394  ELSE IF( ldp.LT.max( 1, n ) ) THEN
395  info = -8
396  END IF
397  IF( info.NE.0 ) THEN
398  CALL xerbla( 'DTGEVC', -info )
399  RETURN
400  END IF
401 *
402 * Count the number of eigenvectors to be computed
403 *
404  IF( .NOT.ilall ) THEN
405  im = 0
406  ilcplx = .false.
407  DO 10 j = 1, n
408  IF( ilcplx ) THEN
409  ilcplx = .false.
410  GO TO 10
411  END IF
412  IF( j.LT.n ) THEN
413  IF( s( j+1, j ).NE.zero )
414  $ ilcplx = .true.
415  END IF
416  IF( ilcplx ) THEN
417  IF( SELECT( j ) .OR. SELECT( j+1 ) )
418  $ im = im + 2
419  ELSE
420  IF( SELECT( j ) )
421  $ im = im + 1
422  END IF
423  10 CONTINUE
424  ELSE
425  im = n
426  END IF
427 *
428 * Check 2-by-2 diagonal blocks of A, B
429 *
430  ilabad = .false.
431  ilbbad = .false.
432  DO 20 j = 1, n - 1
433  IF( s( j+1, j ).NE.zero ) THEN
434  IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
435  $ p( j, j+1 ).NE.zero )ilbbad = .true.
436  IF( j.LT.n-1 ) THEN
437  IF( s( j+2, j+1 ).NE.zero )
438  $ ilabad = .true.
439  END IF
440  END IF
441  20 CONTINUE
442 *
443  IF( ilabad ) THEN
444  info = -5
445  ELSE IF( ilbbad ) THEN
446  info = -7
447  ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
448  info = -10
449  ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
450  info = -12
451  ELSE IF( mm.LT.im ) THEN
452  info = -13
453  END IF
454  IF( info.NE.0 ) THEN
455  CALL xerbla( 'DTGEVC', -info )
456  RETURN
457  END IF
458 *
459 * Quick return if possible
460 *
461  m = im
462  IF( n.EQ.0 )
463  $ RETURN
464 *
465 * Machine Constants
466 *
467  safmin = dlamch( 'Safe minimum' )
468  big = one / safmin
469  CALL dlabad( safmin, big )
470  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
471  small = safmin*n / ulp
472  big = one / small
473  bignum = one / ( safmin*n )
474 *
475 * Compute the 1-norm of each column of the strictly upper triangular
476 * part (i.e., excluding all elements belonging to the diagonal
477 * blocks) of A and B to check for possible overflow in the
478 * triangular solver.
479 *
480  anorm = abs( s( 1, 1 ) )
481  IF( n.GT.1 )
482  $ anorm = anorm + abs( s( 2, 1 ) )
483  bnorm = abs( p( 1, 1 ) )
484  work( 1 ) = zero
485  work( n+1 ) = zero
486 *
487  DO 50 j = 2, n
488  temp = zero
489  temp2 = zero
490  IF( s( j, j-1 ).EQ.zero ) THEN
491  iend = j - 1
492  ELSE
493  iend = j - 2
494  END IF
495  DO 30 i = 1, iend
496  temp = temp + abs( s( i, j ) )
497  temp2 = temp2 + abs( p( i, j ) )
498  30 CONTINUE
499  work( j ) = temp
500  work( n+j ) = temp2
501  DO 40 i = iend + 1, min( j+1, n )
502  temp = temp + abs( s( i, j ) )
503  temp2 = temp2 + abs( p( i, j ) )
504  40 CONTINUE
505  anorm = max( anorm, temp )
506  bnorm = max( bnorm, temp2 )
507  50 CONTINUE
508 *
509  ascale = one / max( anorm, safmin )
510  bscale = one / max( bnorm, safmin )
511 *
512 * Left eigenvectors
513 *
514  IF( compl ) THEN
515  ieig = 0
516 *
517 * Main loop over eigenvalues
518 *
519  ilcplx = .false.
520  DO 220 je = 1, n
521 *
522 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
523 * (b) this would be the second of a complex pair.
524 * Check for complex eigenvalue, so as to be sure of which
525 * entry(-ies) of SELECT to look at.
526 *
527  IF( ilcplx ) THEN
528  ilcplx = .false.
529  GO TO 220
530  END IF
531  nw = 1
532  IF( je.LT.n ) THEN
533  IF( s( je+1, je ).NE.zero ) THEN
534  ilcplx = .true.
535  nw = 2
536  END IF
537  END IF
538  IF( ilall ) THEN
539  ilcomp = .true.
540  ELSE IF( ilcplx ) THEN
541  ilcomp = SELECT( je ) .OR. SELECT( je+1 )
542  ELSE
543  ilcomp = SELECT( je )
544  END IF
545  IF( .NOT.ilcomp )
546  $ GO TO 220
547 *
548 * Decide if (a) singular pencil, (b) real eigenvalue, or
549 * (c) complex eigenvalue.
550 *
551  IF( .NOT.ilcplx ) THEN
552  IF( abs( s( je, je ) ).LE.safmin .AND.
553  $ abs( p( je, je ) ).LE.safmin ) THEN
554 *
555 * Singular matrix pencil -- return unit eigenvector
556 *
557  ieig = ieig + 1
558  DO 60 jr = 1, n
559  vl( jr, ieig ) = zero
560  60 CONTINUE
561  vl( ieig, ieig ) = one
562  GO TO 220
563  END IF
564  END IF
565 *
566 * Clear vector
567 *
568  DO 70 jr = 1, nw*n
569  work( 2*n+jr ) = zero
570  70 CONTINUE
571 * T
572 * Compute coefficients in ( a A - b B ) y = 0
573 * a is ACOEF
574 * b is BCOEFR + i*BCOEFI
575 *
576  IF( .NOT.ilcplx ) THEN
577 *
578 * Real eigenvalue
579 *
580  temp = one / max( abs( s( je, je ) )*ascale,
581  $ abs( p( je, je ) )*bscale, safmin )
582  salfar = ( temp*s( je, je ) )*ascale
583  sbeta = ( temp*p( je, je ) )*bscale
584  acoef = sbeta*ascale
585  bcoefr = salfar*bscale
586  bcoefi = zero
587 *
588 * Scale to avoid underflow
589 *
590  scale = one
591  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
592  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
593  $ small
594  IF( lsa )
595  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
596  IF( lsb )
597  $ scale = max( scale, ( small / abs( salfar ) )*
598  $ min( bnorm, big ) )
599  IF( lsa .OR. lsb ) THEN
600  scale = min( scale, one /
601  $ ( safmin*max( one, abs( acoef ),
602  $ abs( bcoefr ) ) ) )
603  IF( lsa ) THEN
604  acoef = ascale*( scale*sbeta )
605  ELSE
606  acoef = scale*acoef
607  END IF
608  IF( lsb ) THEN
609  bcoefr = bscale*( scale*salfar )
610  ELSE
611  bcoefr = scale*bcoefr
612  END IF
613  END IF
614  acoefa = abs( acoef )
615  bcoefa = abs( bcoefr )
616 *
617 * First component is 1
618 *
619  work( 2*n+je ) = one
620  xmax = one
621  ELSE
622 *
623 * Complex eigenvalue
624 *
625  CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
626  $ safmin*safety, acoef, temp, bcoefr, temp2,
627  $ bcoefi )
628  bcoefi = -bcoefi
629  IF( bcoefi.EQ.zero ) THEN
630  info = je
631  RETURN
632  END IF
633 *
634 * Scale to avoid over/underflow
635 *
636  acoefa = abs( acoef )
637  bcoefa = abs( bcoefr ) + abs( bcoefi )
638  scale = one
639  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
640  $ scale = ( safmin / ulp ) / acoefa
641  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
642  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
643  IF( safmin*acoefa.GT.ascale )
644  $ scale = ascale / ( safmin*acoefa )
645  IF( safmin*bcoefa.GT.bscale )
646  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
647  IF( scale.NE.one ) THEN
648  acoef = scale*acoef
649  acoefa = abs( acoef )
650  bcoefr = scale*bcoefr
651  bcoefi = scale*bcoefi
652  bcoefa = abs( bcoefr ) + abs( bcoefi )
653  END IF
654 *
655 * Compute first two components of eigenvector
656 *
657  temp = acoef*s( je+1, je )
658  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
659  temp2i = -bcoefi*p( je, je )
660  IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
661  work( 2*n+je ) = one
662  work( 3*n+je ) = zero
663  work( 2*n+je+1 ) = -temp2r / temp
664  work( 3*n+je+1 ) = -temp2i / temp
665  ELSE
666  work( 2*n+je+1 ) = one
667  work( 3*n+je+1 ) = zero
668  temp = acoef*s( je, je+1 )
669  work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
670  $ s( je+1, je+1 ) ) / temp
671  work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
672  END IF
673  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
674  $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
675  END IF
676 *
677  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
678 *
679 * T
680 * Triangular solve of (a A - b B) y = 0
681 *
682 * T
683 * (rowwise in (a A - b B) , or columnwise in (a A - b B) )
684 *
685  il2by2 = .false.
686 *
687  DO 160 j = je + nw, n
688  IF( il2by2 ) THEN
689  il2by2 = .false.
690  GO TO 160
691  END IF
692 *
693  na = 1
694  bdiag( 1 ) = p( j, j )
695  IF( j.LT.n ) THEN
696  IF( s( j+1, j ).NE.zero ) THEN
697  il2by2 = .true.
698  bdiag( 2 ) = p( j+1, j+1 )
699  na = 2
700  END IF
701  END IF
702 *
703 * Check whether scaling is necessary for dot products
704 *
705  xscale = one / max( one, xmax )
706  temp = max( work( j ), work( n+j ),
707  $ acoefa*work( j )+bcoefa*work( n+j ) )
708  IF( il2by2 )
709  $ temp = max( temp, work( j+1 ), work( n+j+1 ),
710  $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
711  IF( temp.GT.bignum*xscale ) THEN
712  DO 90 jw = 0, nw - 1
713  DO 80 jr = je, j - 1
714  work( ( jw+2 )*n+jr ) = xscale*
715  $ work( ( jw+2 )*n+jr )
716  80 CONTINUE
717  90 CONTINUE
718  xmax = xmax*xscale
719  END IF
720 *
721 * Compute dot products
722 *
723 * j-1
724 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
725 * k=je
726 *
727 * To reduce the op count, this is done as
728 *
729 * _ j-1 _ j-1
730 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
731 * k=je k=je
732 *
733 * which may cause underflow problems if A or B are close
734 * to underflow. (E.g., less than SMALL.)
735 *
736 *
737  DO 120 jw = 1, nw
738  DO 110 ja = 1, na
739  sums( ja, jw ) = zero
740  sump( ja, jw ) = zero
741 *
742  DO 100 jr = je, j - 1
743  sums( ja, jw ) = sums( ja, jw ) +
744  $ s( jr, j+ja-1 )*
745  $ work( ( jw+1 )*n+jr )
746  sump( ja, jw ) = sump( ja, jw ) +
747  $ p( jr, j+ja-1 )*
748  $ work( ( jw+1 )*n+jr )
749  100 CONTINUE
750  110 CONTINUE
751  120 CONTINUE
752 *
753  DO 130 ja = 1, na
754  IF( ilcplx ) THEN
755  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
756  $ bcoefr*sump( ja, 1 ) -
757  $ bcoefi*sump( ja, 2 )
758  sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
759  $ bcoefr*sump( ja, 2 ) +
760  $ bcoefi*sump( ja, 1 )
761  ELSE
762  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
763  $ bcoefr*sump( ja, 1 )
764  END IF
765  130 CONTINUE
766 *
767 * T
768 * Solve ( a A - b B ) y = SUM(,)
769 * with scaling and perturbation of the denominator
770 *
771  CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
772  $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
773  $ bcoefi, work( 2*n+j ), n, scale, temp,
774  $ iinfo )
775  IF( scale.LT.one ) THEN
776  DO 150 jw = 0, nw - 1
777  DO 140 jr = je, j - 1
778  work( ( jw+2 )*n+jr ) = scale*
779  $ work( ( jw+2 )*n+jr )
780  140 CONTINUE
781  150 CONTINUE
782  xmax = scale*xmax
783  END IF
784  xmax = max( xmax, temp )
785  160 CONTINUE
786 *
787 * Copy eigenvector to VL, back transforming if
788 * HOWMNY='B'.
789 *
790  ieig = ieig + 1
791  IF( ilback ) THEN
792  DO 170 jw = 0, nw - 1
793  CALL dgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
794  $ work( ( jw+2 )*n+je ), 1, zero,
795  $ work( ( jw+4 )*n+1 ), 1 )
796  170 CONTINUE
797  CALL dlacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
798  $ ldvl )
799  ibeg = 1
800  ELSE
801  CALL dlacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
802  $ ldvl )
803  ibeg = je
804  END IF
805 *
806 * Scale eigenvector
807 *
808  xmax = zero
809  IF( ilcplx ) THEN
810  DO 180 j = ibeg, n
811  xmax = max( xmax, abs( vl( j, ieig ) )+
812  $ abs( vl( j, ieig+1 ) ) )
813  180 CONTINUE
814  ELSE
815  DO 190 j = ibeg, n
816  xmax = max( xmax, abs( vl( j, ieig ) ) )
817  190 CONTINUE
818  END IF
819 *
820  IF( xmax.GT.safmin ) THEN
821  xscale = one / xmax
822 *
823  DO 210 jw = 0, nw - 1
824  DO 200 jr = ibeg, n
825  vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
826  200 CONTINUE
827  210 CONTINUE
828  END IF
829  ieig = ieig + nw - 1
830 *
831  220 CONTINUE
832  END IF
833 *
834 * Right eigenvectors
835 *
836  IF( compr ) THEN
837  ieig = im + 1
838 *
839 * Main loop over eigenvalues
840 *
841  ilcplx = .false.
842  DO 500 je = n, 1, -1
843 *
844 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
845 * (b) this would be the second of a complex pair.
846 * Check for complex eigenvalue, so as to be sure of which
847 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
848 * or SELECT(JE-1).
849 * If this is a complex pair, the 2-by-2 diagonal block
850 * corresponding to the eigenvalue is in rows/columns JE-1:JE
851 *
852  IF( ilcplx ) THEN
853  ilcplx = .false.
854  GO TO 500
855  END IF
856  nw = 1
857  IF( je.GT.1 ) THEN
858  IF( s( je, je-1 ).NE.zero ) THEN
859  ilcplx = .true.
860  nw = 2
861  END IF
862  END IF
863  IF( ilall ) THEN
864  ilcomp = .true.
865  ELSE IF( ilcplx ) THEN
866  ilcomp = SELECT( je ) .OR. SELECT( je-1 )
867  ELSE
868  ilcomp = SELECT( je )
869  END IF
870  IF( .NOT.ilcomp )
871  $ GO TO 500
872 *
873 * Decide if (a) singular pencil, (b) real eigenvalue, or
874 * (c) complex eigenvalue.
875 *
876  IF( .NOT.ilcplx ) THEN
877  IF( abs( s( je, je ) ).LE.safmin .AND.
878  $ abs( p( je, je ) ).LE.safmin ) THEN
879 *
880 * Singular matrix pencil -- unit eigenvector
881 *
882  ieig = ieig - 1
883  DO 230 jr = 1, n
884  vr( jr, ieig ) = zero
885  230 CONTINUE
886  vr( ieig, ieig ) = one
887  GO TO 500
888  END IF
889  END IF
890 *
891 * Clear vector
892 *
893  DO 250 jw = 0, nw - 1
894  DO 240 jr = 1, n
895  work( ( jw+2 )*n+jr ) = zero
896  240 CONTINUE
897  250 CONTINUE
898 *
899 * Compute coefficients in ( a A - b B ) x = 0
900 * a is ACOEF
901 * b is BCOEFR + i*BCOEFI
902 *
903  IF( .NOT.ilcplx ) THEN
904 *
905 * Real eigenvalue
906 *
907  temp = one / max( abs( s( je, je ) )*ascale,
908  $ abs( p( je, je ) )*bscale, safmin )
909  salfar = ( temp*s( je, je ) )*ascale
910  sbeta = ( temp*p( je, je ) )*bscale
911  acoef = sbeta*ascale
912  bcoefr = salfar*bscale
913  bcoefi = zero
914 *
915 * Scale to avoid underflow
916 *
917  scale = one
918  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
919  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
920  $ small
921  IF( lsa )
922  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
923  IF( lsb )
924  $ scale = max( scale, ( small / abs( salfar ) )*
925  $ min( bnorm, big ) )
926  IF( lsa .OR. lsb ) THEN
927  scale = min( scale, one /
928  $ ( safmin*max( one, abs( acoef ),
929  $ abs( bcoefr ) ) ) )
930  IF( lsa ) THEN
931  acoef = ascale*( scale*sbeta )
932  ELSE
933  acoef = scale*acoef
934  END IF
935  IF( lsb ) THEN
936  bcoefr = bscale*( scale*salfar )
937  ELSE
938  bcoefr = scale*bcoefr
939  END IF
940  END IF
941  acoefa = abs( acoef )
942  bcoefa = abs( bcoefr )
943 *
944 * First component is 1
945 *
946  work( 2*n+je ) = one
947  xmax = one
948 *
949 * Compute contribution from column JE of A and B to sum
950 * (See "Further Details", above.)
951 *
952  DO 260 jr = 1, je - 1
953  work( 2*n+jr ) = bcoefr*p( jr, je ) -
954  $ acoef*s( jr, je )
955  260 CONTINUE
956  ELSE
957 *
958 * Complex eigenvalue
959 *
960  CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
961  $ safmin*safety, acoef, temp, bcoefr, temp2,
962  $ bcoefi )
963  IF( bcoefi.EQ.zero ) THEN
964  info = je - 1
965  RETURN
966  END IF
967 *
968 * Scale to avoid over/underflow
969 *
970  acoefa = abs( acoef )
971  bcoefa = abs( bcoefr ) + abs( bcoefi )
972  scale = one
973  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
974  $ scale = ( safmin / ulp ) / acoefa
975  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
976  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
977  IF( safmin*acoefa.GT.ascale )
978  $ scale = ascale / ( safmin*acoefa )
979  IF( safmin*bcoefa.GT.bscale )
980  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
981  IF( scale.NE.one ) THEN
982  acoef = scale*acoef
983  acoefa = abs( acoef )
984  bcoefr = scale*bcoefr
985  bcoefi = scale*bcoefi
986  bcoefa = abs( bcoefr ) + abs( bcoefi )
987  END IF
988 *
989 * Compute first two components of eigenvector
990 * and contribution to sums
991 *
992  temp = acoef*s( je, je-1 )
993  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
994  temp2i = -bcoefi*p( je, je )
995  IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
996  work( 2*n+je ) = one
997  work( 3*n+je ) = zero
998  work( 2*n+je-1 ) = -temp2r / temp
999  work( 3*n+je-1 ) = -temp2i / temp
1000  ELSE
1001  work( 2*n+je-1 ) = one
1002  work( 3*n+je-1 ) = zero
1003  temp = acoef*s( je-1, je )
1004  work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1005  $ s( je-1, je-1 ) ) / temp
1006  work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1007  END IF
1008 *
1009  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1010  $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1011 *
1012 * Compute contribution from columns JE and JE-1
1013 * of A and B to the sums.
1014 *
1015  creala = acoef*work( 2*n+je-1 )
1016  cimaga = acoef*work( 3*n+je-1 )
1017  crealb = bcoefr*work( 2*n+je-1 ) -
1018  $ bcoefi*work( 3*n+je-1 )
1019  cimagb = bcoefi*work( 2*n+je-1 ) +
1020  $ bcoefr*work( 3*n+je-1 )
1021  cre2a = acoef*work( 2*n+je )
1022  cim2a = acoef*work( 3*n+je )
1023  cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1024  cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1025  DO 270 jr = 1, je - 2
1026  work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1027  $ crealb*p( jr, je-1 ) -
1028  $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1029  work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1030  $ cimagb*p( jr, je-1 ) -
1031  $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1032  270 CONTINUE
1033  END IF
1034 *
1035  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1036 *
1037 * Columnwise triangular solve of (a A - b B) x = 0
1038 *
1039  il2by2 = .false.
1040  DO 370 j = je - nw, 1, -1
1041 *
1042 * If a 2-by-2 block, is in position j-1:j, wait until
1043 * next iteration to process it (when it will be j:j+1)
1044 *
1045  IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1046  IF( s( j, j-1 ).NE.zero ) THEN
1047  il2by2 = .true.
1048  GO TO 370
1049  END IF
1050  END IF
1051  bdiag( 1 ) = p( j, j )
1052  IF( il2by2 ) THEN
1053  na = 2
1054  bdiag( 2 ) = p( j+1, j+1 )
1055  ELSE
1056  na = 1
1057  END IF
1058 *
1059 * Compute x(j) (and x(j+1), if 2-by-2 block)
1060 *
1061  CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1062  $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1063  $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1064  $ iinfo )
1065  IF( scale.LT.one ) THEN
1066 *
1067  DO 290 jw = 0, nw - 1
1068  DO 280 jr = 1, je
1069  work( ( jw+2 )*n+jr ) = scale*
1070  $ work( ( jw+2 )*n+jr )
1071  280 CONTINUE
1072  290 CONTINUE
1073  END IF
1074  xmax = max( scale*xmax, temp )
1075 *
1076  DO 310 jw = 1, nw
1077  DO 300 ja = 1, na
1078  work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1079  300 CONTINUE
1080  310 CONTINUE
1081 *
1082 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1083 *
1084  IF( j.GT.1 ) THEN
1085 *
1086 * Check whether scaling is necessary for sum.
1087 *
1088  xscale = one / max( one, xmax )
1089  temp = acoefa*work( j ) + bcoefa*work( n+j )
1090  IF( il2by2 )
1091  $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1092  $ work( n+j+1 ) )
1093  temp = max( temp, acoefa, bcoefa )
1094  IF( temp.GT.bignum*xscale ) THEN
1095 *
1096  DO 330 jw = 0, nw - 1
1097  DO 320 jr = 1, je
1098  work( ( jw+2 )*n+jr ) = xscale*
1099  $ work( ( jw+2 )*n+jr )
1100  320 CONTINUE
1101  330 CONTINUE
1102  xmax = xmax*xscale
1103  END IF
1104 *
1105 * Compute the contributions of the off-diagonals of
1106 * column j (and j+1, if 2-by-2 block) of A and B to the
1107 * sums.
1108 *
1109 *
1110  DO 360 ja = 1, na
1111  IF( ilcplx ) THEN
1112  creala = acoef*work( 2*n+j+ja-1 )
1113  cimaga = acoef*work( 3*n+j+ja-1 )
1114  crealb = bcoefr*work( 2*n+j+ja-1 ) -
1115  $ bcoefi*work( 3*n+j+ja-1 )
1116  cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1117  $ bcoefr*work( 3*n+j+ja-1 )
1118  DO 340 jr = 1, j - 1
1119  work( 2*n+jr ) = work( 2*n+jr ) -
1120  $ creala*s( jr, j+ja-1 ) +
1121  $ crealb*p( jr, j+ja-1 )
1122  work( 3*n+jr ) = work( 3*n+jr ) -
1123  $ cimaga*s( jr, j+ja-1 ) +
1124  $ cimagb*p( jr, j+ja-1 )
1125  340 CONTINUE
1126  ELSE
1127  creala = acoef*work( 2*n+j+ja-1 )
1128  crealb = bcoefr*work( 2*n+j+ja-1 )
1129  DO 350 jr = 1, j - 1
1130  work( 2*n+jr ) = work( 2*n+jr ) -
1131  $ creala*s( jr, j+ja-1 ) +
1132  $ crealb*p( jr, j+ja-1 )
1133  350 CONTINUE
1134  END IF
1135  360 CONTINUE
1136  END IF
1137 *
1138  il2by2 = .false.
1139  370 CONTINUE
1140 *
1141 * Copy eigenvector to VR, back transforming if
1142 * HOWMNY='B'.
1143 *
1144  ieig = ieig - nw
1145  IF( ilback ) THEN
1146 *
1147  DO 410 jw = 0, nw - 1
1148  DO 380 jr = 1, n
1149  work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1150  $ vr( jr, 1 )
1151  380 CONTINUE
1152 *
1153 * A series of compiler directives to defeat
1154 * vectorization for the next loop
1155 *
1156 *
1157  DO 400 jc = 2, je
1158  DO 390 jr = 1, n
1159  work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1160  $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1161  390 CONTINUE
1162  400 CONTINUE
1163  410 CONTINUE
1164 *
1165  DO 430 jw = 0, nw - 1
1166  DO 420 jr = 1, n
1167  vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1168  420 CONTINUE
1169  430 CONTINUE
1170 *
1171  iend = n
1172  ELSE
1173  DO 450 jw = 0, nw - 1
1174  DO 440 jr = 1, n
1175  vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1176  440 CONTINUE
1177  450 CONTINUE
1178 *
1179  iend = je
1180  END IF
1181 *
1182 * Scale eigenvector
1183 *
1184  xmax = zero
1185  IF( ilcplx ) THEN
1186  DO 460 j = 1, iend
1187  xmax = max( xmax, abs( vr( j, ieig ) )+
1188  $ abs( vr( j, ieig+1 ) ) )
1189  460 CONTINUE
1190  ELSE
1191  DO 470 j = 1, iend
1192  xmax = max( xmax, abs( vr( j, ieig ) ) )
1193  470 CONTINUE
1194  END IF
1195 *
1196  IF( xmax.GT.safmin ) THEN
1197  xscale = one / xmax
1198  DO 490 jw = 0, nw - 1
1199  DO 480 jr = 1, iend
1200  vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1201  480 CONTINUE
1202  490 CONTINUE
1203  END IF
1204  500 CONTINUE
1205  END IF
1206 *
1207  RETURN
1208 *
1209 * End of DTGEVC
1210 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
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 dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: dlag2.f:158
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: