![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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 (version 3.2.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * June 2010 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: * 19: * Purpose 20: * ======= 21: * 22: * DLAEIN 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) DOUBLE PRECISION 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) DOUBLE PRECISION 47: * WI (input) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N) 52: * VI (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDB,N) 69: * 70: * LDB (input) INTEGER 71: * The leading dimension of the array B. LDB >= N+1. 72: * 73: * WORK (workspace) DOUBLE PRECISION array, dimension (N) 74: * 75: * EPS3 (input) DOUBLE PRECISION 76: * A small machine-dependent value which is used to perturb 77: * close eigenvalues, and to replace zero pivots. 78: * 79: * SMLNUM (input) DOUBLE PRECISION 80: * A machine-dependent value close to the underflow threshold. 81: * 82: * BIGNUM (input) DOUBLE PRECISION 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: DOUBLE PRECISION ZERO, ONE, TENTH 94: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 ) 95: * .. 96: * .. Local Scalars .. 97: CHARACTER NORMIN, TRANS 98: INTEGER I, I1, I2, I3, IERR, ITS, J 99: DOUBLE PRECISION 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 IDAMAX 105: DOUBLE PRECISION DASUM, DLAPY2, DNRM2 106: EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2 107: * .. 108: * .. External Subroutines .. 109: EXTERNAL DLADIV, DLATRS, DSCAL 110: * .. 111: * .. Intrinsic Functions .. 112: INTRINSIC ABS, DBLE, MAX, 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( DBLE( 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 = DNRM2( N, VR, 1 ) 151: CALL DSCAL( 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 DLATRS( '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 = DASUM( 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 = IDAMAX( N, VR, 1 ) 267: CALL DSCAL( 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 = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) ) 285: REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML ) 286: CALL DSCAL( N, REC, VR, 1 ) 287: CALL DSCAL( 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 = DLAPY2( 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 ) = DASUM( N-I, B( I, I+1 ), LDB ) + 347: $ DASUM( 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 = DLAPY2( 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 ) = DASUM( J-1, B( 1, J ), 1 ) + 413: $ DASUM( 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 DSCAL( N, REC, VR, 1 ) 438: CALL DSCAL( 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 DSCAL( N, REC, VR, 1 ) 465: CALL DSCAL( 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 DLADIV( 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 = DASUM( N, VR, 1 ) + DASUM( 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 DSCAL( N, ONE / VNORM, VR, 1 ) 524: CALL DSCAL( N, ONE / VNORM, VI, 1 ) 525: * 526: END IF 527: * 528: RETURN 529: * 530: * End of DLAEIN 531: * 532: END