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

[Tobias Burnus <burnus at net-b dot de>]

 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
Index: gcc/testsuite/ChangeLog
===================================================================
--- gcc/testsuite/ChangeLog	(Revision 163453)
+++ gcc/testsuite/ChangeLog	(Arbeitskopie)
@@ -1,3 +1,11 @@
+2010-08-22  Tobias Burnus  <burnus@net-b.de>
+	    Dominique d'Humieres <dominiq@lps.ens.fr>
+
+	PR fortran/45367
+	PR fortran/36158
+	* gfortran.dg/bessel_6.f90: Increase numeric tolerence.
+	* gfortran.dg/bessel_7.f90: Increase numeric tolerence.
+
 2010-08-21  Janus Weil  <janus@gcc.gnu.org>
 
 	PR fortran/44863
Index: gcc/testsuite/gfortran.dg/bessel_7.f90
===================================================================
--- gcc/testsuite/gfortran.dg/bessel_7.f90	(Revision 163453)
+++ gcc/testsuite/gfortran.dg/bessel_7.f90	(Arbeitskopie)
@@ -8,12 +8,12 @@
 implicit none
 real,parameter :: values(*) = [0.0, 0.5, 1.0, 0.9, 1.8,2.0,3.0,4.0,4.25,8.0,34.53, 475.78] 
 real,parameter :: myeps(size(values)) = epsilon(0.0) &
-                  * [2, 2, 2, 5, 5, 2, 12, 2, 4, 3, 30, 130 ]
+                  * [2, 3, 3, 5, 7, 2, 12, 4, 7, 3, 30, 168 ]
 ! The following is sufficient for me - the values above are a bit
 ! more tolerant
 !                  * [0, 0, 0, 3, 3, 0, 9, 0, 2, 1, 22, 130 ]
 integer,parameter :: nit(size(values)) =  &
-                 [100, 100, 100, 100, 100, 100, 10, 100, 100, 100, 10, 25 ]
+                 [100, 100, 100, 25, 15, 100, 10, 31, 7, 100, 7, 25 ]
 integer, parameter :: Nmax = 100
 real :: rec(0:Nmax), lib(0:Nmax)
 integer :: i
@@ -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 ()

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 ()