使用 SciPy 接口和 Cython 直接调用 BLAS/LAPACK

Question

这里有一个 post：https://gist.github.com/JonathanRaiman/f2ce5331750da7b2d4e9 通过调用 Fortran 库（BLAS / LAPACK / Intel MKL / OpenBLAS / 无论你用 NumPy 安装什么），它显示了巨大的速度提升.经过数小时的工作（因为 SciPy 库被弃用），我终于让它编译了，但没有结果。它比 NumPy 快 2 倍。不幸的是，正如另一位用户指出的那样，Fortran 例程总是将输出矩阵添加到计算的新结果中，因此它仅在第一个运行上匹配 NumPy。 IE。 A := alpha*x*y.T + A。因此，这仍然需要通过快速解决方案来解决。

[更新：对于那些希望使用 SCIPY 界面的人，只需转到此处 https://github.com/scipy/scipy/blob/master/scipy/linalg/cython_blas.pyx 因为他们已经在 CPDEF 语句中优化了对 BLAS/LAPACK 的调用，只需复制/ 粘贴到您的 CYTHON 脚本中 # Python-accessible wrappers for testing: 同样在 link 上面 cython_lapack.pyx 可用但没有 Cython 测试脚本]

测试脚本

import numpy as np;
from cyblas import outer_prod;
a=np.random.randint(0,100, 1000);
b=np.random.randint(0,100, 1000);
a=a.astype(np.float64)
b=b.astype(np.float64)
cy_outer=np.zeros((a.shape[0],b.shape[0]));
np_outer=np.zeros((a.shape[0],b.shape[0]));

%timeit outer_prod(a,b,cy_outer)
#%timeit outer_prod(a,b) #use with fixed version instead of above line, results will automatically update cy_outer
%timeit np.outer(a,b, np_outer)
100 loops, best of 3: 2.83 ms per loop
100 loops, best of 3: 6.58 ms per loop

# 结束测试脚本

要编译的 PYX 文件 cyblas.pyx（基本上是一个 np.ndarray 版本）

import cython
import numpy as np
cimport numpy as np

from cpython cimport PyCapsule_GetPointer 
cimport scipy.linalg.cython_blas
cimport scipy.linalg.cython_lapack
import scipy.linalg as LA

REAL = np.float64
ctypedef np.float64_t REAL_t
ctypedef np.uint64_t  INT_t

cdef int ONE = 1
cdef REAL_t ONEF = <REAL_t>1.0

ctypedef void (*dger_ptr) (const int *M, const int *N, const double *alpha, const double *X, const int *incX, double *Y, const int *incY, double *A, const int * LDA) nogil
cdef dger_ptr dger=<dger_ptr>PyCapsule_GetPointer(LA.blas.dger._cpointer, NULL)  # A := alpha*x*y.T + A

cpdef outer_prod(_x, _y, _output):
#cpdef outer_prod(_x, _y): #comment above line & use this to use the reset output matrix to zeros
    cdef REAL_t *x = <REAL_t *>(np.PyArray_DATA(_x))
    cdef int M = _y.shape[0]
    cdef int N = _x.shape[0]
    #cdef np.ndarray[np.float64_t, ndim=2, order='c'] _output = np.zeros((M,N)) #slow fix to uncomment to reset output matrix to zeros
    cdef REAL_t *y = <REAL_t *>(np.PyArray_DATA(_y))
    cdef REAL_t *output = <REAL_t *>(np.PyArray_DATA(_output))
    with nogil:
        dger(&M, &N, &ONEF, y, &ONE, x, &ONE, output, &M)

非常感谢。希望这能为其他人节省一些时间（它几乎可以工作）——实际上正如我评论的那样，它工作 1x 并匹配 NumPy，然后每个后续调用再次添加到结果矩阵。如果我将输出矩阵重置为 0 并重新运行结果匹配 NumPy。奇怪......虽然如果取消注释上面的几行它会工作，但只能以 NumPy 的速度。已找到替代方案 memset 并将在另一个 post 中...我只是还没有弄清楚如何调用它。

Answer 1

根据netlibdger(M, N, ALPHA, X INCX, Y, INCY, A, LDA)执行A := alpha*x*y**T + A。所以 A 应该全为零才能得到 X 和 Y.

的外积

使用 SciPy 接口和 Cython 直接调用 BLAS/LAPACK

Calling BLAS / LAPACK directly using the SciPy interface and Cython

python

numpy

cython

blas

intel-mkl

测试脚本

要编译的 PYX 文件 cyblas.pyx（基本上是一个 np.ndarray 版本）