From patchwork Sun Aug 22 09:55:40 2010
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Subject: [Fortran,
 committed] PR 45367/36158 - Increase numeric tolerance in bessel_{6, 7}.f90
Date: Sat, 21 Aug 2010 23:55:40 -0000
From: Tobias Burnus <burnus@net-b.de>
X-Patchwork-Id: 62368
Message-Id: <4C70F41C.8050801@net-b.de>
To: gcc patches <gcc-patches@gcc.gnu.org>, gfortran <fortran@gcc.gnu.org>

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