Max K Max K - 9 months ago 33
Python Question

Code runs much faster in C than in NumPy

I wrote physics simulation code in python using numpy and than rewrote it to C++. in C++ it takes only 0.5 seconds while in python around 40s. Can someone please help my find what I did horribly wrong?

import numpy as np

def myFunc(i):
uH = np.copy(u)
for j in range(1, xmax-1):
u[i][j] = a*uH[i][j-1]+(1-2*a)*uH[i][j]+a*uH[i][j+1]
u[i][0] = u[i][0]/b
for x in range(1, xmax):
u[i][x] = (u[i][x]+a*u[i][x-1])/(b+a*c[x-1])
for x in range(xmax-2,-1,-1):
u[i][x]=u[i][x]-c[x]*u[i][x+1]
xmax = 101
tmax = 2000
#All other variables are defined here but I removed that for visibility
uH = np.zeros((xmax,xmax))
u = np.zeros((xmax,xmax))
c = np.full(xmax,-a)
uH[50][50] = 10000
for t in range(1, tmax):
if t % 2 == 0:
for i in range(0,xmax):
myFunc(i)
else:
for i in range(0, xmax):
myFunc(i)


In case someones wants to run it here is whole code: http://pastebin.com/20ZSpBqQ
EDIT: all variables are defined in the whole code which can be found on pastebin. Sorry for confusion, I thought removing all the clutter will make the code easier to understand

Answer Source

So when I ran your numpy python code it took four minutes to run, once I removed the numpy code and replaced it with standard python code it only took one minute! (I have a not so fast computer)

Here's that code:

#import numpy as np

def impl(i,row):
    if row:
        uH = u[:][:] # this copys the array 'u'
        for j in range(1, xmax-1):
            u[i][j] = a*uH[i][j-1]+(1-2*a)*uH[i][j]+a*uH[i][j+1]
        u[i][0] = u[i][0]/b
        for x in range(1, xmax):
            u[i][x] = (u[i][x]+a*u[i][x-1])/(b+a*c[x-1])
        for x in range(xmax-2,-1,-1):
            u[i][x]=u[i][x]-c[x]*u[i][x+1]
    else:
        uH = u[:][:] # this copys the array 'u'
        for j in range(1, xmax-1):
            u[j][i]= a*uH[j-1][i]+(1-2*a)*uH[j][i]+a*uH[j+1][i]
        u[0][i] = u[0][i]/b
        for y in range(1, xmax):
            u[y][i] = (u[y][i]+a*u[y-1][i])/(b+a*c[y-1])
        for y in range(xmax-2,-1,-1):
            u[y][i]=u[y][i]-c[y]*u[y+1][i]

#Init
xmax = 101
tmax = 2000
D = 0.5
l = 1
tSec = 0.1
uH = [[0.0]*xmax]*xmax #np.zeros((xmax,xmax))
u = [[0.0]*xmax]*xmax #np.zeros((xmax,xmax))
dx = l / xmax
dt = tSec / tmax
a = (D*dt)/(dx*dx);
b=1+2*a
print("dx=="+str(dx))
print("dt=="+str(dt))
print(" a=="+str(a))
#koeficient c v trojdiagonalnej matici
c = [-a]*xmax #np.full(xmax,-a)
c[0]=c[0]/b
for i in range(1, xmax):
    c[i]=c[i]/(b+a*c[i-1])
uH[50][50] = 10000
u = uH
for t in range(1, tmax):
    if t % 2 == 0:
        for i in range(0,xmax):
            impl(i,False)
    else:
        for i in range(0, xmax):
            impl(i,True)

I believe that this could be much faster if you were to have used numpy the correct way rather than as a substitute for arrays, however, not using numpy arrays cut the time to 1/4th of the original.