return 浮点到 return 数组的 Fortran 函数
Fortran function that returns float to return arrays
我有一个来自 scipy 的 Fortran 代码，看起来像这样：
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
我正在 python 中制作一个模块，我希望此函数也能与 numpy 数组一起使用，就像 scipy 那样。
我发现使这项工作起作用的唯一方法是在代码上方创建一个子例程，该子例程采用数组并将每个元素传递给 erf 函数，然后使用 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
用f2py编译后
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> array(0.52049988)
>>> array(0.52049988)
应该是这样的：
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> 0.5204998778130465
>>> array(0.52049988)
但现在我在传递浮点数时失去了精度，结果是一个数字较少的数组。当我传递一个浮点数时，Scipy 以某种方式设法 return 一个浮点数，当我传递一个 numpy 数组时（第二种情况），Scipy 以某种方式设法 return 一个浮点数。我怎样才能得到相同的结果？
这里没有精度问题。你没有失去任何精度。 Python 只是决定在文本表示 array(0.52049988) 中为您打印更少的数字，仅此而已。值相同。
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
我有一个来自 scipy 的 Fortran 代码，看起来像这样：
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
我正在 python 中制作一个模块，我希望此函数也能与 numpy 数组一起使用，就像 scipy 那样。 我发现使这项工作起作用的唯一方法是在代码上方创建一个子例程，该子例程采用数组并将每个元素传递给 erf 函数，然后使用 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
用f2py编译后
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> array(0.52049988)
>>> array(0.52049988)
应该是这样的：
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> 0.5204998778130465
>>> array(0.52049988)
但现在我在传递浮点数时失去了精度，结果是一个数字较少的数组。当我传递一个浮点数时，Scipy 以某种方式设法 return 一个浮点数，当我传递一个 numpy 数组时（第二种情况），Scipy 以某种方式设法 return 一个浮点数。我怎样才能得到相同的结果？
这里没有精度问题。你没有失去任何精度。 Python 只是决定在文本表示 array(0.52049988) 中为您打印更少的数字，仅此而已。值相同。
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