我开始研究 Fortran 库 (BLAS/LAPACK) 的 SciPy 接口,可以在这里看到:Calling BLAS / LAPACK directly using the SciPy interface and Cython并想出了一个解决方案,但不得不求助于使用numpy.zeros
它实际上杀死了任何速度直接调用 Fortran 代码的好处。问题是 Fortran 代码需要一个 0 值的输出矩阵(它在内存中就地处理矩阵)以匹配 Numpy 版本(np.outer
)。
所以我有点困惑,因为 Python 中的 1000x1000 零矩阵只需要 8us(使用 %timeit 或 0.008ms)那么为什么添加 Cython 代码会杀死运行时,注意到我也在 memoryview 上创建它?(基本上 1000 x 1000 矩阵乘法从 3ms 到 8ms 左右)。然后在 SO 上搜索后,我在其他地方找到了一个使用memset
最快的数组更改机制的解决方案。因此,我将引用帖子中的整个代码重写为更新的memoryview
格式,并得到了类似这样的内容(Cythoncyblas.pyx
文件):
import cython
import numpy as np
cimport numpy as np
from libc.string cimport memset #faster than np.zeros
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
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cpdef outer_prod(double[::1] _x, double[::1] _y, double[:, ::1] _output):
cdef int M = _y.shape[0]
cdef int N = _x.shape[0]
memset(&_output[0,0], 0, M*N)
with nogil:
dger(&M, &N, &ONEF, &_y[0], &ONE, &_x[0], &ONE, &_output[0,0], &M)
测试脚本
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 np.outer(a,b, np_outer)
因此,这解决了我将输出矩阵值重置为第 125 行的问题,该问题仍然存在(已解决见下文)。我认为将memset
长度参数设置为 M*N 会清除内存中的 1000*1000 ,但显然不是。
有人知道如何使用 将整个输出矩阵重置为 0memset
吗?非常感激。