LAPACK 3.3.1
Linear Algebra PACKage

cgeesx.f

Go to the documentation of this file.
00001       SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
00002      $                   VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
00003      $                   BWORK, INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.2.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     June 2010
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          JOBVS, SENSE, SORT
00012       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
00013       REAL               RCONDE, RCONDV
00014 *     ..
00015 *     .. Array Arguments ..
00016       LOGICAL            BWORK( * )
00017       REAL               RWORK( * )
00018       COMPLEX            A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
00019 *     ..
00020 *     .. Function Arguments ..
00021       LOGICAL            SELECT
00022       EXTERNAL           SELECT
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  CGEESX computes for an N-by-N complex nonsymmetric matrix A, the
00029 *  eigenvalues, the Schur form T, and, optionally, the matrix of Schur
00030 *  vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
00031 *
00032 *  Optionally, it also orders the eigenvalues on the diagonal of the
00033 *  Schur form so that selected eigenvalues are at the top left;
00034 *  computes a reciprocal condition number for the average of the
00035 *  selected eigenvalues (RCONDE); and computes a reciprocal condition
00036 *  number for the right invariant subspace corresponding to the
00037 *  selected eigenvalues (RCONDV).  The leading columns of Z form an
00038 *  orthonormal basis for this invariant subspace.
00039 *
00040 *  For further explanation of the reciprocal condition numbers RCONDE
00041 *  and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
00042 *  these quantities are called s and sep respectively).
00043 *
00044 *  A complex matrix is in Schur form if it is upper triangular.
00045 *
00046 *  Arguments
00047 *  =========
00048 *
00049 *  JOBVS   (input) CHARACTER*1
00050 *          = 'N': Schur vectors are not computed;
00051 *          = 'V': Schur vectors are computed.
00052 *
00053 *  SORT    (input) CHARACTER*1
00054 *          Specifies whether or not to order the eigenvalues on the
00055 *          diagonal of the Schur form.
00056 *          = 'N': Eigenvalues are not ordered;
00057 *          = 'S': Eigenvalues are ordered (see SELECT).
00058 *
00059 *  SELECT  (external procedure) LOGICAL FUNCTION of one COMPLEX argument
00060 *          SELECT must be declared EXTERNAL in the calling subroutine.
00061 *          If SORT = 'S', SELECT is used to select eigenvalues to order
00062 *          to the top left of the Schur form.
00063 *          If SORT = 'N', SELECT is not referenced.
00064 *          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
00065 *
00066 *  SENSE   (input) CHARACTER*1
00067 *          Determines which reciprocal condition numbers are computed.
00068 *          = 'N': None are computed;
00069 *          = 'E': Computed for average of selected eigenvalues only;
00070 *          = 'V': Computed for selected right invariant subspace only;
00071 *          = 'B': Computed for both.
00072 *          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
00073 *
00074 *  N       (input) INTEGER
00075 *          The order of the matrix A. N >= 0.
00076 *
00077 *  A       (input/output) COMPLEX array, dimension (LDA, N)
00078 *          On entry, the N-by-N matrix A.
00079 *          On exit, A is overwritten by its Schur form T.
00080 *
00081 *  LDA     (input) INTEGER
00082 *          The leading dimension of the array A.  LDA >= max(1,N).
00083 *
00084 *  SDIM    (output) INTEGER
00085 *          If SORT = 'N', SDIM = 0.
00086 *          If SORT = 'S', SDIM = number of eigenvalues for which
00087 *                         SELECT is true.
00088 *
00089 *  W       (output) COMPLEX array, dimension (N)
00090 *          W contains the computed eigenvalues, in the same order
00091 *          that they appear on the diagonal of the output Schur form T.
00092 *
00093 *  VS      (output) COMPLEX array, dimension (LDVS,N)
00094 *          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
00095 *          vectors.
00096 *          If JOBVS = 'N', VS is not referenced.
00097 *
00098 *  LDVS    (input) INTEGER
00099 *          The leading dimension of the array VS.  LDVS >= 1, and if
00100 *          JOBVS = 'V', LDVS >= N.
00101 *
00102 *  RCONDE  (output) REAL
00103 *          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
00104 *          condition number for the average of the selected eigenvalues.
00105 *          Not referenced if SENSE = 'N' or 'V'.
00106 *
00107 *  RCONDV  (output) REAL
00108 *          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
00109 *          condition number for the selected right invariant subspace.
00110 *          Not referenced if SENSE = 'N' or 'E'.
00111 *
00112 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00113 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00114 *
00115 *  LWORK   (input) INTEGER
00116 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00117 *          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
00118 *          where SDIM is the number of selected eigenvalues computed by
00119 *          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
00120 *          that an error is only returned if LWORK < max(1,2*N), but if
00121 *          SENSE = 'E' or 'V' or 'B' this may not be large enough.
00122 *          For good performance, LWORK must generally be larger.
00123 *
00124 *          If LWORK = -1, then a workspace query is assumed; the routine
00125 *          only calculates upper bound on the optimal size of the
00126 *          array WORK, returns this value as the first entry of the WORK
00127 *          array, and no error message related to LWORK is issued by
00128 *          XERBLA.
00129 *
00130 *  RWORK   (workspace) REAL array, dimension (N)
00131 *
00132 *  BWORK   (workspace) LOGICAL array, dimension (N)
00133 *          Not referenced if SORT = 'N'.
00134 *
00135 *  INFO    (output) INTEGER
00136 *          = 0: successful exit
00137 *          < 0: if INFO = -i, the i-th argument had an illegal value.
00138 *          > 0: if INFO = i, and i is
00139 *             <= N: the QR algorithm failed to compute all the
00140 *                   eigenvalues; elements 1:ILO-1 and i+1:N of W
00141 *                   contain those eigenvalues which have converged; if
00142 *                   JOBVS = 'V', VS contains the transformation which
00143 *                   reduces A to its partially converged Schur form.
00144 *             = N+1: the eigenvalues could not be reordered because some
00145 *                   eigenvalues were too close to separate (the problem
00146 *                   is very ill-conditioned);
00147 *             = N+2: after reordering, roundoff changed values of some
00148 *                   complex eigenvalues so that leading eigenvalues in
00149 *                   the Schur form no longer satisfy SELECT=.TRUE.  This
00150 *                   could also be caused by underflow due to scaling.
00151 *
00152 *  =====================================================================
00153 *
00154 *     .. Parameters ..
00155       REAL               ZERO, ONE
00156       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00157 *     ..
00158 *     .. Local Scalars ..
00159       LOGICAL            LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
00160      $                   WANTSV, WANTVS
00161       INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
00162      $                   ITAU, IWRK, LWRK, MAXWRK, MINWRK
00163       REAL               ANRM, BIGNUM, CSCALE, EPS, SMLNUM
00164 *     ..
00165 *     .. Local Arrays ..
00166       REAL               DUM( 1 )
00167 *     ..
00168 *     .. External Subroutines ..
00169       EXTERNAL           CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY,
00170      $                   CLASCL, CTRSEN, CUNGHR, SLABAD, SLASCL, XERBLA
00171 *     ..
00172 *     .. External Functions ..
00173       LOGICAL            LSAME
00174       INTEGER            ILAENV
00175       REAL               CLANGE, SLAMCH
00176       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
00177 *     ..
00178 *     .. Intrinsic Functions ..
00179       INTRINSIC          MAX, SQRT
00180 *     ..
00181 *     .. Executable Statements ..
00182 *
00183 *     Test the input arguments
00184 *
00185       INFO = 0
00186       WANTVS = LSAME( JOBVS, 'V' )
00187       WANTST = LSAME( SORT, 'S' )
00188       WANTSN = LSAME( SENSE, 'N' )
00189       WANTSE = LSAME( SENSE, 'E' )
00190       WANTSV = LSAME( SENSE, 'V' )
00191       WANTSB = LSAME( SENSE, 'B' )
00192       LQUERY = ( LWORK.EQ.-1 )
00193 *
00194       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
00195          INFO = -1
00196       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00197          INFO = -2
00198       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
00199      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
00200          INFO = -4
00201       ELSE IF( N.LT.0 ) THEN
00202          INFO = -5
00203       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00204          INFO = -7
00205       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
00206          INFO = -11
00207       END IF
00208 *
00209 *     Compute workspace
00210 *      (Note: Comments in the code beginning "Workspace:" describe the
00211 *       minimal amount of real workspace needed at that point in the
00212 *       code, as well as the preferred amount for good performance.
00213 *       CWorkspace refers to complex workspace, and RWorkspace to real
00214 *       workspace. NB refers to the optimal block size for the
00215 *       immediately following subroutine, as returned by ILAENV.
00216 *       HSWORK refers to the workspace preferred by CHSEQR, as
00217 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00218 *       the worst case.
00219 *       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
00220 *       depends on SDIM, which is computed by the routine CTRSEN later
00221 *       in the code.)
00222 *
00223       IF( INFO.EQ.0 ) THEN
00224          IF( N.EQ.0 ) THEN
00225             MINWRK = 1
00226             LWRK = 1
00227          ELSE
00228             MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
00229             MINWRK = 2*N
00230 *
00231             CALL CHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
00232      $             WORK, -1, IEVAL )
00233             HSWORK = WORK( 1 )
00234 *
00235             IF( .NOT.WANTVS ) THEN
00236                MAXWRK = MAX( MAXWRK, HSWORK )
00237             ELSE
00238                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00239      $                       ' ', N, 1, N, -1 ) )
00240                MAXWRK = MAX( MAXWRK, HSWORK )
00241             END IF
00242             LWRK = MAXWRK
00243             IF( .NOT.WANTSN )
00244      $         LWRK = MAX( LWRK, ( N*N )/2 )
00245          END IF
00246          WORK( 1 ) = LWRK
00247 *
00248          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00249             INFO = -15
00250          END IF
00251       END IF
00252 *
00253       IF( INFO.NE.0 ) THEN
00254          CALL XERBLA( 'CGEESX', -INFO )
00255          RETURN
00256       ELSE IF( LQUERY ) THEN
00257          RETURN
00258       END IF
00259 *
00260 *     Quick return if possible
00261 *
00262       IF( N.EQ.0 ) THEN
00263          SDIM = 0
00264          RETURN
00265       END IF
00266 *
00267 *     Get machine constants
00268 *
00269       EPS = SLAMCH( 'P' )
00270       SMLNUM = SLAMCH( 'S' )
00271       BIGNUM = ONE / SMLNUM
00272       CALL SLABAD( SMLNUM, BIGNUM )
00273       SMLNUM = SQRT( SMLNUM ) / EPS
00274       BIGNUM = ONE / SMLNUM
00275 *
00276 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00277 *
00278       ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
00279       SCALEA = .FALSE.
00280       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00281          SCALEA = .TRUE.
00282          CSCALE = SMLNUM
00283       ELSE IF( ANRM.GT.BIGNUM ) THEN
00284          SCALEA = .TRUE.
00285          CSCALE = BIGNUM
00286       END IF
00287       IF( SCALEA )
00288      $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00289 *
00290 *
00291 *     Permute the matrix to make it more nearly triangular
00292 *     (CWorkspace: none)
00293 *     (RWorkspace: need N)
00294 *
00295       IBAL = 1
00296       CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00297 *
00298 *     Reduce to upper Hessenberg form
00299 *     (CWorkspace: need 2*N, prefer N+N*NB)
00300 *     (RWorkspace: none)
00301 *
00302       ITAU = 1
00303       IWRK = N + ITAU
00304       CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00305      $             LWORK-IWRK+1, IERR )
00306 *
00307       IF( WANTVS ) THEN
00308 *
00309 *        Copy Householder vectors to VS
00310 *
00311          CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS )
00312 *
00313 *        Generate unitary matrix in VS
00314 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00315 *        (RWorkspace: none)
00316 *
00317          CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
00318      $                LWORK-IWRK+1, IERR )
00319       END IF
00320 *
00321       SDIM = 0
00322 *
00323 *     Perform QR iteration, accumulating Schur vectors in VS if desired
00324 *     (CWorkspace: need 1, prefer HSWORK (see comments) )
00325 *     (RWorkspace: none)
00326 *
00327       IWRK = ITAU
00328       CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
00329      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
00330       IF( IEVAL.GT.0 )
00331      $   INFO = IEVAL
00332 *
00333 *     Sort eigenvalues if desired
00334 *
00335       IF( WANTST .AND. INFO.EQ.0 ) THEN
00336          IF( SCALEA )
00337      $      CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
00338          DO 10 I = 1, N
00339             BWORK( I ) = SELECT( W( I ) )
00340    10    CONTINUE
00341 *
00342 *        Reorder eigenvalues, transform Schur vectors, and compute
00343 *        reciprocal condition numbers
00344 *        (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
00345 *                     otherwise, need none )
00346 *        (RWorkspace: none)
00347 *
00348          CALL CTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
00349      $                RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
00350      $                ICOND )
00351          IF( .NOT.WANTSN )
00352      $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
00353          IF( ICOND.EQ.-14 ) THEN
00354 *
00355 *           Not enough complex workspace
00356 *
00357             INFO = -15
00358          END IF
00359       END IF
00360 *
00361       IF( WANTVS ) THEN
00362 *
00363 *        Undo balancing
00364 *        (CWorkspace: none)
00365 *        (RWorkspace: need N)
00366 *
00367          CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
00368      $                IERR )
00369       END IF
00370 *
00371       IF( SCALEA ) THEN
00372 *
00373 *        Undo scaling for the Schur form of A
00374 *
00375          CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
00376          CALL CCOPY( N, A, LDA+1, W, 1 )
00377          IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
00378             DUM( 1 ) = RCONDV
00379             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
00380             RCONDV = DUM( 1 )
00381          END IF
00382       END IF
00383 *
00384       WORK( 1 ) = MAXWRK
00385       RETURN
00386 *
00387 *     End of CGEESX
00388 *
00389       END
 All Files Functions