LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ dgegv()

subroutine dgegv ( character  JOBVL,
character  JOBVR,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
 This routine is deprecated and has been replaced by routine DGGEV.

 DGEGV computes the eigenvalues and, optionally, the left and/or right
 eigenvectors of a real matrix pair (A,B).
 Given two square matrices A and B,
 the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
 eigenvalues lambda and corresponding (non-zero) eigenvectors x such
 that

    A*x = lambda*B*x.

 An alternate form is to find the eigenvalues mu and corresponding
 eigenvectors y such that

    mu*A*y = B*y.

 These two forms are equivalent with mu = 1/lambda and x = y if
 neither lambda nor mu is zero.  In order to deal with the case that
 lambda or mu is zero or small, two values alpha and beta are returned
 for each eigenvalue, such that lambda = alpha/beta and
 mu = beta/alpha.

 The vectors x and y in the above equations are right eigenvectors of
 the matrix pair (A,B).  Vectors u and v satisfying

    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B

 are left eigenvectors of (A,B).

 Note: this routine performs "full balancing" on A and B
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors (returned
                  in VL).
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors (returned
                  in VR).
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the matrix A.
          If JOBVL = 'V' or JOBVR = 'V', then on exit A
          contains the real Schur form of A from the generalized Schur
          factorization of the pair (A,B) after balancing.
          If no eigenvectors were computed, then only the diagonal
          blocks from the Schur form will be correct.  See DGGHRD and
          DHGEQZ for details.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB, N)
          On entry, the matrix B.
          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
          upper triangular matrix obtained from B in the generalized
          Schur factorization of the pair (A,B) after balancing.
          If no eigenvectors were computed, then only those elements of
          B corresponding to the diagonal blocks from the Schur form of
          A will be correct.  See DGGHRD and DHGEQZ for details.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
          The real parts of each scalar alpha defining an eigenvalue of
          GNEP.
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
          The imaginary parts of each scalar alpha defining an
          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
          eigenvalue is real; if positive, then the j-th and
          (j+1)-st eigenvalues are a complex conjugate pair, with
          ALPHAI(j+1) = -ALPHAI(j).
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          The scalars beta that define the eigenvalues of GNEP.

          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
          beta = BETA(j) represent the j-th eigenvalue of the matrix
          pair (A,B), in one of the forms lambda = alpha/beta or
          mu = beta/alpha.  Since either lambda or mu may overflow,
          they should not, in general, be computed.
