![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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