LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 
)

DGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B).

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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
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
subroutine dgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
DGGHRD
Definition dgghrd.f:207
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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 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: