Return to dggev.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

    1:       SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
    2:      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
    3: *
    4: *  -- LAPACK driver routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       CHARACTER          JOBVL, JOBVR
   11:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
   15:      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
   16:      $                   VR( LDVR, * ), WORK( * )
   17: *     ..
   18: *
   19: *  Purpose
   20: *  =======
   21: *
   22: *  DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)
   23: *  the generalized eigenvalues, and optionally, the left and/or right
   24: *  generalized eigenvectors.
   25: *
   26: *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar
   27: *  lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
   28: *  singular. It is usually represented as the pair (alpha,beta), as
   29: *  there is a reasonable interpretation for beta=0, and even for both
   30: *  being zero.
   31: *
   32: *  The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
   33: *  of (A,B) satisfies
   34: *
   35: *                   A * v(j) = lambda(j) * B * v(j).
   36: *
   37: *  The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
   38: *  of (A,B) satisfies
   39: *
   40: *                   u(j)**H * A  = lambda(j) * u(j)**H * B .
   41: *
   42: *  where u(j)**H is the conjugate-transpose of u(j).
   43: *
   44: *
   45: *  Arguments
   46: *  =========
   47: *
   48: *  JOBVL   (input) CHARACTER*1
   49: *          = 'N':  do not compute the left generalized eigenvectors;
   50: *          = 'V':  compute the left generalized eigenvectors.
   51: *
   52: *  JOBVR   (input) CHARACTER*1
   53: *          = 'N':  do not compute the right generalized eigenvectors;
   54: *          = 'V':  compute the right generalized eigenvectors.
   55: *
   56: *  N       (input) INTEGER
   57: *          The order of the matrices A, B, VL, and VR.  N >= 0.
   58: *
   59: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
   60: *          On entry, the matrix A in the pair (A,B).
   61: *          On exit, A has been overwritten.
   62: *
   63: *  LDA     (input) INTEGER
   64: *          The leading dimension of A.  LDA >= max(1,N).
   65: *
   66: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
   67: *          On entry, the matrix B in the pair (A,B).
   68: *          On exit, B has been overwritten.
   69: *
   70: *  LDB     (input) INTEGER
   71: *          The leading dimension of B.  LDB >= max(1,N).
   72: *
   73: *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
   74: *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
   75: *  BETA    (output) DOUBLE PRECISION array, dimension (N)
   76: *          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
   77: *          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
   78: *          the j-th eigenvalue is real; if positive, then the j-th and
   79: *          (j+1)-st eigenvalues are a complex conjugate pair, with
   80: *          ALPHAI(j+1) negative.
   81: *
   82: *          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
   83: *          may easily over- or underflow, and BETA(j) may even be zero.
   84: *          Thus, the user should avoid naively computing the ratio
   85: *          alpha/beta.  However, ALPHAR and ALPHAI will be always less
   86: *          than and usually comparable with norm(A) in magnitude, and
   87: *          BETA always less than and usually comparable with norm(B).
   88: *
   89: *  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
   90: *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
   91: *          after another in the columns of VL, in the same order as
   92: *          their eigenvalues. If the j-th eigenvalue is real, then
   93: *          u(j) = VL(:,j), the j-th column of VL. If the j-th and
   94: *          (j+1)-th eigenvalues form a complex conjugate pair, then
   95: *          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
   96: *          Each eigenvector is scaled so the largest component has
   97: *          abs(real part)+abs(imag. part)=1.
   98: *          Not referenced if JOBVL = 'N'.
   99: *
  100: *  LDVL    (input) INTEGER
  101: *          The leading dimension of the matrix VL. LDVL >= 1, and
  102: *          if JOBVL = 'V', LDVL >= N.
  103: *
  104: *  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
  105: *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
  106: *          after another in the columns of VR, in the same order as
  107: *          their eigenvalues. If the j-th eigenvalue is real, then
  108: *          v(j) = VR(:,j), the j-th column of VR. If the j-th and
  109: *          (j+1)-th eigenvalues form a complex conjugate pair, then
  110: *          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
  111: *          Each eigenvector is scaled so the largest component has
  112: *          abs(real part)+abs(imag. part)=1.
  113: *          Not referenced if JOBVR = 'N'.
  114: *
  115: *  LDVR    (input) INTEGER
  116: *          The leading dimension of the matrix VR. LDVR >= 1, and
  117: *          if JOBVR = 'V', LDVR >= N.
  118: *
  119: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  120: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  121: *
  122: *  LWORK   (input) INTEGER
  123: *          The dimension of the array WORK.  LWORK >= max(1,8*N).
  124: *          For good performance, LWORK must generally be larger.
  125: *
  126: *          If LWORK = -1, then a workspace query is assumed; the routine
  127: *          only calculates the optimal size of the WORK array, returns
  128: *          this value as the first entry of the WORK array, and no error
  129: *          message related to LWORK is issued by XERBLA.
  130: *
  131: *  INFO    (output) INTEGER
  132: *          = 0:  successful exit
  133: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  134: *          = 1,...,N:
  135: *                The QZ iteration failed.  No eigenvectors have been
  136: *                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
  137: *                should be correct for j=INFO+1,...,N.
  138: *          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
  139: *                =N+2: error return from DTGEVC.
  140: *
  141: *  =====================================================================
  142: *
  143: *     .. Parameters ..
  144:       DOUBLE PRECISION   ZERO, ONE
  145:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  146: *     ..
  147: *     .. Local Scalars ..
  148:       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
  149:       CHARACTER          CHTEMP
  150:       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
  151:      $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
  152:      $                   MINWRK
  153:       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
  154:      $                   SMLNUM, TEMP
  155: *     ..
  156: *     .. Local Arrays ..
  157:       LOGICAL            LDUMMA( 1 )
  158: *     ..
  159: *     .. External Subroutines ..
  160:       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
  161:      $                   DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
  162:      $                   XERBLA
  163: *     ..
  164: *     .. External Functions ..
  165:       LOGICAL            LSAME
  166:       INTEGER            ILAENV
  167:       DOUBLE PRECISION   DLAMCH, DLANGE
  168:       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
  169: *     ..
  170: *     .. Intrinsic Functions ..
  171:       INTRINSIC          ABS, MAX, SQRT
  172: *     ..
  173: *     .. Executable Statements ..
  174: *
  175: *     Decode the input arguments
  176: *
  177:       IF( LSAME( JOBVL, 'N' ) ) THEN
  178:          IJOBVL = 1
  179:          ILVL = .FALSE.
  180:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
  181:          IJOBVL = 2
  182:          ILVL = .TRUE.
  183:       ELSE
  184:          IJOBVL = -1
  185:          ILVL = .FALSE.
  186:       END IF
  187: *
  188:       IF( LSAME( JOBVR, 'N' ) ) THEN
  189:          IJOBVR = 1
  190:          ILVR = .FALSE.
  191:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
  192:          IJOBVR = 2
  193:          ILVR = .TRUE.
  194:       ELSE
  195:          IJOBVR = -1
  196:          ILVR = .FALSE.
  197:       END IF
  198:       ILV = ILVL .OR. ILVR
  199: *
  200: *     Test the input arguments
  201: *
  202:       INFO = 0
  203:       LQUERY = ( LWORK.EQ.-1 )
  204:       IF( IJOBVL.LE.0 ) THEN
  205:          INFO = -1
  206:       ELSE IF( IJOBVR.LE.0 ) THEN
  207:          INFO = -2
  208:       ELSE IF( N.LT.0 ) THEN
  209:          INFO = -3
  210:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  211:          INFO = -5
  212:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  213:          INFO = -7
  214:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
  215:          INFO = -12
  216:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
  217:          INFO = -14
  218:       END IF
  219: *
  220: *     Compute workspace
  221: *      (Note: Comments in the code beginning "Workspace:" describe the
  222: *       minimal amount of workspace needed at that point in the code,
  223: *       as well as the preferred amount for good performance.
  224: *       NB refers to the optimal block size for the immediately
  225: *       following subroutine, as returned by ILAENV. The workspace is
  226: *       computed assuming ILO = 1 and IHI = N, the worst case.)
  227: *
  228:       IF( INFO.EQ.0 ) THEN
  229:          MINWRK = MAX( 1, 8*N )
  230:          MAXWRK = MAX( 1, N*( 7 +
  231:      $                 ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) )
  232:          MAXWRK = MAX( MAXWRK, N*( 7 +
  233:      $                 ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) )
  234:          IF( ILVL ) THEN
  235:             MAXWRK = MAX( MAXWRK, N*( 7 +
  236:      $                 ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) )
  237:          END IF
  238:          WORK( 1 ) = MAXWRK
  239: *
  240:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
  241:      $      INFO = -16
  242:       END IF
  243: *
  244:       IF( INFO.NE.0 ) THEN
  245:          CALL XERBLA( 'DGGEV ', -INFO )
  246:          RETURN
  247:       ELSE IF( LQUERY ) THEN
  248:          RETURN
  249:       END IF
  250: *
  251: *     Quick return if possible
  252: *
  253:       IF( N.EQ.0 )
  254:      $   RETURN
  255: *
  256: *     Get machine constants
  257: *
  258:       EPS = DLAMCH( 'P' )
  259:       SMLNUM = DLAMCH( 'S' )
  260:       BIGNUM = ONE / SMLNUM
  261:       CALL DLABAD( SMLNUM, BIGNUM )
  262:       SMLNUM = SQRT( SMLNUM ) / EPS
  263:       BIGNUM = ONE / SMLNUM
  264: *
  265: *     Scale A if max element outside range [SMLNUM,BIGNUM]
  266: *
  267:       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
  268:       ILASCL = .FALSE.
  269:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  270:          ANRMTO = SMLNUM
  271:          ILASCL = .TRUE.
  272:       ELSE IF( ANRM.GT.BIGNUM ) THEN
  273:          ANRMTO = BIGNUM
  274:          ILASCL = .TRUE.
  275:       END IF
  276:       IF( ILASCL )
  277:      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
  278: *
  279: *     Scale B if max element outside range [SMLNUM,BIGNUM]
  280: *
  281:       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
  282:       ILBSCL = .FALSE.
  283:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  284:          BNRMTO = SMLNUM
  285:          ILBSCL = .TRUE.
  286:       ELSE IF( BNRM.GT.BIGNUM ) THEN
  287:          BNRMTO = BIGNUM
  288:          ILBSCL = .TRUE.
  289:       END IF
  290:       IF( ILBSCL )
  291:      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
  292: *
  293: *     Permute the matrices A, B to isolate eigenvalues if possible
  294: *     (Workspace: need 6*N)
  295: *
  296:       ILEFT = 1
  297:       IRIGHT = N + 1
  298:       IWRK = IRIGHT + N
  299:       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
  300:      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
  301: *
  302: *     Reduce B to triangular form (QR decomposition of B)
  303: *     (Workspace: need N, prefer N*NB)
  304: *
  305:       IROWS = IHI + 1 - ILO
  306:       IF( ILV ) THEN
  307:          ICOLS = N + 1 - ILO
  308:       ELSE
  309:          ICOLS = IROWS
  310:       END IF
  311:       ITAU = IWRK
  312:       IWRK = ITAU + IROWS
  313:       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
  314:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
  315: *
  316: *     Apply the orthogonal transformation to matrix A
  317: *     (Workspace: need N, prefer N*NB)
  318: *
  319:       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
  320:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
  321:      $             LWORK+1-IWRK, IERR )
  322: *
  323: *     Initialize VL
  324: *     (Workspace: need N, prefer N*NB)
  325: *
  326:       IF( ILVL ) THEN
  327:          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
  328:          IF( IROWS.GT.1 ) THEN
  329:             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
  330:      $                   VL( ILO+1, ILO ), LDVL )
  331:          END IF
  332:          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
  333:      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
  334:       END IF
  335: *
  336: *     Initialize VR
  337: *
  338:       IF( ILVR )
  339:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
  340: *
  341: *     Reduce to generalized Hessenberg form
  342: *     (Workspace: none needed)
  343: *
  344:       IF( ILV ) THEN
  345: *
  346: *        Eigenvectors requested -- work on whole matrix.
  347: *
  348:          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
  349:      $                LDVL, VR, LDVR, IERR )
  350:       ELSE
  351:          CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
  352:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
  353:       END IF
  354: *
  355: *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
  356: *     Schur forms and Schur vectors)
  357: *     (Workspace: need N)
  358: *
  359:       IWRK = ITAU
  360:       IF( ILV ) THEN
  361:          CHTEMP = 'S'
  362:       ELSE
  363:          CHTEMP = 'E'
  364:       END IF
  365:       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
  366:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
  367:      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
  368:       IF( IERR.NE.0 ) THEN
  369:          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
  370:             INFO = IERR
  371:          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
  372:             INFO = IERR - N
  373:          ELSE
  374:             INFO = N + 1
  375:          END IF
  376:          GO TO 110
  377:       END IF
  378: *
  379: *     Compute Eigenvectors
  380: *     (Workspace: need 6*N)
  381: *
  382:       IF( ILV ) THEN
  383:          IF( ILVL ) THEN
  384:             IF( ILVR ) THEN
  385:                CHTEMP = 'B'
  386:             ELSE
  387:                CHTEMP = 'L'
  388:             END IF
  389:          ELSE
  390:             CHTEMP = 'R'
  391:          END IF
  392:          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
  393:      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
  394:          IF( IERR.NE.0 ) THEN
  395:             INFO = N + 2
  396:             GO TO 110
  397:          END IF
  398: *
  399: *        Undo balancing on VL and VR and normalization
  400: *        (Workspace: none needed)
  401: *
  402:          IF( ILVL ) THEN
  403:             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
  404:      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
  405:             DO 50 JC = 1, N
  406:                IF( ALPHAI( JC ).LT.ZERO )
  407:      $            GO TO 50
  408:                TEMP = ZERO
  409:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
  410:                   DO 10 JR = 1, N
  411:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
  412:    10             CONTINUE
  413:                ELSE
  414:                   DO 20 JR = 1, N
  415:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
  416:      $                      ABS( VL( JR, JC+1 ) ) )
  417:    20             CONTINUE
  418:                END IF
  419:                IF( TEMP.LT.SMLNUM )
  420:      $            GO TO 50
  421:                TEMP = ONE / TEMP
  422:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
  423:                   DO 30 JR = 1, N
  424:                      VL( JR, JC ) = VL( JR, JC )*TEMP
  425:    30             CONTINUE
  426:                ELSE
  427:                   DO 40 JR = 1, N
  428:                      VL( JR, JC ) = VL( JR, JC )*TEMP
  429:                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
  430:    40             CONTINUE
  431:                END IF
  432:    50       CONTINUE
  433:          END IF
  434:          IF( ILVR ) THEN
  435:             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
  436:      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
  437:             DO 100 JC = 1, N
  438:                IF( ALPHAI( JC ).LT.ZERO )
  439:      $            GO TO 100
  440:                TEMP = ZERO
  441:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
  442:                   DO 60 JR = 1, N
  443:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
  444:    60             CONTINUE
  445:                ELSE
  446:                   DO 70 JR = 1, N
  447:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
  448:      $                      ABS( VR( JR, JC+1 ) ) )
  449:    70             CONTINUE
  450:                END IF
  451:                IF( TEMP.LT.SMLNUM )
  452:      $            GO TO 100
  453:                TEMP = ONE / TEMP
  454:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
  455:                   DO 80 JR = 1, N
  456:                      VR( JR, JC ) = VR( JR, JC )*TEMP
  457:    80             CONTINUE
  458:                ELSE
  459:                   DO 90 JR = 1, N
  460:                      VR( JR, JC ) = VR( JR, JC )*TEMP
  461:                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
  462:    90             CONTINUE
  463:                END IF
  464:   100       CONTINUE
  465:          END IF
  466: *
  467: *        End of eigenvector calculation
  468: *
  469:       END IF
  470: *
  471: *     Undo scaling if necessary
  472: *
  473:       IF( ILASCL ) THEN
  474:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
  475:          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
  476:       END IF
  477: *
  478:       IF( ILBSCL ) THEN
  479:          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
  480:       END IF
  481: *
  482:   110 CONTINUE
  483: *
  484:       WORK( 1 ) = MAXWRK
  485: *
  486:       RETURN
  487: *
  488: *     End of DGGEV
  489: *
  490:       END