Mercurial > octave-libgccjit
changeset 3848:562c1b1fa5f4
[project @ 2001-08-13 17:26:42 by jwe]
| author | jwe |
|---|---|
| date | Mon, 13 Aug 2001 17:26:42 +0000 |
| parents | 92fb162eba24 |
| children | 5266e351a19c |
| files | libcruft/ChangeLog libcruft/lapack/dlasq3.f libcruft/lapack/dlasq5.f |
| diffstat | 3 files changed, 65 insertions(+), 63 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Fri Aug 10 17:20:24 2001 +0000 +++ b/libcruft/ChangeLog Mon Aug 13 17:26:42 2001 +0000 @@ -1,3 +1,8 @@ +2001-08-13 John W. Eaton <jwe@bevo.che.wisc.edu> + + * lapack/dlasq3.f: Update from netlib. + * lapack/dlasq5.f: Ditto. + 2001-05-02 Mumit Khan <khan@nanotech.wisc.edu> * Makefile.in (CRUFT_DIRS): Substitute @FFT_DIR@.
--- a/libcruft/lapack/dlasq3.f Fri Aug 10 17:20:24 2001 +0000 +++ b/libcruft/lapack/dlasq3.f Mon Aug 13 17:26:42 2001 +0000 @@ -4,7 +4,7 @@ * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* October 31, 1999 +* May 17, 2000 * * .. Scalar Arguments .. LOGICAL IEEE @@ -32,7 +32,7 @@ * Last index. * * Z (input) DOUBLE PRECISION array, dimension ( 4*N ) -* Z holds the qd array. +* Z holds the qd array. * * PP (input) INTEGER * PP=0 for ping, PP=1 for pong. @@ -47,7 +47,7 @@ * Lower order part of SIGMA * * QMAX (input) DOUBLE PRECISION -* Maximum value of q. +* Maximum value of q. * * NFAIL (output) INTEGER * Number of times shift was too big. @@ -75,7 +75,7 @@ * .. * .. Local Scalars .. INTEGER IPN4, J4, N0IN, NN, TTYPE - DOUBLE PRECISION DMIN1, DMIN2, DN, DN1, DN2, EPS, S, SAFMIN, T, + DOUBLE PRECISION DMIN1, DMIN2, DN, DN1, DN2, EPS, S, SAFMIN, T, $ TAU, TEMP, TOL, TOL2 * .. * .. External Subroutines .. @@ -105,16 +105,16 @@ TOL = EPS*HUNDRD TOL2 = TOL**2 * -* Check for deflation. +* Check for deflation. * 10 CONTINUE * IF( N0.LT.I0 ) $ RETURN - IF( N0.EQ.I0 ) + IF( N0.EQ.I0 ) $ GO TO 20 NN = 4*N0 + PP - IF( N0.EQ.( I0+1 ) ) + IF( N0.EQ.( I0+1 ) ) $ GO TO 40 * * Check whether E(N0-1) is negligible, 1 eigenvalue. @@ -217,57 +217,39 @@ NDIV = NDIV + ( N0-I0+2 ) ITER = ITER + 1 * - IF( DMIN.NE.DMIN ) THEN -* -* Check for NaN: "DMIN.NE.DMIN" +* Check status. * - Z( 4*N0+PP-1 ) = ZERO - GO TO 70 - ELSE IF( Z( 4*N0+PP ).LE.ZERO ) THEN + IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN * -* Possible unnecessary underflow in the e's. -* Call safe dqd. +* Success. * - Z( 4*N0+PP-1 ) = ZERO - DMIN = ZERO - GO TO 70 - ELSE IF( DMIN.EQ.ZERO .AND. DMIN1.LE.ZERO ) THEN + GO TO 100 * -* Possible unnecessary underflow in the d's. -* Call safe dqd. + ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. + $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. + $ ABS( DN ).LT.TOL*SIGMA ) THEN * - Z( 4*N0+PP-1 ) = ZERO - GO TO 70 - END IF -* -* Check for convergence hidden by negative DN. +* Convergence hidden by negative DN. * - IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. Z( 4*( N0-1 )-PP ) - $ .LT.TOL*( SIGMA+DN1 ) .AND. ABS( DN ).LT.TOL*SIGMA ) THEN Z( 4*( N0-1 )-PP+2 ) = ZERO - DMIN = ABS( DMIN ) - END IF + DMIN = ZERO + GO TO 100 + ELSE IF( DMIN.LT.ZERO ) THEN * - IF( DMIN.LT.ZERO ) THEN -* -* Failure. Select new TAU and try again. +* TAU too big. Select new TAU and try again. * NFAIL = NFAIL + 1 -* -* Failed twice. Play it safe. + IF( TTYPE.LT.-22 ) THEN * - IF( TTYPE.LT.-22 ) THEN - Z( 4*N0+PP-1 ) = ZERO - DMIN = ZERO - GO TO 70 - END IF +* Failed twice. Play it safe. * - IF( DMIN1.GT.ZERO ) THEN + TAU = ZERO + ELSE IF( DMIN1.GT.ZERO ) THEN * * Late failure. Gives excellent shift. * - TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) - TTYPE = TTYPE - 11 + TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) + TTYPE = TTYPE - 11 ELSE * * Early failure. Divide by 4. @@ -276,17 +258,32 @@ TTYPE = TTYPE - 12 END IF GO TO 80 + ELSE IF( DMIN.NE.DMIN ) THEN +* +* NaN. +* + TAU = ZERO + GO TO 80 + ELSE +* +* Possible underflow. Play it safe. +* + GO TO 90 END IF - ELSE - CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) - NDIV = NDIV + ( N0-I0 ) - ITER = ITER + 1 - TAU = ZERO END IF * +* Risk of underflow. +* + 90 CONTINUE + CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) + NDIV = NDIV + ( N0-I0+2 ) + ITER = ITER + 1 + TAU = ZERO +* + 100 CONTINUE IF( TAU.LT.SIGMA ) THEN DESIG = DESIG + TAU - T = SIGMA + DESIG + T = SIGMA + DESIG DESIG = DESIG - ( T-SIGMA ) ELSE T = SIGMA + TAU
--- a/libcruft/lapack/dlasq5.f Fri Aug 10 17:20:24 2001 +0000 +++ b/libcruft/lapack/dlasq5.f Mon Aug 13 17:26:42 2001 +0000 @@ -4,7 +4,7 @@ * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* October 31, 1999 +* May 17, 2000 * * .. Scalar Arguments .. LOGICAL IEEE @@ -80,7 +80,7 @@ $ RETURN * J4 = 4*I0 + PP - 3 - EMIN = Z( J4+4 ) + EMIN = Z( J4+4 ) D = Z( J4 ) - TAU DMIN = D DMIN1 = -Z( J4 ) @@ -91,7 +91,7 @@ * IF( PP.EQ.0 ) THEN DO 10 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-2 ) = D + Z( J4-1 ) + Z( J4-2 ) = D + Z( J4-1 ) TEMP = Z( J4+1 ) / Z( J4-2 ) D = D*TEMP - TAU DMIN = MIN( DMIN, D ) @@ -100,7 +100,7 @@ 10 CONTINUE ELSE DO 20 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-3 ) = D + Z( J4 ) + Z( J4-3 ) = D + Z( J4 ) TEMP = Z( J4+2 ) / Z( J4-3 ) D = D*TEMP - TAU DMIN = MIN( DMIN, D ) @@ -109,7 +109,7 @@ 20 CONTINUE END IF * -* Unroll last two steps. +* Unroll last two steps. * DNM2 = D DMIN2 = DMIN @@ -134,10 +134,10 @@ * IF( PP.EQ.0 ) THEN DO 30 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-2 ) = D + Z( J4-1 ) - IF( D.LE.ZERO ) THEN + Z( J4-2 ) = D + Z( J4-1 ) + IF( D.LT.ZERO ) THEN RETURN - ELSE + ELSE Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU END IF @@ -146,10 +146,10 @@ 30 CONTINUE ELSE DO 40 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-3 ) = D + Z( J4 ) - IF( D.LE.ZERO ) THEN + Z( J4-3 ) = D + Z( J4 ) + IF( D.LT.ZERO ) THEN RETURN - ELSE + ELSE Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU END IF @@ -158,14 +158,14 @@ 40 CONTINUE END IF * -* Unroll last two steps. +* Unroll last two steps. * DNM2 = D DMIN2 = DMIN J4 = 4*( N0-2 ) - PP J4P2 = J4 + 2*PP - 1 Z( J4-2 ) = DNM2 + Z( J4P2 ) - IF( DNM2.LE.ZERO ) THEN + IF( DNM2.LT.ZERO ) THEN RETURN ELSE Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) @@ -177,7 +177,7 @@ J4 = J4 + 4 J4P2 = J4 + 2*PP - 1 Z( J4-2 ) = DNM1 + Z( J4P2 ) - IF( DNM1.LE.ZERO ) THEN + IF( DNM1.LT.ZERO ) THEN RETURN ELSE Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )