SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO ) * * -- LAPACK auxiliary routine (instrumented to count operations) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. LOGICAL NOINIT, RIGHTV INTEGER INFO, LDB, LDH, N REAL BIGNUM, EPS3, SMLNUM, WI, WR * .. * .. Array Arguments .. REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ), $ WORK( * ) * .. * Common block to return operation count. * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. REAL ITCNT, OPS * .. * * Purpose * ======= * * SLAEIN 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) REAL 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) REAL * WI (input) REAL * The real and imaginary parts of the eigenvalue of H whose * corresponding right or left eigenvector is to be computed. * * VR (input/output) REAL array, dimension (N) * VI (input/output) REAL 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) REAL array, dimension (LDB,N) * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= N+1. * * WORK (workspace) REAL array, dimension (N) * * EPS3 (input) REAL * A small machine-dependent value which is used to perturb * close eigenvalues, and to replace zero pivots. * * SMLNUM (input) REAL * A machine-dependent value close to the underflow threshold. * * BIGNUM (input) REAL * 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 .. REAL ZERO, ONE, TENTH PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TENTH = 1.0E-1 ) * .. * .. Local Scalars .. CHARACTER NORMIN, TRANS INTEGER I, I1, I2, I3, IERR, ITS, J REAL ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML, $ OPST, REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, $ VNORM, W, W1, X, XI, XR, Y * .. * .. External Functions .. INTEGER ISAMAX REAL SASUM, SLAPY2, SNRM2 EXTERNAL ISAMAX, SASUM, SLAPY2, SNRM2 * .. * .. External Subroutines .. EXTERNAL SLADIV, SLATRS, SSCAL * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, REAL, SQRT * .. * .. Executable Statements .. * INFO = 0 *** * Initialize OPST = 0 *** * * GROWTO is the threshold used in the acceptance test for an * eigenvector. * ROOTN = SQRT( REAL( N ) ) GROWTO = TENTH / ROOTN NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM *** * Increment op count for computing ROOTN, GROWTO and NRMSML OPST = OPST + 4 *** * * 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 *** OPST = OPST + N *** * 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 = SNRM2( N, VR, 1 ) CALL SSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR, $ 1 ) *** OPST = OPST + ( 3*N+2 ) *** 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 *** * Increment op count for LU decomposition OPS = OPS + ( N-1 )*( N+1 ) *** * 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 *** * Increment op count for UL decomposition OPS = OPS + ( N-1 )*( N+1 ) *** * 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 SLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, $ VR, SCALE, WORK, IERR ) *** * Increment opcount for triangular solver, assuming that * ops SLATRS = ops STRSV, with no scaling in SLATRS. OPS = OPS + N*N *** NORMIN = 'Y' * * Test for sufficient growth in the norm of v. * VNORM = SASUM( N, VR, 1 ) *** OPST = OPST + N *** 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 *** OPST = OPST + 4 *** 110 CONTINUE * * Failure to find eigenvector in N iterations. * INFO = 1 * 120 CONTINUE * * Normalize eigenvector. * I = ISAMAX( N, VR, 1 ) CALL SSCAL( N, ONE / ABS( VR( I ) ), VR, 1 ) *** OPST = OPST + ( 2*N+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 = SLAPY2( SNRM2( N, VR, 1 ), SNRM2( N, VI, 1 ) ) REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML ) CALL SSCAL( N, REC, VR, 1 ) CALL SSCAL( N, REC, VI, 1 ) *** OPST = OPST + ( 6*N+5 ) *** 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 = SLAPY2( 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 *** OPST = OPST + ( 4*( N-I )+6 ) *** 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 *** OPST = OPST + ( 7*( N-I )+4 ) *** END IF * * Compute 1-norm of offdiagonal elements of i-th row. * WORK( I ) = SASUM( N-I, B( I, I+1 ), LDB ) + $ SASUM( N-I, B( I+2, I ), 1 ) *** OPST = OPST + ( 2*( N-I )+4 ) *** 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 = SLAPY2( 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 *** OPST = OPST + ( 4*( J-1 )+6 ) *** 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 *** OPST = OPST + ( 7*( J-1 )+4 ) *** END IF * * Compute 1-norm of offdiagonal elements of j-th column. * WORK( J ) = SASUM( J-1, B( 1, J ), 1 ) + $ SASUM( J-1, B( J+1, 1 ), LDB ) *** OPST = OPST + ( 2*( J-1 )+4 ) *** 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 SSCAL( N, REC, VR, 1 ) CALL SSCAL( 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 SSCAL( N, REC, VR, 1 ) CALL SSCAL( 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 SLADIV( 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 *** OPST = OPST + 9 *** 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 *** * Increment op count for loop 260, assuming no scaling OPS = OPS + 4*N*( N-1 ) *** * * Test for sufficient growth in the norm of (VR,VI). * VNORM = SASUM( N, VR, 1 ) + SASUM( N, VI, 1 ) *** OPST = OPST + 2*N *** 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 *** OPST = OPST + 4 *** 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 SSCAL( N, ONE / VNORM, VR, 1 ) CALL SSCAL( N, ONE / VNORM, VI, 1 ) *** OPST = OPST + ( 4*N+1 ) *** * END IF * *** * Compute final op count OPS = OPS + OPST *** RETURN * * End of SLAEIN * END