Python C Api: transfer dynamically allocated array to Python
<p>I am using Python and I want to do numerical calculations (Runge-Kutta for a scalar ODE, <em>dy/dt = f(t,y)</em> ) 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).</p>

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

<p>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]</p>

<pre><code># include &lt;stdlib.h&gt;
# include &lt;stdio.h&gt;
# include &lt;math.h&gt;
# include &lt;Python.h&gt;
# include &lt;numpy\arrayobject.h&gt;

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", &amp;U, &amp;dt, &amp;rho0, &amp;m, &amp;eta1, &amp;eta2,&amp;rho2decay, &amp;Nt, &amp;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&lt;Nr; k++ )
    {
     u0 = rho0;
     for (i=0; i&lt;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);
}
</code></pre>

<p>Thank you very much for any help you can give me:)</p>

<p>In your allocation, :</p>

<pre><code>mag size = Nt*Nr,
</code></pre>

<p>In your loop, when i is max and k is max:</p>

<pre><code>k=Nr-1, i=Nt-1, 
</code></pre>

<p>Hence:</p>

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

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

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

python c python-c-api