[out]VL
          VL is DOUBLE PRECISION array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored
          in the columns of VL, in the same order as their eigenvalues.
          If the j-th eigenvalue is real, then u(j) = VL(:,j).
          If the j-th and (j+1)-st eigenvalues form a complex conjugate
          pair, then
             u(j) = VL(:,j) + i*VL(:,j+1)
          and
            u(j+1) = VL(:,j) - i*VL(:,j+1).

          Each eigenvector is scaled so that its largest component has
          abs(real part) + abs(imag. part) = 1, except for eigenvectors
          corresponding to an eigenvalue with alpha = beta = 0, which
          are set to zero.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is DOUBLE PRECISION array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors x(j) are stored
          in the columns of VR, in the same order as their eigenvalues.
          If the j-th eigenvalue is real, then x(j) = VR(:,j).
          If the j-th and (j+1)-st eigenvalues form a complex conjugate
          pair, then
            x(j) = VR(:,j) + i*VR(:,j+1)
          and
            x(j+1) = VR(:,j) - i*VR(:,j+1).

          Each eigenvector is scaled so that its largest component has
          abs(real part) + abs(imag. part) = 1, except for eigenvalues
          corresponding to an eigenvalue with alpha = beta = 0, which
          are set to zero.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,8*N).
          For good performance, LWORK must generally be larger.
          To compute the optimal value of LWORK, call ILAENV to get
          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
          The optimal LWORK is:
              2*N + MAX( 6*N, N*(NB+1) ).

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          = 1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
                should be correct for j=INFO+1,...,N.
          > N:  errors that usually indicate LAPACK problems:
                =N+1: error return from DGGBAL
                =N+2: error return from DGEQRF
                =N+3: error return from DORMQR
                =N+4: error return from DORGQR
                =N+5: error return from DGGHRD
                =N+6: error return from DHGEQZ (other than failed
                                                iteration)
                =N+7: error return from DTGEVC
                =N+8: error return from DGGBAK (computing VL)
                =N+9: error return from DGGBAK (computing VR)
                =N+10: error return from DLASCL (various calls)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  Balancing
  ---------

  This driver calls DGGBAL to both permute and scale rows and columns
  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
  and PL*B*R will be upper triangular except for the diagonal blocks
  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
  possible.  The diagonal scaling matrices DL and DR are chosen so
  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
  one (except for the elements that start out zero.)

  After the eigenvalues and eigenvectors of the balanced matrices
  have been computed, DGGBAK transforms the eigenvectors back to what
  they would have been (in perfect arithmetic) if they had not been
  balanced.

  Contents of A and B on Exit
  -------- -- - --- - -- ----

  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
  both), then on exit the arrays A and B will contain the real Schur
  form[*] of the "balanced" versions of A and B.  If no eigenvectors
  are computed, then only the diagonal blocks will be correct.

  [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
      by Golub & van Loan, pub. by Johns Hopkins U. Press.

Definition at line 304 of file dgegv.f.

306 *
307 * -- LAPACK driver routine --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310 *
311 * .. Scalar Arguments ..
312  CHARACTER JOBVL, JOBVR
313  INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314 * ..
315 * .. Array Arguments ..
316  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317  $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
318  $ VR( LDVR, * ), WORK( * )
319 * ..
320 *
321 * =====================================================================
322 *
323 * .. Parameters ..
324  DOUBLE PRECISION ZERO, ONE
325  parameter( zero = 0.0d0, one = 1.0d0 )
326 * ..
327 * .. Local Scalars ..
328  LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329  CHARACTER CHTEMP
330  INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331  $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
332  $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
333  DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334  $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
335  $ SALFAI, SALFAR, SBETA, SCALE, TEMP
336 * ..
337 * .. Local Arrays ..
338  LOGICAL LDUMMA( 1 )
339 * ..
340 * .. External Subroutines ..
341  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
343 * ..
344 * .. External Functions ..
345  LOGICAL LSAME
346  INTEGER ILAENV
347  DOUBLE PRECISION DLAMCH, DLANGE
348  EXTERNAL lsame, ilaenv, dlamch, dlange
349 * ..
350 * .. Intrinsic Functions ..
351  INTRINSIC abs, int, max
352 * ..
353 * .. Executable Statements ..
354 *
355 * Decode the input arguments
356 *
357  IF( lsame( jobvl, 'N' ) ) THEN
358  ijobvl = 1
359  ilvl = .false.
360  ELSE IF( lsame( jobvl, 'V' ) ) THEN
361  ijobvl = 2
362  ilvl = .true.
363  ELSE
364  ijobvl = -1
365  ilvl = .false.
366  END IF
367 *
368  IF( lsame( jobvr, 'N' ) ) THEN
369  ijobvr = 1
370  ilvr = .false.
371  ELSE IF( lsame( jobvr, 'V' ) ) THEN
372  ijobvr = 2
373  ilvr = .true.
374  ELSE
375  ijobvr = -1
376  ilvr = .false.
377  END IF
378  ilv = ilvl .OR. ilvr
379 *
380 * Test the input arguments
381 *
382  lwkmin = max( 8*n, 1 )
383  lwkopt = lwkmin
384  work( 1 ) = lwkopt
385  lquery = ( lwork.EQ.-1 )
386  info = 0
387  IF( ijobvl.LE.0 ) THEN
388  info = -1
389  ELSE IF( ijobvr.LE.0 ) THEN
390  info = -2
391  ELSE IF( n.LT.0 ) THEN
392  info = -3
393  ELSE IF( lda.LT.max( 1, n ) ) THEN
394  info = -5
395  ELSE IF( ldb.LT.max( 1, n ) ) THEN
396  info = -7
397  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
398  info = -12
399  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
400  info = -14
401  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
402  info = -16
403  END IF
404 *
405  IF( info.EQ.0 ) THEN
406  nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
407  nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
408  nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
409  nb = max( nb1, nb2, nb3 )
410  lopt = 2*n + max( 6*n, n*( nb+1 ) )
411  work( 1 ) = lopt
412  END IF
413 *
414  IF( info.NE.0 ) THEN
415  CALL xerbla( 'DGEGV ', -info )
416  RETURN
417  ELSE IF( lquery ) THEN
418  RETURN
419  END IF
420 *
421 * Quick return if possible
422 *
423  IF( n.EQ.0 )
424  $ RETURN
425 *
426 * Get machine constants
427 *
428  eps = dlamch( 'E' )*dlamch( 'B' )
429  safmin = dlamch( 'S' )
430  safmin = safmin + safmin
431  safmax = one / safmin
432  onepls = one + ( 4*eps )
433 *
434 * Scale A
435 *
436  anrm = dlange( 'M', n, n, a, lda, work )
437  anrm1 = anrm
438  anrm2 = one
439  IF( anrm.LT.one ) THEN
440  IF( safmax*anrm.LT.one ) THEN
441  anrm1 = safmin
442  anrm2 = safmax*anrm
443  END IF
444  END IF
445 *
446  IF( anrm.GT.zero ) THEN
447  CALL dlascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
448  IF( iinfo.NE.0 ) THEN
449  info = n + 10
450  RETURN
451  END IF
452  END IF
453 *
454 * Scale B
455 *
456  bnrm = dlange( 'M', n, n, b, ldb, work )
457  bnrm1 = bnrm
458  bnrm2 = one
459  IF( bnrm.LT.one ) THEN
460  IF( safmax*bnrm.LT.one ) THEN
461  bnrm1 = safmin
462  bnrm2 = safmax*bnrm
463  END IF
464  END IF
465 *
466  IF( bnrm.GT.zero ) THEN
467  CALL dlascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
468  IF( iinfo.NE.0 ) THEN
469  info = n + 10
470  RETURN
471  END IF
472  END IF
473 *
474 * Permute the matrix to make it more nearly triangular
475 * Workspace layout: (8*N words -- "work" requires 6*N words)
476 * left_permutation, right_permutation, work...
477 *
478  ileft = 1
479  iright = n + 1
480  iwork = iright + n
481  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
482  $ work( iright ), work( iwork ), iinfo )
483  IF( iinfo.NE.0 ) THEN
484  info = n + 1
485  GO TO 120
486  END IF
487 *
488 * Reduce B to triangular form, and initialize VL and/or VR
489 * Workspace layout: ("work..." must have at least N words)
490 * left_permutation, right_permutation, tau, work...
491 *
492  irows = ihi + 1 - ilo
493  IF( ilv ) THEN
494  icols = n + 1 - ilo
495  ELSE
496  icols = irows
497  END IF
498  itau = iwork
499  iwork = itau + irows
500  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
501  $ work( iwork ), lwork+1-iwork, iinfo )
502  IF( iinfo.GE.0 )
503  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
504  IF( iinfo.NE.0 ) THEN
505  info = n + 2
506  GO TO 120
507  END IF
508 *
509  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
510  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
511  $ lwork+1-iwork, iinfo )
512  IF( iinfo.GE.0 )
513  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514  IF( iinfo.NE.0 ) THEN
515  info = n + 3
516  GO TO 120
517  END IF
518 *
519  IF( ilvl ) THEN
520  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
521  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
522  $ vl( ilo+1, ilo ), ldvl )
523  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
524  $ work( itau ), work( iwork ), lwork+1-iwork,
525  $ iinfo )
526  IF( iinfo.GE.0 )
527  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
528  IF( iinfo.NE.0 ) THEN
529  info = n + 4
530  GO TO 120
531  END IF
532  END IF
533 *
534  IF( ilvr )
535  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
536 *
537 * Reduce to generalized Hessenberg form
538 *
539  IF( ilv ) THEN
540 *
541 * Eigenvectors requested -- work on whole matrix.
542 *
543  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
544  $ ldvl, vr, ldvr, iinfo )
545  ELSE
546  CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
547  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548  END IF
549  IF( iinfo.NE.0 ) THEN
550  info = n + 5
551  GO TO 120
552  END IF
553 *
554 * Perform QZ algorithm
555 * Workspace layout: ("work..." must have at least 1 word)
556 * left_permutation, right_permutation, work...
557 *
558  iwork = itau
559  IF( ilv ) THEN
560  chtemp = 'S'
561  ELSE
562  chtemp = 'E'
563  END IF
564  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
565  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
566  $ work( iwork ), lwork+1-iwork, iinfo )
567  IF( iinfo.GE.0 )
568  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
569  IF( iinfo.NE.0 ) THEN
570  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
571  info = iinfo
572  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
573  info = iinfo - n
574  ELSE
575  info = n + 6
576  END IF
577  GO TO 120
578  END IF
579 *
580  IF( ilv ) THEN
581 *
582 * Compute Eigenvectors (DTGEVC requires 6*N words of workspace)
583 *
584  IF( ilvl ) THEN
585  IF( ilvr ) THEN
586  chtemp = 'B'
587  ELSE
588  chtemp = 'L'
589  END IF
590  ELSE
591  chtemp = 'R'
592  END IF
593 *
594  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
595  $ vr, ldvr, n, in, work( iwork ), iinfo )
596  IF( iinfo.NE.0 ) THEN
597  info = n + 7
598  GO TO 120
599  END IF
600 *
601 * Undo balancing on VL and VR, rescale
602 *
603  IF( ilvl ) THEN
604  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
605  $ work( iright ), n, vl, ldvl, iinfo )
606  IF( iinfo.NE.0 ) THEN
607  info = n + 8
608  GO TO 120
609  END IF
610  DO 50 jc = 1, n
611  IF( alphai( jc ).LT.zero )
612  $ GO TO 50
613  temp = zero
614  IF( alphai( jc ).EQ.zero ) THEN
615  DO 10 jr = 1, n
616  temp = max( temp, abs( vl( jr, jc ) ) )
617  10 CONTINUE
618  ELSE
619  DO 20 jr = 1, n
620  temp = max( temp, abs( vl( jr, jc ) )+
621  $ abs( vl( jr, jc+1 ) ) )
622  20 CONTINUE
623  END IF
624  IF( temp.LT.safmin )
625  $ GO TO 50
626  temp = one / temp
627  IF( alphai( jc ).EQ.zero ) THEN
628  DO 30 jr = 1, n
629  vl( jr, jc ) = vl( jr, jc )*temp
630  30 CONTINUE
631  ELSE
632  DO 40 jr = 1, n
633  vl( jr, jc ) = vl( jr, jc )*temp
634  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
635  40 CONTINUE
636  END IF
637  50 CONTINUE
638  END IF
639  IF( ilvr ) THEN
640  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
641  $ work( iright ), n, vr, ldvr, iinfo )
642  IF( iinfo.NE.0 ) THEN
643  info = n + 9
644  GO TO 120
645  END IF
646  DO 100 jc = 1, n
647  IF( alphai( jc ).LT.zero )
648  $ GO TO 100
649  temp = zero
650  IF( alphai( jc ).EQ.zero ) THEN
651  DO 60 jr = 1, n
652  temp = max( temp, abs( vr( jr, jc ) ) )
653  60 CONTINUE
654  ELSE
655  DO 70 jr = 1, n
656  temp = max( temp, abs( vr( jr, jc ) )+
657  $ abs( vr( jr, jc+1 ) ) )
658  70 CONTINUE
659  END IF
660  IF( temp.LT.safmin )
661  $ GO TO 100
662  temp = one / temp
663  IF( alphai( jc ).EQ.zero ) THEN
664  DO 80 jr = 1, n
665  vr( jr, jc ) = vr( jr, jc )*temp
666  80 CONTINUE
667  ELSE
668  DO 90 jr = 1, n
669  vr( jr, jc ) = vr( jr, jc )*temp
670  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
671  90 CONTINUE
672  END IF
673  100 CONTINUE
674  END IF
675 *
676 * End of eigenvector calculation
677 *
678  END IF
679 *
680 * Undo scaling in alpha, beta
681 *
682 * Note: this does not give the alpha and beta for the unscaled
683 * problem.
684 *
685 * Un-scaling is limited to avoid underflow in alpha and beta
686 * if they are significant.
687 *
688  DO 110 jc = 1, n
689  absar = abs( alphar( jc ) )
690  absai = abs( alphai( jc ) )
691  absb = abs( beta( jc ) )
692  salfar = anrm*alphar( jc )
693  salfai = anrm*alphai( jc )
694  sbeta = bnrm*beta( jc )
695  ilimit = .false.
696  scale = one
697 *
698 * Check for significant underflow in ALPHAI
699 *
700  IF( abs( salfai ).LT.safmin .AND. absai.GE.
701  $ max( safmin, eps*absar, eps*absb ) ) THEN
702  ilimit = .true.
703  scale = ( onepls*safmin / anrm1 ) /
704  $ max( onepls*safmin, anrm2*absai )
705 *
706  ELSE IF( salfai.EQ.zero ) THEN
707 *
708 * If insignificant underflow in ALPHAI, then make the
709 * conjugate eigenvalue real.
710 *
711  IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
712  alphai( jc-1 ) = zero
713  ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
714  alphai( jc+1 ) = zero
715  END IF
716  END IF
717 *
718 * Check for significant underflow in ALPHAR
719 *
720  IF( abs( salfar ).LT.safmin .AND. absar.GE.
721  $ max( safmin, eps*absai, eps*absb ) ) THEN
722  ilimit = .true.
723  scale = max( scale, ( onepls*safmin / anrm1 ) /
724  $ max( onepls*safmin, anrm2*absar ) )
725  END IF
726 *
727 * Check for significant underflow in BETA
728 *
729  IF( abs( sbeta ).LT.safmin .AND. absb.GE.
730  $ max( safmin, eps*absar, eps*absai ) ) THEN
731  ilimit = .true.
732  scale = max( scale, ( onepls*safmin / bnrm1 ) /
733  $ max( onepls*safmin, bnrm2*absb ) )
734  END IF
735 *
736 * Check for possible overflow when limiting scaling
737 *
738  IF( ilimit ) THEN
739  temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
740  $ abs( sbeta ) )
741  IF( temp.GT.one )
742  $ scale = scale / temp
743  IF( scale.LT.one )
744  $ ilimit = .false.
745  END IF
746 *
747 * Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748 *
749  IF( ilimit ) THEN
750  salfar = ( scale*alphar( jc ) )*anrm
751  salfai = ( scale*alphai( jc ) )*anrm
752  sbeta = ( scale*beta( jc ) )*bnrm
753  END IF
754  alphar( jc ) = salfar
755  alphai( jc ) = salfai
756  beta( jc ) = sbeta
757  110 CONTINUE
758 *
759  120 CONTINUE
760  work( 1 ) = lwkopt
761 *
762  RETURN
763 *
764 * End of DGEGV
765 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
Definition: dggbak.f:147
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
Definition: dggbal.f:177
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:304
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:295
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:128
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
Definition: dgghrd.f:207
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
Here is the call graph for this function: