gcc/testsuite/gfortran.dg/bessel_6.f90

   1 ! { dg-do run { xfail spu-*-* } }
   2 ! { dg-add-options ieee }
   3 !
   4 ! PR fortran/36158
   5 ! PR fortran/33197
   6 !
   7 ! XFAILed for SPU targets since we don't have an accurate library
   8 ! implementation of the single-precision Bessel functions.
   9 !
  10 ! Run-time tests for transformations BESSEL_JN
  11 !
  12 implicit none
  13 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]
  14 real,parameter :: myeps(size(values)) = epsilon(0.0) &
  15                   * [2, 7, 5, 6, 9, 12, 12, 7, 7, 8, 75, 15 ]
  16 ! The following is sufficient for me - the values above are a bit
  17 ! more tolerant
  18 !                  * [0, 5, 3, 4, 6, 7, 7, 5, 5, 6, 66, 4 ]
  19 integer,parameter :: mymax(size(values)) =  &
  20                  [100, 17, 23, 21, 27, 28, 32, 35, 31, 41, 47, 37 ]
  21 integer, parameter :: Nmax = 100
  22 real :: rec(0:Nmax), lib(0:Nmax)
  23 integer :: i
  24
  25 do i = 1, ubound(values,dim=1)
  26   call compare(mymax(i), values(i), myeps(i))
  27 end do
  28
  29 contains
  30
  31 subroutine compare(mymax, X, myeps)
  32
  33 integer :: i, nit, mymax
  34 real X, myeps, myeps2
  35
  36 rec(0:mymax) = BESSEL_JN(0, mymax, X)
  37 lib(0:mymax) = [ (BESSEL_JN(i, X), i=0,mymax) ]
  38
  39 !print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x)
  40 do i = 0, mymax
  41 !  print '(i2,2e17.9,e12.2,f18.10,2l3)', i, rec(i), lib(i), &
  42 !        rec(i)-lib(i),           ((rec(i)-lib(i))/rec(i))/epsilon(x), &
  43 !        rec(i) == lib(i), abs((rec(i)-lib(i))/rec(i)) < myeps
  44 if (rec(i) == lib(i)) CYCLE
  45 if (abs((rec(i)-lib(i))/rec(i)) > myeps) &
  46   call abort()
  47 end do
  48
  49 end
  50 end