Annotation of rpl/lapack/lapack/dlaein.f, revision 1.1
1.1 ! bertrand 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) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack