Fw: [SciPy-dev] scipy.weave versus simple C++.

[eric]

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