2012-10-25 69 views
3

我想加快我的代码,目前需要一个多小时才能在Python/Numpy中运行。大部分计算时间出现在粘贴在下面的函数中。向Python化/ Numpy中的循环向量化可能吗?

我想向量化Z,但是我发现它对于三重循环很难。我能否在某处实施numpy.diff功能?请看:

def MyFESolver(KK,D,r,Z): 
    global tdim 
    global xdim 
    global q1 
    global q2 
    for k in range(1,tdim): 
     for i in range(1,xdim-1): 
      for j in range (1,xdim-1): 
       Z[k,i,j]=Z[k-1,i,j]+r*q1*Z[k-1,i,j]*(KK-Z[k-1,i,j])+D*q2*(Z[k-1,i-1,j]-4*Z[k-1,i,j]+Z[k-1,i+1,j]+Z[k-1,i,j-1]+Z[k-1,i,j+1]) 
    return Z 

tdim = 75xdim = 25

+0

尝试行'Z [1 :, 1:-1,1:-1] = Z [: - 1,1:-1,1:-1 ] + R * q1 * Z [:-1,1:-1,1:-1] *(KK-Z [:-1,1:-1,1:-1])+ D * q2 *(Z [:-1,... -2,1:-1] -4 * Z [:-1,1:-1,1:-1] + Z [:-1,2 :, 1:-1] + Z [:-1,1:-1,:-2] + Z [: - 1,1:-1,2:])'而不是你的三元组。 – halex

+1

不相关,但很重要:您可能想重新检查使用'global'关键字。你的情况没用。 – Simon

回答

2

这看起来像用Cython理想的情况下。我建议在Cython中编写这个函数,它可能会快几百倍。

4

我同意,这很棘手,因为所有四边的BCs都破坏了刚度矩阵的简单结构。你可以摆脱空间的循环这样:

from pylab import * 
from scipy.sparse.lil import lil_matrix 
tdim = 3;  xdim = 4; r = 1.0; q1, q2 = .05, .05; KK= 1.0; D = .5 #random values 
Z = ones((tdim, xdim, xdim)) 
#Iterate in time 
for k in range(1,tdim): 
    Z_prev = Z[k-1,:,:] #may need to flatten 
    Z_up = Z_prev[1:-1,2:] 
    Z_down = Z_prev[1:-1,:-2] 

    Z_left = Z_prev[:-2,1:-1] 
    Z_right = Z_prev[2:,1:-1] 

    centre_term = (q1*r*(Z_prev[1:-1,1:-1] + KK) - 4*D*q2)* Z_prev[1:-1,1:-1] 

    Z[k,1:-1,1:-1]= Z_prev[1:-1,1:-1]+ centre_term + q2*(Z_up+Z_left+Z_right+Z_down) 

但我不认为你可以摆脱时间循环的...

我想表达:

Z_up = Z_prev[1:-1,2:] 

在numpy中制作副本,而你想要的是一个视图 - 如果你能想出如何做到这一点 - 它应该更快(多少?)

最后,我同意其余的回答者 - 来自经验e,这种循环最好在C中完成,然后包装成numpy。但上面应该比原来快...

+2

是的,除了时间循环外,其他所有的东西都可以矢量化。 Numpy总是为切片操作创建视图。 – seberg