[SciPy-user] Performance Python with Weave article updated

[Francesc Alted]

Hi Prabhu,

Very interesting comparison indeed.

By the way, I've managed to create a Pyrex version of your small benchmark
that is a little simpler than yours (it does not need to expose the Numeric
array structure) and that works for both Numeric and numarray. The times
taken by this version are very similar to your version (for both numarray
and Numeric).

I'm attaching the code at the end of this message.

By the way, here you have the figures for my own platform (Pentium IV @
2GHz):

$ python laplace.py
Doing 100 iterations on a 500x500 grid
numeric took 8.49 seconds
blitz took 2.96 seconds
inline took 0.9 seconds
fastinline took 0.6 seconds
fortran took 0.83 seconds
pyrex took 0.6 seconds
slow (1 iteration) took 3.13 seconds
100 iterations should take about 313.000000 seconds
slow with Psyco (1 iteration) took 2.19 seconds
100 iterations should take about 219.000000 seconds

Cheers,

-- 
Francesc Alted


# Pyrex sources to compute the laplacian.
#
# Author: Prabhu Ramachandran <prabhu_r at users dot sf dot net>
# Modified by F. Alted for yet another way to access Numeric/numarray objects

# Some helper routines from the Python API
cdef extern from "Python.h":
  int PyObject_AsReadBuffer(object, void **rbuf, int *len)
  int PyObject_AsWriteBuffer(object, void **rbuf, int *len)
    
cdef extern from "numarray/libnumarray.h":
    # The numarray initialization funtion
    void import_libnumarray()

# The Numeric API requires this function to be called before
# using any Numeric facilities in an extension module.
import_libnumarray()

cdef extern from "math.h":
    double sqrt(double x)

def pyrexTimeStep(object u, double dx, double dy):
    cdef int nx, ny
    cdef double dx2, dy2, dnr_inv, err
    cdef double *elem
    cdef int i, j
    cdef double *uc, *uu, *ud, *ul, *ur
    cdef double diff, tmp
    cdef int buflen
    cdef void *data
    
    if u.typecode() <> "d":
        raise TypeError("Double array required")
    if len(u.shape) <> 2:
        raise ValueError("2 dimensional array required")

    nx = u.shape[0]
    ny = u.shape[1]
    dx2, dy2 = dx**2, dy**2
    dnr_inv = 0.5/(dx2 + dy2)

    # Get the pointer to the buffer data area
    if hasattr(u, "__class__"):
        # numarray case
        if PyObject_AsReadBuffer(u._data, &data, &buflen) <> 0:
            raise RuntimeError("Error getting the array data buffer")
    else:
        # Numeric case
        if PyObject_AsReadBuffer(u, &data, &buflen) <> 0:
            raise RuntimeError("Error getting the array data buffer")
    elem = <double *>data
    
    err = 0.0
    for i from 1 <= i < nx-1:
        uc = elem + i*ny + 1
        ur = elem + i*ny + 2
        ul = elem + i*ny
        uu = elem + (i+1)*ny + 1
        ud = elem + (i-1)*ny + 1
        
        for j from 1 <= j < ny-1:
            tmp = uc[0]
            uc[0] = ((ul[0] + ur[0])*dy2 +
                     (uu[0] + ud[0])*dx2)*dnr_inv
            diff = uc[0] - tmp
            err = err + diff*diff
            uc = uc + 1; ur = ur + 1;  ul = ul + 1
            uu = uu + 1; ud = ud + 1

    return sqrt(err)