parallelize (not symmetric) loops in python
<p>The following code is written in python and it works, i.e. returns the expected result. However, it is very slow and I believe that can be optimized.</p>

<pre><code>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
</code></pre>

<p>Do you know how can I parallelize it? You should note that the two loops are not symmetric.</p>

<p>Edit one: A numpythonic solution was presented above and I made a comparison between the c++ implementation, my loop version in python and thr numpythonic. Results are the following:
- c++ = 0.14seg
- numpythonic version = 1.39seg
- python loop version = 46.56seg
Probably results can get better if we use the intel version of numpy.</p>