2016-02-11 76 views
0

我正在尝试为cython和numpy获取FLOPS基准。为此我写了一个cython程序。这里是:Cython性能基准

cimport numpy as np 
import numpy as np 
import time 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark(): 

    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.random.rand(3) 
    cdef np.ndarray[np.float64_t, ndim=1] res 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     res = np.dot(m1, m2) 
    t2 = time.time() 
    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

在我的机器上我测量“计算MFLops是:7.42390102416”。我的机器具有Intel Core i7-6700HQ CPU @ 2.6 GHz并运行Windows 10.

如果要在机器上运行它,请将代码保存在名为“benchmark.pyx”的文件中。然后创建一个名为“setup.py”与文件,内容如下:

from distutils.core import setup                 
from Cython.Build import cythonize                
import numpy                      

setup(                       
    ext_modules = cythonize("benchmark.pyx"), 
    include_dirs=[numpy.get_include()]               
) 

那么你应该能够“蟒蛇setup.py build_ext --inplace”进行编译。在Windows中,遇到可怕的“无法找到vcvarsall.bat”错误并且不得不花费大量精力来解决这个问题时,它可能会更加困难。

这个表现对我来说似乎很差。我想知道是否有人可以在他们的平台上运行它并告诉我你得到了什么?或者指出我的代码中对性能有不利影响的任何明显错误?

谢谢!

+1

这似乎是它会更适合codereview.stackexchange.com –

回答

1

Cython实际上并没有消除np.dot上的任何Python调用开销。这包括(请注意,名单并不详尽,并可能地方稍有不妥,但它给人的要点):

  1. 寻找np.dot打电话:

    • 在字典查找全局名字空间为np
    • np的字典查询为dot的名字空间。 (请注意,以上所有内容都可以通过在功能中执行dot = np.dot来消除,然后调用dot
    • 关于dot__call__的字典查询。 (这可通过更快的机制来完成,如果点是一个C/Fortran的编译函数)
  2. 包装制备的论据np.dot

    • 创建包含传递给np.dot
    • 两个参数的元组
    • 增加每个参数的引用计数。
  3. np.dot接着处理所述参数...

    • 解包元组的每个参数
    • 检查在numpy的阵列。
    • 检查每个numpy数组的dtype是否相同,并基于dtype选择哪个BLAS函数来调用。
    • 检查数组尺寸并确保它们匹配。
  4. ...的输出参数分配空间...

    • 分配新对象np.ndarray
    • 递增该
    • 分配空间的引用计数为ndarray内的物理阵列
  5. ...调用BLAS操作,给出您在浮点运算...

  6. ...而递减的它传递的输入参数引用计数(检查,看是否有应被释放,但没有人会是)

  7. 您的通话功能然后必须:

    • 检查是否异常,通过np.dot
    • 凸起接收输出阵列(可能还有一些的refcount这里杂耍)
    • 递减引用计数以前的内容res
    • 释放以前的内容res,记住它至少需要2步处理,因为阵列与ndarray持有者分开持有。

如果你想使这个最(可能除了分配)微不足道相比,矩阵向量乘法,那么你需要做显著更大的阵列的测量。您可以使用np.dot中的out可选参数摆脱分配。如果您想使其全部消失,则可以使用scipy Cython BLAS interface直接调用BLAS功能。

+0

谢谢大卫!很高兴知道。 – jungleman007

1

在仔细阅读DavidW的文章并做了一些实验之后,我找到了避免所有numpy开销的方法。它涉及使用指针,特别是不会将numpy数组传递给循环内的函数。

这里是全码:

cimport numpy as np 
import numpy as np 
import time 


cdef matrixdotvector(double* mat, int numrows, int numcols, double* vec, double* outputvec): 
    outputvec[0] = mat[0+0*numcols] * vec[0] + mat[1+0*numcols] * vec[1] + mat[2+0*numcols] * vec[2] 
    outputvec[1] = mat[0+1*numcols] * vec[0] + mat[1+1*numcols] * vec[1] + mat[2+1*numcols] * vec[2] 
    outputvec[2] = mat[0+2*numcols] * vec[0] + mat[1+2*numcols] * vec[1] + mat[2+2*numcols] * vec[2] 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark(): 

    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.transpose(np.random.rand(3)) 
    cdef np.ndarray[np.float64_t, ndim=1] res 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     res = np.dot(m1, m2) 
    t2 = time.time() 
    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark2(): 

    cdef int numrows = 3 
    cdef int numcols = 3 
    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.transpose(np.random.rand(3)) 
    cdef np.ndarray[np.float64_t, ndim=1] res = np.zeros(3) 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     matrixdotvector(&m1[0,0], numrows, numcols, &m2[0], &res[0]) 
    t2 = time.time() 

    assert (np.linalg.norm(np.dot(m1,m2) - res) < 1.0e-6), "Arrays do not match" 

    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

numpybenchmark之间的大的差异()和numpybenchmark2()是,我通过使指针指向numpy的数据阵列中numpybenchmark2避免所有numpy的开销()(而不是传递键入的numpy对象,这也很慢)。我通过展开并在代码中重新实现它来避免np.dot计算的开销。

所以,现在我得到的基准测试结果是:

在[13]:benchmark.numpybenchmark() 已计算MFLOPS是:7.3412268815

在[14]:benchmark.numpybenchmark2() 计算MFLops是:1521.81908107

所以这是一个相当大的增加。老实说,这不是一种“pythonic”的方式,但它是快速尖叫,所以在某些情况下可能有用。由于matrixdotvector()中的代码看起来像C代码,所以我们可以说这应该全部用C编写。就我个人而言,我发现使用类C的Cython代码而不是直接跳入C的原型会更快。

无论如何,对于正在学习cython的人来说,这篇文章可能有用。