Fortran、Numpy、Cython 和 Numexpr 的性能比较
Performance comparison Fortran, Numpy,Cython and Numexpr
我有以下功能：
def get_denom(n_comp,qs,x,cp,cs):
'''
len(n_comp) = 1 # number of proteins
len(cp) = n_comp # protein concentration
len(qp) = n_comp # protein capacity
len(x) = 3*n_comp + 1 # fit parameters
len(cs) = 1
'''
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
我将它与 Fortran 实现（我的第一个 Fortran 函数）进行比较：
subroutine get_denom (qs,x,cp,cs,n_comp,denom)
! Calculates the denominator in the SMA model (Brooks and Cramer 1992)
! The function is called at a specific salt concentration and isotherm point
! I loops over the number of components
implicit none
! declaration of input variables
integer, intent(in) :: n_comp ! number of components
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration
double precision, dimension(n_comp), INTENT(IN) ::cp ! protein concentration
double precision, dimension(3*n_comp + 1), INTENT(IN) :: x ! parameters
! declaration of local variables
double precision, dimension(n_comp) :: k,sigma,z
double precision :: a
integer :: i
! declaration of outpur variables
double precision, intent(out) :: denom
k = x(1:n_comp) ! equlibrium constant
sigma = x(n_comp+1:2*n_comp) ! steric hindrance factor
z = x(2*n_comp+1:3*n_comp) ! charge of protein
a = 0.
do i = 1,n_comp
a = a + (sigma(i) + z(i))*(k(i)*(qs/cs)**(z(i)-1.))*cp(i)
end do
denom = a + cs
end subroutine get_denom
我使用以下方法编译了 .f95 文件：
1) f2py -c -m get_denom get_denom.f95 --fcompiler=gfortran
2) f2py -c -m get_denom_vec get_denom.f95 --fcompiler=gfortran --f90flags='-msse2'（最后一个选项应该打开自动矢量化）
我通过以下方式测试功能：
import numpy as np
import get_denom as fort_denom
import get_denom_vec as fort_denom_vec
from matplotlib import pyplot as plt
%matplotlib inline
def get_denom(n_comp,qs,x,cp,cs):
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
n_comp = 100
cp = np.tile(1.243,n_comp)
cs = 100.
qs = np.tile(1100.,n_comp)
x= np.random.rand(3*n_comp+1)
denom = np.empty(1)
%timeit get_denom(n_comp,qs,x,cp,cs)
%timeit fort_denom.get_denom(qs,x,cp,cs,n_comp)
%timeit fort_denom_vec.get_denom(qs,x,cp,cs,n_comp)
我添加了以下 Cython 代码：
import cython
# import both numpy and the Cython declarations for numpy
import numpy as np
cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
def get_denom(int n_comp,np.ndarray[double, ndim=1, mode='c'] qs, np.ndarray[double, ndim=1, mode='c'] x,np.ndarray[double, ndim=1, mode='c'] cp, double cs):
cdef int i
cdef double a
cdef double denom
cdef double[:] k = x[0:n_comp]
cdef double[:] sigma = x[n_comp:2*n_comp]
cdef double[:] z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = 0.
for i in range(n_comp):
#a += (sigma[i] + z[i])*( pow( k[i]*(qs[i]/cs), (z[i]-1) ) )*cp[i]
a += (sigma[i] + z[i])*( k[i]*(qs[i]/cs)**(z[i]-1) )*cp[i]
denom = a + cs
return denom
编辑：
添加了 Numexpr，使用一个线程：
def get_denom_numexp(n_comp,qs,x,cp,cs):
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = ne.evaluate('(sigma + z)*( k*(qs/cs)**(z-1) )*cp' )
return cs + np.sum(a)
ne.set_num_threads(1) # using just 1 thread
%timeit get_denom_numexp(n_comp,qs,x,cp,cs)
结果是（越小越好）：
为什么随着数组大小的增加，Fortran 的速度越来越接近 Numpy？我怎样才能加快 Cython 的速度？使用指针？
注释中的信息不足，但以下一些内容可能会有所帮助：
1) Fortran 优化了内部函数，例如 "Sum()" 和 "Dot_Product"，您可能希望使用它们代替 Do 循环进行求和等
在某些情况下（这里不一定），可能 "better" 使用 ForAll 或其他创建要求和的 "meta" 数组，然后将求和应用于 "meta"数组。
但是，Fortran 允许数组部分，因此您不需要创建 automatic/intermediate 数组 sigma、k 和 z，以及相关的开销。相反可以有类似
的东西
n_compP1 = n_comp+1
n_compT2 = n_comp*2
a = Sum( x(1:n_comp)+2*x(n_compP1,n_compT2) ) ! ... just for example
2) 有时（取决于编译器、机器等），如果数组大小不是特定的 "binary intervals"（例如 1024 与 1000）等，则可能存在 "memory contentions" 等
您可能希望在图表中的更多点（即其他 "n_comps"）重复您的实验，尤其是在 "boundaries".
附近
3) 无法确定您是否正在使用完整的编译器优化（标志）来编译您的 Fortran 代码。您可能希望查找各种“-o”标志等。
4) 您可能希望包含 OpemMP 指令（或者至少在您的标志中包含 openmp 等）。这有时可以改善某些开销问题，即使没有明确依赖循环中的 OpenMP 指令等。
5) 一般：这可能适用于您使用循环的每个方法
a) "summation" 公式中的 "constant operations" 可以在循环外执行，例如。创建类似 qsDcs = qs/cs 的东西，并在循环中使用 qsDcs。
b) 同样，有时创建类似 zM1(:) = z(:) - 1 的东西并在循环中使用 zM1(:) 是有用的。
根据我之前的回答和 Vladimir 的薄弱推测，我设置了两个 s/r：一个作为原始给定，一个使用数组部分和 Sum()。我还想证明 Vladimir 对 Do 循环优化的评论是薄弱的。
此外，我通常会在基准测试中提出一点，上面示例中 n_comp 的大小太小了。下面的结果将每个 "original" 和 "better" SumArraySection (SAS) 变体放入在计时调用中重复 1,000 次的循环中，因此结果是每个 s/r 的 1000 次计算。如果您的时间是几分之一秒，则它们可能不可靠。
有许多变体值得考虑，none 带有显式指针。此插图使用的一种变体是
subroutine get_denomSAS (qs,x,cp,cs,n_comp,denom)
! Calculates the denominator in the SMA model (Brooks and Cramer 1992)
! The function is called at a specific salt concentration and isotherm point
! I loops over the number of components
implicit none
! declaration of input variables
integer, intent(in) :: n_comp ! number of components
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration
double precision, Intent(In) :: cp(:)
double precision, Intent(In) :: x(:)
! declaration of local variables
integer :: i
! declaration of outpur variables
double precision, intent(out) :: denom
!
!
double precision :: qsDcs
!
!
qsDcs = qs/cs
!
denom = Sum( (x(n_comp+1:2*n_comp) + x(2*n_comp+1:3*n_comp))*(x(1:n_comp)*(qsDcs) &
**(x(2*n_comp+1:3*n_comp)-1))*cp(1:n_comp) ) + cs
!
!
end subroutine get_denomSAS
主要区别是：
a) 传递的数组是 (:)
b) s/r 中没有数组赋值，而是使用数组部分（相当于 "effective " 指针）。
c) 使用 Sum() 而不是 Do
然后还尝试两种不同的编译器优化来证明其含义。
如两个图表所示，原始代码（蓝色菱形）要慢得多 c.f。低优化的 SAS（红色方块）。 SAS 在高度优化的情况下仍然更好，但它们正在接近。当使用低编译器优化时，Sum() 是 "better optimised" 的部分原因。
黄线表示两个 s/r 的时间比率。忽略顶部图像中“0”处的黄线值（n_comp太小导致其中一个时间不稳定）
由于我没有用户的原始数据与 Numpy 的比率，我只能声明他图表上的 SAS 曲线应该低于他当前的 Fortran 结果，并且可能更平坦甚至下降趋势。
换句话说，可能实际上不存在原始帖子中看到的分歧，或者至少没有到那个程度。
...虽然更多的实验可能有助于证明其他 comments/answers 已经提供。
亲爱的莫里茨：哎呀，我忘了提，关于你关于指针的问题。如前所述，SAS 变体得到改进的一个关键原因是它更好地利用了 "effective pointers"，因为它避免了将数组 x() 重新分配给三个新的本地数组的需要（即因为 x 由参考，使用数组部分是 Fortran 内置的一种指针方法，因此不需要显式指针），但随后需要 Sum() 或 Dot_Product() 或其他任何东西。
相反，您可以保留 Do 并通过将 x 更改为 n_compx3 二维数组或直接传递顺序为 n_comp 的三个显式一维数组来实现类似的效果。这个决定很可能是由代码的大小和复杂性驱动的，因为它需要重写 calling/sr 语句等，以及使用 x() 的任何其他地方。我们的一些项目超过 300,000 行代码，因此在这种情况下，在本地更改代码（例如 SAS 等）的成本要低得多。
我仍在等待获得在我们的一台机器上安装 Numpy 的许可。如前所述，有趣的是为什么您的相对时间暗示 Numpy 随着 n_comp ...
的增加而改进
当然，关于 "proper" 基准测试等的评论，以及您对 fpy 的使用暗示了哪些编译器开关的问题，仍然适用，因为这些可能会极大地改变您的结果的特征。
如果您的结果针对这些排列进行了更新，我很想看看您的结果。
被其他答案点名，必须回复
我知道这并不能真正回答最初的问题，但发帖者在他的评论中鼓励朝着这个方向努力。
我的观点是：
1. 我不相信数组内在函数在任何方面都更好地优化。如果幸运的话，它们会被翻译成与手动循环相同的循环代码。否则，可能会出现性能问题。因此，必须要小心。有可能触发临时数组。
我将提供的 SAS 数组转换为通常的 do 循环。我称之为 DOS。我证明了 DO 循环一点也不慢，在这种情况下，两个子例程或多或少会产生相同的代码。
qsDcs = qs/cs
denom = 0
do j = 1, n_comp
denom = denom + (x(n_comp+j) + x(2*n_comp+j)) * (x(j)*(qsDcs)**(x(2*n_comp+j)-1))*cp(j)
end do
denom = denom + cs
重要的是要说，我不认为这个版本的可读性差只是因为它多了一两行。其实也很简单，看看这里发生了什么。
现在是这些的时间安排
f2py -c -m sas sas.f90 --opt='-Ofast'
f2py -c -m dos dos.f90 --opt='-Ofast'
In [24]: %timeit test_sas(10000)
1000 loops, best of 3: 796 µs per loop
In [25]: %timeit test_sas(10000)
1000 loops, best of 3: 793 µs per loop
In [26]: %timeit test_dos(10000)
1000 loops, best of 3: 795 µs per loop
In [27]: %timeit test_dos(10000)
1000 loops, best of 3: 797 µs per loop
他们是一样的。数组内在函数和数组表达式算术中没有隐藏的性能魔法。在这种情况下，它们只是被转换为引擎盖下的循环。
如果你检查生成的 GIMPLE 代码，SAS 和 DOS 都被翻译成相同的优化代码主块，这里没有调用 SUM() 的神奇版本：
<bb 8>:
# val.8_59 = PHI <val.8_49(9), 0.0(7)>
# ivtmp.18_123 = PHI <ivtmp.18_122(9), 0(7)>
# ivtmp.25_121 = PHI <ivtmp.25_120(9), ivtmp.25_117(7)>
# ivtmp.28_116 = PHI <ivtmp.28_115(9), ivtmp.28_112(7)>
_111 = (void *) ivtmp.25_121;
_32 = MEM[base: _111, index: _106, step: 8, offset: 0B];
_36 = MEM[base: _111, index: _99, step: 8, offset: 0B];
_37 = _36 + _32;
_40 = MEM[base: _111, offset: 0B];
_41 = _36 - 1.0e+0;
_42 = __builtin_pow (qsdcs_18, _41);
_97 = (void *) ivtmp.28_116;
_47 = MEM[base: _97, offset: 0B];
_43 = _40 * _47;
_44 = _43 * _42;
_48 = _44 * _37;
val.8_49 = val.8_59 + _48;
ivtmp.18_122 = ivtmp.18_123 + 1;
ivtmp.25_120 = ivtmp.25_121 + _118;
ivtmp.28_115 = ivtmp.28_116 + _113;
if (ivtmp.18_122 == _96)
goto <bb 10>;
else
goto <bb 9>;
<bb 9>:
goto <bb 8>;
<bb 10>:
# val.8_13 = PHI <val.8_49(8), 0.0(6)>
_51 = val.8_13 + _17;
*denom_52(D) = _51;
代码在功能上与 do 循环版本相同，只是变量的名称不同。
2. 他们假设形状数组参数 (:) 有可能降低性能。尽管在假定大小参数 (*) 或显式大小 (n) 中接收到的参数总是简单连续的，但理论上不必是假定形状，编译器必须为此做好准备。因此，我总是建议将 contiguous 属性用于您假定的形状参数，只要您知道您将使用连续数组调用它。
在另一个答案中，更改毫无意义，因为它没有使用假定形状参数的任何优点。也就是说，您不必传递带有数组大小的参数，您可以使用 size() 和 shape().
等内部函数
综合比较结果如下。我让它尽可能公平。 Fortran文件用-Ofast编译如上图：
import numpy as np
import dos as dos
import sas as sas
from matplotlib import pyplot as plt
import timeit
import numexpr as ne
#%matplotlib inline
ne.set_num_threads(1)
def test_n(n_comp):
cp = np.tile(1.243,n_comp)
cs = 100.
qs = np.tile(1100.,n_comp)
x= np.random.rand(3*n_comp+1)
def test_dos():
denom = np.empty(1)
dos.get_denomsas(qs,x,cp,cs,n_comp)
def test_sas():
denom = np.empty(1)
sas.get_denomsas(qs,x,cp,cs,n_comp)
def get_denom():
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
def get_denom_numexp():
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
loc_cp = cp
loc_cs = cs
loc_qs = qs
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = ne.evaluate('(sigma + z)*( k*(loc_qs/loc_cs)**(z-1) )*loc_cp' )
return cs + np.sum(a)
print 'py', timeit.Timer(get_denom).timeit(1000000/n_comp)
print 'dos', timeit.Timer(test_dos).timeit(1000000/n_comp)
print 'sas', timeit.Timer(test_sas).timeit(1000000/n_comp)
print 'ne', timeit.Timer(get_denom_numexp).timeit(1000000/n_comp)
def test():
for n in [10,100,1000,10000,100000,1000000]:
print "-----"
print n
test_n(n)
结果：
py dos sas numexpr
10 11.2188110352 1.8704519272 1.8659651279 28.6881871223
100 1.6688809395 0.6675260067 0.667083025 3.4943861961
1000 0.7014708519 0.5406000614 0.5441288948 0.9069931507
10000 0.5825948715 0.5269498825 0.5309231281 0.6178650856
100000 0.5736029148 0.526198864 0.5304090977 0.5886831284
1000000 0.6355218887 0.5294830799 0.5366530418 0.5983200073
10000000 0.7903120518 0.5301260948 0.5367569923 0.6030929089
您可以看到两个 Fortran 版本之间的差异非常小。数组语法稍微慢一点，但真的没什么好说的。
结论 1：在这个比较中，所有的开销应该是公平的，你会看到对于理想的向量长度 Numpy 和 Numexpr CAN 几乎达到 Fortran 的性能，但是当向量太小甚至可能太大时，Python 解决方案的开销占上风。
结论2：另一个比较中SAS版本较高是与原OP版本不等价造成的。我的答案上面包含等效的优化 DO 循环版本。
怀疑它。
好的，最后，我们获准在我们的一台机器上安装 Numpy 等，这可能是对您的原始内容的全面解释 post。
简短的回答是，从某种意义上说，您最初的问题是 "the wrong question"。此外，其中一位贡献者有很多无理取闹的辱骂和错误信息，这些错误和捏造值得引起注意，以免有人误以为他们，并为此付出代价。
此外，我决定将此答案作为一个单独的答案提交，而不是编辑我 4 月 14 日的答案，原因见下文，并且是正当的。
A 部分：OP 的答案
首先，处理原文中的问题 post：您可能还记得我只能对 Fortran 方面发表评论，因为我们的政策对可以安装什么软件以及我们的位置有严格的规定机器，我们手头没有 Python 等（直到刚才）。我还反复声明，就我们所说的曲线特征或 "concavity".
而言，你的结果的特征很有趣
此外，纯粹使用 "relative" 结果（因为你没有 post 绝对时间，而且我当时手头没有 Numpy），我已经指出了几次一些重要的信息可能潜伏在其中。
正是如此
首先，我想确保我可以重现你的结果，因为我们通常不使用 Python/F2py，你的结果中隐含的编译器设置等并不明显，所以我进行了各种的测试以确保我们正在谈论苹果对苹果（我 4 月 14 日的回答表明 Debug 与 Release/O2 有很大的不同）。
图 1 显示了我的 Python 三种情况的结果：你原来的 Python/Numpy 内部子程序（称之为 P，我只是 cut/pasted 你的原始程序），你的基于 Fortran 的原始 Do s/r 你已经 posted（称之为 FDo，我只是 copy/pasted 你的原始，和以前一样），以及我之前建议的依赖数组部分的变体之一，因此需要 Sum()（调用此 FSAS，通过编辑您的原始 FDo 创建）。图 1 显示了通过 timeit 的绝对时序。
图 2 显示了根据您的相对策略除以 Python/Numpy (P) 时序的相对结果。仅显示了两个（相对的）Fortran 变体。
很明显，那些重现了您原始情节中的角色，我们可以确信我的测试似乎与您的测试一致。
现在，您原来的问题是 "Why is is the speed of Fortran getting closer to Numpy with increasing size of the arrays?"。
其实不然。这纯粹是 artefact/distortion 纯粹依靠 "relative" 时间安排可能给人的印象。
查看图 1，从绝对时间来看，很明显 Numpy 和 Fortran 存在分歧。因此，事实上，如果您愿意，Fortran 结果正在远离 Numpy，反之亦然。
一个更好的问题，也是我之前评论中反复出现的一个问题，为什么这些首先是向上弯曲的（例如 c.f. 线性）？我之前的 Fortran-only 结果显示 "mostly" 持平的相对性能比（我的 4 月 14 日 chart/answer 中的黄线），所以我推测 Python 方面发生了一些有趣的事情并建议检查一下。
显示这一点的一种方法是使用另一种相对度量。我将每个（绝对）系列与其自身的最高值（即 n_comp = 10k）分开，以查看此 "internal relative" 性能如何展开（这些被称为 ??10k 值，代表时间n_comp = 10,000）。图 3 将 P、FDo 和 FSAS 的这些结果显示为 P/P10k、FDo/FDo10k 和 FSAS/FSAS10k。为清楚起见，y 轴采用对数刻度。很明显，随着 n_comp c.f 的降低，Fortran 变体的表现相对更好。 P 结果（例如，红色圆圈注释部分）。
换句话说，Fortran 在减小数组大小方面非常非常（非线性）更有效。不确定为什么 Python 会随着 n_comp 的减少而变得更糟......但它确实存在，并且可能是内部 overhead/set-up 等的问题，以及口译员与口译员之间的差异编译器等
因此，Fortran 并不是 "catching up" 与 Python，恰恰相反，它正在继续与 Python 保持距离（见图 1）。然而，Python 和 Fortran 之间的非线性差异随着 n_comp 的减少产生了 "relative" 具有明显反直觉和非线性特征的时序。
因此，随着 n_comp 的增加和每个方法 "stabilises" 或多或少的线性模式，曲线变平表明它们的差异呈线性增长，并且相对比率趋于近似常量（忽略内存争用、smp 问题等）...如果允许 n_comp > 10k，这更容易看出，但我 4 月 14 日的回答中的黄线已经显示了 Fortran-only s/r的。
旁白：我的偏好是创建自己的时间 routines/functions。 timeit 看起来很方便，但是"black box"里面有很多事情要做。设置您自己的循环和结构，并确定您的计时函数的属性s/r解决方案对于正确评估很重要。
我有以下功能：
def get_denom(n_comp,qs,x,cp,cs):
'''
len(n_comp) = 1 # number of proteins
len(cp) = n_comp # protein concentration
len(qp) = n_comp # protein capacity
len(x) = 3*n_comp + 1 # fit parameters
len(cs) = 1
'''
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
我将它与 Fortran 实现（我的第一个 Fortran 函数）进行比较：
subroutine get_denom (qs,x,cp,cs,n_comp,denom)
! Calculates the denominator in the SMA model (Brooks and Cramer 1992)
! The function is called at a specific salt concentration and isotherm point
! I loops over the number of components
implicit none
! declaration of input variables
integer, intent(in) :: n_comp ! number of components
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration
double precision, dimension(n_comp), INTENT(IN) ::cp ! protein concentration
double precision, dimension(3*n_comp + 1), INTENT(IN) :: x ! parameters
! declaration of local variables
double precision, dimension(n_comp) :: k,sigma,z
double precision :: a
integer :: i
! declaration of outpur variables
double precision, intent(out) :: denom
k = x(1:n_comp) ! equlibrium constant
sigma = x(n_comp+1:2*n_comp) ! steric hindrance factor
z = x(2*n_comp+1:3*n_comp) ! charge of protein
a = 0.
do i = 1,n_comp
a = a + (sigma(i) + z(i))*(k(i)*(qs/cs)**(z(i)-1.))*cp(i)
end do
denom = a + cs
end subroutine get_denom
我使用以下方法编译了 .f95 文件：
1) f2py -c -m get_denom get_denom.f95 --fcompiler=gfortran
2) f2py -c -m get_denom_vec get_denom.f95 --fcompiler=gfortran --f90flags='-msse2'（最后一个选项应该打开自动矢量化）
我通过以下方式测试功能：
import numpy as np
import get_denom as fort_denom
import get_denom_vec as fort_denom_vec
from matplotlib import pyplot as plt
%matplotlib inline
def get_denom(n_comp,qs,x,cp,cs):
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
n_comp = 100
cp = np.tile(1.243,n_comp)
cs = 100.
qs = np.tile(1100.,n_comp)
x= np.random.rand(3*n_comp+1)
denom = np.empty(1)
%timeit get_denom(n_comp,qs,x,cp,cs)
%timeit fort_denom.get_denom(qs,x,cp,cs,n_comp)
%timeit fort_denom_vec.get_denom(qs,x,cp,cs,n_comp)
我添加了以下 Cython 代码：
import cython
# import both numpy and the Cython declarations for numpy
import numpy as np
cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
def get_denom(int n_comp,np.ndarray[double, ndim=1, mode='c'] qs, np.ndarray[double, ndim=1, mode='c'] x,np.ndarray[double, ndim=1, mode='c'] cp, double cs):
cdef int i
cdef double a
cdef double denom
cdef double[:] k = x[0:n_comp]
cdef double[:] sigma = x[n_comp:2*n_comp]
cdef double[:] z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = 0.
for i in range(n_comp):
#a += (sigma[i] + z[i])*( pow( k[i]*(qs[i]/cs), (z[i]-1) ) )*cp[i]
a += (sigma[i] + z[i])*( k[i]*(qs[i]/cs)**(z[i]-1) )*cp[i]
denom = a + cs
return denom
编辑：
添加了 Numexpr，使用一个线程：
def get_denom_numexp(n_comp,qs,x,cp,cs):
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = ne.evaluate('(sigma + z)*( k*(qs/cs)**(z-1) )*cp' )
return cs + np.sum(a)
ne.set_num_threads(1) # using just 1 thread
%timeit get_denom_numexp(n_comp,qs,x,cp,cs)
结果是（越小越好）：
为什么随着数组大小的增加，Fortran 的速度越来越接近 Numpy？我怎样才能加快 Cython 的速度？使用指针？
注释中的信息不足，但以下一些内容可能会有所帮助：
1) Fortran 优化了内部函数，例如 "Sum()" 和 "Dot_Product"，您可能希望使用它们代替 Do 循环进行求和等
在某些情况下（这里不一定），可能 "better" 使用 ForAll 或其他创建要求和的 "meta" 数组，然后将求和应用于 "meta"数组。
但是，Fortran 允许数组部分，因此您不需要创建 automatic/intermediate 数组 sigma、k 和 z，以及相关的开销。相反可以有类似
的东西n_compP1 = n_comp+1
n_compT2 = n_comp*2
a = Sum( x(1:n_comp)+2*x(n_compP1,n_compT2) ) ! ... just for example
2) 有时（取决于编译器、机器等），如果数组大小不是特定的 "binary intervals"（例如 1024 与 1000）等，则可能存在 "memory contentions" 等
您可能希望在图表中的更多点（即其他 "n_comps"）重复您的实验，尤其是在 "boundaries".
附近3) 无法确定您是否正在使用完整的编译器优化（标志）来编译您的 Fortran 代码。您可能希望查找各种“-o”标志等。
4) 您可能希望包含 OpemMP 指令（或者至少在您的标志中包含 openmp 等）。这有时可以改善某些开销问题，即使没有明确依赖循环中的 OpenMP 指令等。
5) 一般：这可能适用于您使用循环的每个方法
a) "summation" 公式中的 "constant operations" 可以在循环外执行，例如。创建类似 qsDcs = qs/cs 的东西，并在循环中使用 qsDcs。
b) 同样，有时创建类似 zM1(:) = z(:) - 1 的东西并在循环中使用 zM1(:) 是有用的。
根据我之前的回答和 Vladimir 的薄弱推测，我设置了两个 s/r：一个作为原始给定，一个使用数组部分和 Sum()。我还想证明 Vladimir 对 Do 循环优化的评论是薄弱的。
此外，我通常会在基准测试中提出一点，上面示例中 n_comp 的大小太小了。下面的结果将每个 "original" 和 "better" SumArraySection (SAS) 变体放入在计时调用中重复 1,000 次的循环中，因此结果是每个 s/r 的 1000 次计算。如果您的时间是几分之一秒，则它们可能不可靠。
有许多变体值得考虑，none 带有显式指针。此插图使用的一种变体是
subroutine get_denomSAS (qs,x,cp,cs,n_comp,denom)
! Calculates the denominator in the SMA model (Brooks and Cramer 1992)
! The function is called at a specific salt concentration and isotherm point
! I loops over the number of components
implicit none
! declaration of input variables
integer, intent(in) :: n_comp ! number of components
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration
double precision, Intent(In) :: cp(:)
double precision, Intent(In) :: x(:)
! declaration of local variables
integer :: i
! declaration of outpur variables
double precision, intent(out) :: denom
!
!
double precision :: qsDcs
!
!
qsDcs = qs/cs
!
denom = Sum( (x(n_comp+1:2*n_comp) + x(2*n_comp+1:3*n_comp))*(x(1:n_comp)*(qsDcs) &
**(x(2*n_comp+1:3*n_comp)-1))*cp(1:n_comp) ) + cs
!
!
end subroutine get_denomSAS
主要区别是：
a) 传递的数组是 (:) b) s/r 中没有数组赋值，而是使用数组部分（相当于 "effective " 指针）。 c) 使用 Sum() 而不是 Do
然后还尝试两种不同的编译器优化来证明其含义。
如两个图表所示，原始代码（蓝色菱形）要慢得多 c.f。低优化的 SAS（红色方块）。 SAS 在高度优化的情况下仍然更好，但它们正在接近。当使用低编译器优化时，Sum() 是 "better optimised" 的部分原因。
黄线表示两个 s/r 的时间比率。忽略顶部图像中“0”处的黄线值（n_comp太小导致其中一个时间不稳定）
由于我没有用户的原始数据与 Numpy 的比率，我只能声明他图表上的 SAS 曲线应该低于他当前的 Fortran 结果，并且可能更平坦甚至下降趋势。
换句话说，可能实际上不存在原始帖子中看到的分歧，或者至少没有到那个程度。
...虽然更多的实验可能有助于证明其他 comments/answers 已经提供。
亲爱的莫里茨：哎呀，我忘了提，关于你关于指针的问题。如前所述，SAS 变体得到改进的一个关键原因是它更好地利用了 "effective pointers"，因为它避免了将数组 x() 重新分配给三个新的本地数组的需要（即因为 x 由参考，使用数组部分是 Fortran 内置的一种指针方法，因此不需要显式指针），但随后需要 Sum() 或 Dot_Product() 或其他任何东西。
相反，您可以保留 Do 并通过将 x 更改为 n_compx3 二维数组或直接传递顺序为 n_comp 的三个显式一维数组来实现类似的效果。这个决定很可能是由代码的大小和复杂性驱动的，因为它需要重写 calling/sr 语句等，以及使用 x() 的任何其他地方。我们的一些项目超过 300,000 行代码，因此在这种情况下，在本地更改代码（例如 SAS 等）的成本要低得多。
我仍在等待获得在我们的一台机器上安装 Numpy 的许可。如前所述，有趣的是为什么您的相对时间暗示 Numpy 随着 n_comp ...
的增加而改进当然，关于 "proper" 基准测试等的评论，以及您对 fpy 的使用暗示了哪些编译器开关的问题，仍然适用，因为这些可能会极大地改变您的结果的特征。
如果您的结果针对这些排列进行了更新，我很想看看您的结果。
被其他答案点名，必须回复
我知道这并不能真正回答最初的问题，但发帖者在他的评论中鼓励朝着这个方向努力。
我的观点是：
1. 我不相信数组内在函数在任何方面都更好地优化。如果幸运的话，它们会被翻译成与手动循环相同的循环代码。否则，可能会出现性能问题。因此，必须要小心。有可能触发临时数组。
我将提供的 SAS 数组转换为通常的 do 循环。我称之为 DOS。我证明了 DO 循环一点也不慢，在这种情况下，两个子例程或多或少会产生相同的代码。
qsDcs = qs/cs
denom = 0
do j = 1, n_comp
denom = denom + (x(n_comp+j) + x(2*n_comp+j)) * (x(j)*(qsDcs)**(x(2*n_comp+j)-1))*cp(j)
end do
denom = denom + cs
重要的是要说，我不认为这个版本的可读性差只是因为它多了一两行。其实也很简单，看看这里发生了什么。
现在是这些的时间安排
f2py -c -m sas sas.f90 --opt='-Ofast'
f2py -c -m dos dos.f90 --opt='-Ofast'
In [24]: %timeit test_sas(10000)
1000 loops, best of 3: 796 µs per loop
In [25]: %timeit test_sas(10000)
1000 loops, best of 3: 793 µs per loop
In [26]: %timeit test_dos(10000)
1000 loops, best of 3: 795 µs per loop
In [27]: %timeit test_dos(10000)
1000 loops, best of 3: 797 µs per loop
他们是一样的。数组内在函数和数组表达式算术中没有隐藏的性能魔法。在这种情况下，它们只是被转换为引擎盖下的循环。
如果你检查生成的 GIMPLE 代码，SAS 和 DOS 都被翻译成相同的优化代码主块，这里没有调用 SUM() 的神奇版本：
<bb 8>:
# val.8_59 = PHI <val.8_49(9), 0.0(7)>
# ivtmp.18_123 = PHI <ivtmp.18_122(9), 0(7)>
# ivtmp.25_121 = PHI <ivtmp.25_120(9), ivtmp.25_117(7)>
# ivtmp.28_116 = PHI <ivtmp.28_115(9), ivtmp.28_112(7)>
_111 = (void *) ivtmp.25_121;
_32 = MEM[base: _111, index: _106, step: 8, offset: 0B];
_36 = MEM[base: _111, index: _99, step: 8, offset: 0B];
_37 = _36 + _32;
_40 = MEM[base: _111, offset: 0B];
_41 = _36 - 1.0e+0;
_42 = __builtin_pow (qsdcs_18, _41);
_97 = (void *) ivtmp.28_116;
_47 = MEM[base: _97, offset: 0B];
_43 = _40 * _47;
_44 = _43 * _42;
_48 = _44 * _37;
val.8_49 = val.8_59 + _48;
ivtmp.18_122 = ivtmp.18_123 + 1;
ivtmp.25_120 = ivtmp.25_121 + _118;
ivtmp.28_115 = ivtmp.28_116 + _113;
if (ivtmp.18_122 == _96)
goto <bb 10>;
else
goto <bb 9>;
<bb 9>:
goto <bb 8>;
<bb 10>:
# val.8_13 = PHI <val.8_49(8), 0.0(6)>
_51 = val.8_13 + _17;
*denom_52(D) = _51;
代码在功能上与 do 循环版本相同，只是变量的名称不同。
2. 他们假设形状数组参数 (:) 有可能降低性能。尽管在假定大小参数 (*) 或显式大小 (n) 中接收到的参数总是简单连续的，但理论上不必是假定形状，编译器必须为此做好准备。因此，我总是建议将 contiguous 属性用于您假定的形状参数，只要您知道您将使用连续数组调用它。
在另一个答案中，更改毫无意义，因为它没有使用假定形状参数的任何优点。也就是说，您不必传递带有数组大小的参数，您可以使用 size() 和 shape().
综合比较结果如下。我让它尽可能公平。 Fortran文件用-Ofast编译如上图：
import numpy as np
import dos as dos
import sas as sas
from matplotlib import pyplot as plt
import timeit
import numexpr as ne
#%matplotlib inline
ne.set_num_threads(1)
def test_n(n_comp):
cp = np.tile(1.243,n_comp)
cs = 100.
qs = np.tile(1100.,n_comp)
x= np.random.rand(3*n_comp+1)
def test_dos():
denom = np.empty(1)
dos.get_denomsas(qs,x,cp,cs,n_comp)
def test_sas():
denom = np.empty(1)
sas.get_denomsas(qs,x,cp,cs,n_comp)
def get_denom():
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
denom = np.sum(a) + cs
return denom
def get_denom_numexp():
k = x[0:n_comp]
sigma = x[n_comp:2*n_comp]
z = x[2*n_comp:3*n_comp]
loc_cp = cp
loc_cs = cs
loc_qs = qs
# calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
a = ne.evaluate('(sigma + z)*( k*(loc_qs/loc_cs)**(z-1) )*loc_cp' )
return cs + np.sum(a)
print 'py', timeit.Timer(get_denom).timeit(1000000/n_comp)
print 'dos', timeit.Timer(test_dos).timeit(1000000/n_comp)
print 'sas', timeit.Timer(test_sas).timeit(1000000/n_comp)
print 'ne', timeit.Timer(get_denom_numexp).timeit(1000000/n_comp)
def test():
for n in [10,100,1000,10000,100000,1000000]:
print "-----"
print n
test_n(n)
结果：
py dos sas numexpr
10 11.2188110352 1.8704519272 1.8659651279 28.6881871223
100 1.6688809395 0.6675260067 0.667083025 3.4943861961
1000 0.7014708519 0.5406000614 0.5441288948 0.9069931507
10000 0.5825948715 0.5269498825 0.5309231281 0.6178650856
100000 0.5736029148 0.526198864 0.5304090977 0.5886831284
1000000 0.6355218887 0.5294830799 0.5366530418 0.5983200073
10000000 0.7903120518 0.5301260948 0.5367569923 0.6030929089
您可以看到两个 Fortran 版本之间的差异非常小。数组语法稍微慢一点，但真的没什么好说的。
结论 1：在这个比较中，所有的开销应该是公平的，你会看到对于理想的向量长度 Numpy 和 Numexpr CAN 几乎达到 Fortran 的性能，但是当向量太小甚至可能太大时，Python 解决方案的开销占上风。
结论2：另一个比较中SAS版本较高是与原OP版本不等价造成的。我的答案上面包含等效的优化 DO 循环版本。
怀疑它。
好的，最后，我们获准在我们的一台机器上安装 Numpy 等，这可能是对您的原始内容的全面解释 post。
简短的回答是，从某种意义上说，您最初的问题是 "the wrong question"。此外，其中一位贡献者有很多无理取闹的辱骂和错误信息，这些错误和捏造值得引起注意，以免有人误以为他们，并为此付出代价。
此外，我决定将此答案作为一个单独的答案提交，而不是编辑我 4 月 14 日的答案，原因见下文，并且是正当的。
A 部分：OP 的答案
首先，处理原文中的问题 post：您可能还记得我只能对 Fortran 方面发表评论，因为我们的政策对可以安装什么软件以及我们的位置有严格的规定机器，我们手头没有 Python 等（直到刚才）。我还反复声明，就我们所说的曲线特征或 "concavity".
而言，你的结果的特征很有趣此外，纯粹使用 "relative" 结果（因为你没有 post 绝对时间，而且我当时手头没有 Numpy），我已经指出了几次一些重要的信息可能潜伏在其中。
正是如此
首先，我想确保我可以重现你的结果，因为我们通常不使用 Python/F2py，你的结果中隐含的编译器设置等并不明显，所以我进行了各种的测试以确保我们正在谈论苹果对苹果（我 4 月 14 日的回答表明 Debug 与 Release/O2 有很大的不同）。
图 1 显示了我的 Python 三种情况的结果：你原来的 Python/Numpy 内部子程序（称之为 P，我只是 cut/pasted 你的原始程序），你的基于 Fortran 的原始 Do s/r 你已经 posted（称之为 FDo，我只是 copy/pasted 你的原始，和以前一样），以及我之前建议的依赖数组部分的变体之一，因此需要 Sum()（调用此 FSAS，通过编辑您的原始 FDo 创建）。图 1 显示了通过 timeit 的绝对时序。
图 2 显示了根据您的相对策略除以 Python/Numpy (P) 时序的相对结果。仅显示了两个（相对的）Fortran 变体。
很明显，那些重现了您原始情节中的角色，我们可以确信我的测试似乎与您的测试一致。
现在，您原来的问题是 "Why is is the speed of Fortran getting closer to Numpy with increasing size of the arrays?"。
其实不然。这纯粹是 artefact/distortion 纯粹依靠 "relative" 时间安排可能给人的印象。
查看图 1，从绝对时间来看，很明显 Numpy 和 Fortran 存在分歧。因此，事实上，如果您愿意，Fortran 结果正在远离 Numpy，反之亦然。
一个更好的问题，也是我之前评论中反复出现的一个问题，为什么这些首先是向上弯曲的（例如 c.f. 线性）？我之前的 Fortran-only 结果显示 "mostly" 持平的相对性能比（我的 4 月 14 日 chart/answer 中的黄线），所以我推测 Python 方面发生了一些有趣的事情并建议检查一下。
显示这一点的一种方法是使用另一种相对度量。我将每个（绝对）系列与其自身的最高值（即 n_comp = 10k）分开，以查看此 "internal relative" 性能如何展开（这些被称为 ??10k 值，代表时间n_comp = 10,000）。图 3 将 P、FDo 和 FSAS 的这些结果显示为 P/P10k、FDo/FDo10k 和 FSAS/FSAS10k。为清楚起见，y 轴采用对数刻度。很明显，随着 n_comp c.f 的降低，Fortran 变体的表现相对更好。 P 结果（例如，红色圆圈注释部分）。
换句话说，Fortran 在减小数组大小方面非常非常（非线性）更有效。不确定为什么 Python 会随着 n_comp 的减少而变得更糟......但它确实存在，并且可能是内部 overhead/set-up 等的问题，以及口译员与口译员之间的差异编译器等
因此，Fortran 并不是 "catching up" 与 Python，恰恰相反，它正在继续与 Python 保持距离（见图 1）。然而，Python 和 Fortran 之间的非线性差异随着 n_comp 的减少产生了 "relative" 具有明显反直觉和非线性特征的时序。
因此，随着 n_comp 的增加和每个方法 "stabilises" 或多或少的线性模式，曲线变平表明它们的差异呈线性增长，并且相对比率趋于近似常量（忽略内存争用、smp 问题等）...如果允许 n_comp > 10k，这更容易看出，但我 4 月 14 日的回答中的黄线已经显示了 Fortran-only s/r的。
旁白：我的偏好是创建自己的时间 routines/functions。 timeit 看起来很方便，但是"black box"里面有很多事情要做。设置您自己的循环和结构，并确定您的计时函数的属性s/r解决方案对于正确评估很重要。