PGI User Forum
 SearchSearch   MemberlistMemberlist     RegisterRegister   ProfileProfile    Log inLog in 

Free OpenACC Webinar

Something horrbily wrong with ERF(x)

 
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming
View previous topic :: View next topic  
Author Message
alfvenwave



Joined: 08 Apr 2010
Posts: 79

PostPosted: Wed Apr 27, 2011 5:52 am    Post subject: Something horrbily wrong with ERF(x) Reply with quote

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:

Code:

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:

Code:
 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:

Code:

 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.
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 6213
Location: The Portland Group Inc.

PostPosted: Wed Apr 27, 2011 11:30 am    Post subject: Reply with quote

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.
Quote:
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
Back to top
View user's profile
alfvenwave



Joined: 08 Apr 2010
Posts: 79

PostPosted: Thu Apr 28, 2011 1:35 am    Post subject: Reply with quote

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.
Back to top
View user's profile
TheMatt



Joined: 06 Jul 2009
Posts: 322
Location: Greenbelt, MD

PostPosted: Thu Apr 28, 2011 4:44 am    Post subject: Reply with quote

mkcolg wrote:
Other compilers may implicitly declare these functions for you but this would be an extension and not part of the Fortran Standard.

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...
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 6213
Location: The Portland Group Inc.

PostPosted: Thu Apr 28, 2011 3:55 pm    Post subject: Reply with quote

Quote:
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
Back to top
View user's profile
Display posts from previous:   
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming All times are GMT - 7 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © phpBB Group