python中的并行(不对称)循环

Question

下面的代码是用python编写的，它可以工作，即返回预期的结果。然而，它非常缓慢，我相信这是可以优化的。

G_tensor = numpy.matlib.identity(N_particles*3,dtype=complex)

for i in range(N_particles):
    for j in range(i, N_particles):
        if i != j:

            #Do lots of things, here is shown an example.
            # However you should not be scared because 
            #it only fills the G_tensor
            R = numpy.linalg.norm(numpy.array(positions[i])-numpy.array(positions[j]))
            rx = numpy.array(positions[i][0])-numpy.array(positions[j][0])
            ry = numpy.array(positions[i][1])-numpy.array(positions[j][1])
            rz = numpy.array(positions[i][2])-numpy.array(positions[j][2])
            krq = (k*R)**2
            pf = -k**2*alpha*numpy.exp(1j*k*R)/(4*math.pi*R)
            a = 1.+(1j*k*R-1.)/(krq)
            b = (3.-3.*1j*k*R-krq)/(krq) 
            G_tensor[3*i+0,3*j+0] = pf*(a + b * (rx*rx)/(R**2))  #Gxx
            G_tensor[3*i+1,3*j+1] = pf*(a + b * (ry*ry)/(R**2))  #Gyy
            G_tensor[3*i+2,3*j+2] = pf*(a + b * (rz*rz)/(R**2))  #Gzz
            G_tensor[3*i+0,3*j+1] = pf*(b * (rx*ry)/(R**2))      #Gxy
            G_tensor[3*i+0,3*j+2] = pf*(b * (rx*rz)/(R**2))      #Gxz
            G_tensor[3*i+1,3*j+0] = pf*(b * (ry*rx)/(R**2))      #Gyx
            G_tensor[3*i+1,3*j+2] = pf*(b * (ry*rz)/(R**2))      #Gyz
            G_tensor[3*i+2,3*j+0] = pf*(b * (rz*rx)/(R**2))      #Gzx
            G_tensor[3*i+2,3*j+1] = pf*(b * (rz*ry)/(R**2))      #Gzy

            G_tensor[3*j+0,3*i+0] = pf*(a + b * (rx*rx)/(R**2))  #Gxx
            G_tensor[3*j+1,3*i+1] = pf*(a + b * (ry*ry)/(R**2))  #Gyy
            G_tensor[3*j+2,3*i+2] = pf*(a + b * (rz*rz)/(R**2))  #Gzz
            G_tensor[3*j+0,3*i+1] = pf*(b * (rx*ry)/(R**2))      #Gxy
            G_tensor[3*j+0,3*i+2] = pf*(b * (rx*rz)/(R**2))      #Gxz
            G_tensor[3*j+1,3*i+0] = pf*(b * (ry*rx)/(R**2))      #Gyx
            G_tensor[3*j+1,3*i+2] = pf*(b * (ry*rz)/(R**2))      #Gyz
            G_tensor[3*j+2,3*i+0] = pf*(b * (rz*rx)/(R**2))      #Gzx
            G_tensor[3*j+2,3*i+1] = pf*(b * (rz*ry)/(R**2))      #Gzy

你知道我怎么把它并行化吗？您应该注意到，这两个循环不对称。

编辑一:上面介绍了一个numpythonic解决方案，我比较了c++实现、python循环版本和thr实现。结果如下：- c++ =0.14seg-numpythonic版本=1.39seg-python循环版本= 46.56seg，如果我们使用intel版本的numpy，结果可能会更好。

Answer 1

Python不是一种快速语言。使用python进行数字处理时，应该始终使用用编译语言编写的关键部分代码。将编译降到CPU级别，您可以将代码的速度提高到100倍，然后继续进行并行化。因此，我不会把目光放在使用更多的核心来做低效的事情上，而是更有效率地工作。我看到了以下加速代码的方法：

1)更好地利用numpy:您能否直接在向量/矩阵级别上进行标量级的计算？例如：rx =0,:：，0-位置为，但类似于这些位置。

如果这在你的计算中是不可能的，那么你可以选择2或3。

2)使用cython。Cython将Python代码编译为C，然后将其编译到CPU中。通过在正确的位置使用静态输入，您可以使代码更快，请参见cython教程(例如：http://cython.readthedocs.io/en/latest/src/quickstart/cythonize.html )。

3)如果您熟悉FORTRAN，那么最好用FORTRAN编写这个部分，然后使用f2py从Python调用它。实际上，您的代码看起来很像FORTRAN。对于C和C++来说，SWIG是一个很好的工具，可以让编译的代码在Python中可用，但是还有很多其他技术(cython：：Python，numba等等)。

当您已经这样做了，而且它仍然慢，使用GPU电源与pyCUDA或并行化与mpi4py或多处理可能是一种选择。

Answer 2

下面是一个建议，即现在应该工作(我更正了几个错误)，但是它会给您提供如何将版本化应用到代码中以便有效地使用numpy数组的一般概念。所有东西都是以“一次”(即没有任何for-循环)的方式构建的，这是"numpythonic“方式：

import numpy as np
import math
N=2
k,alpha=1,1
G = np.zeros((N,3,N,3),dtype=complex)

# np.mgrid gives convenient arrays of indices that 
# can be used to write readable code
i,x_i,j,x_j = np.ogrid[0:N,0:3,0:N,0:3]
# A quick demo on how we can make the identity tensor with it
G[np.where((i == j) & (x_i == x_j))] = 1
#print(G.reshape(N*3,N*3))

positions=np.random.rand(N,3)
# Here I assumed position has shape [N,3] 
# I build arr[i,j]=position[i] - position[j] using broadcasting 
# by turning position into a column and a row 
R = np.linalg.norm(positions[None,:,:]-positions[:,None,:],axis=-1)
# R is now a N,N matrix of all the distances
#we reshape R to N,1,N,1 so that it can be broadcated to N,3,N,3 
R=R.reshape(N,1,N,1)

r=positions[None,:,:]-positions[:,None,:]

krq = (k*R)**2
pf = -k**2*alpha*np.exp(1j*k*R)/(4*math.pi*R)
a = 1.+(1j*k*R-1.)/(krq)
b = (3.-3.*1j*k*R-krq)/(krq)

#print(np.isnan(pf[:,0,:,0]))

# here we build all the combination rx*rx rx*ry etc...
comb_r=(r[:,:,:,None]*r[:,:,None,:]).transpose([0,2,1,3])
#we compute G without the pf*A term
G = pf*(b * comb_r/(R**2))
#we add pf*a term where it is due
G[np.where(x_i == x_j)] = (G + pf*a)[np.where(x_i == x_j)]
# we didn't bother with the identity or condition i!=j so we enforce it here
G[np.where(i == j)] = 0
G[np.where((i == j) & (x_i == x_j))] = 1

print(G.reshape(N*3,N*3))