Something horrbily wrong with ERF(x)

Accelerated Computing HPC Compilers Legacy PGI Compilers
1

I’ve been tearing my hair out over the last few days and have finally homed in on something that seems very wrong indeed. Here is some test code - all it should do is to return the real error function for the value 1.0:

function erf_stegun( x ) result( erf_x )
!------------------------------------------------------------------------------
!
! *** Description:
!
!     Real error function ERF(x).
!
! *** Details:
!
! Handbook of Mathematical Functions. Edited by Milton Abramowitz and 
! Irene A. Stegun, Dover Publications, Inc., New York, 1965.
!
! Rational approximation for 0 <= x <= Inf.
!
! Error function and Fresnel Integrals, EQN. 7.1.28.
! Valid to |E(x)| <= 3e-7. Calculation in double precision, result returned 
! in gpu precision.
!
!------------------------------------------------------------------------------

implicit none

real*8              :: x
real*8              :: erf_x

real*8              :: a1 = 0.0705230784d0
real*8              :: a2 = 0.0422820123d0
real*8              :: a3 = 0.0092705272d0
real*8              :: a4 = 0.0001520143d0
real*8              :: a5 = 0.0002765672d0
real*8              :: a6 = 0.0000430638d0

erf_x  = 1.0d0 - 1.0d0/(1.0d0 + a1*x    + a2*x**2 +                           &
                                a3*x**3 + a4*x**4 +                           &
                                a5*x**5 + a6*x**6)**16

print*, 'erf_stegun (function)   : ', x, erf_x

end function erf_stegun

subroutine erf_stegun_s( x, erf_x )
!------------------------------------------------------------------------------
!
! *** Description:
!
!     Real error function ERF(x).
!
! *** Details:
!
! Handbook of Mathematical Functions. Edited by Milton Abramowitz and 
! Irene A. Stegun, Dover Publications, Inc., New York, 1965.
!
! Rational approximation for 0 <= x <= Inf.
!
! Error function and Fresnel Integrals, EQN. 7.1.28.
! Valid to |E(x)| <= 3e-7. Calculation in double precision, result returned 
! in gpu precision.
!
!------------------------------------------------------------------------------

implicit none

real*8,  intent(in) :: x
real*8, intent(out) :: erf_x

real*8              :: a1 = 0.0705230784d0
real*8              :: a2 = 0.0422820123d0
real*8              :: a3 = 0.0092705272d0
real*8              :: a4 = 0.0001520143d0
real*8              :: a5 = 0.0002765672d0
real*8              :: a6 = 0.0000430638d0

erf_x  = 1.0d0 - 1.0d0/(1.0d0 + a1*x    + a2*x**2 +                           &
                                a3*x**3 + a4*x**4 +                           &
                                a5*x**5 + a6*x**6)**16

print*, 'erf_stegun (subroutine) : ', x, erf_x

end subroutine erf_stegun_s

program erf_test

integer :: i

real*8  :: x, result

x = 1.0d0

! Call my ERF(x) "subroutine":

call erf_stegun_s( x, result )

print*, 'Returned result         : ', result

! Call my function version of ERF(x):

result = erf_stegun(x)

print*, 'Returned result         : ', result

! Call PGI intrinsic ERF(x):

result = derf(x)

print*, 'Returned result         : ', result

end program erf_test

If I compile this with intel fortran: f95 -o a.out erf_test.f90, I get the following output:

 erf_stegun (subroutine) :    1.00000000000000       0.842701046333892
 Returned result         :   0.842701046333892
 erf_stegun (function)   :    1.00000000000000       0.842701046333892
 Returned result         :   0.842701046333892
 Returned result         :   0.842700792949715

Note, my own versions of ERF(1.0d0) are the same for the subroutine and the function version I have included, and they are the same as the intel intrinsic ERF(1.0d0) value to 5 d.p.

Now, if I compile with PGI fortran: pgfortran -o a.out erf_test.f90, I get this output:

 erf_stegun (subroutine) :     1.000000000000000        0.8427010463338918
 Returned result         :    0.8427010463338918
 erf_stegun (function)   :     1.000000000000000        0.8427010463338918
 Returned result         :    1.8361739906325170E-010
 Returned result         :   -2.6788761621721014E-015

There are two things that are horribly wrong. Firstly, my function version of ERF is not returning its value correctly (the print statement inside the function prints the right value but the returned value is wrong). The second thing that’s wrong is the double precision intrinsic DERF function, which is clearly producing junk. I’m not sure how “new” the DERF function is in PGI fortran. About a year ago I had to write my own function. Periodic updates must have added in the ERF and DERF intrinsics which at some point ended up replacing my own functions - resulting in a mamoth bug hunt.

Apologies if I’m doing something silly - any advice/info would be greatly appreciated…

Rob.

2

Hi Rob,

The problem is that you haven’t declared derf or erf_stegun so they are being implicitly typed as REAL. If you add ‘implicit none’ in the erf_test program all compilers should let you know that they aren’t declared.

To fix, declare these functions as REAL*8 or compile with “-r8” to change the default kind.

real*8 :: erf_stegun,derf

Note that if you moved the erf_stegun function to a separate file, f95 will most likely fail in the same way. Also, derf is an external lib3f function, not an intrinsic, and why you need to declare it. Other compilers may implicitly declare these functions for you but this would be an extension and not part of the Fortran Standard.

Hope this helps,
Mat

3

D’OH! I’m an eejut. I usually never forget to put implicit none in. Basically, I recently changed from s.p. to d.p. (hence code used to work) - I did change my function from ERF to DERF and assumed that was all I needed to do. I didn’t realise that these functions aren’t intrinsics. Bit annoying that intel declares DERF for you as a non standard feature.

Anyway, thanks again as always Mat…

Rob.

4

Mat,

Does PGI support the F2008* erf yet?

Matt

*Ahh…standards bodies. Never understood why it isn’t renamed for the year it’s actually approved…

5

Does PGI support the F2008* erf yet?

No, we need to finish-up F2003 first. We were suppose to be done with F2003 by now but unfortunately got behind. Almost there though.

Once complete, I think engineering will start prioritizing F2008 features and at least get some of these easy ones in. No time frame yet, though. Co-arrays will be a while.

  • Mat