|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
