Fw: [SciPy-dev] scipy.weave versus simple C++.
eric
eric at scipy.org
Fri Jan 11 11:50:43 CST 2002
Hey Prabhu,
Here is a near verbatim conversion of your laplace to use Python/Inline. I
didn't spend much time, and there are errors. The results don't match.
You'll get the general idea from it though.
I used blitz arrays because they're indexing is kinda like double** indexing
so the conversion was easy. Standard Numeric arrays don't have the pointer
arrays to the rows that C uses, so your code wouldn't translate easily to it
without macros or calculating offsets -- you could do this though. I'm not
sure that indexing blitz arrays is as fast as indexing normal arrays -- this
example would suggest not:
500x500, 100 iterations, PII 850 laptop, W2K, gcc-2.95.2
laplace.py: 3.47 seconds
laplace.cxx: 2.34 seconds
The Python calling overhead should be close to nil, so there is no reason
why the python version should be slower, other than the fact that I used
blitz arrays and indexing. But the current implementation still pays a 50%
penalty for using Python.
eric
---------------------------------------------------------
# laplace.py
import time
try:
from scipy import weave
from scipy import *
except:
import weave
from Numeric import *
from weave.blitz_tools import blitz_type_factories
def BC(x,y):
return x**2-y**2
class Grid:
def __init__(self,shape,typecode=Float64):
self.shape = shape
# should really handle typecode here.
self.dx = 1.0 / shape[0] - 1
self.dy = 1.0 / shape[1] - 1
self.u = zeros(shape,typecode)
def setBCFunc(self,f):
xmin, ymin, xmax, ymax = 0.0, 0.0, 1.0, 1.0
x = arange(self.shape[0])*self.dy
y = arange(self.shape[1])*self.dy
self.u[0 ,:] = f(xmin,y)
self.u[-1,:] = f(xmax,y)
self.u[:, 0] = f(x,ymin)
self.u[:,-1] = f(x,ymax)
class LaplaceSolver:
def __init__(self,grid):
self.g = grid
def timeStep(self, dt = 0.0):
dx2 = self.g.dx**2
dy2 = self.g.dy**2
dnr_inv = .5/(dx2+dy2)
nx, ny = self.g.shape
u = self.g.u
code = """
#line 39 "laplace.py"
double tmp, err, diff;
for (int i=1; i<nx-1; ++i) {
for (int j=1; j<ny-1; ++j) {
tmp = u(i,j);
u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2 +
(u(i,j-1) + u(i,j+1))*dx2)*dnr_inv;
diff = u(i,j) - tmp;
err += diff*diff;
}
}
return_val = Py::new_reference_to(Py::Float(sqrt(err)));
"""
# compiler keyword only needed on windows with MSVC installed
err = weave.inline(code, ['u', 'dx2', 'dy2', 'dnr_inv', 'nx','ny'],
type_factories = blitz_type_factories,
compiler = 'gcc')
return err
def solve(self, n_iter=0, eps=1e-16):
err = self.timeStep()
#print err
count = 1
while err > eps:
if n_iter and (count > n_iter):
return err;
err = self.timeStep()
#print err
count += 1
return count
if __name__ == "__main__":
grid = Grid((500,500))
grid.setBCFunc(BC)
s = LaplaceSolver(grid)
t1 = time.time()
err = s.solve(20)
t2 = time.time()
print "Iterations took ", t2 - t1, " seconds."
print "Error: ", err
More information about the Scipy-dev
mailing list