我有以下功能：
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？使用指针？
4 回答
抱了它 .
好的，最后，我们被允许在我们的一个盒子上安装Numpy等，这允许对原始帖子进行全面的解释 .
简短的回答是，从某种意义上说，你的原始问题是“错误的问题” . 此外，其中一个贡献者有很多无理取闹的滥用和错误信息，这些错误和捏造值得关注，以免任何人犯错误相信他们和他们的成本 .
另外，我已经决定将这个答案作为一个单独的答案提交，而不是编辑我4月14日的答案，原因如下所示，并且是适当的 .
A部分：OP的答案
首先，处理原始帖子中的问题：您可能还记得我只能向Fortran方面发表评论，因为我们的政策严格规定可能安装的软件和我们机器的位置，而且我们没有Python等交出（直到现在） . 我也一再声明，你的结果的特征在我们可以称之为弯曲的角色或者可能是“凹陷”方面是有趣的 .
另外，纯粹用“相对”结果工作（因为你没有公布绝对时间，而我当时没有Numpy提交），我曾多次指出一些重要信息可能潜伏在其中 .
情况确实如此 .
首先，我想确保我可以重现你的结果，因为我们不正常使用Python / F2py，结果中隐含的是什么编译器设置并不明显，所以我进行了各种测试以确保我们是谈论苹果对苹果（我的4月14日的回答表明Debug vs Release / O2产生了很大的不同） .
图1显示了我的Python结果，仅用于以下三种情况：您的原始Python / Numpy内部子程序（称为P，我只是剪切/粘贴您的原始文件），您原来基于Do的Fortran s / r您已发布（称为此FDo，我只是像以前一样复制/粘贴你的原文），以及我之前建议的其中一个依赖于数组截面的变体，因此需要Sum（）（调用这个FSAS，通过编辑你原来的FDo创建） . 图1显示了通过timeit的绝对时序 .
图2显示了基于您按Python / Numpy（P）时序划分的相对策略的相对结果 . 仅显示了两个（相对）Fortran变体 .
显然，那些在原始情节中重现了这个角色，我们可能确信我的测试似乎与您的测试一致 .
现在，您最初的问题是“为什么随着阵列数量的增加，Fortran的速度越来越接近Numpy？” .
事实上，事实并非如此 . 纯粹是纯粹依赖于“相对”时间的人工制品/歪曲可能会给人留下这样的印象 .
从图1可以看出，绝对时间很明显，Numpy和Fortran正在发生分歧 . 所以，事实上，如果你愿意的话，Fortran的结果正在从Numpy转移，反之亦然 .
一个更好的问题，以及我在之前的评论中反复出现的一个问题，是为什么这些问题首先是向上弯曲的（例如，线性的）？我以前的Fortran结果显示了“大部分”平坦的相对性能比（我在4月14日的图表/答案中的黄线），所以我推测Python方面有一些有趣的事情，并建议检查一下 .
证明这一点的一种方法是采用不同的相对衡量标准 . 我将每个（绝对）系列除以其自己的最高值（即在n_comp = 10k），以查看这个"internal relative"性能如何展开（这些被称为?? 10k值，表示n_comp = 10,000的时序） . 图3显示P，FDo和FSAS的这些结果为P / P10k，FDo / FDo10k和FSAS / FSAS10k . 为清楚起见，y轴具有对数标度 . 很明显，Fortran变体的预制形式相对非常好，随着下降n_comp c.f. P结果（例如红色圆圈注释部分） .
换句话说，Fortran非常（非线性）更有效地减少阵列大小 . 不确定为什么Python会因为减少n_comp而做得更糟......但是它存在，并且可能是内部开销/设置等问题，以及解释器与编译器之间的差异等 .
因此，并不是说Fortran正在“追赶”Python，恰恰相反，它继续与Python保持距离（见图1） . 然而，随着n_comp的减少，Python和Fortran之间的非线性差异产生了“相对”时间，显然是反直觉和非线性的 .
因此，随着n_comp增加并且每种方法“稳定”到或多或少线性模式，曲线变平以显示它们的差异线性增长，并且相对比率稳定到近似常数（忽略存储器争用，smp问题等 . ）...如果允许n_comp大于10k，则更容易看到，但我4月14日答案中的黄线已经显示了仅适用于Fortran的s / r .
旁白：我的偏好是创建自己的计时例程/函数 . timeit似乎很方便，但内部有很多“黑匣子” . 设置自己的循环和结构，以及确定计时功能的属性/分辨率对于正确评估非常重要 .
在另一个答案中被命名，我必须回应 .
我知道这并没有真正回答原来的问题，但原始海报鼓励在他的评论中追求这个方向 .
我的观点是这些：
1. 我不相信数组内在因素可以更好地进行优化 . 如果幸运的话，它们被转换为与手动循环相同的循环代码 . 如果不是，则可能出现性能问题 . 因此，必须要小心 . 有可能触发临时数组 .
我将提供的SAS数组转换为通常的do循环 . 我称之为DOS . 我证明DO循环决不会慢，在这种情况下，两个子程序都会产生或多或少相同的代码 .
重要的是要说我不相信这个版本的可读性较差，因为它有一两行 . 这实际上也非常简单，看看这里发生了什么 .
现在是这些的时机
他们是一样的 . 数组内在函数和数组表达式算术中没有隐藏的性能魔法 . 在这种情况下，它们只是被转换为引擎盖下的循环 .
如果检查生成的GIMPLE代码，SAS和DOS都会转换为优化代码的相同主块，这里不会调用
SUM()的神奇版本：代码在功能上与do循环版本相同，只是变量的名称不同 .
2. 他们假设形状数组参数
(:)有可能降低性能 . 在假定大小参数(*)或显式大小(n)中接收的参数始终是简单连续的，理论上假设的形状不必是必须的，编译器必须为此做好准备 . 因此，我始终建议将contiguous属性用于假定的形状参数，只要您知道将使用连续数组调用它 .在另一个答案中，这种变化毫无意义，因为它没有使用假设形状参数的任何优点 . 也就是说，您不必传递具有数组大小的参数，并且可以使用内部函数，例如
size()和shape().以下是综合比较的结果 . 我让它尽可能公平 . Fortran文件使用
-Ofast编译，如上所示：结果：
您可以看到两个Fortran版本之间存在非常小的差异 . 数组语法稍慢，但真的没什么可说的 .
Conclusion 1 ：在这种比较中，所有人的开销应该是合理的，你会看到理想的向量长度Numpy和Numexpr CAN几乎达到Fortran的性能，但是当向量太小或者甚至太大时，Python解决方案的开销就会占上风 .
Conclusion 2 ：在另一个比较中，性能更高的SAS版本是通过与原始OP的版本进行比较而产生的，该版本不相同 . 我的答案中包含了等效的优化DO循环版本 .
继我之前的回答，弗拉基米尔的弱点推测，我设置了两个s / r：一个作为原始给定，一个使用数组部分和Sum（） . 我还希望证明弗拉基米尔关于Do循环优化的评论很弱 .
此外，我通常用于基准测试的点，上面示例中的n_comp的大小太小 . 下面的结果将每个“原始”和“更好”的SumArraySection（SAS）变量放入在定时调用中重复1,000次的循环中，因此结果是每个s / r的1000个计算结果 . 如果你的时间是几分之一秒，它们可能不可靠 .
有许多变化值得考虑，没有明确的指针 . 用于此插图的一种变体是
关键的区别是：
a）传递的数组是（:) b）s / r中没有数组赋值，而是使用数组部分（相当于“有效”指针） . c）使用Sum（）代替Do
然后还尝试两种不同的编译器优化来演示含义 .
正如两张图所示，原始代码（蓝色钻石）的速度要慢得多 . SAS（红色方块）具有低优化 . SAS在高优化方面仍然更好，但它们越来越接近 . 当使用低编译器优化时，Sum（）被“更好地优化”，这部分地解释了这一点 .
黄线表示两个s / r时序之间的比率 . 忽略顶部图像中“0”处的黄线值（n_comp太小导致其中一个时间变得不稳定）
由于我没有用户的原始数据与Numpy的比率，我只能声明他的图表上的SAS曲线应该低于他当前的Fortran结果，并且可能更平坦甚至是趋势向下 .
换句话说，实际上可能并不存在原始发布中看到的分歧，或者至少不存在这种差异 .
......虽然更多的实验可能有助于证明已经提供的其他评论/答案 .
亲爱的莫里茨：哎呀，我忘了提及，关于指针的问题 . 如前所述，SAS变化改进的一个关键原因是它更好地利用了“有效指针”，因为它不需要将数组x（）重新分配到三个新的本地数组中（即，因为x通过ref传递，使用数组部分是一种内置于Fortran中的指针方法，因此不需要显式指针），但是需要Sum（）或Dot_Product（）等等 .
相反，您可以通过将x更改为n_compx3 2D数组来保持Do并实现类似的功能，或者直接传递n_comp的三个显式1D数组 . 这个决定很可能是由代码的大小和复杂性驱动的，因为它需要重写调用/ sr语句等，以及其他任何地方使用x（） . 我们的一些项目是> 300,000行代码，因此在这种情况下，在本地更改代码（例如SAS等）要便宜得多 .
我还在等待获得在我们的一个盒子上安装Numpy的许可 . 如前所述，为什么你的相对时间暗示Numpy随着n_comp的增加而改善，这是有趣的 .
当然，关于“正确”基准测试等的评论，以及使用fpy暗示的编译器切换的问题仍然适用，因为这些可能会极大地改变您的结果的特征 .
如果他们针对这些排列进行了更新，我会很有兴趣看到您的结果 .
注释中没有足够的信息，但以下某些内容可能会有所帮助：
1）Fortran优化了内部函数，例如“Sum（）”和“Dot_Product”，您可能希望使用它们代替Do循环进行求和等 .
在某些情况下（这里不是neccessiraly），使用ForAll或其他任何东西创建“meta”数组进行求和可能会“更好”，然后在“meta”数组上应用求和 .
但是，Fortran允许数组部分，因此您不需要创建自动/中间数组sigma，k和z，以及相应的开销 . 相反可以有类似的东西
2）有时（取决于编译器，机器等），如果数组大小不是特定的“二进制间隔”（例如1024对1000）等，则可能存在“存储器争用” .
您可能希望在图表中的更多点（即各种其他“n_comps”）重复您的实验，特别是在这样的“边界”附近 .
3）无法判断您是否正在使用完整的编译器optimsation（flags）来编译fortran代码 . 您可能希望查找各种“-o”标志等 .
4）您可能希望包括OpemMP指令（或者至少在你的标志中包含openmp等） . 这有时可以改善某些开销问题，即使不依赖于循环中的OpenMP指令等 .
5）概述：这可能适用于使用循环的每个方法
a）“求和”公式中的“常数运算”可以在循环外执行，例如 . 创建类似qsDcs = qs / cs的东西，并在循环中使用qsDcs .
b）同样，有时创建类似zM1（:) = z（:) - 1的东西很有用，而在循环中使用zM1（:) .