Python C Api: transfer dynamically allocated array to Python

Question

I am using Python and I want to do numerical calculations (Runge-Kutta for a scalar ODE, dy/dt = f(t,y) ) on huge arrays (size can be up to 8000 x 500 ). To speed it up, I implemented the routine in C. I am using the Python C API to pass the data to C, do the calculation and then send it back. Since such a huge array is involved, I use dynamic memory allocation to store the calculated result (labelled R, see code).

Now I have problems sending this array back to Python. At the moment I am using Py_BuildValue (see code). Building the .pyd extension works, but when I execute the Python code, I get the message "python.exe stops working".

I would like to ask, if you could help me implementing it the right way:) [I omitted some variables in the code, to make it a little shorter]

# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <Python.h>
# include <numpy\arrayobject.h>

static PyObject* integrateODE(PyObject* self, PyObject* args);
double func ( double t, double u , double a );

int Nt, Nr, m;
double dt, rho0, eta1, eta2, rho2decay;

static PyObject* integrateODE(PyObject* self, PyObject* args)
  {  
  int i,k, len;
  double *mag;
  double u, u0, u1,u2,u3 ;  
  double a0, a1;  
  double f0, f1, f2, f3;
  PyObject *U;
  double *R;

  PyArg_ParseTuple(args, "Offifffii", &U, &dt, &rho0, &m, &eta1, &eta2,&rho2decay, &Nt, &Nr);  
  R = malloc(sizeof(double) * Nt*Nr);  
  len = PySequence_Length(U); 
  mag= malloc(sizeof(double) * Nt*Nr);

  while (len--) {  
  *(mag + len) = (double) PyFloat_AsDouble(PySequence_GetItem(U, len)); 
  }

  for (k=0; k<Nr; k++ )
    {
     u0 = rho0;
     for (i=0; i<Nt; i++ )
         {
          a0=mag[k*Nt + i];
          a1=mag[k*Nt + i+1];
          f0 = func( t0, u0, a0 );
          u1 = u0 + dt * f0 / 2.0;            
          f1 = func( t1, u1, a1 );
          u2 = u0 + dt * f1 / 2.0;  
          f2 = func( t2, u2, a1 );    
          u3 = u0 + dt * f2;
          f3 = func( t3, u3, a0 );
          u = u0 + dt * ( f0 + 2.0 * f1 + 2.0 * f2 + f3 ) / 6.0;    
          R[k*Nt + i] = u;
          u0 = u;       
          }           
     }

     return Py_BuildValue("O", R);
}

double func ( double t, double u, double a ){
  dydt = -u*u  +  a * u + a;  
  return dydt;
}
static PyMethodDef ODE_Methods[] =
{
 {"integrateODE", integrateODE, METH_VARARGS, NULL}, 
 {NULL, NULL, 0, NULL} 
};

PyMODINIT_FUNC initodeNRP(void) {
(void) Py_InitModule("odeNRP", ODE_Methods);
}

Thank you very much for any help you can give me:)

Hi David, thanks for the link. Interesting idea to modify the numpy array directly, without passing arrays at all. I will try this :) — Impaulator, Oct 19 '16 at 07:17
You can also create the numpy array yourself https://docs.scipy.org/doc/numpy/user/c-info.how-to-extend.html#creating-a-brand-new-ndarray. Either way numpy is the way to go — DavidW, Oct 19 '16 at 07:49

Anatoly Strashkevich · Answer 1 · 2016-10-19T07:21:38.773

0

In your allocation, :

mag size = Nt*Nr,

In your loop, when i is max and k is max:

k=Nr-1, i=Nt-1,

Hence:

a0=mag[k*Nt + i]; -> (Nt-1)+Nt*(Nr-1)=Nt-1+NtNr-Nt=NtNr-1, and you call

a1=mag[k*Nt + i+1];, (Nt-1)+Nt*(Nr-1)+1=Nt-1+NtNr-Nt+1=NtNr which is NtNr -> out of bounds

Also you may try to use PyList_New(length); or PyArray_SimpleNewFromData(nd, dims, NPY_DOUBLE, a); to return your array.

edited Oct 19 '16 at 07:21

answered Oct 19 '16 at 07:09

Anatoly Strashkevich

1,834
4
16
32

Hi Anatoly, thank you very much - Davids and your method both work. I feel stupid for asking such easy questions, but sometimes a little guidance does wonders:) – Impaulator Oct 19 '16 at 12:35
No problem! Good luck! – Anatoly Strashkevich Oct 19 '16 at 12:39

Python C Api: transfer dynamically allocated array to Python

1 Answers1