Python C扩展

问题描述

我在将2D数组从C扩展名返回到Python时遇到问题。当我使用malloc分配内存时,返回的数据是垃圾。当我仅初始化sol_matrix [nt] [nvar]之类的数组时,返回的数据就如预期的那样。

#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));
}

您能提供有关我哪里出问题的任何指导吗?

谢谢。

解决方法

sol_matrix中的数据不在连续的内存中,而是在nt 单独分配的数组中。因此,线

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

不起作用。

我不是指针对指针数组的忠实拥护者,因此,您最好的选择是将sol_matrix分配为一个大块:

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

这确实意味着您无法进行2D索引编制,因此需要这样做

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

相反

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

是一个很大的内存块,因此副本可以正常工作。