Return to zgees.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

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