Submitter | Tobias Burnus |
---|---|

Date | Aug. 22, 2010, 2:06 p.m. |

Message ID | <4C712EF2.4020404@net-b.de> |

Download | mbox | patch |

Permalink | /patch/62388/ |

State | New |

Headers | show |

## Comments

On Sun, Aug 22, 2010 at 04:06:42PM +0200, Tobias Burnus wrote: > Tobias Burnus wrote: >> 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 > > Dominique found additional failures for bessel_6.f90. Let's see how many > systems still fail with the numbers below. I have still the hope that > there are reasonable numbers which work on all systems with working > jn/yn functions ... > > Committed as Rev. 163460 > > Tobias Tobias, No regressions in the gfortran testsuite at r163460 using -m32 or -m64 on x86_64-apple-darwin10. Jack

On Sun, Aug 22, 2010 at 7:06 AM, Tobias Burnus <burnus@net-b.de> wrote: > Tobias Burnus wrote: >> >> 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 > > Dominique found additional failures for bessel_6.f90. Let's see how many > systems still fail with the numbers below. I have still the hope that there > are reasonable numbers which work on all systems with working jn/yn > functions ... > > Committed as Rev. 163460 > gfortran.dg/bessel_7.f90 still fails with -m32 -mfpmath=sse: http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02282.html http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02284.html

H.J. Lu wrote: > gfortran.dg/bessel_7.f90 still fails with -m32 -mfpmath=sse: > > http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02282.html > http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02284.html I think that's because we have different libm versions. And because the recurrence algorithm is somehow not very stable for certain combinations of N and X. $ gfortran -m32 -mfpmath=sse bessel_7.f90 && ./a.out $ Can you uncomment the "print" lines in bessel_7 and comment the if and abort line - and send me the output? Tobias

On Sun, Aug 22, 2010 at 2:17 PM, Tobias Burnus <burnus@net-b.de> wrote: > H.J. Lu wrote: >> >> gfortran.dg/bessel_7.f90 still fails with -m32 -mfpmath=sse: >> >> http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02282.html >> http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02284.html > > I think that's because we have different libm versions. And because the > recurrence algorithm is somehow not very stable for certain combinations of > N and X. > > $ gfortran -m32 -mfpmath=sse bessel_7.f90 && ./a.out > $ > > Can you uncomment the "print" lines in bessel_7 and comment the if and abort > line - and send me the output? > Here it is.

On Sun, Aug 22, 2010 at 11:17:19PM +0200, Tobias Burnus wrote: > H.J. Lu wrote: > >gfortran.dg/bessel_7.f90 still fails with -m32 -mfpmath=sse: > > > >http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02282.html > >http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02284.html > > I think that's because we have different libm versions. And because the > recurrence algorithm is somehow not very stable for certain combinations > of N and X. > The recurrence algorithm is stable. What you are see is a cancellation issue. If j(n-1) is a small number then the most significant bits in (2n/x)*j(n) and j(n+1) are the same, so in j(n-1) = 2n/x * j(n) - j(n+1) one loses precision. When one refers to the stability of the recurrence algorithm, one is normally interested in the propagation of error. For bessel functions, backwards recurrence is stable for all x and n with the caveat that one needs j(n) and j(n+1) to start the recurrence. If n < x, then forward recurrence is stable but once n exceeds x the error starts to grow. ! ! Forward recurrence for bessel fcns. ! program testing real j2, j0, j1 j0 = besjn(0,2.) j1 = besjn(1,2.) do i = 1, 10 j2 = (2 * i / 2.) * j1 - j0 print *, j2, besjn(i+1,2.) j0 = j1 j1 = j2 end do end program laptop:kargl[214] gfc4x -o z test.f90 && ./z 0.35283405 0.35283405 0.12894326 0.12894325 3.39957476E-02 3.39957215E-02 7.03972578E-03 7.03962985E-03 1.20288134E-03 1.20242906E-03 1.77562237E-04 1.74944071E-04 4.00543213E-05 2.21795508E-05 1.42872334E-04 2.49234358E-06 1.24579668E-03 2.51538637E-07 1.23150945E-02 2.30428494E-08 There is a paper in an obscure applied math journal that shows the error analysis and shows that the stability of the recurrence. I don't have that paper handy, but I could probably dig it up if interested.

## Patch

Index: gcc/testsuite/ChangeLog =================================================================== --- gcc/testsuite/ChangeLog (Revision 163458) +++ gcc/testsuite/ChangeLog (Arbeitskopie) @@ -1,6 +1,12 @@ 2010-08-22 Tobias Burnus <burnus@net-b.de> + Dominique d'Humieres <dominiq@lps.ens.fr> PR fortran/45367 + * gfortran.dg/bessel_6.f90: Further reduce required accuracy. + +2010-08-22 Tobias Burnus <burnus@net-b.de> + + PR fortran/45367 * gfortran.dg/bessel_6.f90: Fix numeric tolerence. 2010-08-22 Tobias Burnus <burnus@net-b.de> Index: gcc/testsuite/gfortran.dg/bessel_6.f90 =================================================================== --- gcc/testsuite/gfortran.dg/bessel_6.f90 (Revision 163458) +++ gcc/testsuite/gfortran.dg/bessel_6.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, 7, 5, 6, 9, 12, 12, 7, 7, 8, 75, 6 ] + * [2, 7, 5, 6, 9, 12, 12, 7, 7, 8, 75, 15 ] ! The following is sufficient for me - the values above are a bit ! more tolerant ! * [0, 5, 3, 4, 6, 7, 7, 5, 5, 6, 66, 4 ] integer,parameter :: mymax(size(values)) = & - [100, 17, 23, 21, 27, 28, 32, 35, 36, 41, 49, 50 ] + [100, 17, 23, 21, 27, 28, 32, 35, 36, 41, 47, 37 ] integer, parameter :: Nmax = 100 real :: rec(0:Nmax), lib(0:Nmax) integer :: i