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