1 SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 2 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO ) 3* 4* -- LAPACK routine (instrumented to count operations, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER HOWMNY, SIDE 11 INTEGER INFO, LDA, LDB, LDVL, LDVR, M, MM, N 12* .. 13* .. Array Arguments .. 14 LOGICAL SELECT( * ) 15 REAL RWORK( * ) 16 COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 17 $ VR( LDVR, * ), WORK( * ) 18* .. 19* 20* ----------------------- Begin Timing Code ------------------------ 21* Common block to return operation count and iteration count 22* ITCNT is initialized to 0, OPS is only incremented 23* OPST is used to accumulate small contributions to OPS 24* to avoid roundoff error 25* .. Common blocks .. 26 COMMON / LATIME / OPS, ITCNT 27* .. 28* .. Scalars in Common .. 29 REAL ITCNT, OPS 30* .. 31* ------------------------ End Timing Code ------------------------- 32* 33* 34* Purpose 35* ======= 36* 37* CTGEVC computes some or all of the right and/or left generalized 38* eigenvectors of a pair of complex upper triangular matrices (A,B). 39* 40* The right generalized eigenvector x and the left generalized 41* eigenvector y of (A,B) corresponding to a generalized eigenvalue 42* w are defined by: 43* 44* (A - wB) * x = 0 and y**H * (A - wB) = 0 45* 46* where y**H denotes the conjugate tranpose of y. 47* 48* If an eigenvalue w is determined by zero diagonal elements of both A 49* and B, a unit vector is returned as the corresponding eigenvector. 50* 51* If all eigenvectors are requested, the routine may either return 52* the matrices X and/or Y of right or left eigenvectors of (A,B), or 53* the products Z*X and/or Q*Y, where Z and Q are input unitary 54* matrices. If (A,B) was obtained from the generalized Schur 55* factorization of an original pair of matrices 56* (A0,B0) = (Q*A*Z**H,Q*B*Z**H), 57* then Z*X and Q*Y are the matrices of right or left eigenvectors of 58* A. 59* 60* Arguments 61* ========= 62* 63* SIDE (input) CHARACTER*1 64* = 'R': compute right eigenvectors only; 65* = 'L': compute left eigenvectors only; 66* = 'B': compute both right and left eigenvectors. 67* 68* HOWMNY (input) CHARACTER*1 69* = 'A': compute all right and/or left eigenvectors; 70* = 'B': compute all right and/or left eigenvectors, and 71* backtransform them using the input matrices supplied 72* in VR and/or VL; 73* = 'S': compute selected right and/or left eigenvectors, 74* specified by the logical array SELECT. 75* 76* SELECT (input) LOGICAL array, dimension (N) 77* If HOWMNY='S', SELECT specifies the eigenvectors to be 78* computed. 79* If HOWMNY='A' or 'B', SELECT is not referenced. 80* To select the eigenvector corresponding to the j-th 81* eigenvalue, SELECT(j) must be set to .TRUE.. 82* 83* N (input) INTEGER 84* The order of the matrices A and B. N >= 0. 85* 86* A (input) COMPLEX array, dimension (LDA,N) 87* The upper triangular matrix A. 88* 89* LDA (input) INTEGER 90* The leading dimension of array A. LDA >= max(1,N). 91* 92* B (input) COMPLEX array, dimension (LDB,N) 93* The upper triangular matrix B. B must have real diagonal 94* elements. 95* 96* LDB (input) INTEGER 97* The leading dimension of array B. LDB >= max(1,N). 98* 99* VL (input/output) COMPLEX array, dimension (LDVL,MM) 100* On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 101* contain an N-by-N matrix Q (usually the unitary matrix Q 102* of left Schur vectors returned by CHGEQZ). 103* On exit, if SIDE = 'L' or 'B', VL contains: 104* if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B); 105* if HOWMNY = 'B', the matrix Q*Y; 106* if HOWMNY = 'S', the left eigenvectors of (A,B) specified by 107* SELECT, stored consecutively in the columns of 108* VL, in the same order as their eigenvalues. 109* If SIDE = 'R', VL is not referenced. 110* 111* LDVL (input) INTEGER 112* The leading dimension of array VL. 113* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. 114* 115* VR (input/output) COMPLEX array, dimension (LDVR,MM) 116* On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 117* contain an N-by-N matrix Q (usually the unitary matrix Z 118* of right Schur vectors returned by CHGEQZ). 119* On exit, if SIDE = 'R' or 'B', VR contains: 120* if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B); 121* if HOWMNY = 'B', the matrix Z*X; 122* if HOWMNY = 'S', the right eigenvectors of (A,B) specified by 123* SELECT, stored consecutively in the columns of 124* VR, in the same order as their eigenvalues. 125* If SIDE = 'L', VR is not referenced. 126* 127* LDVR (input) INTEGER 128* The leading dimension of the array VR. 129* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 130* 131* MM (input) INTEGER 132* The leading dimension of the array VR. 133* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 134* 135* M (output) INTEGER 136* The number of columns in the arrays VL and/or VR actually 137* used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 138* is set to N. Each selected eigenvector occupies one column. 139* 140* WORK (workspace) COMPLEX array, dimension (2*N) 141* 142* RWORK (workspace) REAL array, dimension (2*N) 143* 144* INFO (output) INTEGER 145* = 0: successful exit. 146* < 0: if INFO = -i, the i-th argument had an illegal value. 147* 148* ===================================================================== 149* 150* .. Parameters .. 151 REAL ZERO, ONE 152 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 153 COMPLEX CZERO, CONE 154 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 155 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 156* .. 157* .. Local Scalars .. 158 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP, 159 $ LSA, LSB 160 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, IOPST, ISIDE, 161 $ ISRC, J, JE, JR 162 REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG, 163 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA, 164 $ SCALE, SMALL, TEMP, ULP, XMAX 165 COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X 166* .. 167* .. External Functions .. 168 LOGICAL LSAME 169 REAL SLAMCH 170 COMPLEX CLADIV 171 EXTERNAL CLADIV, LSAME, SLAMCH 172* .. 173* .. External Subroutines .. 174 EXTERNAL CGEMV, SLABAD, XERBLA 175* .. 176* .. Intrinsic Functions .. 177 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL 178* .. 179* .. Statement Functions .. 180 REAL ABS1 181* .. 182* .. Statement Function definitions .. 183 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 184* .. 185* .. Executable Statements .. 186* 187* Decode and Test the input parameters 188* 189 IF( LSAME( HOWMNY, 'A' ) ) THEN 190 IHWMNY = 1 191 ILALL = .TRUE. 192 ILBACK = .FALSE. 193 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 194 IHWMNY = 2 195 ILALL = .FALSE. 196 ILBACK = .FALSE. 197 ELSE IF( LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'T' ) ) THEN 198 IHWMNY = 3 199 ILALL = .TRUE. 200 ILBACK = .TRUE. 201 ELSE 202 IHWMNY = -1 203 END IF 204* 205 IF( LSAME( SIDE, 'R' ) ) THEN 206 ISIDE = 1 207 COMPL = .FALSE. 208 COMPR = .TRUE. 209 ELSE IF( LSAME( SIDE, 'L' ) ) THEN 210 ISIDE = 2 211 COMPL = .TRUE. 212 COMPR = .FALSE. 213 ELSE IF( LSAME( SIDE, 'B' ) ) THEN 214 ISIDE = 3 215 COMPL = .TRUE. 216 COMPR = .TRUE. 217 ELSE 218 ISIDE = -1 219 END IF 220* 221 INFO = 0 222 IF( ISIDE.LT.0 ) THEN 223 INFO = -1 224 ELSE IF( IHWMNY.LT.0 ) THEN 225 INFO = -2 226 ELSE IF( N.LT.0 ) THEN 227 INFO = -4 228 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 229 INFO = -6 230 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 231 INFO = -8 232 END IF 233 IF( INFO.NE.0 ) THEN 234 CALL XERBLA( 'CTGEVC', -INFO ) 235 RETURN 236 END IF 237* 238* Count the number of eigenvectors 239* 240 IF( .NOT.ILALL ) THEN 241 IM = 0 242 DO 10 J = 1, N 243 IF( SELECT( J ) ) 244 $ IM = IM + 1 245 10 CONTINUE 246 ELSE 247 IM = N 248 END IF 249* 250* Check diagonal of B 251* 252 ILBBAD = .FALSE. 253 DO 20 J = 1, N 254 IF( AIMAG( B( J, J ) ).NE.ZERO ) 255 $ ILBBAD = .TRUE. 256 20 CONTINUE 257* 258 IF( ILBBAD ) THEN 259 INFO = -7 260 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 261 INFO = -10 262 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 263 INFO = -12 264 ELSE IF( MM.LT.IM ) THEN 265 INFO = -13 266 END IF 267 IF( INFO.NE.0 ) THEN 268 CALL XERBLA( 'CTGEVC', -INFO ) 269 RETURN 270 END IF 271* 272* Quick return if possible 273* 274 M = IM 275 IF( N.EQ.0 ) 276 $ RETURN 277* 278* Machine Constants 279* 280 SAFMIN = SLAMCH( 'Safe minimum' ) 281 BIG = ONE / SAFMIN 282 CALL SLABAD( SAFMIN, BIG ) 283 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 284 SMALL = SAFMIN*N / ULP 285 BIG = ONE / SMALL 286 BIGNUM = ONE / ( SAFMIN*N ) 287* 288* Compute the 1-norm of each column of the strictly upper triangular 289* part of A and B to check for possible overflow in the triangular 290* solver. 291* 292 ANORM = ABS1( A( 1, 1 ) ) 293 BNORM = ABS1( B( 1, 1 ) ) 294 RWORK( 1 ) = ZERO 295 RWORK( N+1 ) = ZERO 296 DO 40 J = 2, N 297 RWORK( J ) = ZERO 298 RWORK( N+J ) = ZERO 299 DO 30 I = 1, J - 1 300 RWORK( J ) = RWORK( J ) + ABS1( A( I, J ) ) 301 RWORK( N+J ) = RWORK( N+J ) + ABS1( B( I, J ) ) 302 30 CONTINUE 303 ANORM = MAX( ANORM, RWORK( J )+ABS1( A( J, J ) ) ) 304 BNORM = MAX( BNORM, RWORK( N+J )+ABS1( B( J, J ) ) ) 305 40 CONTINUE 306* 307 ASCALE = ONE / MAX( ANORM, SAFMIN ) 308 BSCALE = ONE / MAX( BNORM, SAFMIN ) 309* ---------------------- Begin Timing Code ------------------------- 310 OPS = OPS + REAL( 2*N**2+2*N+6 ) 311* ----------------------- End Timing Code -------------------------- 312* 313* Left eigenvectors 314* 315 IF( COMPL ) THEN 316 IEIG = 0 317* 318* Main loop over eigenvalues 319* 320 DO 140 JE = 1, N 321 IF( ILALL ) THEN 322 ILCOMP = .TRUE. 323 ELSE 324 ILCOMP = SELECT( JE ) 325 END IF 326 IF( ILCOMP ) THEN 327 IEIG = IEIG + 1 328* 329 IF( ABS1( A( JE, JE ) ).LE.SAFMIN .AND. 330 $ ABS( REAL( B( JE, JE ) ) ).LE.SAFMIN ) THEN 331* 332* Singular matrix pencil -- return unit eigenvector 333* 334 DO 50 JR = 1, N 335 VL( JR, IEIG ) = CZERO 336 50 CONTINUE 337 VL( IEIG, IEIG ) = CONE 338 GO TO 140 339 END IF 340* 341* Non-singular eigenvalue: 342* Compute coefficients a and b in 343* H 344* y ( a A - b B ) = 0 345* 346 TEMP = ONE / MAX( ABS1( A( JE, JE ) )*ASCALE, 347 $ ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN ) 348 SALPHA = ( TEMP*A( JE, JE ) )*ASCALE 349 SBETA = ( TEMP*REAL( B( JE, JE ) ) )*BSCALE 350 ACOEFF = SBETA*ASCALE 351 BCOEFF = SALPHA*BSCALE 352* 353* Scale to avoid underflow 354* 355 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL 356 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. 357 $ SMALL 358* 359 SCALE = ONE 360 IF( LSA ) 361 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 362 IF( LSB ) 363 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* 364 $ MIN( BNORM, BIG ) ) 365 IF( LSA .OR. LSB ) THEN 366 SCALE = MIN( SCALE, ONE / 367 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), 368 $ ABS1( BCOEFF ) ) ) ) 369 IF( LSA ) THEN 370 ACOEFF = ASCALE*( SCALE*SBETA ) 371 ELSE 372 ACOEFF = SCALE*ACOEFF 373 END IF 374 IF( LSB ) THEN 375 BCOEFF = BSCALE*( SCALE*SALPHA ) 376 ELSE 377 BCOEFF = SCALE*BCOEFF 378 END IF 379* ----------------- Begin Timing Code ------------------ 380* Calculation of SALPHA through DMIN 381 IOPST = 34 382 ELSE 383 IOPST = 20 384* ------------------ End Timing Code ------------------- 385 END IF 386* 387 ACOEFA = ABS( ACOEFF ) 388 BCOEFA = ABS1( BCOEFF ) 389 XMAX = ONE 390 DO 60 JR = 1, N 391 WORK( JR ) = CZERO 392 60 CONTINUE 393 WORK( JE ) = CONE 394 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 395* 396* H 397* Triangular solve of (a A - b B) y = 0 398* 399* H 400* (rowwise in (a A - b B) , or columnwise in a A - b B) 401* 402 DO 100 J = JE + 1, N 403* 404* Compute 405* j-1 406* SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k) 407* k=je 408* (Scale if necessary) 409* 410 TEMP = ONE / XMAX 411 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM* 412 $ TEMP ) THEN 413 DO 70 JR = JE, J - 1 414 WORK( JR ) = TEMP*WORK( JR ) 415 70 CONTINUE 416 XMAX = ONE 417* ---------------- Begin Timing Code ---------------- 418 IOPST = IOPST + 2*( J-JE ) 419* ----------------- End Timing Code ----------------- 420 END IF 421 SUMA = CZERO 422 SUMB = CZERO 423* 424 DO 80 JR = JE, J - 1 425 SUMA = SUMA + CONJG( A( JR, J ) )*WORK( JR ) 426 SUMB = SUMB + CONJG( B( JR, J ) )*WORK( JR ) 427 80 CONTINUE 428 SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB 429* 430* Form x(j) = - SUM / conjg( a*A(j,j) - b*B(j,j) ) 431* 432* with scaling and perturbation of the denominator 433* 434 D = CONJG( ACOEFF*A( J, J )-BCOEFF*B( J, J ) ) 435 IF( ABS1( D ).LE.DMIN ) 436 $ D = CMPLX( DMIN ) 437* 438 IF( ABS1( D ).LT.ONE ) THEN 439 IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN 440 TEMP = ONE / ABS1( SUM ) 441 DO 90 JR = JE, J - 1 442 WORK( JR ) = TEMP*WORK( JR ) 443 90 CONTINUE 444 XMAX = TEMP*XMAX 445 SUM = TEMP*SUM 446* -------------- Begin Timing Code --------------- 447 IOPST = IOPST + 2*( J-JE ) + 5 448* --------------- End Timing Code ---------------- 449 END IF 450 END IF 451 WORK( J ) = CLADIV( -SUM, D ) 452 XMAX = MAX( XMAX, ABS1( WORK( J ) ) ) 453 100 CONTINUE 454* 455* Back transform eigenvector if HOWMNY='B'. 456* 457 IF( ILBACK ) THEN 458 CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL, 459 $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 ) 460 ISRC = 2 461 IBEG = 1 462* ----------------- Begin Timing Code ------------------ 463 IOPST = IOPST + 8*N*( N+1-JE ) 464* ------------------ End Timing Code ------------------- 465 ELSE 466 ISRC = 1 467 IBEG = JE 468 END IF 469* 470* Copy and scale eigenvector into column of VL 471* 472 XMAX = ZERO 473 DO 110 JR = IBEG, N 474 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) 475 110 CONTINUE 476* 477 IF( XMAX.GT.SAFMIN ) THEN 478 TEMP = ONE / XMAX 479 DO 120 JR = IBEG, N 480 VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) 481 120 CONTINUE 482 ELSE 483 IBEG = N + 1 484 END IF 485* 486 DO 130 JR = 1, IBEG - 1 487 VL( JR, IEIG ) = CZERO 488 130 CONTINUE 489* 490* ------------------- Begin Timing Code ------------------- 491 OPS = OPS + ( REAL( 36*( N-JE )+8*( N-JE )*( N+1-JE )+3* 492 $ ( N+1-IBEG )+1 )+REAL( IOPST ) ) 493* -------------------- End Timing Code -------------------- 494 END IF 495 140 CONTINUE 496 END IF 497* 498* Right eigenvectors 499* 500 IF( COMPR ) THEN 501 IEIG = IM + 1 502* 503* Main loop over eigenvalues 504* 505 DO 250 JE = N, 1, -1 506 IF( ILALL ) THEN 507 ILCOMP = .TRUE. 508 ELSE 509 ILCOMP = SELECT( JE ) 510 END IF 511 IF( ILCOMP ) THEN 512 IEIG = IEIG - 1 513* 514 IF( ABS1( A( JE, JE ) ).LE.SAFMIN .AND. 515 $ ABS( REAL( B( JE, JE ) ) ).LE.SAFMIN ) THEN 516* 517* Singular matrix pencil -- return unit eigenvector 518* 519 DO 150 JR = 1, N 520 VR( JR, IEIG ) = CZERO 521 150 CONTINUE 522 VR( IEIG, IEIG ) = CONE 523 GO TO 250 524 END IF 525* 526* Non-singular eigenvalue: 527* Compute coefficients a and b in 528* 529* ( a A - b B ) x = 0 530* 531 TEMP = ONE / MAX( ABS1( A( JE, JE ) )*ASCALE, 532 $ ABS( REAL( B( JE, JE ) ) )*BSCALE, SAFMIN ) 533 SALPHA = ( TEMP*A( JE, JE ) )*ASCALE 534 SBETA = ( TEMP*REAL( B( JE, JE ) ) )*BSCALE 535 ACOEFF = SBETA*ASCALE 536 BCOEFF = SALPHA*BSCALE 537* 538* Scale to avoid underflow 539* 540 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL 541 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. 542 $ SMALL 543* 544 SCALE = ONE 545 IF( LSA ) 546 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 547 IF( LSB ) 548 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* 549 $ MIN( BNORM, BIG ) ) 550 IF( LSA .OR. LSB ) THEN 551 SCALE = MIN( SCALE, ONE / 552 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), 553 $ ABS1( BCOEFF ) ) ) ) 554 IF( LSA ) THEN 555 ACOEFF = ASCALE*( SCALE*SBETA ) 556 ELSE 557 ACOEFF = SCALE*ACOEFF 558 END IF 559 IF( LSB ) THEN 560 BCOEFF = BSCALE*( SCALE*SALPHA ) 561 ELSE 562 BCOEFF = SCALE*BCOEFF 563 END IF 564* ----------------- Begin Timing Code ------------------ 565* Calculation of SALPHA through DMIN 566 IOPST = 34 567 ELSE 568 IOPST = 20 569* ------------------ End Timing Code ------------------- 570 END IF 571* 572 ACOEFA = ABS( ACOEFF ) 573 BCOEFA = ABS1( BCOEFF ) 574 XMAX = ONE 575 DO 160 JR = 1, N 576 WORK( JR ) = CZERO 577 160 CONTINUE 578 WORK( JE ) = CONE 579 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 580* 581* Triangular solve of (a A - b B) x = 0 (columnwise) 582* 583* WORK(1:j-1) contains sums w, 584* WORK(j+1:JE) contains x 585* 586 DO 170 JR = 1, JE - 1 587 WORK( JR ) = ACOEFF*A( JR, JE ) - BCOEFF*B( JR, JE ) 588 170 CONTINUE 589 WORK( JE ) = CONE 590* 591 DO 210 J = JE - 1, 1, -1 592* 593* Form x(j) := - w(j) / d 594* with scaling and perturbation of the denominator 595* 596 D = ACOEFF*A( J, J ) - BCOEFF*B( J, J ) 597 IF( ABS1( D ).LE.DMIN ) 598 $ D = CMPLX( DMIN ) 599* 600 IF( ABS1( D ).LT.ONE ) THEN 601 IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN 602 TEMP = ONE / ABS1( WORK( J ) ) 603 DO 180 JR = 1, JE 604 WORK( JR ) = TEMP*WORK( JR ) 605 180 CONTINUE 606* -------------- Begin Timing Code --------------- 607 IOPST = IOPST + 2*JE + 5 608 ELSE 609 IOPST = IOPST + 3 610* --------------- End Timing Code ---------------- 611 END IF 612 END IF 613* 614 WORK( J ) = CLADIV( -WORK( J ), D ) 615* 616 IF( J.GT.1 ) THEN 617* 618* w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling 619* 620 IF( ABS1( WORK( J ) ).GT.ONE ) THEN 621 TEMP = ONE / ABS1( WORK( J ) ) 622 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE. 623 $ BIGNUM*TEMP ) THEN 624 DO 190 JR = 1, JE 625 WORK( JR ) = TEMP*WORK( JR ) 626 190 CONTINUE 627* ------------- Begin Timing Code ------------- 628 IOPST = IOPST + 2*JE + 6 629 ELSE 630 IOPST = IOPST + 6 631* -------------- End Timing Code -------------- 632 END IF 633 END IF 634* 635 CA = ACOEFF*WORK( J ) 636 CB = BCOEFF*WORK( J ) 637 DO 200 JR = 1, J - 1 638 WORK( JR ) = WORK( JR ) + CA*A( JR, J ) - 639 $ CB*B( JR, J ) 640 200 CONTINUE 641 END IF 642 210 CONTINUE 643* 644* Back transform eigenvector if HOWMNY='B'. 645* 646 IF( ILBACK ) THEN 647 CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1, 648 $ CZERO, WORK( N+1 ), 1 ) 649 ISRC = 2 650 IEND = N 651* ----------------- Begin Timing Code ------------------ 652 IOPST = IOPST + 8*N*JE 653* ------------------ End Timing Code ------------------- 654 ELSE 655 ISRC = 1 656 IEND = JE 657 END IF 658* 659* Copy and scale eigenvector into column of VR 660* 661 XMAX = ZERO 662 DO 220 JR = 1, IEND 663 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) 664 220 CONTINUE 665* 666 IF( XMAX.GT.SAFMIN ) THEN 667 TEMP = ONE / XMAX 668 DO 230 JR = 1, IEND 669 VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) 670 230 CONTINUE 671 ELSE 672 IEND = 0 673 END IF 674* 675 DO 240 JR = IEND + 1, N 676 VR( JR, IEIG ) = CZERO 677 240 CONTINUE 678* 679* ------------------- Begin Timing Code ------------------- 680 OPS = OPS + ( REAL( 30*( JE-2 )+8*( JE-1 )*( JE-2 )+3* 681 $ IEND+22 )+REAL( IOPST ) ) 682* -------------------- End Timing Code -------------------- 683 END IF 684 250 CONTINUE 685 END IF 686* 687 RETURN 688* 689* End of CTGEVC 690* 691 END 692