實現2D拉普拉斯在用Cython週期性邊界counditions

Question

我有實現二維拉普拉斯對於偏微分方程finite differences集成方法的代碼，使用roll method of Numpy：

def lapOp(u): 
    """ 
    This is the laplacian operator on 2D array 
    of stencil of 4th accuracy terms 
    """ 
    lap = ((4.0/3.0)*np.roll(u,1,axis=0) + (4.0/3.0)*np.roll(u,-1,axis=0) + (4.0/3.0)*np.roll(u,1,axis=1) + (4.0/3.0)*np.roll(u,-1,axis=1) -5.0*u) 
    lap -= ((1.0/12.0)*np.roll(u,2,axis=0) + (1.0/12.0)*np.roll(u,-2,axis=0) + (1.0/12.0)*np.roll(u,2,axis=1) + (1.0/12.0)*np.roll(u,-2,axis=1)) 
    lap = lap/hh 
    return lap 

我想cythonize我的代碼 - 將滾動方法在pyx代碼中工作，還是應該使用C實現滾動方法？

Answer 1

簡短的回答是：roll將在Cython中工作，但它不會更快（任何？）。

如果你想加快你應該使用像roll東西完全可能避免（這是緩慢的，因爲它創建每次被稱爲時間的完整副本），而使用索引來獲得享有numpy的陣列u的大塊。你不應該需要Cython，並且可能不會從中受益。

一個不完整的例子如下（示出足以給主旨）：

def lapOp(u): 
    lap = np.empty_like(u) 
    # this bit is equivalent to (4.0/3)*np.roll(u,1,axis=0) 
    lap[1:,:] = (4.0/3.0)*u[:-1,:] 
    lap[0,:] = (4.0/3.0)*u[-1,:] 

    # add (4.0/3)*np.roll(u,-1,axis=0) 
    lap[:-1,:] += (4.0/3.0)*u[1:,:] 
    lap[-1,:] += (4.0/3.0)*u[0,:] 

    # add (4.0/3)*np.roll(u,1,axis=1) 
    lap[:,1:] += (4.0/3.0)*u[:,:-1] 
    lap[:,0] += (4.0/3.0)*u[:,-1] 

    # the remainder is left as a rather tedious exercise for the reader 

    return lap/hh