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