我開始研究一個到Fortran庫（BLAS/LAPACK）的SciPy接口，如下所示：Calling BLAS/LAPACK directly using the SciPy interface and Cython並提出了一個解決方案，但不得不求助於使用numpy.zeros，這實際上是直接通過調用Fortran代碼來殺死任何速度增益。問題是Fortran代碼需要一個0值輸出矩陣（它在內存中的矩陣上運行）以匹配Numpy版本（np.outer）。在Cython腳本中使用memset代替np.zeros以提高速度

所以我有點困惑，因爲Python中的1000x1000零點矩陣只需要8us（使用％timeit或0.008ms），所以爲什麼添加Cython代碼會殺死運行時，注意到我也在memoryview上創建它？（基本上它在1000×1000矩陣乘法上從3ms到8ms左右）。然後，搜索後，我發現一個解決方案使用memset作爲最快的數組更改機制。所以，我從被引用後爲較新的memoryview格式重寫整個代碼，並得到了這樣的事情（用Cython cyblas.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，但顯然不是。

有沒有人有一個想法如何重置整個輸出矩陣爲0使用memset？非常感激。

來源

2017-07-01 Matt

不......任何大小的'memset' _will_都可以工作。你是否100％肯定M和N是你想要的？你打印出來並檢查了嗎？ –

@Coldspeed我剛剛添加了打印輸出，它顯示M和N都是1000，M * N是1000000，正如預期的那樣。（1000,1000,1000000），因爲它等於輸出數組中的＃ – Matt

好的，還有一個問題。無論如何，你將一個0的數組傳遞給該函數，那麼如何知道它只能工作到第125行？ –

[更新 - 修復是：它需要#bytes而不僅僅是M*N的數組大小，即M*N*variable_bytes作爲長度輸入。通過前面的結果可以看出邏輯：第125行是停止的地方* 8字節= 1000行，因此回答了問題。仍然指標不是很好：100循環，最好的3：每循環5.41毫秒（循環）100循環，最好的3：每循環3.95毫秒（numpy），但仍然解決。對上述代碼的更改是：cdef variable_bytes = np.dtype(REAL).itemsize #added to get bytes for memset, after REAL is defined, in this case 8 bytes 然後當你調用memset：memset(&_output[0,0], 0, M*N*variable_bytes) # gives the same output as np.zeros function現在唯一加快速度的地方是我可以看到在大型矩陣上使用了OpenMP語句，但還沒有被看到。

來源

2017-07-02 00:54:41 Matt

在Cython腳本中使用memset代替np.zeros以提高速度

測試腳本

回答

相關問題