1 SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, 2 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO ) 3* 4* -- LAPACK auxiliary routine (instrumented to count operations) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* September 30, 1994 8* 9* .. Scalar Arguments .. 10 LOGICAL NOINIT, RIGHTV 11 INTEGER INFO, LDB, LDH, N 12 DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR 13* .. 14* .. Array Arguments .. 15 DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ), 16 $ WORK( * ) 17* .. 18* Common block to return operation count. 19* .. Common blocks .. 20 COMMON / LATIME / OPS, ITCNT 21* .. 22* .. Scalars in Common .. 23 DOUBLE PRECISION ITCNT, OPS 24* .. 25* 26* Purpose 27* ======= 28* 29* DLAEIN uses inverse iteration to find a right or left eigenvector 30* corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg 31* matrix H. 32* 33* Arguments 34* ========= 35* 36* RIGHTV (input) LOGICAL 37* = .TRUE. : compute right eigenvector; 38* = .FALSE.: compute left eigenvector. 39* 40* NOINIT (input) LOGICAL 41* = .TRUE. : no initial vector supplied in (VR,VI). 42* = .FALSE.: initial vector supplied in (VR,VI). 43* 44* N (input) INTEGER 45* The order of the matrix H. N >= 0. 46* 47* H (input) DOUBLE PRECISION array, dimension (LDH,N) 48* The upper Hessenberg matrix H. 49* 50* LDH (input) INTEGER 51* The leading dimension of the array H. LDH >= max(1,N). 52* 53* WR (input) DOUBLE PRECISION 54* WI (input) DOUBLE PRECISION 55* The real and imaginary parts of the eigenvalue of H whose 56* corresponding right or left eigenvector is to be computed. 57* 58* VR (input/output) DOUBLE PRECISION array, dimension (N) 59* VI (input/output) DOUBLE PRECISION array, dimension (N) 60* On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain 61* a real starting vector for inverse iteration using the real 62* eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI 63* must contain the real and imaginary parts of a complex 64* starting vector for inverse iteration using the complex 65* eigenvalue (WR,WI); otherwise VR and VI need not be set. 66* On exit, if WI = 0.0 (real eigenvalue), VR contains the 67* computed real eigenvector; if WI.ne.0.0 (complex eigenvalue), 68* VR and VI contain the real and imaginary parts of the 69* computed complex eigenvector. The eigenvector is normalized 70* so that the component of largest magnitude has magnitude 1; 71* here the magnitude of a complex number (x,y) is taken to be 72* |x| + |y|. 73* VI is not referenced if WI = 0.0. 74* 75* B (workspace) DOUBLE PRECISION array, dimension (LDB,N) 76* 77* LDB (input) INTEGER 78* The leading dimension of the array B. LDB >= N+1. 79* 80* WORK (workspace) DOUBLE PRECISION array, dimension (N) 81* 82* EPS3 (input) DOUBLE PRECISION 83* A small machine-dependent value which is used to perturb 84* close eigenvalues, and to replace zero pivots. 85* 86* SMLNUM (input) DOUBLE PRECISION 87* A machine-dependent value close to the underflow threshold. 88* 89* BIGNUM (input) DOUBLE PRECISION 90* A machine-dependent value close to the overflow threshold. 91* 92* INFO (output) INTEGER 93* = 0: successful exit 94* = 1: inverse iteration did not converge; VR is set to the 95* last iterate, and so is VI if WI.ne.0.0. 96* 97* ===================================================================== 98* 99* .. Parameters .. 100 DOUBLE PRECISION ZERO, ONE, TENTH 101 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 ) 102* .. 103* .. Local Scalars .. 104 CHARACTER NORMIN, TRANS 105 INTEGER I, I1, I2, I3, IERR, ITS, J 106 DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML, 107 $ OPST, REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, 108 $ VNORM, W, W1, X, XI, XR, Y 109* .. 110* .. External Functions .. 111 INTEGER IDAMAX 112 DOUBLE PRECISION DASUM, DLAPY2, DNRM2 113 EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2 114* .. 115* .. External Subroutines .. 116 EXTERNAL DLADIV, DLATRS, DSCAL 117* .. 118* .. Intrinsic Functions .. 119 INTRINSIC ABS, DBLE, MAX, SQRT 120* .. 121* .. Executable Statements .. 122* 123 INFO = 0 124*** 125* Initialize 126 OPST = 0 127*** 128* 129* GROWTO is the threshold used in the acceptance test for an 130* eigenvector. 131* 132 ROOTN = SQRT( DBLE( N ) ) 133 GROWTO = TENTH / ROOTN 134 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM 135*** 136* Increment op count for computing ROOTN, GROWTO and NRMSML 137 OPST = OPST + 4 138*** 139* 140* Form B = H - (WR,WI)*I (except that the subdiagonal elements and 141* the imaginary parts of the diagonal elements are not stored). 142* 143 DO 20 J = 1, N 144 DO 10 I = 1, J - 1 145 B( I, J ) = H( I, J ) 146 10 CONTINUE 147 B( J, J ) = H( J, J ) - WR 148 20 CONTINUE 149*** 150 OPST = OPST + N 151*** 152* 153 IF( WI.EQ.ZERO ) THEN 154* 155* Real eigenvalue. 156* 157 IF( NOINIT ) THEN 158* 159* Set initial vector. 160* 161 DO 30 I = 1, N 162 VR( I ) = EPS3 163 30 CONTINUE 164 ELSE 165* 166* Scale supplied initial vector. 167* 168 VNORM = DNRM2( N, VR, 1 ) 169 CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR, 170 $ 1 ) 171*** 172 OPST = OPST + ( 3*N+2 ) 173*** 174 END IF 175* 176 IF( RIGHTV ) THEN 177* 178* LU decomposition with partial pivoting of B, replacing zero 179* pivots by EPS3. 180* 181 DO 60 I = 1, N - 1 182 EI = H( I+1, I ) 183 IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN 184* 185* Interchange rows and eliminate. 186* 187 X = B( I, I ) / EI 188 B( I, I ) = EI 189 DO 40 J = I + 1, N 190 TEMP = B( I+1, J ) 191 B( I+1, J ) = B( I, J ) - X*TEMP 192 B( I, J ) = TEMP 193 40 CONTINUE 194 ELSE 195* 196* Eliminate without interchange. 197* 198 IF( B( I, I ).EQ.ZERO ) 199 $ B( I, I ) = EPS3 200 X = EI / B( I, I ) 201 IF( X.NE.ZERO ) THEN 202 DO 50 J = I + 1, N 203 B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 204 50 CONTINUE 205 END IF 206 END IF 207 60 CONTINUE 208 IF( B( N, N ).EQ.ZERO ) 209 $ B( N, N ) = EPS3 210*** 211* Increment op count for LU decomposition 212 OPS = OPS + ( N-1 )*( N+1 ) 213*** 214* 215 TRANS = 'N' 216* 217 ELSE 218* 219* UL decomposition with partial pivoting of B, replacing zero 220* pivots by EPS3. 221* 222 DO 90 J = N, 2, -1 223 EJ = H( J, J-1 ) 224 IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN 225* 226* Interchange columns and eliminate. 227* 228 X = B( J, J ) / EJ 229 B( J, J ) = EJ 230 DO 70 I = 1, J - 1 231 TEMP = B( I, J-1 ) 232 B( I, J-1 ) = B( I, J ) - X*TEMP 233 B( I, J ) = TEMP 234 70 CONTINUE 235 ELSE 236* 237* Eliminate without interchange. 238* 239 IF( B( J, J ).EQ.ZERO ) 240 $ B( J, J ) = EPS3 241 X = EJ / B( J, J ) 242 IF( X.NE.ZERO ) THEN 243 DO 80 I = 1, J - 1 244 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 245 80 CONTINUE 246 END IF 247 END IF 248 90 CONTINUE 249 IF( B( 1, 1 ).EQ.ZERO ) 250 $ B( 1, 1 ) = EPS3 251*** 252* Increment op count for UL decomposition 253 OPS = OPS + ( N-1 )*( N+1 ) 254*** 255* 256 TRANS = 'T' 257* 258 END IF 259* 260 NORMIN = 'N' 261 DO 110 ITS = 1, N 262* 263* Solve U*x = scale*v for a right eigenvector 264* or U'*x = scale*v for a left eigenvector, 265* overwriting x on v. 266* 267 CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, 268 $ VR, SCALE, WORK, IERR ) 269*** 270* Increment opcount for triangular solver, assuming that 271* ops DLATRS = ops DTRSV, with no scaling in DLATRS. 272 OPS = OPS + N*N 273*** 274 NORMIN = 'Y' 275* 276* Test for sufficient growth in the norm of v. 277* 278 VNORM = DASUM( N, VR, 1 ) 279*** 280 OPST = OPST + N 281*** 282 IF( VNORM.GE.GROWTO*SCALE ) 283 $ GO TO 120 284* 285* Choose new orthogonal starting vector and try again. 286* 287 TEMP = EPS3 / ( ROOTN+ONE ) 288 VR( 1 ) = EPS3 289 DO 100 I = 2, N 290 VR( I ) = TEMP 291 100 CONTINUE 292 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN 293*** 294 OPST = OPST + 4 295*** 296 110 CONTINUE 297* 298* Failure to find eigenvector in N iterations. 299* 300 INFO = 1 301* 302 120 CONTINUE 303* 304* Normalize eigenvector. 305* 306 I = IDAMAX( N, VR, 1 ) 307 CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 ) 308*** 309 OPST = OPST + ( 2*N+1 ) 310*** 311 ELSE 312* 313* Complex eigenvalue. 314* 315 IF( NOINIT ) THEN 316* 317* Set initial vector. 318* 319 DO 130 I = 1, N 320 VR( I ) = EPS3 321 VI( I ) = ZERO 322 130 CONTINUE 323 ELSE 324* 325* Scale supplied initial vector. 326* 327 NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) ) 328 REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML ) 329 CALL DSCAL( N, REC, VR, 1 ) 330 CALL DSCAL( N, REC, VI, 1 ) 331*** 332 OPST = OPST + ( 6*N+5 ) 333*** 334 END IF 335* 336 IF( RIGHTV ) THEN 337* 338* LU decomposition with partial pivoting of B, replacing zero 339* pivots by EPS3. 340* 341* The imaginary part of the (i,j)-th element of U is stored in 342* B(j+1,i). 343* 344 B( 2, 1 ) = -WI 345 DO 140 I = 2, N 346 B( I+1, 1 ) = ZERO 347 140 CONTINUE 348* 349 DO 170 I = 1, N - 1 350 ABSBII = DLAPY2( B( I, I ), B( I+1, I ) ) 351 EI = H( I+1, I ) 352 IF( ABSBII.LT.ABS( EI ) ) THEN 353* 354* Interchange rows and eliminate. 355* 356 XR = B( I, I ) / EI 357 XI = B( I+1, I ) / EI 358 B( I, I ) = EI 359 B( I+1, I ) = ZERO 360 DO 150 J = I + 1, N 361 TEMP = B( I+1, J ) 362 B( I+1, J ) = B( I, J ) - XR*TEMP 363 B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP 364 B( I, J ) = TEMP 365 B( J+1, I ) = ZERO 366 150 CONTINUE 367 B( I+2, I ) = -WI 368 B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI 369 B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI 370*** 371 OPST = OPST + ( 4*( N-I )+6 ) 372*** 373 ELSE 374* 375* Eliminate without interchanging rows. 376* 377 IF( ABSBII.EQ.ZERO ) THEN 378 B( I, I ) = EPS3 379 B( I+1, I ) = ZERO 380 ABSBII = EPS3 381 END IF 382 EI = ( EI / ABSBII ) / ABSBII 383 XR = B( I, I )*EI 384 XI = -B( I+1, I )*EI 385 DO 160 J = I + 1, N 386 B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) + 387 $ XI*B( J+1, I ) 388 B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J ) 389 160 CONTINUE 390 B( I+2, I+1 ) = B( I+2, I+1 ) - WI 391*** 392 OPST = OPST + ( 7*( N-I )+4 ) 393*** 394 END IF 395* 396* Compute 1-norm of offdiagonal elements of i-th row. 397* 398 WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) + 399 $ DASUM( N-I, B( I+2, I ), 1 ) 400*** 401 OPST = OPST + ( 2*( N-I )+4 ) 402*** 403 170 CONTINUE 404 IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO ) 405 $ B( N, N ) = EPS3 406 WORK( N ) = ZERO 407* 408 I1 = N 409 I2 = 1 410 I3 = -1 411 ELSE 412* 413* UL decomposition with partial pivoting of conjg(B), 414* replacing zero pivots by EPS3. 415* 416* The imaginary part of the (i,j)-th element of U is stored in 417* B(j+1,i). 418* 419 B( N+1, N ) = WI 420 DO 180 J = 1, N - 1 421 B( N+1, J ) = ZERO 422 180 CONTINUE 423* 424 DO 210 J = N, 2, -1 425 EJ = H( J, J-1 ) 426 ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) ) 427 IF( ABSBJJ.LT.ABS( EJ ) ) THEN 428* 429* Interchange columns and eliminate 430* 431 XR = B( J, J ) / EJ 432 XI = B( J+1, J ) / EJ 433 B( J, J ) = EJ 434 B( J+1, J ) = ZERO 435 DO 190 I = 1, J - 1 436 TEMP = B( I, J-1 ) 437 B( I, J-1 ) = B( I, J ) - XR*TEMP 438 B( J, I ) = B( J+1, I ) - XI*TEMP 439 B( I, J ) = TEMP 440 B( J+1, I ) = ZERO 441 190 CONTINUE 442 B( J+1, J-1 ) = WI 443 B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI 444 B( J, J-1 ) = B( J, J-1 ) - XR*WI 445*** 446 OPST = OPST + ( 4*( J-1 )+6 ) 447*** 448 ELSE 449* 450* Eliminate without interchange. 451* 452 IF( ABSBJJ.EQ.ZERO ) THEN 453 B( J, J ) = EPS3 454 B( J+1, J ) = ZERO 455 ABSBJJ = EPS3 456 END IF 457 EJ = ( EJ / ABSBJJ ) / ABSBJJ 458 XR = B( J, J )*EJ 459 XI = -B( J+1, J )*EJ 460 DO 200 I = 1, J - 1 461 B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) + 462 $ XI*B( J+1, I ) 463 B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J ) 464 200 CONTINUE 465 B( J, J-1 ) = B( J, J-1 ) + WI 466*** 467 OPST = OPST + ( 7*( J-1 )+4 ) 468*** 469 END IF 470* 471* Compute 1-norm of offdiagonal elements of j-th column. 472* 473 WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) + 474 $ DASUM( J-1, B( J+1, 1 ), LDB ) 475*** 476 OPST = OPST + ( 2*( J-1 )+4 ) 477*** 478 210 CONTINUE 479 IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO ) 480 $ B( 1, 1 ) = EPS3 481 WORK( 1 ) = ZERO 482* 483 I1 = 1 484 I2 = N 485 I3 = 1 486 END IF 487* 488 DO 270 ITS = 1, N 489 SCALE = ONE 490 VMAX = ONE 491 VCRIT = BIGNUM 492* 493* Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector, 494* or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector, 495* overwriting (xr,xi) on (vr,vi). 496* 497 DO 250 I = I1, I2, I3 498* 499 IF( WORK( I ).GT.VCRIT ) THEN 500 REC = ONE / VMAX 501 CALL DSCAL( N, REC, VR, 1 ) 502 CALL DSCAL( N, REC, VI, 1 ) 503 SCALE = SCALE*REC 504 VMAX = ONE 505 VCRIT = BIGNUM 506 END IF 507* 508 XR = VR( I ) 509 XI = VI( I ) 510 IF( RIGHTV ) THEN 511 DO 220 J = I + 1, N 512 XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J ) 513 XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J ) 514 220 CONTINUE 515 ELSE 516 DO 230 J = 1, I - 1 517 XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J ) 518 XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J ) 519 230 CONTINUE 520 END IF 521* 522 W = ABS( B( I, I ) ) + ABS( B( I+1, I ) ) 523 IF( W.GT.SMLNUM ) THEN 524 IF( W.LT.ONE ) THEN 525 W1 = ABS( XR ) + ABS( XI ) 526 IF( W1.GT.W*BIGNUM ) THEN 527 REC = ONE / W1 528 CALL DSCAL( N, REC, VR, 1 ) 529 CALL DSCAL( N, REC, VI, 1 ) 530 XR = VR( I ) 531 XI = VI( I ) 532 SCALE = SCALE*REC 533 VMAX = VMAX*REC 534 END IF 535 END IF 536* 537* Divide by diagonal element of B. 538* 539 CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ), 540 $ VI( I ) ) 541 VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX ) 542 VCRIT = BIGNUM / VMAX 543*** 544 OPST = OPST + 9 545*** 546 ELSE 547 DO 240 J = 1, N 548 VR( J ) = ZERO 549 VI( J ) = ZERO 550 240 CONTINUE 551 VR( I ) = ONE 552 VI( I ) = ONE 553 SCALE = ZERO 554 VMAX = ONE 555 VCRIT = BIGNUM 556 END IF 557 250 CONTINUE 558*** 559* Increment op count for loop 260, assuming no scaling 560 OPS = OPS + 4*N*( N-1 ) 561*** 562* 563* Test for sufficient growth in the norm of (VR,VI). 564* 565 VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 ) 566*** 567 OPST = OPST + 2*N 568*** 569 IF( VNORM.GE.GROWTO*SCALE ) 570 $ GO TO 280 571* 572* Choose a new orthogonal starting vector and try again. 573* 574 Y = EPS3 / ( ROOTN+ONE ) 575 VR( 1 ) = EPS3 576 VI( 1 ) = ZERO 577* 578 DO 260 I = 2, N 579 VR( I ) = Y 580 VI( I ) = ZERO 581 260 CONTINUE 582 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN 583*** 584 OPST = OPST + 4 585*** 586 270 CONTINUE 587* 588* Failure to find eigenvector in N iterations 589* 590 INFO = 1 591* 592 280 CONTINUE 593* 594* Normalize eigenvector. 595* 596 VNORM = ZERO 597 DO 290 I = 1, N 598 VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) ) 599 290 CONTINUE 600 CALL DSCAL( N, ONE / VNORM, VR, 1 ) 601 CALL DSCAL( N, ONE / VNORM, VI, 1 ) 602*** 603 OPST = OPST + ( 4*N+1 ) 604*** 605* 606 END IF 607* 608*** 609* Compute final op count 610 OPS = OPS + OPST 611*** 612 RETURN 613* 614* End of DLAEIN 615* 616 END 617