diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/dlaein.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/dlaein.f')
-rw-r--r-- | SRC/dlaein.f | 531 |
1 files changed, 531 insertions, 0 deletions
diff --git a/SRC/dlaein.f b/SRC/dlaein.f new file mode 100644 index 00000000..9f9b5fa5 --- /dev/null +++ b/SRC/dlaein.f @@ -0,0 +1,531 @@ + SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, + $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + LOGICAL NOINIT, RIGHTV + INTEGER INFO, LDB, LDH, N + DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR +* .. +* .. Array Arguments .. + DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* DLAEIN uses inverse iteration to find a right or left eigenvector +* corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg +* matrix H. +* +* Arguments +* ========= +* +* RIGHTV (input) LOGICAL +* = .TRUE. : compute right eigenvector; +* = .FALSE.: compute left eigenvector. +* +* NOINIT (input) LOGICAL +* = .TRUE. : no initial vector supplied in (VR,VI). +* = .FALSE.: initial vector supplied in (VR,VI). +* +* N (input) INTEGER +* The order of the matrix H. N >= 0. +* +* H (input) DOUBLE PRECISION array, dimension (LDH,N) +* The upper Hessenberg matrix H. +* +* LDH (input) INTEGER +* The leading dimension of the array H. LDH >= max(1,N). +* +* WR (input) DOUBLE PRECISION +* WI (input) DOUBLE PRECISION +* The real and imaginary parts of the eigenvalue of H whose +* corresponding right or left eigenvector is to be computed. +* +* VR (input/output) DOUBLE PRECISION array, dimension (N) +* VI (input/output) DOUBLE PRECISION array, dimension (N) +* On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain +* a real starting vector for inverse iteration using the real +* eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI +* must contain the real and imaginary parts of a complex +* starting vector for inverse iteration using the complex +* eigenvalue (WR,WI); otherwise VR and VI need not be set. +* On exit, if WI = 0.0 (real eigenvalue), VR contains the +* computed real eigenvector; if WI.ne.0.0 (complex eigenvalue), +* VR and VI contain the real and imaginary parts of the +* computed complex eigenvector. The eigenvector is normalized +* so that the component of largest magnitude has magnitude 1; +* here the magnitude of a complex number (x,y) is taken to be +* |x| + |y|. +* VI is not referenced if WI = 0.0. +* +* B (workspace) DOUBLE PRECISION array, dimension (LDB,N) +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= N+1. +* +* WORK (workspace) DOUBLE PRECISION array, dimension (N) +* +* EPS3 (input) DOUBLE PRECISION +* A small machine-dependent value which is used to perturb +* close eigenvalues, and to replace zero pivots. +* +* SMLNUM (input) DOUBLE PRECISION +* A machine-dependent value close to the underflow threshold. +* +* BIGNUM (input) DOUBLE PRECISION +* A machine-dependent value close to the overflow threshold. +* +* INFO (output) INTEGER +* = 0: successful exit +* = 1: inverse iteration did not converge; VR is set to the +* last iterate, and so is VI if WI.ne.0.0. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TENTH + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 ) +* .. +* .. Local Scalars .. + CHARACTER NORMIN, TRANS + INTEGER I, I1, I2, I3, IERR, ITS, J + DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML, + $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W, + $ W1, X, XI, XR, Y +* .. +* .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DASUM, DLAPY2, DNRM2 + EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2 +* .. +* .. External Subroutines .. + EXTERNAL DLADIV, DLATRS, DSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, MAX, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 +* +* GROWTO is the threshold used in the acceptance test for an +* eigenvector. +* + ROOTN = SQRT( DBLE( N ) ) + GROWTO = TENTH / ROOTN + NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM +* +* Form B = H - (WR,WI)*I (except that the subdiagonal elements and +* the imaginary parts of the diagonal elements are not stored). +* + DO 20 J = 1, N + DO 10 I = 1, J - 1 + B( I, J ) = H( I, J ) + 10 CONTINUE + B( J, J ) = H( J, J ) - WR + 20 CONTINUE +* + IF( WI.EQ.ZERO ) THEN +* +* Real eigenvalue. +* + IF( NOINIT ) THEN +* +* Set initial vector. +* + DO 30 I = 1, N + VR( I ) = EPS3 + 30 CONTINUE + ELSE +* +* Scale supplied initial vector. +* + VNORM = DNRM2( N, VR, 1 ) + CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR, + $ 1 ) + END IF +* + IF( RIGHTV ) THEN +* +* LU decomposition with partial pivoting of B, replacing zero +* pivots by EPS3. +* + DO 60 I = 1, N - 1 + EI = H( I+1, I ) + IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN +* +* Interchange rows and eliminate. +* + X = B( I, I ) / EI + B( I, I ) = EI + DO 40 J = I + 1, N + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - X*TEMP + B( I, J ) = TEMP + 40 CONTINUE + ELSE +* +* Eliminate without interchange. +* + IF( B( I, I ).EQ.ZERO ) + $ B( I, I ) = EPS3 + X = EI / B( I, I ) + IF( X.NE.ZERO ) THEN + DO 50 J = I + 1, N + B( I+1, J ) = B( I+1, J ) - X*B( I, J ) + 50 CONTINUE + END IF + END IF + 60 CONTINUE + IF( B( N, N ).EQ.ZERO ) + $ B( N, N ) = EPS3 +* + TRANS = 'N' +* + ELSE +* +* UL decomposition with partial pivoting of B, replacing zero +* pivots by EPS3. +* + DO 90 J = N, 2, -1 + EJ = H( J, J-1 ) + IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN +* +* Interchange columns and eliminate. +* + X = B( J, J ) / EJ + B( J, J ) = EJ + DO 70 I = 1, J - 1 + TEMP = B( I, J-1 ) + B( I, J-1 ) = B( I, J ) - X*TEMP + B( I, J ) = TEMP + 70 CONTINUE + ELSE +* +* Eliminate without interchange. +* + IF( B( J, J ).EQ.ZERO ) + $ B( J, J ) = EPS3 + X = EJ / B( J, J ) + IF( X.NE.ZERO ) THEN + DO 80 I = 1, J - 1 + B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) + 80 CONTINUE + END IF + END IF + 90 CONTINUE + IF( B( 1, 1 ).EQ.ZERO ) + $ B( 1, 1 ) = EPS3 +* + TRANS = 'T' +* + END IF +* + NORMIN = 'N' + DO 110 ITS = 1, N +* +* Solve U*x = scale*v for a right eigenvector +* or U'*x = scale*v for a left eigenvector, +* overwriting x on v. +* + CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, + $ VR, SCALE, WORK, IERR ) + NORMIN = 'Y' +* +* Test for sufficient growth in the norm of v. +* + VNORM = DASUM( N, VR, 1 ) + IF( VNORM.GE.GROWTO*SCALE ) + $ GO TO 120 +* +* Choose new orthogonal starting vector and try again. +* + TEMP = EPS3 / ( ROOTN+ONE ) + VR( 1 ) = EPS3 + DO 100 I = 2, N + VR( I ) = TEMP + 100 CONTINUE + VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN + 110 CONTINUE +* +* Failure to find eigenvector in N iterations. +* + INFO = 1 +* + 120 CONTINUE +* +* Normalize eigenvector. +* + I = IDAMAX( N, VR, 1 ) + CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 ) + ELSE +* +* Complex eigenvalue. +* + IF( NOINIT ) THEN +* +* Set initial vector. +* + DO 130 I = 1, N + VR( I ) = EPS3 + VI( I ) = ZERO + 130 CONTINUE + ELSE +* +* Scale supplied initial vector. +* + NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) ) + REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML ) + CALL DSCAL( N, REC, VR, 1 ) + CALL DSCAL( N, REC, VI, 1 ) + END IF +* + IF( RIGHTV ) THEN +* +* LU decomposition with partial pivoting of B, replacing zero +* pivots by EPS3. +* +* The imaginary part of the (i,j)-th element of U is stored in +* B(j+1,i). +* + B( 2, 1 ) = -WI + DO 140 I = 2, N + B( I+1, 1 ) = ZERO + 140 CONTINUE +* + DO 170 I = 1, N - 1 + ABSBII = DLAPY2( B( I, I ), B( I+1, I ) ) + EI = H( I+1, I ) + IF( ABSBII.LT.ABS( EI ) ) THEN +* +* Interchange rows and eliminate. +* + XR = B( I, I ) / EI + XI = B( I+1, I ) / EI + B( I, I ) = EI + B( I+1, I ) = ZERO + DO 150 J = I + 1, N + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - XR*TEMP + B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP + B( I, J ) = TEMP + B( J+1, I ) = ZERO + 150 CONTINUE + B( I+2, I ) = -WI + B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI + B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI + ELSE +* +* Eliminate without interchanging rows. +* + IF( ABSBII.EQ.ZERO ) THEN + B( I, I ) = EPS3 + B( I+1, I ) = ZERO + ABSBII = EPS3 + END IF + EI = ( EI / ABSBII ) / ABSBII + XR = B( I, I )*EI + XI = -B( I+1, I )*EI + DO 160 J = I + 1, N + B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) + + $ XI*B( J+1, I ) + B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J ) + 160 CONTINUE + B( I+2, I+1 ) = B( I+2, I+1 ) - WI + END IF +* +* Compute 1-norm of offdiagonal elements of i-th row. +* + WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) + + $ DASUM( N-I, B( I+2, I ), 1 ) + 170 CONTINUE + IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO ) + $ B( N, N ) = EPS3 + WORK( N ) = ZERO +* + I1 = N + I2 = 1 + I3 = -1 + ELSE +* +* UL decomposition with partial pivoting of conjg(B), +* replacing zero pivots by EPS3. +* +* The imaginary part of the (i,j)-th element of U is stored in +* B(j+1,i). +* + B( N+1, N ) = WI + DO 180 J = 1, N - 1 + B( N+1, J ) = ZERO + 180 CONTINUE +* + DO 210 J = N, 2, -1 + EJ = H( J, J-1 ) + ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) ) + IF( ABSBJJ.LT.ABS( EJ ) ) THEN +* +* Interchange columns and eliminate +* + XR = B( J, J ) / EJ + XI = B( J+1, J ) / EJ + B( J, J ) = EJ + B( J+1, J ) = ZERO + DO 190 I = 1, J - 1 + TEMP = B( I, J-1 ) + B( I, J-1 ) = B( I, J ) - XR*TEMP + B( J, I ) = B( J+1, I ) - XI*TEMP + B( I, J ) = TEMP + B( J+1, I ) = ZERO + 190 CONTINUE + B( J+1, J-1 ) = WI + B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI + B( J, J-1 ) = B( J, J-1 ) - XR*WI + ELSE +* +* Eliminate without interchange. +* + IF( ABSBJJ.EQ.ZERO ) THEN + B( J, J ) = EPS3 + B( J+1, J ) = ZERO + ABSBJJ = EPS3 + END IF + EJ = ( EJ / ABSBJJ ) / ABSBJJ + XR = B( J, J )*EJ + XI = -B( J+1, J )*EJ + DO 200 I = 1, J - 1 + B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) + + $ XI*B( J+1, I ) + B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J ) + 200 CONTINUE + B( J, J-1 ) = B( J, J-1 ) + WI + END IF +* +* Compute 1-norm of offdiagonal elements of j-th column. +* + WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) + + $ DASUM( J-1, B( J+1, 1 ), LDB ) + 210 CONTINUE + IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO ) + $ B( 1, 1 ) = EPS3 + WORK( 1 ) = ZERO +* + I1 = 1 + I2 = N + I3 = 1 + END IF +* + DO 270 ITS = 1, N + SCALE = ONE + VMAX = ONE + VCRIT = BIGNUM +* +* Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector, +* or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector, +* overwriting (xr,xi) on (vr,vi). +* + DO 250 I = I1, I2, I3 +* + IF( WORK( I ).GT.VCRIT ) THEN + REC = ONE / VMAX + CALL DSCAL( N, REC, VR, 1 ) + CALL DSCAL( N, REC, VI, 1 ) + SCALE = SCALE*REC + VMAX = ONE + VCRIT = BIGNUM + END IF +* + XR = VR( I ) + XI = VI( I ) + IF( RIGHTV ) THEN + DO 220 J = I + 1, N + XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J ) + XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J ) + 220 CONTINUE + ELSE + DO 230 J = 1, I - 1 + XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J ) + XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J ) + 230 CONTINUE + END IF +* + W = ABS( B( I, I ) ) + ABS( B( I+1, I ) ) + IF( W.GT.SMLNUM ) THEN + IF( W.LT.ONE ) THEN + W1 = ABS( XR ) + ABS( XI ) + IF( W1.GT.W*BIGNUM ) THEN + REC = ONE / W1 + CALL DSCAL( N, REC, VR, 1 ) + CALL DSCAL( N, REC, VI, 1 ) + XR = VR( I ) + XI = VI( I ) + SCALE = SCALE*REC + VMAX = VMAX*REC + END IF + END IF +* +* Divide by diagonal element of B. +* + CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ), + $ VI( I ) ) + VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX ) + VCRIT = BIGNUM / VMAX + ELSE + DO 240 J = 1, N + VR( J ) = ZERO + VI( J ) = ZERO + 240 CONTINUE + VR( I ) = ONE + VI( I ) = ONE + SCALE = ZERO + VMAX = ONE + VCRIT = BIGNUM + END IF + 250 CONTINUE +* +* Test for sufficient growth in the norm of (VR,VI). +* + VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 ) + IF( VNORM.GE.GROWTO*SCALE ) + $ GO TO 280 +* +* Choose a new orthogonal starting vector and try again. +* + Y = EPS3 / ( ROOTN+ONE ) + VR( 1 ) = EPS3 + VI( 1 ) = ZERO +* + DO 260 I = 2, N + VR( I ) = Y + VI( I ) = ZERO + 260 CONTINUE + VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN + 270 CONTINUE +* +* Failure to find eigenvector in N iterations +* + INFO = 1 +* + 280 CONTINUE +* +* Normalize eigenvector. +* + VNORM = ZERO + DO 290 I = 1, N + VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) ) + 290 CONTINUE + CALL DSCAL( N, ONE / VNORM, VR, 1 ) + CALL DSCAL( N, ONE / VNORM, VI, 1 ) +* + END IF +* + RETURN +* +* End of DLAEIN +* + END |