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

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

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.