[Pyrex] Newbie questions: Numeric vs Numpy

Daniel daniel at rozengardt.net
Mon Oct 1 18:11:16 CEST 2007


Hi folks,
I’m newby in Pyrex. I’m trying to optimize a python application. I’m
working on Windows XP enviroment, Python 2.5 and the last binary pirex
version 0.9.3.1
First I tried, following an example from pirex´s guide using Numeric,
and it works exactly like my old python version. The problem is I used
numpy.ndarray so I have to transform into Numeric.array, use the rutin
and back to numpy. It’s not good.
Looking for a solution, I followed an example in pirex’s code at
http://www.scipy.org/PerformancePython
It works with ndarray but the results are not the same. I can not find
where the error is.
 
First my Numeric code: (c1, c2, c3 are arrays of dimension 1, a is a
matrix (3,7))
 
cdef extern from "Numeric/arrayobject.h":
 
  struct PyArray_Descr:
    int type_num, elsize
    char type
 
  ctypedef class Numeric.ArrayType [object PyArrayObject]:
    cdef char *data
    cdef int nd
    cdef int *dimensions, *strides
    cdef object base
    cdef PyArray_Descr *descr
    cdef int flags
 
def calculo1(ArrayType c2, ArrayType c3, ArrayType a, ArrayType c1):
    cdef int k, j , ncola
    cdef double pivote 
    cdef double *elema, *elemc2, *elemc3, *elemc1 
    elema = <double *>a.data
    elemc2 = <double *>c2.data
    elemc3 = <double *>c3.data
    elemc1 = <double *>c1.data
    ncola = a.dimensions[1]
    for k from 0 <= k < 2:
        pivote = elema[k * ncola] / elemc2[0]
        for j from 0 <= j < 2:
            elema[k * ncola + j + 1] = elema[k * ncola + j + 1] - pivote
* elemc2[j+1]
        elema[k * ncola + 6] = elema[k * ncola + 6] - pivote * elemc2[3]
    for k from 0 <= k < 3:
        pivote = elema[k * ncola + 1] / elemc3[0]
        elema[k * ncola + 2] = elema[k * ncola + 2] - pivote * elemc3[1]
        elema[k * ncola + 6] = elema[k * ncola + 6] - pivote * elemc3[2]
    for k from 2 <= k < 7:
        elemc1[k-2] = elema[k]
    return a, c1
 
Now my Numpy approach:
 
cdef extern from "Python.h":
  int PyObject_AsReadBuffer(object, void **rbuf, int *len)
  int PyObject_AsWriteBuffer(object, void **rbuf, int *len)
 
def calculo1(object c2, object c3, object a, object c1):
    cdef int k, j , ncola
    cdef double pivote 
    cdef double *elema, *elemc2, *elemc3, *elemc1 
    elema = <double *>a.data
    #elema = a.ravel()
    elemc2 = <double *>c2.data
    elemc3 = <double *>c3.data
    elemc1 = <double *>c1.data
    ncola = a.shape[1]
    for k from 0 <= k < 2:
        pivote = elema[k * ncola] / elemc2[0]
        for j from 0 <= j < 2:
            elema[k * ncola + j + 1] = elema[k * ncola + j + 1] - pivote
* elemc2[j+1]
        elema[k * ncola + 6] = elema[k * ncola + 6] - pivote * elemc2[3]
    for k from 0 <= k < 3:
        pivote = elema[k * ncola + 1] / elemc3[0]
        elema[k * ncola + 2] = elema[k * ncola + 2] - pivote * elemc3[1]
        elema[k * ncola + 6] = elema[k * ncola + 6] - pivote * elemc3[2]
    for k from 2 <= k < 7:
        elemc1[k-2] = elema[k]
    return a, c1
 
Thanks!
Daniel Rozengardt
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.copyleft.no/pipermail/pyrex/attachments/20071001/d926497d/attachment.html 


More information about the Pyrex mailing list