Integración Numérica - ¿Cómo se puede paralelizar?

Question

Estoy empezando con OpenCL, pude ver el ejemplo de agregar vector y entenderlo. Pero estaba pensando en el método del trapecio. Este es el código (C) para el cálculo integral para x^2 en [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; 
} 

La pregunta es ¿cómo obtener un kernel para el cálculo integral usando el método trapecio?

---

Por lo tanto, yo estaba pensando en la idea: parciales [i] = integran (a, a + offset), y luego hacer un kernel para calcular la suma de los parciales que se han mencionado Patrick87.

Pero, esta es la mejor manera?

Answer 1

Esto es lo que ocurrió. No pude hacer una prueba de extremo a extremo de este kernel. Haré una actualización cuando tenga un poco más de tiempo.

comp_trap es el método de conquista división & básica basada en el código que tiene encima. comp_trap_multi aumenta la precisión haciendo que cada elemento de trabajo divida su subsección

Solo necesita asignar una matriz de dobles en el host para que cada grupo de trabajo tenga un doble para devolver el resultado. Esto debería ayudar a reducir la asignación de vectores que quería evitar.

Por favor, hágamelo saber si hay algún problema con esto.

Actualizado:

1) cambió todas las referencias dobles a flotar, ya que el doble es opcional en OpenCL

2) codificado el tamaño del grupo de trabajo a 64. Este valor es óptima en mi sistema, y ​​debe ser determinado experimentalmente Prefiero codificar duro este valor sobre pasar en una matriz local de flotantes para usar, porque el programa host eventualmente usará solo el valor óptimo en el sistema objetivo de todos modos.

3) fijado un cálculo incorrecto (a1 estaba mal, debe ser mejor ahora)

/* 
numerical-integration.cl 
*/ 

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

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

__kernel void comp_trap(
    float a, 
    float b, 
    __global float* sums) 
{ 
/* 
- assumes 1D global and local work dimensions 
- each work unit will calculate 1/get_global_size of the total sum 
- the 0th work unit of each group then accumulates the sum for the 
group and stores it in __global * sums 
- memory allocation: sizeof(sums) = get_num_groups(0) * sizeof(float) 
- assumes local scratchpad size is at lease 8 bytes per work unit in the group 
ie sizeof(wiSums) = get_local_size(0) * sizeof(float) 
*/ 
    __local float wiSums[64]; 
    int l_id = get_local_id(0); 

    //cumpute range for this work item is: a1, b1 
    float a1 = a+((b-a)/get_global_size(0))*get_global_id(0); 
    float b1 = a1+(b-a)/get_global_size(0); 

    wiSums[l_id] = simple_trap(a1,b1); 

    barrier(CLK_LOCAL_MEM_FENCE); 

    int i; 
    if(l_id == 0){ 
     for(i=1;i<get_local_size(0);i++){ 
      wiSums[0] += wiSums[i]; 
     } 
     sums[get_group_id(0)] = wiSums[0]; 
    } 
} 

__kernel void comp_trap_multi(
    float a, 
    float b, 
    __global float* sums, 
    int divisions) 
{ 
/* 
- same as above, but each work unit further divides its range into 
'divisions' equal parts, yielding a more accurate result 
- work units still store only one sum in the local array, which is 
used later for the final group accumulation 
*/ 
    __local float wiSums[64]; 
    int l_id = get_local_id(0); 

    float a1 = a+((b-a)/get_global_size(0))*get_global_id(0); 
    float b1 = a1+(b-a)/get_global_size(0); 
    float range; 
    if(divisions > 0){ 
     range = (b1-a1)/divisions; 
    }else{ 
     range = (b1-a1); 
    } 

    int i; 
    wiSums[l_id] = 0; 
    for(i=0;i<divisions;i++){ 
     wiSums[l_id] += simple_trap(a1+range*i,a1+range*(i+1)); 
    } 

    barrier(CLK_LOCAL_MEM_FENCE); 

    if(l_id == 0){ 
     for(i=1;i<get_local_size(0);i++){ 
      wiSums[0] += wiSums[i]; 
     } 
     sums[get_group_id(0)] = wiSums[0]; 
    } 
} 

Answer 2

método del trapecio es sólo un ligero refinamiento de hacer las sumas de Riemann. Para hacer esto en paralelo, querrás dividir el intervalo en tantos subintervalos como quieras que haya hilos; luego, haga que cada subproceso integre la función sobre su subintervalo. Finalmente, realice una reducción de suma global sobre todas las integrales calculadas en la fase anterior. Puedes experimentar con cuántos hilos usar para cada etapa.