1*> \brief <b> CGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CGEESX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeesx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeesx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeesx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, 22* VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 23* BWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER JOBVS, SENSE, SORT 27* INTEGER INFO, LDA, LDVS, LWORK, N, SDIM 28* REAL RCONDE, RCONDV 29* .. 30* .. Array Arguments .. 31* LOGICAL BWORK( * ) 32* REAL RWORK( * ) 33* COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) 34* .. 35* .. Function Arguments .. 36* LOGICAL SELECT 37* EXTERNAL SELECT 38* .. 39* 40* 41*> \par Purpose: 42* ============= 43*> 44*> \verbatim 45*> 46*> CGEESX computes for an N-by-N complex nonsymmetric matrix A, the 47*> eigenvalues, the Schur form T, and, optionally, the matrix of Schur 48*> vectors Z. This gives the Schur factorization A = Z*T*(Z**H). 49*> 50*> Optionally, it also orders the eigenvalues on the diagonal of the 51*> Schur form so that selected eigenvalues are at the top left; 52*> computes a reciprocal condition number for the average of the 53*> selected eigenvalues (RCONDE); and computes a reciprocal condition 54*> number for the right invariant subspace corresponding to the 55*> selected eigenvalues (RCONDV). The leading columns of Z form an 56*> orthonormal basis for this invariant subspace. 57*> 58*> For further explanation of the reciprocal condition numbers RCONDE 59*> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where 60*> these quantities are called s and sep respectively). 61*> 62*> A complex matrix is in Schur form if it is upper triangular. 63*> \endverbatim 64* 65* Arguments: 66* ========== 67* 68*> \param[in] JOBVS 69*> \verbatim 70*> JOBVS is CHARACTER*1 71*> = 'N': Schur vectors are not computed; 72*> = 'V': Schur vectors are computed. 73*> \endverbatim 74*> 75*> \param[in] SORT 76*> \verbatim 77*> SORT is CHARACTER*1 78*> Specifies whether or not to order the eigenvalues on the 79*> diagonal of the Schur form. 80*> = 'N': Eigenvalues are not ordered; 81*> = 'S': Eigenvalues are ordered (see SELECT). 82*> \endverbatim 83*> 84*> \param[in] SELECT 85*> \verbatim 86*> SELECT is a LOGICAL FUNCTION of one COMPLEX argument 87*> SELECT must be declared EXTERNAL in the calling subroutine. 88*> If SORT = 'S', SELECT is used to select eigenvalues to order 89*> to the top left of the Schur form. 90*> If SORT = 'N', SELECT is not referenced. 91*> An eigenvalue W(j) is selected if SELECT(W(j)) is true. 92*> \endverbatim 93*> 94*> \param[in] SENSE 95*> \verbatim 96*> SENSE is CHARACTER*1 97*> Determines which reciprocal condition numbers are computed. 98*> = 'N': None are computed; 99*> = 'E': Computed for average of selected eigenvalues only; 100*> = 'V': Computed for selected right invariant subspace only; 101*> = 'B': Computed for both. 102*> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'. 103*> \endverbatim 104*> 105*> \param[in] N 106*> \verbatim 107*> N is INTEGER 108*> The order of the matrix A. N >= 0. 109*> \endverbatim 110*> 111*> \param[in,out] A 112*> \verbatim 113*> A is COMPLEX array, dimension (LDA, N) 114*> On entry, the N-by-N matrix A. 115*> On exit, A is overwritten by its Schur form T. 116*> \endverbatim 117*> 118*> \param[in] LDA 119*> \verbatim 120*> LDA is INTEGER 121*> The leading dimension of the array A. LDA >= max(1,N). 122*> \endverbatim 123*> 124*> \param[out] SDIM 125*> \verbatim 126*> SDIM is INTEGER 127*> If SORT = 'N', SDIM = 0. 128*> If SORT = 'S', SDIM = number of eigenvalues for which 129*> SELECT is true. 130*> \endverbatim 131*> 132*> \param[out] W 133*> \verbatim 134*> W is COMPLEX array, dimension (N) 135*> W contains the computed eigenvalues, in the same order 136*> that they appear on the diagonal of the output Schur form T. 137*> \endverbatim 138*> 139*> \param[out] VS 140*> \verbatim 141*> VS is COMPLEX array, dimension (LDVS,N) 142*> If JOBVS = 'V', VS contains the unitary matrix Z of Schur 143*> vectors. 144*> If JOBVS = 'N', VS is not referenced. 145*> \endverbatim 146*> 147*> \param[in] LDVS 148*> \verbatim 149*> LDVS is INTEGER 150*> The leading dimension of the array VS. LDVS >= 1, and if 151*> JOBVS = 'V', LDVS >= N. 152*> \endverbatim 153*> 154*> \param[out] RCONDE 155*> \verbatim 156*> RCONDE is REAL 157*> If SENSE = 'E' or 'B', RCONDE contains the reciprocal 158*> condition number for the average of the selected eigenvalues. 159*> Not referenced if SENSE = 'N' or 'V'. 160*> \endverbatim 161*> 162*> \param[out] RCONDV 163*> \verbatim 164*> RCONDV is REAL 165*> If SENSE = 'V' or 'B', RCONDV contains the reciprocal 166*> condition number for the selected right invariant subspace. 167*> Not referenced if SENSE = 'N' or 'E'. 168*> \endverbatim 169*> 170*> \param[out] WORK 171*> \verbatim 172*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 173*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 174*> \endverbatim 175*> 176*> \param[in] LWORK 177*> \verbatim 178*> LWORK is INTEGER 179*> The dimension of the array WORK. LWORK >= max(1,2*N). 180*> Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM), 181*> where SDIM is the number of selected eigenvalues computed by 182*> this routine. Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also 183*> that an error is only returned if LWORK < max(1,2*N), but if 184*> SENSE = 'E' or 'V' or 'B' this may not be large enough. 185*> For good performance, LWORK must generally be larger. 186*> 187*> If LWORK = -1, then a workspace query is assumed; the routine 188*> only calculates upper bound on the optimal size of the 189*> array WORK, returns this value as the first entry of the WORK 190*> array, and no error message related to LWORK is issued by 191*> XERBLA. 192*> \endverbatim 193*> 194*> \param[out] RWORK 195*> \verbatim 196*> RWORK is REAL array, dimension (N) 197*> \endverbatim 198*> 199*> \param[out] BWORK 200*> \verbatim 201*> BWORK is LOGICAL array, dimension (N) 202*> Not referenced if SORT = 'N'. 203*> \endverbatim 204*> 205*> \param[out] INFO 206*> \verbatim 207*> INFO is INTEGER 208*> = 0: successful exit 209*> < 0: if INFO = -i, the i-th argument had an illegal value. 210*> > 0: if INFO = i, and i is 211*> <= N: the QR algorithm failed to compute all the 212*> eigenvalues; elements 1:ILO-1 and i+1:N of W 213*> contain those eigenvalues which have converged; if 214*> JOBVS = 'V', VS contains the transformation which 215*> reduces A to its partially converged Schur form. 216*> = N+1: the eigenvalues could not be reordered because some 217*> eigenvalues were too close to separate (the problem 218*> is very ill-conditioned); 219*> = N+2: after reordering, roundoff changed values of some 220*> complex eigenvalues so that leading eigenvalues in 221*> the Schur form no longer satisfy SELECT=.TRUE. This 222*> could also be caused by underflow due to scaling. 223*> \endverbatim 224* 225* Authors: 226* ======== 227* 228*> \author Univ. of Tennessee 229*> \author Univ. of California Berkeley 230*> \author Univ. of Colorado Denver 231*> \author NAG Ltd. 232* 233*> \ingroup complexGEeigen 234* 235* ===================================================================== 236 SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, 237 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 238 $ BWORK, INFO ) 239* 240* -- LAPACK driver routine -- 241* -- LAPACK is a software package provided by Univ. of Tennessee, -- 242* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 243* 244* .. Scalar Arguments .. 245 CHARACTER JOBVS, SENSE, SORT 246 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM 247 REAL RCONDE, RCONDV 248* .. 249* .. Array Arguments .. 250 LOGICAL BWORK( * ) 251 REAL RWORK( * ) 252 COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) 253* .. 254* .. Function Arguments .. 255 LOGICAL SELECT 256 EXTERNAL SELECT 257* .. 258* 259* ===================================================================== 260* 261* .. Parameters .. 262 REAL ZERO, ONE 263 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 264* .. 265* .. Local Scalars .. 266 LOGICAL LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST, 267 $ WANTSV, WANTVS 268 INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO, 269 $ ITAU, IWRK, LWRK, MAXWRK, MINWRK 270 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM 271* .. 272* .. Local Arrays .. 273 REAL DUM( 1 ) 274* .. 275* .. External Subroutines .. 276 EXTERNAL CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, 277 $ CLASCL, CTRSEN, CUNGHR, SLABAD, SLASCL, XERBLA 278* .. 279* .. External Functions .. 280 LOGICAL LSAME 281 INTEGER ILAENV 282 REAL CLANGE, SLAMCH 283 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH 284* .. 285* .. Intrinsic Functions .. 286 INTRINSIC MAX, SQRT 287* .. 288* .. Executable Statements .. 289* 290* Test the input arguments 291* 292 INFO = 0 293 WANTVS = LSAME( JOBVS, 'V' ) 294 WANTST = LSAME( SORT, 'S' ) 295 WANTSN = LSAME( SENSE, 'N' ) 296 WANTSE = LSAME( SENSE, 'E' ) 297 WANTSV = LSAME( SENSE, 'V' ) 298 WANTSB = LSAME( SENSE, 'B' ) 299 LQUERY = ( LWORK.EQ.-1 ) 300* 301 IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN 302 INFO = -1 303 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 304 INFO = -2 305 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 306 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 307 INFO = -4 308 ELSE IF( N.LT.0 ) THEN 309 INFO = -5 310 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 311 INFO = -7 312 ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN 313 INFO = -11 314 END IF 315* 316* Compute workspace 317* (Note: Comments in the code beginning "Workspace:" describe the 318* minimal amount of real workspace needed at that point in the 319* code, as well as the preferred amount for good performance. 320* CWorkspace refers to complex workspace, and RWorkspace to real 321* workspace. NB refers to the optimal block size for the 322* immediately following subroutine, as returned by ILAENV. 323* HSWORK refers to the workspace preferred by CHSEQR, as 324* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 325* the worst case. 326* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed 327* depends on SDIM, which is computed by the routine CTRSEN later 328* in the code.) 329* 330 IF( INFO.EQ.0 ) THEN 331 IF( N.EQ.0 ) THEN 332 MINWRK = 1 333 LWRK = 1 334 ELSE 335 MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 ) 336 MINWRK = 2*N 337* 338 CALL CHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS, 339 $ WORK, -1, IEVAL ) 340 HSWORK = REAL( WORK( 1 ) ) 341* 342 IF( .NOT.WANTVS ) THEN 343 MAXWRK = MAX( MAXWRK, HSWORK ) 344 ELSE 345 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR', 346 $ ' ', N, 1, N, -1 ) ) 347 MAXWRK = MAX( MAXWRK, HSWORK ) 348 END IF 349 LWRK = MAXWRK 350 IF( .NOT.WANTSN ) 351 $ LWRK = MAX( LWRK, ( N*N )/2 ) 352 END IF 353 WORK( 1 ) = LWRK 354* 355 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 356 INFO = -15 357 END IF 358 END IF 359* 360 IF( INFO.NE.0 ) THEN 361 CALL XERBLA( 'CGEESX', -INFO ) 362 RETURN 363 ELSE IF( LQUERY ) THEN 364 RETURN 365 END IF 366* 367* Quick return if possible 368* 369 IF( N.EQ.0 ) THEN 370 SDIM = 0 371 RETURN 372 END IF 373* 374* Get machine constants 375* 376 EPS = SLAMCH( 'P' ) 377 SMLNUM = SLAMCH( 'S' ) 378 BIGNUM = ONE / SMLNUM 379 CALL SLABAD( SMLNUM, BIGNUM ) 380 SMLNUM = SQRT( SMLNUM ) / EPS 381 BIGNUM = ONE / SMLNUM 382* 383* Scale A if max element outside range [SMLNUM,BIGNUM] 384* 385 ANRM = CLANGE( 'M', N, N, A, LDA, DUM ) 386 SCALEA = .FALSE. 387 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 388 SCALEA = .TRUE. 389 CSCALE = SMLNUM 390 ELSE IF( ANRM.GT.BIGNUM ) THEN 391 SCALEA = .TRUE. 392 CSCALE = BIGNUM 393 END IF 394 IF( SCALEA ) 395 $ CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 396* 397* 398* Permute the matrix to make it more nearly triangular 399* (CWorkspace: none) 400* (RWorkspace: need N) 401* 402 IBAL = 1 403 CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR ) 404* 405* Reduce to upper Hessenberg form 406* (CWorkspace: need 2*N, prefer N+N*NB) 407* (RWorkspace: none) 408* 409 ITAU = 1 410 IWRK = N + ITAU 411 CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 412 $ LWORK-IWRK+1, IERR ) 413* 414 IF( WANTVS ) THEN 415* 416* Copy Householder vectors to VS 417* 418 CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS ) 419* 420* Generate unitary matrix in VS 421* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) 422* (RWorkspace: none) 423* 424 CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ), 425 $ LWORK-IWRK+1, IERR ) 426 END IF 427* 428 SDIM = 0 429* 430* Perform QR iteration, accumulating Schur vectors in VS if desired 431* (CWorkspace: need 1, prefer HSWORK (see comments) ) 432* (RWorkspace: none) 433* 434 IWRK = ITAU 435 CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS, 436 $ WORK( IWRK ), LWORK-IWRK+1, IEVAL ) 437 IF( IEVAL.GT.0 ) 438 $ INFO = IEVAL 439* 440* Sort eigenvalues if desired 441* 442 IF( WANTST .AND. INFO.EQ.0 ) THEN 443 IF( SCALEA ) 444 $ CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR ) 445 DO 10 I = 1, N 446 BWORK( I ) = SELECT( W( I ) ) 447 10 CONTINUE 448* 449* Reorder eigenvalues, transform Schur vectors, and compute 450* reciprocal condition numbers 451* (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM) 452* otherwise, need none ) 453* (RWorkspace: none) 454* 455 CALL CTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM, 456 $ RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1, 457 $ ICOND ) 458 IF( .NOT.WANTSN ) 459 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 460 IF( ICOND.EQ.-14 ) THEN 461* 462* Not enough complex workspace 463* 464 INFO = -15 465 END IF 466 END IF 467* 468 IF( WANTVS ) THEN 469* 470* Undo balancing 471* (CWorkspace: none) 472* (RWorkspace: need N) 473* 474 CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS, 475 $ IERR ) 476 END IF 477* 478 IF( SCALEA ) THEN 479* 480* Undo scaling for the Schur form of A 481* 482 CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR ) 483 CALL CCOPY( N, A, LDA+1, W, 1 ) 484 IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN 485 DUM( 1 ) = RCONDV 486 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 487 RCONDV = DUM( 1 ) 488 END IF 489 END IF 490* 491 WORK( 1 ) = MAXWRK 492 RETURN 493* 494* End of CGEESX 495* 496 END 497