Return to zgees.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

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