Return to zgeesx.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief <b> ZGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZGEESX + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeesx.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeesx.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeesx.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
   22: *                          VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
   23: *                          BWORK, INFO )
   24: * 
   25: *       .. Scalar Arguments ..
   26: *       CHARACTER          JOBVS, SENSE, SORT
   27: *       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
   28: *       DOUBLE PRECISION   RCONDE, RCONDV
   29: *       ..
   30: *       .. Array Arguments ..
   31: *       LOGICAL            BWORK( * )
   32: *       DOUBLE PRECISION   RWORK( * )
   33: *       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
   34: *       ..
   35: *       .. Function Arguments ..
   36: *       LOGICAL            SELECT
   37: *       EXTERNAL           SELECT
   38: *       ..
   39: *  
   40: *
   41: *> \par Purpose:
   42: *  =============
   43: *>
   44: *> \verbatim
   45: *>
   46: *> ZGEESX computes for an N-by-N complex nonsymmetric matrix A, the
   47: *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
   48: *> vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
   49: *>
   50: *> Optionally, it also orders the eigenvalues on the diagonal of the
   51: *> Schur form so that selected eigenvalues are at the top left;
   52: *> computes a reciprocal condition number for the average of the
   53: *> selected eigenvalues (RCONDE); and computes a reciprocal condition
   54: *> number for the right invariant subspace corresponding to the
   55: *> selected eigenvalues (RCONDV).  The leading columns of Z form an
   56: *> orthonormal basis for this invariant subspace.
   57: *>
   58: *> For further explanation of the reciprocal condition numbers RCONDE
   59: *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
   60: *> these quantities are called s and sep respectively).
   61: *>
   62: *> A complex matrix is in Schur form if it is upper triangular.
   63: *> \endverbatim
   64: *
   65: *  Arguments:
   66: *  ==========
   67: *
   68: *> \param[in] JOBVS
   69: *> \verbatim
   70: *>          JOBVS is CHARACTER*1
   71: *>          = 'N': Schur vectors are not computed;
   72: *>          = 'V': Schur vectors are computed.
   73: *> \endverbatim
   74: *>
   75: *> \param[in] SORT
   76: *> \verbatim
   77: *>          SORT is CHARACTER*1
   78: *>          Specifies whether or not to order the eigenvalues on the
   79: *>          diagonal of the Schur form.
   80: *>          = 'N': Eigenvalues are not ordered;
   81: *>          = 'S': Eigenvalues are ordered (see SELECT).
   82: *> \endverbatim
   83: *>
   84: *> \param[in] SELECT
   85: *> \verbatim
   86: *>          SELECT is procedure) LOGICAL FUNCTION of one COMPLEX*16 argument
   87: *>          SELECT must be declared EXTERNAL in the calling subroutine.
   88: *>          If SORT = 'S', SELECT is used to select eigenvalues to order
   89: *>          to the top left of the Schur form.
   90: *>          If SORT = 'N', SELECT is not referenced.
   91: *>          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
   92: *> \endverbatim
   93: *>
   94: *> \param[in] SENSE
   95: *> \verbatim
   96: *>          SENSE is CHARACTER*1
   97: *>          Determines which reciprocal condition numbers are computed.
   98: *>          = 'N': None are computed;
   99: *>          = 'E': Computed for average of selected eigenvalues only;
  100: *>          = 'V': Computed for selected right invariant subspace only;
  101: *>          = 'B': Computed for both.
  102: *>          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
  103: *> \endverbatim
  104: *>
  105: *> \param[in] N
  106: *> \verbatim
  107: *>          N is INTEGER
  108: *>          The order of the matrix A. N >= 0.
  109: *> \endverbatim
  110: *>
  111: *> \param[in,out] A
  112: *> \verbatim
  113: *>          A is COMPLEX*16 array, dimension (LDA, N)
  114: *>          On entry, the N-by-N matrix A.
  115: *>          On exit, A is overwritten by its Schur form T.
  116: *> \endverbatim
  117: *>
  118: *> \param[in] LDA
  119: *> \verbatim
  120: *>          LDA is INTEGER
  121: *>          The leading dimension of the array A.  LDA >= max(1,N).
  122: *> \endverbatim
  123: *>
  124: *> \param[out] SDIM
  125: *> \verbatim
  126: *>          SDIM is INTEGER
  127: *>          If SORT = 'N', SDIM = 0.
  128: *>          If SORT = 'S', SDIM = number of eigenvalues for which
  129: *>                         SELECT is true.
  130: *> \endverbatim
  131: *>
  132: *> \param[out] W
  133: *> \verbatim
  134: *>          W is COMPLEX*16 array, dimension (N)
  135: *>          W contains the computed eigenvalues, in the same order
  136: *>          that they appear on the diagonal of the output Schur form T.
  137: *> \endverbatim
  138: *>
  139: *> \param[out] VS
  140: *> \verbatim
  141: *>          VS is COMPLEX*16 array, dimension (LDVS,N)
  142: *>          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
  143: *>          vectors.
  144: *>          If JOBVS = 'N', VS is not referenced.
  145: *> \endverbatim
  146: *>
  147: *> \param[in] LDVS
  148: *> \verbatim
  149: *>          LDVS is INTEGER
  150: *>          The leading dimension of the array VS.  LDVS >= 1, and if
  151: *>          JOBVS = 'V', LDVS >= N.
  152: *> \endverbatim
  153: *>
  154: *> \param[out] RCONDE
  155: *> \verbatim
  156: *>          RCONDE is DOUBLE PRECISION
  157: *>          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
  158: *>          condition number for the average of the selected eigenvalues.
  159: *>          Not referenced if SENSE = 'N' or 'V'.
  160: *> \endverbatim
  161: *>
  162: *> \param[out] RCONDV
  163: *> \verbatim
  164: *>          RCONDV is DOUBLE PRECISION
  165: *>          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
  166: *>          condition number for the selected right invariant subspace.
  167: *>          Not referenced if SENSE = 'N' or 'E'.
  168: *> \endverbatim
  169: *>
  170: *> \param[out] WORK
  171: *> \verbatim
  172: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
  173: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  174: *> \endverbatim
  175: *>
  176: *> \param[in] LWORK
  177: *> \verbatim
  178: *>          LWORK is INTEGER
  179: *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
  180: *>          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
  181: *>          where SDIM is the number of selected eigenvalues computed by
  182: *>          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
  183: *>          that an error is only returned if LWORK < max(1,2*N), but if
  184: *>          SENSE = 'E' or 'V' or 'B' this may not be large enough.
  185: *>          For good performance, LWORK must generally be larger.
  186: *>
  187: *>          If LWORK = -1, then a workspace query is assumed; the routine
  188: *>          only calculates upper bound on the optimal size of the
  189: *>          array WORK, returns this value as the first entry of the WORK
  190: *>          array, and no error message related to LWORK is issued by
  191: *>          XERBLA.
  192: *> \endverbatim
  193: *>
  194: *> \param[out] RWORK
  195: *> \verbatim
  196: *>          RWORK is DOUBLE PRECISION array, dimension (N)
  197: *> \endverbatim
  198: *>
  199: *> \param[out] BWORK
  200: *> \verbatim
  201: *>          BWORK is LOGICAL array, dimension (N)
  202: *>          Not referenced if SORT = 'N'.
  203: *> \endverbatim
  204: *>
  205: *> \param[out] INFO
  206: *> \verbatim
  207: *>          INFO is INTEGER
  208: *>          = 0: successful exit
  209: *>          < 0: if INFO = -i, the i-th argument had an illegal value.
  210: *>          > 0: if INFO = i, and i is
  211: *>             <= N: the QR algorithm failed to compute all the
  212: *>                   eigenvalues; elements 1:ILO-1 and i+1:N of W
  213: *>                   contain those eigenvalues which have converged; if
  214: *>                   JOBVS = 'V', VS contains the transformation which
  215: *>                   reduces A to its partially converged Schur form.
  216: *>             = N+1: the eigenvalues could not be reordered because some
  217: *>                   eigenvalues were too close to separate (the problem
  218: *>                   is very ill-conditioned);
  219: *>             = N+2: after reordering, roundoff changed values of some
  220: *>                   complex eigenvalues so that leading eigenvalues in
  221: *>                   the Schur form no longer satisfy SELECT=.TRUE.  This
  222: *>                   could also be caused by underflow due to scaling.
  223: *> \endverbatim
  224: *
  225: *  Authors:
  226: *  ========
  227: *
  228: *> \author Univ. of Tennessee 
  229: *> \author Univ. of California Berkeley 
  230: *> \author Univ. of Colorado Denver 
  231: *> \author NAG Ltd. 
  232: *
  233: *> \date November 2011
  234: *
  235: *> \ingroup complex16GEeigen
  236: *
  237: *  =====================================================================
  238:       SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
  239:      $                   VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
  240:      $                   BWORK, INFO )
  241: *
  242: *  -- LAPACK driver routine (version 3.4.0) --
  243: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  244: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  245: *     November 2011
  246: *
  247: *     .. Scalar Arguments ..
  248:       CHARACTER          JOBVS, SENSE, SORT
  249:       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
  250:       DOUBLE PRECISION   RCONDE, RCONDV
  251: *     ..
  252: *     .. Array Arguments ..
  253:       LOGICAL            BWORK( * )
  254:       DOUBLE PRECISION   RWORK( * )
  255:       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
  256: *     ..
  257: *     .. Function Arguments ..
  258:       LOGICAL            SELECT
  259:       EXTERNAL           SELECT
  260: *     ..
  261: *
  262: *  =====================================================================
  263: *
  264: *     .. Parameters ..
  265:       DOUBLE PRECISION   ZERO, ONE
  266:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  267: *     ..
  268: *     .. Local Scalars ..
  269:       LOGICAL            LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
  270:      $                   WANTSV, WANTVS
  271:       INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
  272:      $                   ITAU, IWRK, LWRK, MAXWRK, MINWRK
  273:       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SMLNUM
  274: *     ..
  275: *     .. Local Arrays ..
  276:       DOUBLE PRECISION   DUM( 1 )
  277: *     ..
  278: *     .. External Subroutines ..
  279:       EXTERNAL           DLABAD, DLASCL, XERBLA, ZCOPY, ZGEBAK, ZGEBAL,
  280:      $                   ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
  281: *     ..
  282: *     .. External Functions ..
  283:       LOGICAL            LSAME
  284:       INTEGER            ILAENV
  285:       DOUBLE PRECISION   DLAMCH, ZLANGE
  286:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
  287: *     ..
  288: *     .. Intrinsic Functions ..
  289:       INTRINSIC          MAX, SQRT
  290: *     ..
  291: *     .. Executable Statements ..
  292: *
  293: *     Test the input arguments
  294: *
  295:       INFO = 0
  296:       WANTVS = LSAME( JOBVS, 'V' )
  297:       WANTST = LSAME( SORT, 'S' )
  298:       WANTSN = LSAME( SENSE, 'N' )
  299:       WANTSE = LSAME( SENSE, 'E' )
  300:       WANTSV = LSAME( SENSE, 'V' )
  301:       WANTSB = LSAME( SENSE, 'B' )
  302:       LQUERY = ( LWORK.EQ.-1 )
  303: *
  304:       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
  305:          INFO = -1
  306:       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
  307:          INFO = -2
  308:       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
  309:      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
  310:          INFO = -4
  311:       ELSE IF( N.LT.0 ) THEN
  312:          INFO = -5
  313:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  314:          INFO = -7
  315:       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
  316:          INFO = -11
  317:       END IF
  318: *
  319: *     Compute workspace
  320: *      (Note: Comments in the code beginning "Workspace:" describe the
  321: *       minimal amount of real workspace needed at that point in the
  322: *       code, as well as the preferred amount for good performance.
  323: *       CWorkspace refers to complex workspace, and RWorkspace to real
  324: *       workspace. NB refers to the optimal block size for the
  325: *       immediately following subroutine, as returned by ILAENV.
  326: *       HSWORK refers to the workspace preferred by ZHSEQR, as
  327: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
  328: *       the worst case.
  329: *       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
  330: *       depends on SDIM, which is computed by the routine ZTRSEN later
  331: *       in the code.)
  332: *
  333:       IF( INFO.EQ.0 ) THEN
  334:          IF( N.EQ.0 ) THEN
  335:             MINWRK = 1
  336:             LWRK = 1
  337:          ELSE
  338:             MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
  339:             MINWRK = 2*N
  340: *
  341:             CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
  342:      $             WORK, -1, IEVAL )
  343:             HSWORK = WORK( 1 )
  344: *
  345:             IF( .NOT.WANTVS ) THEN
  346:                MAXWRK = MAX( MAXWRK, HSWORK )
  347:             ELSE
  348:                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
  349:      $                       ' ', N, 1, N, -1 ) )
  350:                MAXWRK = MAX( MAXWRK, HSWORK )
  351:             END IF
  352:             LWRK = MAXWRK
  353:             IF( .NOT.WANTSN )
  354:      $         LWRK = MAX( LWRK, ( N*N )/2 )
  355:          END IF
  356:          WORK( 1 ) = LWRK
  357: *
  358:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
  359:             INFO = -15
  360:          END IF
  361:       END IF
  362: *
  363:       IF( INFO.NE.0 ) THEN
  364:          CALL XERBLA( 'ZGEESX', -INFO )
  365:          RETURN
  366:       ELSE IF( LQUERY ) THEN
  367:          RETURN
  368:       END IF
  369: *
  370: *     Quick return if possible
  371: *
  372:       IF( N.EQ.0 ) THEN
  373:          SDIM = 0
  374:          RETURN
  375:       END IF
  376: *
  377: *     Get machine constants
  378: *
  379:       EPS = DLAMCH( 'P' )
  380:       SMLNUM = DLAMCH( 'S' )
  381:       BIGNUM = ONE / SMLNUM
  382:       CALL DLABAD( SMLNUM, BIGNUM )
  383:       SMLNUM = SQRT( SMLNUM ) / EPS
  384:       BIGNUM = ONE / SMLNUM
  385: *
  386: *     Scale A if max element outside range [SMLNUM,BIGNUM]
  387: *
  388:       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
  389:       SCALEA = .FALSE.
  390:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  391:          SCALEA = .TRUE.
  392:          CSCALE = SMLNUM
  393:       ELSE IF( ANRM.GT.BIGNUM ) THEN
  394:          SCALEA = .TRUE.
  395:          CSCALE = BIGNUM
  396:       END IF
  397:       IF( SCALEA )
  398:      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
  399: *
  400: *
  401: *     Permute the matrix to make it more nearly triangular
  402: *     (CWorkspace: none)
  403: *     (RWorkspace: need N)
  404: *
  405:       IBAL = 1
  406:       CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
  407: *
  408: *     Reduce to upper Hessenberg form
  409: *     (CWorkspace: need 2*N, prefer N+N*NB)
  410: *     (RWorkspace: none)
  411: *
  412:       ITAU = 1
  413:       IWRK = N + ITAU
  414:       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
  415:      $             LWORK-IWRK+1, IERR )
  416: *
  417:       IF( WANTVS ) THEN
  418: *
  419: *        Copy Householder vectors to VS
  420: *
  421:          CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
  422: *
  423: *        Generate unitary matrix in VS
  424: *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
  425: *        (RWorkspace: none)
  426: *
  427:          CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
  428:      $                LWORK-IWRK+1, IERR )
  429:       END IF
  430: *
  431:       SDIM = 0
  432: *
  433: *     Perform QR iteration, accumulating Schur vectors in VS if desired
  434: *     (CWorkspace: need 1, prefer HSWORK (see comments) )
  435: *     (RWorkspace: none)
  436: *
  437:       IWRK = ITAU
  438:       CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
  439:      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
  440:       IF( IEVAL.GT.0 )
  441:      $   INFO = IEVAL
  442: *
  443: *     Sort eigenvalues if desired
  444: *
  445:       IF( WANTST .AND. INFO.EQ.0 ) THEN
  446:          IF( SCALEA )
  447:      $      CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
  448:          DO 10 I = 1, N
  449:             BWORK( I ) = SELECT( W( I ) )
  450:    10    CONTINUE
  451: *
  452: *        Reorder eigenvalues, transform Schur vectors, and compute
  453: *        reciprocal condition numbers
  454: *        (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
  455: *                     otherwise, need none )
  456: *        (RWorkspace: none)
  457: *
  458:          CALL ZTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
  459:      $                RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
  460:      $                ICOND )
  461:          IF( .NOT.WANTSN )
  462:      $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
  463:          IF( ICOND.EQ.-14 ) THEN
  464: *
  465: *           Not enough complex workspace
  466: *
  467:             INFO = -15
  468:          END IF
  469:       END IF
  470: *
  471:       IF( WANTVS ) THEN
  472: *
  473: *        Undo balancing
  474: *        (CWorkspace: none)
  475: *        (RWorkspace: need N)
  476: *
  477:          CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
  478:      $                IERR )
  479:       END IF
  480: *
  481:       IF( SCALEA ) THEN
  482: *
  483: *        Undo scaling for the Schur form of A
  484: *
  485:          CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
  486:          CALL ZCOPY( N, A, LDA+1, W, 1 )
  487:          IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
  488:             DUM( 1 ) = RCONDV
  489:             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
  490:             RCONDV = DUM( 1 )
  491:          END IF
  492:       END IF
  493: *
  494:       WORK( 1 ) = MAXWRK
  495:       RETURN
  496: *
  497: *     End of ZGEESX
  498: *
  499:       END