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