python - C 的 python lambda 函数的等效项(Python 扩展)

Question

我用 C 语言编写了一个 Python 扩展模块来加快计算速度。第一步是函数 f(x,y,k) 的 2D 积分，速度非常快，允许我对 [y1(x),y2(x)] 中的 y 和 [a,b] 中的 x 进行积分同时为 k 分配一个 float 。但我确实需要在 [c,d] 范围内对 k 进行积分。目前，我正在 Python 中做类似的事情

inner = lambda k: calc.kernel(l,k,ki)
I = quad(inner,c,d)[0]

其中 calc 是我的 C 扩展模块，calc.kernel 调用 gauss2 来执行 2D 积分。 l 和 ki 只是其他变量。但根据我的数据，quad 仍然需要很多小时才能完成。我想在 C 扩展模块内进行所有计算，但我真的很难理解如何实现这个外部积分。这是我的 C 代码

#include <Python.h>
#include <math.h>

double A96[96]={        /* abscissas for 96-point Gauss   quadrature */
    };

double W96[96]={            /* weights for 96-point Gauss quadrature */
    };

double Y1(double x){
    return 0;
    }

double Y2(double x){
    return x;
    }

double gauss1(double F(double),double a,double b)
{ /* 96-pt Gauss qaudrature integrates F(x) from a to b */
int i;
double cx,dx,q;
cx=(a+b)/2;
dx=(b-a)/2;
q=0;
for(i=0;i<48;i++)
    q+=W96[i]*(F(cx-dx*A96[i])+F(cx+dx*A96[i]));
return(q*dx);
}

double gauss2(double F(double,double,int,double,double),double Y1(double),double Y2(double),double a,double b,int l,double k, double ki)
{/* 96x96-pt 2-D Gauss qaudrature integrates
           F(x,y) from y=Y1(x) to Y2(x) and x=a to b */
int i,j,h;
double cx,cy,dx,dy,q,w,x,y1,y2;
cx=(a+b)/2;
dx=(b-a)/2;
q=0;
for(i=0;i<48;i++)
    {
    for(h=-1;h<=1;h+=2)
        {
        x=cx+h*dx*A96[i];
        y1=Y1(x);
        y2=Y2(x);
        cy=(y1+y2)/2;
        dy=(y2-y1)/2;
        w=dy*W96[i];
        for(j=0;j<48;j++)
            q+=w*W96[j]*(F(x,cy-dy*A96[j],l,k,ki)+F(x,cy+dy*A96[j],l,k,ki));
        }
    }
return(q*dx);
}

double ps_fact(double z){
    double M = 0.3;
    return 3/2*(M*(1+z)*(1+z)*(1+z) + (1-M))*(M*(1+z)*(1+z)*(1+z) + (1-M))*(M*(1+z)*(1+z)*(1+z) + (1-M))/(1+z)/(1+z);
    }
double drdz(double z){
    double M = 0.3;
    return 3000/sqrt(M*(1+z)*(1+z)*(1+z) + (1-M));
}

double rInt(double z){
    double M = 0.3;
    return 3000/sqrt(M*(1+z)*(1+z)*(1+z) + (1-M));
}

double kernel_func ( double y , double x, int l,double k, double ki) {
    return ps_fact(y)*ki*rInt(x)*sqrt(M_PI/2/rInt(x))*jn(l+0.5,ki*rInt(x))*drdz(x)*(rInt(x)-rInt(y))/rInt(y)*sqrt(M_PI/2/rInt(y))*jn(l+0.5,k*rInt(y))*drdz(y);
}

static PyObject* calc(PyObject* self, PyObject* args)
{
    int l;
    double k, ki;
    if (!PyArg_ParseTuple(args, "idd", &l, &k, &ki))
        return NULL;

    double res;
    res = gauss2(kernel_func,Y1, Y2, 0,10,l, k, ki);

    return Py_BuildValue("d", res);
}

static PyMethodDef CalcMethods[] = {
    {"kernel", calc, METH_VARARGS, "Calculates kernel values."},
    {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC initcalc(void){
    (void) Py_InitModule("calc", CalcMethods);

A96 和 W96 都包含高斯求积的点，所以不用担心它们在这里是空的。我应该补充一点，我对函数 gauss1 和 gauss2 没有任何功劳。

编辑:python 代码错误 - 现已编辑。

Answer 1

如果您还没有看过，也许 scipyintegratequad 的源代码是一个很好的起点:https://github.com/scipy/scipy/blob/v0.17.0/scipy/integrate/quadpack.py#L45-L360

看起来大部分工作已经由 native Fortran 代码完成，通常与 C/C++ 代码一样快或更快。除非您创建/找到 CUDA 实现，否则您将很难对此进行改进。

您可以将 Fortran 代码设置为多线程(如果尚未设置且源代码已开放)。最后，您可以在 C/Fortran 中创建一个线程调度程序(Python 由于 GIL 的原因不支持真正的线程)，并且至少可以相互调用四并行。直接将 calc 与 Fortran 四元接口(interface)连接可能也会为您节省一些可观的开销。