Numerical Integration - How to parallelize it?

Question

I'm getting starting with OpenCL, I could see the add vector example and understand it. But I was thinking about the trapezium method. This is the code ( C ) for the integral calculation for x^2 in [a,b].

double f(double x)
{
    return x*x;
}

double Simple_Trap(double a, double b)
{
    double fA, fB;
    fA = f(a);
    fB = f(b);
    return ((fA + fB) * (b-a)) / 2;
}

double Comp_Trap( double a, double b)
{
    double Suma = 0;
    double i = 0;
    i = a + INC;
    Suma += Simple_Trap(a,i);
    while(i < b)
    {
        i+=INC;
        Suma += Simple_Trap(i,i + INC);
    }
    return Suma;
}

The question is ¿how to obtain a kernel for integral calculation using the trapezium method?

---

So, I was thinking about the idea: partials[i] = integrate(a,a+offset), and then make a kernel to compute the sum of partials as mentioned Patrick87.

But, this is the best way?

Answer 1

The trapezoid method is just a slight refinement of doing Riemann sums. To do this in parallel, you'll want to break the interval into as many subintervals as you want there to be threads; then, have each thread integrate the function over its subinterval. Finally, do a global sum reduction over all the integrals computed in the previous phase. You can experiment with how many threads to use for each stage.