Python C Extension

Question

I am having issues returning a 2D array from a C extension back to Python. When I allocate memory using malloc the returned data is rubbish. When I just initialise an array like sol_matrix[nt][nvar] the returned data is as expected.

#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

// function to be solved by Euler solver
double func (double xt, double y){

    double y_temp = pow(xt, 2); 
    y = y_temp;

    return y;
}

static PyObject* C_Euler(double h, double xn)
{
    double y_temp, dydx;                //temps required for solver
    double y_sav = 0;                   //temp required for solver
    double xt = 0;                      //starting value for xt
    int nvar = 2;                       //number of variables (including time)
    int nt = xn/h;                      //timesteps
    double y = 0;                       //y starting value
            
    //double sol_matrix[nt][nvar];        //works fine
        
    double **sol_matrix = malloc(nt * sizeof(double*));  //doesn't work

    for (int i=0; i<nt; ++i){
        sol_matrix[i] = malloc (nvar * sizeof(double));
    }
    
    int i=0;
    
    //solution loop - Euler method.
    while (i < nt){

        sol_matrix[i][0]=xt;
        sol_matrix[i][1]=y_sav;
        dydx = func(xt, y);          
        y_temp = y_sav + h*dydx;
        xt = xt+h;
        y_sav=y_temp;
        i=i+1;   
    }

    npy_intp dims[2];
    dims[0] = nt;
    dims[1] = 2;
  
    //Create Python object to copy solution array into, get pointer to
    //beginning of array, memcpy the data from the C colution matrix 
    //to the Python object. 
    PyObject *newarray = PyArray_SimpleNew(2, dims, NPY_DOUBLE);
    double *p = (double *) PyArray_DATA(newarray);
    memcpy(p, sol_matrix, sizeof(double)*(nt*nvar));   
    

    // return array to Python
    
    return newarray;      
}

static PyObject* Euler(PyObject* self, PyObject* args)
{
    double h, xn;

    if (!PyArg_ParseTuple(args, "dd", &h, &xn)){
        
        return NULL;
    }
    return Py_BuildValue("O", C_Euler(h,xn));
}

Could you provide any guidance on where I am going wrong?

Thank you.

score 0 · Accepted Answer · answered Oct 10 '20 at 09:15

The data in sol_matrix is not in contiguous memory, it's in nt separately allocated arrays. Therefore the line

memcpy(p, sol_matrix, sizeof(double)*(nt*nvar));

is not going to work.

I'm not a big fan of pointer-to-pointer arrays so believe your best option is to allocate sol_matrix as one big chunk:

double *sol_matrix = malloc(nt*nvar * sizeof(double));

This does mean you can't do 2D indexing so will need to do

// OLD: sol_matrix[i][0]=xt;
sol_matrix[i*nvar + 0] = xt;

In contrast

double sol_matrix[nt][nvar];        //works fine

is a single big chunk of memory so the copy works fine.

Python C Extension

1 Answers1