SUBROUTINE SLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, $ DN1, DN2, TAU, TTYPE ) * * -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * May 17, 2000 * * .. Scalar Arguments .. INTEGER I0, N0, N0IN, PP, TTYPE REAL DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU * .. * .. Array Arguments .. REAL Z( * ) * .. * .. Common block to return operation count .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. REAL ITCNT, OPS * .. * * Purpose * ======= * * SLASQ4 computes an approximation TAU to the smallest eigenvalue * using values of d from the previous transform. * * I0 (input) INTEGER * First index. * * N0 (input) INTEGER * Last index. * * Z (input) REAL array, dimension ( 4*N ) * Z holds the qd array. * * PP (input) INTEGER * PP=0 for ping, PP=1 for pong. * * NOIN (input) INTEGER * The value of N0 at start of EIGTEST. * * DMIN (input) REAL * Minimum value of d. * * DMIN1 (input) REAL * Minimum value of d, excluding D( N0 ). * * DMIN2 (input) REAL * Minimum value of d, excluding D( N0 ) and D( N0-1 ). * * DN (input) REAL * d(N) * * DN1 (input) REAL * d(N-1) * * DN2 (input) REAL * d(N-2) * * TAU (output) REAL * This is the shift. * * TTYPE (output) INTEGER * Shift type. * * Further Details * =============== * CNST1 = 9/16 * * ===================================================================== * * .. Parameters .. REAL CNST1, CNST2, CNST3 PARAMETER ( CNST1 = 0.5630E0, CNST2 = 1.010E0, $ CNST3 = 1.050E0 ) REAL QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD PARAMETER ( QURTR = 0.250E0, THIRD = 0.3330E0, $ HALF = 0.50E0, ZERO = 0.0E0, ONE = 1.0E0, $ TWO = 2.0E0, HUNDRD = 100.0E0 ) * .. * .. Local Scalars .. INTEGER I4, NN, NP REAL A2, B1, B2, G, GAM, GAP1, GAP2, S * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN, REAL, SQRT * .. * .. Save statement .. SAVE G * .. * .. Data statement .. DATA G / ZERO / * .. * .. Executable Statements .. * * A negative DMIN forces the shift to take that absolute value * TTYPE records the type of shift. * IF( DMIN.LE.ZERO ) THEN TAU = -DMIN TTYPE = -1 RETURN END IF * NN = 4*N0 + PP IF( N0IN.EQ.N0 ) THEN * * No eigenvalues deflated. * IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN * OPS = OPS + REAL( 7 ) B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) ) B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) ) A2 = Z( NN-7 ) + Z( NN-5 ) * * Cases 2 and 3. * IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN OPS = OPS + REAL( 3 ) GAP2 = DMIN2 - A2 - DMIN2*QURTR IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN OPS = OPS + REAL( 4 ) GAP1 = A2 - DN - ( B2 / GAP2 )*B2 ELSE OPS = OPS + REAL( 3 ) GAP1 = A2 - DN - ( B1+B2 ) END IF IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN OPS = OPS + REAL( 4 ) S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN ) TTYPE = -2 ELSE OPS = OPS + REAL( 2 ) S = ZERO IF( DN.GT.B1 ) $ S = DN - B1 IF( A2.GT.( B1+B2 ) ) $ S = MIN( S, A2-( B1+B2 ) ) S = MAX( S, THIRD*DMIN ) TTYPE = -3 END IF ELSE * * Case 4. * TTYPE = -4 OPS = OPS + REAL( 1 ) S = QURTR*DMIN IF( DMIN.EQ.DN ) THEN OPS = OPS + REAL( 1 ) GAM = DN A2 = ZERO IF( Z( NN-5 ) .GT. Z( NN-7 ) ) $ RETURN B2 = Z( NN-5 ) / Z( NN-7 ) NP = NN - 9 ELSE OPS = OPS + REAL( 2 ) NP = NN - 2*PP B2 = Z( NP-2 ) GAM = DN1 IF( Z( NP-4 ) .GT. Z( NP-2 ) ) $ RETURN A2 = Z( NP-4 ) / Z( NP-2 ) IF( Z( NN-9 ) .GT. Z( NN-11 ) ) $ RETURN B2 = Z( NN-9 ) / Z( NN-11 ) NP = NN - 13 END IF * * Approximate contribution to norm squared from I < NN-1. * A2 = A2 + B2 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 OPS = OPS + REAL( 5 ) IF( B2.EQ.ZERO ) $ GO TO 20 B1 = B2 IF( Z( I4 ) .GT. Z( I4-2 ) ) $ RETURN B2 = B2*( Z( I4 ) / Z( I4-2 ) ) A2 = A2 + B2 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) $ GO TO 20 10 CONTINUE 20 CONTINUE OPS = OPS + REAL( 1 ) A2 = CNST3*A2 * * Rayleigh quotient residual bound. * OPS = OPS + REAL( 5 ) IF( A2.LT.CNST1 ) $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) END IF ELSE IF( DMIN.EQ.DN2 ) THEN * * Case 5. * TTYPE = -5 OPS = OPS + REAL( 1 ) S = QURTR*DMIN * * Compute contribution to norm squared from I > NN-2. * OPS = OPS + REAL( 4 ) NP = NN - 2*PP B1 = Z( NP-2 ) B2 = Z( NP-6 ) GAM = DN2 IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) $ RETURN A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) * * Approximate contribution to norm squared from I < NN-2. * IF( N0-I0.GT.2 ) THEN OPS = OPS + REAL( 3 ) B2 = Z( NN-13 ) / Z( NN-15 ) A2 = A2 + B2 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 OPS = OPS + REAL( 5 ) IF( B2.EQ.ZERO ) $ GO TO 40 B1 = B2 IF( Z( I4 ) .GT. Z( I4-2 ) ) $ RETURN B2 = B2*( Z( I4 ) / Z( I4-2 ) ) A2 = A2 + B2 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) $ GO TO 40 30 CONTINUE 40 CONTINUE A2 = CNST3*A2 END IF * OPS = OPS + REAL( 5 ) IF( A2.LT.CNST1 ) $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) ELSE * * Case 6, no information to guide us. * IF( TTYPE.EQ.-6 ) THEN OPS = OPS + REAL( 3 ) G = G + THIRD*( ONE-G ) ELSE IF( TTYPE.EQ.-18 ) THEN OPS = OPS + REAL( 1 ) G = QURTR*THIRD ELSE G = QURTR END IF OPS = OPS + REAL( 1 ) S = G*DMIN TTYPE = -6 END IF * ELSE IF( N0IN.EQ.( N0+1 ) ) THEN * * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. * IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN * * Cases 7 and 8. * TTYPE = -7 OPS = OPS + REAL( 2 ) S = THIRD*DMIN1 IF( Z( NN-5 ).GT.Z( NN-7 ) ) $ RETURN B1 = Z( NN-5 ) / Z( NN-7 ) B2 = B1 IF( B2.EQ.ZERO ) $ GO TO 60 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 OPS = OPS + REAL( 4 ) A2 = B1 IF( Z( I4 ).GT.Z( I4-2 ) ) $ RETURN B1 = B1*( Z( I4 ) / Z( I4-2 ) ) B2 = B2 + B1 IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) $ GO TO 60 50 CONTINUE 60 CONTINUE OPS = OPS + REAL( 8 ) B2 = SQRT( CNST3*B2 ) A2 = DMIN1 / ( ONE+B2**2 ) GAP2 = HALF*DMIN2 - A2 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN OPS = OPS + REAL( 7 ) S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) ELSE OPS = OPS + REAL( 4 ) S = MAX( S, A2*( ONE-CNST2*B2 ) ) TTYPE = -8 END IF ELSE * * Case 9. * OPS = OPS + REAL( 2 ) S = QURTR*DMIN1 IF( DMIN1.EQ.DN1 ) $ S = HALF*DMIN1 TTYPE = -9 END IF * ELSE IF( N0IN.EQ.( N0+2 ) ) THEN * * Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. * * Cases 10 and 11. * OPS = OPS + REAL( 1 ) IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN TTYPE = -10 OPS = OPS + REAL( 1 ) S = THIRD*DMIN2 IF( Z( NN-5 ).GT.Z( NN-7 ) ) $ RETURN B1 = Z( NN-5 ) / Z( NN-7 ) B2 = B1 IF( B2.EQ.ZERO ) $ GO TO 80 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 OPS = OPS + REAL( 4 ) IF( Z( I4 ).GT.Z( I4-2 ) ) $ RETURN B1 = B1*( Z( I4 ) / Z( I4-2 ) ) B2 = B2 + B1 IF( HUNDRD*B1.LT.B2 ) $ GO TO 80 70 CONTINUE 80 CONTINUE OPS = OPS + REAL( 12 ) B2 = SQRT( CNST3*B2 ) A2 = DMIN2 / ( ONE+B2**2 ) GAP2 = Z( NN-7 ) + Z( NN-9 ) - $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN OPS = OPS + REAL( 7 ) S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) ELSE OPS = OPS + REAL( 4 ) S = MAX( S, A2*( ONE-CNST2*B2 ) ) END IF ELSE OPS = OPS + REAL( 1 ) S = QURTR*DMIN2 TTYPE = -11 END IF ELSE IF( N0IN.GT.( N0+2 ) ) THEN * * Case 12, more than two eigenvalues deflated. No information. * S = ZERO TTYPE = -12 END IF * TAU = S RETURN * * End of SLASQ4 * END