diff mbox

[Fortran,committed] PR 45367/36158 - Increase numeric tolerance in bessel_{6,7}.f90

Message ID 4C70F41C.8050801@net-b.de
State New
Headers show

Commit Message

Tobias Burnus Aug. 22, 2010, 9:55 a.m. UTC 
As HJ and Dominique reported, for Linux/ia32 and 
powerpc-apple-darwin9, the tolerance values were too tight for the 
comparison of the recurrence algorithm with the library version.

I increased now the tolerance based on Dominique's values and cross 
checked using -m32 on x86-64-linux. I hope it will now work on most 
systems, but we still might need to tweak the values (or use xfail) for 
some architectures as the initial values are slightly different and as 
the architectures might do different rounding or take different code paths.

Thanks to Dominique for testing and adjusting the values - and to HJ for 
keeping an eye on regressions. Committed as Rev. 163454/163455.

  * * *

Note of warning to users: Be careful when applying the recurrence 
algorithm (i.e. when using BESSEL_JN/YN(N1, N2, X): The results can be 
quite different! For instance, using
   BESSEL_JN(0, 100, 475.78)
the deviation is < 6*epsilon(0.0) - except for the following N, where it 
is (much) larger:
N=63: < 9*epsilon(0.0)
N=91: < 13*epsilon(0.0)
N=49: < 326*epsilon(0.0)

Note in particular that checking the last item in the recurrence 
algorithm - for JN that's lowest, here: N=0, does not help: For the 
example above the difference to the value returned by the library for JN 
is just 0.5*epsilon while for N=49 the difference is much larger! (The 
reason is that for N=49 the result is -0.8E-04 while the other values 
are of the order 0.xE-01 and 0.xE-02.)
.
Tobias
diff mbox

Patch

Index: gcc/testsuite/ChangeLog
===================================================================
--- gcc/testsuite/ChangeLog	(Revision 163454)
+++ gcc/testsuite/ChangeLog	(Arbeitskopie)
@@ -1,4 +1,9 @@ 
 2010-08-22  Tobias Burnus  <burnus@net-b.de>
+
+	PR fortran/36158
+	* gfortran.dg/bessel_7.f90: Disable accidently enabled debug output.
+
+2010-08-22  Tobias Burnus  <burnus@net-b.de>
 	    Dominique d'Humieres <dominiq@lps.ens.fr>
 
 	PR fortran/45367
Index: gcc/testsuite/gfortran.dg/bessel_7.f90
===================================================================
--- gcc/testsuite/gfortran.dg/bessel_7.f90	(Revision 163454)
+++ gcc/testsuite/gfortran.dg/bessel_7.f90	(Arbeitskopie)
@@ -32,13 +32,13 @@ 
 rec = BESSEL_YN(0, Nmax, X)
 lib = [ (BESSEL_YN(i, X), i=0,Nmax) ]
 
-print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x)
+!print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x)
 do i = 0, Nmax
-  print '(i2,2e17.9,e12.2,f14.10,2l3)', i, rec(i), lib(i), &
-        rec(i)-lib(i), ((rec(i)-lib(i))/rec(i))/epsilon(x), &
-        i > nit .or. rec(i) == lib(i) &
-                .or. abs((rec(i)-lib(i))/rec(i)) < myeps2, &
-        rec(i) == lib(i) .or. abs((rec(i)-lib(i))/rec(i)) < myeps
+!  print '(i2,2e17.9,e12.2,f14.10,2l3)', i, rec(i), lib(i), &
+!        rec(i)-lib(i), ((rec(i)-lib(i))/rec(i))/epsilon(x), &
+!        i > nit .or. rec(i) == lib(i) &
+!                .or. abs((rec(i)-lib(i))/rec(i)) < myeps2, &
+!        rec(i) == lib(i) .or. abs((rec(i)-lib(i))/rec(i)) < myeps
 if (.not. (i > nit .or. rec(i) == lib(i) &
                    .or. abs((rec(i)-lib(i))/rec(i)) < myeps2)) &
   call abort ()