Annotation of rpl/lapack/lapack/zgeesx.f, revision 1.7

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

CVSweb interface <joel.bertrand@systella.fr>