Äquivalent von Python-Lambda-Funktion für C (Python-Erweiterungen)

Question

Ich habe ein Python-Erweiterungsmodul mit C geschrieben, um die Berechnungszeiten zu beschleunigen. Der erste Schritt ist eine 2D-Integration einer Funktion f (x, y, k), die sehr schnell ist und mir erlaubt, über y in [y1 (x), y2 (x)] und x in [a, b] zu integrieren. während ein Float zu k zugewiesen wird. Aber ich muss k wirklich über den Bereich integrieren [c, d]. Derzeit mache ich so etwas wie dies in Python

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

wo ber meine C-Erweiterungsmodul und calc.kernel sind Anrufe gauss2 2D-Integration durchzuführen. l und ki sind nur andere Variablen. Aber mit meinen Daten dauert Quad noch viele Stunden bis zum Ende. Ich würde gerne alle Berechnungen innerhalb des C-Erweiterungsmoduls durchführen, aber ich bin wirklich ratlos, wie ich dieses äußere Integral implementieren kann. Hier ist meine C-Code

#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 und W96 enthalten sowohl die Punkte für die Gauß-Quadratur, also keine Sorgen machen, dass sie hier sind leer. Ich sollte hinzufügen, dass ich keine Anerkennung für die Funktionen Gauss1 und Gauss2 habe.

EDIT: Python-Code war falsch - jetzt bearbeitet.

Answer 1

Vielleicht ist der Quellcode für scipy integrieren Quad ist ein guter Ort, wenn Sie nicht, dass es ausgesehen haben, zu starten: https://github.com/scipy/scipy/blob/v0.17.0/scipy/integrate/quadpack.py#L45-L360

Sieht aus wie die meiste Arbeit bereits durch native Fortran-Code getan, die normalerweise entweder so schnell oder schneller als C/C++ - Code. Es wird schwierig sein, das zu verbessern, es sei denn, Sie erstellen eine CUDA-Implementierung.

Sie machen den Fortran-Code Multithreading, wenn es nicht bereits und die Quelle geöffnet ist. Zuletzt könnten Sie einen Threading-Dispatcher in C/Fortran machen (Python unterstützt wegen der GIL kein echtes Threading) und machen Sie Ihre Aufrufe von Quads zumindest parallel zueinander. Wenn man direkt mit Fortran Quads rechnet, würde das wahrscheinlich einen ordentlichen Overhead einsparen.