[SciPy-user] Performance Python with Weave article updated

Francesc Alted falted at pytables.org
Wed Sep 22 10:22:32 CDT 2004

Hi again,

I was thinking on why the numarray headers do exist on the code that I've
sent earlier and I've ended with another version. This is simpler, without
any need to importing Numeric or numarray libraries or headers. This one
just recognize the object (Numeric or numarray) and access to its data

Regarding compilation, you don't even need to have Numeric and/or numarray
installed to compile the extension, so the code is very portable to any
platform having Numeric and/or numarray and besides, you don't have to
bother in case of an ABI change on any of these.

As you can see, although I've used PyObject_AsReadBuffer() to get the buffer
data area, nothing forbids me to write in this area as well (in the end,
elem is a C pointer). So you can use this trick to modify data as well.

Maybe this is a hack, but hey! it "just work". And more importantly, it
reveals how powerful Pyrex can be by effectively mixing Python and C calls.
This can easy a lot the transition between Numeric and numarray apps. In
fact, it would be very easy to write apps that support both Numeric and
numarray without any conflict and as I said before without worrying about
library headers or ABI changes.


Francesc Alted

# Pyrex sources to compute the laplacian.
# Some of this code is taken from numeric_demo.pyx
# 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 "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")
        # 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)

More information about the SciPy-user mailing list