<p>按照其他建议,在计时器之前使用<code>p=np.asfortranarray(p)</code>确实可以使性能与我测试时的numpy不相上下。我将二元工作台的范围扩展到<code>n_bi = np.array([2**i for i in xrange(1, 15)])</code>,这样p矩阵将大于我的L3缓存大小。在</p>
<p>为了进一步优化这一点,我认为自动编译器选项不会有多大帮助,因为内部循环有一个依赖关系。只有当您手动展开它时,<code>ifort</code>才会向量化最里面的循环。对于<code>gfortran</code>,需要<code>-O3</code>和{<cd6>}。对于受主内存带宽限制的矩阵大小,这将比numpy提高1到3倍的性能优势。在</p>
<p><strong>Update</strong>:在将此方法应用于单变量代码并使用<code>f2py --opt='-O3 -ffast-math' -c -m polynomial polynomial.f90</code>进行编译之后,我得到以下源代码和结果基准.py公司名称:</p>
<pre><code>subroutine polyval(p, x, pval, nx)
implicit none
real*8, dimension(nx), intent(in) :: p
real*8, intent(in) :: x
real*8, intent(out) :: pval
integer, intent(in) :: nx
integer, parameter :: simd = 8
real*8 :: tmp(simd), vecx(simd), xfactor
integer :: i, j, k
! precompute factors
do i = 1, simd
vecx(i)=x**(i-1)
end do
xfactor = x**simd
tmp = 0.0d0
do i = 1, nx, simd
do k = 1, simd
tmp(k) = tmp(k)*xfactor + p(nx-(i+k-1)+1)*vecx(simd-k+1)
end do
end do
pval = sum(tmp)
end subroutine polyval
subroutine polyval2(p, x, y, pval, nx, ny)
implicit none
real*8, dimension(nx, ny), intent(in) :: p
real*8, intent(in) :: x, y
real*8, intent(out) :: pval
integer, intent(in) :: nx, ny
integer, parameter :: simd = 8
real*8 :: tmp(simd), vecx(simd), xfactor
integer :: i, j, k
! precompute factors
do i = 1, simd
vecx(i)=x**(i-1)
end do
xfactor = x**simd
! horner
pval=0.0d0
do i = 1, ny
tmp = 0.0d0
do j = 1, nx, simd
! inner vectorizable loop
do k = 1, simd
tmp(k) = tmp(k)*xfactor + p(nx-(j+k-1)+1,ny-i+1)*vecx(simd-k+1)
end do
end do
pval = pval*y + sum(tmp)
end do
end subroutine polyval2
</code></pre>
<p><strong>更新2</strong>:如前所述,此代码不正确,至少当大小不能被<code>simd</code>整除时。它只是展示了手动帮助编译器的概念,所以不要像这样使用它。如果大小不是2的幂次方,则一个小的余数循环必须处理悬空的索引。这样做并不难,以下是单变量情况下的正确步骤,应直接将其扩展到二元:</p>
^{pr2}$
<p><img src="https://i.stack.imgur.com/QBN6R.png" alt="univariate"/></p>
<p><img src="https://i.stack.imgur.com/rVYkn.png" alt="bivariate"/></p>
<p>另外,对于非常小的尺寸也要小心,因为时间太短,无法获得准确的性能曲线。而且,关于<code>numpy</code>的相对时间可能是欺骗的,因为numpy的绝对时间可能非常糟糕。以下是最大案件的时间安排:</p>
<p>对于nx=2<strong>20的单变量,numpy的时间为1.21s,而自定义fortran版本的时间为1.69e-3s。对于nx<em>ny=2</em></strong>20的二元变量,numpy的时间是8e-3s,定制版本的时间是1.68e-3s。当nxny总大小相同时,单变量和双变量的时间是相同的这一事实非常重要,因为它支持这样一个事实,即代码的执行速度接近内存带宽限制。在</p>
<p><strong>更新3</strong>:使用新的python脚本来实现更小的大小,<code>simd=4</code>我获得了以下性能:</p>
<p><img src="https://i.stack.imgur.com/0Zz2t.png" alt="enter image description here"/></p>
<p><img src="https://i.stack.imgur.com/Xdxkf.png" alt="enter image description here"/></p>
<p><strong>更新4</strong>:对于正确性,结果在双精度精度精度内是相同的,如果您为单变量示例运行以下python代码,可以看到:</p>
<pre><code>import polynomial as P
import numpy.polynomial.polynomial as PP
import numpy as np
for n in xrange(2,100):
poly1n = np.random.rand(n)
poly1f = np.asfortranarray(poly1n)
x = 2
print "%18.14e" % P.polyval(poly1f, x)
print "%18.14e" % PP.polyval(x, poly1n)
print (P.polyval(poly1f, x) - PP.polyval(x, poly1n))/PP.polyval(x,poly1n), '\n'
</code></pre>