1 SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, 2 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 3 $ BWORK, INFO ) 4* 5* -- LAPACK driver routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER JOBVS, SENSE, SORT 12 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM 13 REAL RCONDE, RCONDV 14* .. 15* .. Array Arguments .. 16 LOGICAL BWORK( * ) 17 REAL RWORK( * ) 18 COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) 19* .. 20* .. Function Arguments .. 21 LOGICAL SELECT 22 EXTERNAL SELECT 23* .. 24* 25* Purpose 26* ======= 27* 28* CGEESX computes for an N-by-N complex nonsymmetric matrix A, the 29* eigenvalues, the Schur form T, and, optionally, the matrix of Schur 30* vectors Z. This gives the Schur factorization A = Z*T*(Z**H). 31* 32* Optionally, it also orders the eigenvalues on the diagonal of the 33* Schur form so that selected eigenvalues are at the top left; 34* computes a reciprocal condition number for the average of the 35* selected eigenvalues (RCONDE); and computes a reciprocal condition 36* number for the right invariant subspace corresponding to the 37* selected eigenvalues (RCONDV). The leading columns of Z form an 38* orthonormal basis for this invariant subspace. 39* 40* For further explanation of the reciprocal condition numbers RCONDE 41* and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where 42* these quantities are called s and sep respectively). 43* 44* A complex matrix is in Schur form if it is upper triangular. 45* 46* Arguments 47* ========= 48* 49* JOBVS (input) CHARACTER*1 50* = 'N': Schur vectors are not computed; 51* = 'V': Schur vectors are computed. 52* 53* SORT (input) CHARACTER*1 54* Specifies whether or not to order the eigenvalues on the 55* diagonal of the Schur form. 56* = 'N': Eigenvalues are not ordered; 57* = 'S': Eigenvalues are ordered (see SELECT). 58* 59* SELECT (input) LOGICAL FUNCTION of one COMPLEX argument 60* SELECT must be declared EXTERNAL in the calling subroutine. 61* If SORT = 'S', SELECT is used to select eigenvalues to order 62* to the top left of the Schur form. 63* If SORT = 'N', SELECT is not referenced. 64* An eigenvalue W(j) is selected if SELECT(W(j)) is true. 65* 66* SENSE (input) CHARACTER*1 67* Determines which reciprocal condition numbers are computed. 68* = 'N': None are computed; 69* = 'E': Computed for average of selected eigenvalues only; 70* = 'V': Computed for selected right invariant subspace only; 71* = 'B': Computed for both. 72* If SENSE = 'E', 'V' or 'B', SORT must equal 'S'. 73* 74* N (input) INTEGER 75* The order of the matrix A. N >= 0. 76* 77* A (input/output) COMPLEX array, dimension (LDA, N) 78* On entry, the N-by-N matrix A. 79* On exit, A is overwritten by its Schur form T. 80* 81* LDA (input) INTEGER 82* The leading dimension of the array A. LDA >= max(1,N). 83* 84* SDIM (output) INTEGER 85* If SORT = 'N', SDIM = 0. 86* If SORT = 'S', SDIM = number of eigenvalues for which 87* SELECT is true. 88* 89* W (output) COMPLEX array, dimension (N) 90* W contains the computed eigenvalues, in the same order 91* that they appear on the diagonal of the output Schur form T. 92* 93* VS (output) COMPLEX array, dimension (LDVS,N) 94* If JOBVS = 'V', VS contains the unitary matrix Z of Schur 95* vectors. 96* If JOBVS = 'N', VS is not referenced. 97* 98* LDVS (input) INTEGER 99* The leading dimension of the array VS. LDVS >= 1, and if 100* JOBVS = 'V', LDVS >= N. 101* 102* RCONDE (output) REAL 103* If SENSE = 'E' or 'B', RCONDE contains the reciprocal 104* condition number for the average of the selected eigenvalues. 105* Not referenced if SENSE = 'N' or 'V'. 106* 107* RCONDV (output) REAL 108* If SENSE = 'V' or 'B', RCONDV contains the reciprocal 109* condition number for the selected right invariant subspace. 110* Not referenced if SENSE = 'N' or 'E'. 111* 112* WORK (workspace/output) COMPLEX array, dimension (LWORK) 113* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 114* 115* LWORK (input) INTEGER 116* The dimension of the array WORK. LWORK >= max(1,2*N). 117* Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM), 118* where SDIM is the number of selected eigenvalues computed by 119* this routine. Note that 2*SDIM*(N-SDIM) <= N*N/2. 120* For good performance, LWORK must generally be larger. 121* 122* RWORK (workspace) REAL array, dimension (N) 123* 124* BWORK (workspace) LOGICAL array, dimension (N) 125* Not referenced if SORT = 'N'. 126* 127* INFO (output) INTEGER 128* = 0: successful exit 129* < 0: if INFO = -i, the i-th argument had an illegal value. 130* > 0: if INFO = i, and i is 131* <= N: the QR algorithm failed to compute all the 132* eigenvalues; elements 1:ILO-1 and i+1:N of W 133* contain those eigenvalues which have converged; if 134* JOBVS = 'V', VS contains the transformation which 135* reduces A to its partially converged Schur form. 136* = N+1: the eigenvalues could not be reordered because some 137* eigenvalues were too close to separate (the problem 138* is very ill-conditioned); 139* = N+2: after reordering, roundoff changed values of some 140* complex eigenvalues so that leading eigenvalues in 141* the Schur form no longer satisfy SELECT=.TRUE. This 142* could also be caused by underflow due to scaling. 143* 144* ===================================================================== 145* 146* .. Parameters .. 147 REAL ZERO, ONE 148 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 149* .. 150* .. Local Scalars .. 151 LOGICAL SCALEA, WANTSB, WANTSE, WANTSN, WANTST, 152 $ WANTSV, WANTVS 153 INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO, 154 $ ITAU, IWRK, K, MAXB, MAXWRK, MINWRK 155 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM 156* .. 157* .. Local Arrays .. 158 REAL DUM( 1 ) 159* .. 160* .. External Subroutines .. 161 EXTERNAL CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, 162 $ CLASCL, CTRSEN, CUNGHR, SLABAD, SLASCL, XERBLA 163* .. 164* .. External Functions .. 165 LOGICAL LSAME 166 INTEGER ILAENV 167 REAL CLANGE, SLAMCH 168 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH 169* .. 170* .. Intrinsic Functions .. 171 INTRINSIC MAX, MIN, SQRT 172* .. 173* .. Executable Statements .. 174* 175* Test the input arguments 176* 177 INFO = 0 178 WANTVS = LSAME( JOBVS, 'V' ) 179 WANTST = LSAME( SORT, 'S' ) 180 WANTSN = LSAME( SENSE, 'N' ) 181 WANTSE = LSAME( SENSE, 'E' ) 182 WANTSV = LSAME( SENSE, 'V' ) 183 WANTSB = LSAME( SENSE, 'B' ) 184 IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN 185 INFO = -1 186 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 187 INFO = -2 188 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 189 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 190 INFO = -4 191 ELSE IF( N.LT.0 ) THEN 192 INFO = -5 193 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 194 INFO = -7 195 ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN 196 INFO = -11 197 END IF 198* 199* Compute workspace 200* (Note: Comments in the code beginning "Workspace:" describe the 201* minimal amount of real workspace needed at that point in the 202* code, as well as the preferred amount for good performance. 203* CWorkspace refers to complex workspace, and RWorkspace to real 204* workspace. NB refers to the optimal block size for the 205* immediately following subroutine, as returned by ILAENV. 206* HSWORK refers to the workspace preferred by CHSEQR, as 207* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 208* the worst case. 209* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed 210* depends on SDIM, which is computed by the routine CTRSEN later 211* in the code.) 212* 213 MINWRK = 1 214 IF( INFO.EQ.0 .AND. ( LWORK.GE.1 ) ) THEN 215 MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 ) 216 MINWRK = MAX( 1, 2*N ) 217 IF( .NOT.WANTVS ) THEN 218 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SN', N, 1, N, -1 ), 2 ) 219 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SN', N, 1, 220 $ N, -1 ) ) ) 221 HSWORK = MAX( K*( K+2 ), 2*N ) 222 MAXWRK = MAX( MAXWRK, HSWORK, 1 ) 223 ELSE 224 MAXWRK = MAX( MAXWRK, N+( N-1 )* 225 $ ILAENV( 1, 'CUNGHR', ' ', N, 1, N, -1 ) ) 226 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SV', N, 1, N, -1 ), 2 ) 227 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SV', N, 1, 228 $ N, -1 ) ) ) 229 HSWORK = MAX( K*( K+2 ), 2*N ) 230 MAXWRK = MAX( MAXWRK, HSWORK, 1 ) 231 END IF 232 WORK( 1 ) = MAXWRK 233 END IF 234 IF( LWORK.LT.MINWRK ) THEN 235 INFO = -15 236 END IF 237 IF( INFO.NE.0 ) THEN 238 CALL XERBLA( 'CGEESX', -INFO ) 239 RETURN 240 END IF 241* 242* Quick return if possible 243* 244 IF( N.EQ.0 ) THEN 245 SDIM = 0 246 RETURN 247 END IF 248* 249* Get machine constants 250* 251 EPS = SLAMCH( 'P' ) 252 SMLNUM = SLAMCH( 'S' ) 253 BIGNUM = ONE / SMLNUM 254 CALL SLABAD( SMLNUM, BIGNUM ) 255 SMLNUM = SQRT( SMLNUM ) / EPS 256 BIGNUM = ONE / SMLNUM 257* 258* Scale A if max element outside range [SMLNUM,BIGNUM] 259* 260 ANRM = CLANGE( 'M', N, N, A, LDA, DUM ) 261 SCALEA = .FALSE. 262 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 263 SCALEA = .TRUE. 264 CSCALE = SMLNUM 265 ELSE IF( ANRM.GT.BIGNUM ) THEN 266 SCALEA = .TRUE. 267 CSCALE = BIGNUM 268 END IF 269 IF( SCALEA ) 270 $ CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 271* 272* 273* Permute the matrix to make it more nearly triangular 274* (CWorkspace: none) 275* (RWorkspace: need N) 276* 277 IBAL = 1 278 CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR ) 279* 280* Reduce to upper Hessenberg form 281* (CWorkspace: need 2*N, prefer N+N*NB) 282* (RWorkspace: none) 283* 284 ITAU = 1 285 IWRK = N + ITAU 286 CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 287 $ LWORK-IWRK+1, IERR ) 288* 289 IF( WANTVS ) THEN 290* 291* Copy Householder vectors to VS 292* 293 CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS ) 294* 295* Generate unitary matrix in VS 296* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) 297* (RWorkspace: none) 298* 299 CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ), 300 $ LWORK-IWRK+1, IERR ) 301 END IF 302* 303 SDIM = 0 304* 305* Perform QR iteration, accumulating Schur vectors in VS if desired 306* (CWorkspace: need 1, prefer HSWORK (see comments) ) 307* (RWorkspace: none) 308* 309 IWRK = ITAU 310 CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS, 311 $ WORK( IWRK ), LWORK-IWRK+1, IEVAL ) 312 IF( IEVAL.GT.0 ) 313 $ INFO = IEVAL 314* 315* Sort eigenvalues if desired 316* 317 IF( WANTST .AND. INFO.EQ.0 ) THEN 318 IF( SCALEA ) 319 $ CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR ) 320 DO 10 I = 1, N 321 BWORK( I ) = SELECT( W( I ) ) 322 10 CONTINUE 323* 324* Reorder eigenvalues, transform Schur vectors, and compute 325* reciprocal condition numbers 326* (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM) 327* otherwise, need none ) 328* (RWorkspace: none) 329* 330 CALL CTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM, 331 $ RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1, 332 $ ICOND ) 333 IF( .NOT.WANTSN ) 334 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 335 IF( ICOND.EQ.-14 ) THEN 336* 337* Not enough complex workspace 338* 339 INFO = -15 340 END IF 341 END IF 342* 343 IF( WANTVS ) THEN 344* 345* Undo balancing 346* (CWorkspace: none) 347* (RWorkspace: need N) 348* 349 CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS, 350 $ IERR ) 351 END IF 352* 353 IF( SCALEA ) THEN 354* 355* Undo scaling for the Schur form of A 356* 357 CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR ) 358 CALL CCOPY( N, A, LDA+1, W, 1 ) 359 IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN 360 DUM( 1 ) = RCONDV 361 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 362 RCONDV = DUM( 1 ) 363 END IF 364 END IF 365* 366 WORK( 1 ) = MAXWRK 367 RETURN 368* 369* End of CGEESX 370* 371 END 372