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