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