Fortran function that returns float to return arrays

Question

I have a fortran code from scipy that looks like this:

erf.f

      DOUBLE PRECISION FUNCTION ERF(x)
C-----------------------------------------------------------------------
C             EVALUATION OF THE REAL ERROR FUNCTION
C-----------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ax,bot,c,t,top,x2
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,exp,sign
C     ..
C     .. Data statements ..
C-------------------------
C-------------------------
C-------------------------
C-------------------------
      DATA c/.564189583547756D0/
      DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/,
     +     a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/,
     +     a(5)/.128379167095513D+00/
      DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/,
     +     b(3)/.375795757275549D+00/
      DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/,
     +     p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/,
     +     p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/,
     +     p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/
      DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/,
     +     q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/,
     +     q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/,
     +     q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/
      DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/,
     +     r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/,
     +     r(5)/2.82094791773523D-01/
      DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/,
     +     s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/
C     ..
C     .. Executable Statements ..
C-------------------------
      ax = abs(x)
      IF (ax.GT.0.5D0) GO TO 10
      t = x*x
      top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0
      bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0
      erf = x* (top/bot)
      RETURN
C
   10 IF (ax.GT.4.0D0) GO TO 20
      top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+
     +      p(7))*ax + p(8)
      bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+
     +      q(7))*ax + q(8)
      erf = 0.5D0 + (0.5D0-exp(-x*x)*top/bot)
      IF (x.LT.0.0D0) erf2 = -erf2
      RETURN
C
   20 IF (ax.GE.5.8D0) GO TO 30
      x2 = x*x
      t = 1.0D0/x2
      top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5)
      bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0
      erf = (c-top/ (x2*bot))/ax
      erf = 0.5D0 + (0.5D0-exp(-x2)*erf)
      IF (x.LT.0.0D0) erf = -erf
      RETURN
C
   30 erf = sign(1.0D0,x)
      RETURN

      END

I'm making a module in python and I want this function to work with numpy arrays too, like scipy does. The only way I found that make this work is creating a subroutine above the code which takes an array and every element is passed to the erf function, and then compile with f2py.

erfmod.f

      subroutine erf(a,n)
      implicit none

      integer :: n,i
      real*8 :: a(n)
Cf2py intent(in,out,copy) :: a
cf2py integer intent(hide),depend(a) :: n=len(a)

      do i=1, n
            a(i) = erf2(a(i))
      end do

      contains
      DOUBLE PRECISION FUNCTION erf2(x)
C-----------------------------------------------------------------------
C             EVALUATION OF THE REAL ERROR FUNCTION
C-----------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ax,bot,c,t,top,x2
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,exp,sign
C     ..
C     .. Data statements ..
C-------------------------
C-------------------------
C-------------------------
C-------------------------
      DATA c/.564189583547756D0/
      DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/,
     +     a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/,
     +     a(5)/.128379167095513D+00/
      DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/,
     +     b(3)/.375795757275549D+00/
      DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/,
     +     p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/,
     +     p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/,
     +     p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/
      DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/,
     +     q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/,
     +     q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/,
     +     q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/
      DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/,
     +     r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/,
     +     r(5)/2.82094791773523D-01/
      DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/,
     +     s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/
C     ..
C     .. Executable Statements ..
C-------------------------
      ax = abs(x)
      IF (ax.GT.0.5D0) GO TO 10
      t = x*x
      top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0
      bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0
      erf2 = x* (top/bot)
      RETURN
C
   10 IF (ax.GT.4.0D0) GO TO 20
      top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+
     +      p(7))*ax + p(8)
      bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+
     +      q(7))*ax + q(8)
      erf2 = 0.5D0 + (0.5D0-exp(-x*x)*top/bot)
      IF (x.LT.0.0D0) erf2 = -erf2
      RETURN
C
   20 IF (ax.GE.5.8D0) GO TO 30
      x2 = x*x
      t = 1.0D0/x2
      top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5)
      bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0
      erf2 = (c-top/ (x2*bot))/ax
      erf2 = 0.5D0 + (0.5D0-exp(-x2)*erf2)
      IF (x.LT.0.0D0) erf2 = -erf2
      RETURN
C
   30 erf2 = sign(1.0D0,x)
      RETURN

      end function erf2

      end subroutine

after compiling with f2py

import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))

>>> array(0.52049988)
>>> array(0.52049988)

shoud look like this:

import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))

>>> 0.5204998778130465
>>> array(0.52049988)

But now I lost the precision when I'm passing a float, the result is an array with less digits. Scipy somehow manage to return a float when I'm passing a float, and an array when I pass a numpy array (the second case). How can I get the same result?

What same result? Please show your result and show your Python code. Python `float` is the same precision as Fortran `DOUBLE PRECISION`. Be aware that there is a Fortran intrinsic named [`ERF()`](https://gcc.gnu.org/onlinedocs/gfortran/ERF.html) for the error function. — Vladimir F Героям слава, Apr 03 '22 at 15:44
I added the python code, the intrinsic erf shouldn't be a problem here — Tito Diego, Apr 03 '22 at 15:59
And how the the results look like? You wrote that you lost some precision. How exactly? The erf intrinsic shouldn't be a problem, but it might be an alternative. — Vladimir F Героям слава, Apr 03 '22 at 17:02
I don't know python at all (...) but searching for "print numpy array precision" gives me https://numpy.org/doc/stable/reference/generated/numpy.set_printoptions.html as the first hit. There it says "Number of digits of precision for floating point output (default 8)." Are 8 figures being reported just because that is the default number to print, rather than any precision issues? https://stackoverflow.com/questions/2891790/how-to-pretty-print-a-numpy-array-without-scientific-notation-and-with-given-pre might also be of use — Ian Bush, Apr 03 '22 at 17:28

Vladimir F Героям слава · Accepted Answer · 2022-04-03T17:44:52.137

There are no precision issues here. You did not lost any precision. Python just decided to print fewer digits for you in the text representation array(0.52049988), that's all. The values are identical.

In [1]: import numpy as np                                                                                                                                                                                

In [2]: np.array(0.5204998778130465)                                                                                                                                                                      
Out[2]: array(0.52049988)

In [3]: a = np.array(0.5204998778130465)                                                                                                                                                                 

In [4]: a                                                                                                                                                                                                
Out[4]: array(0.52049988)

In [5]: a.item()                                                                                                                                                                                         
Out[5]: 0.5204998778130465

you're right, I was confused because in the scipy implementation gives you the whole number, thanks for your help — Tito Diego, Apr 03 '22 at 17:55

Fortran function that returns float to return arrays

1 Answers1