Численная интеграция - Как распараллелить это? - PullRequest
Новая
       7

Численная интеграция - Как распараллелить это?

5 голосов
/ 09 января 2012

Я начинаю с OpenCL, я мог увидеть пример добавления вектора и понять его. Но я думал о методе трапеции. Это код (C) для интегрального вычисления для x ^ 2 в [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;
}

Вопрос в том, to как получить ядро ​​для интегрального вычисления с использованием метода трапеции?

Итак, я подумал об этой идее: partials [i] = integrate (a, a + offset), а затем создать ядро ​​для вычисления суммы партиалов, как упоминалось Patrick87.

Но это лучший способ?

Ответы [ 2 ]

3 голосов
/ 09 января 2012

Трапециевидный метод является лишь небольшим уточнением сумм Римана. Чтобы сделать это параллельно, вам нужно разбить интервал на столько подинтервалов, сколько вы хотите, чтобы были потоки; затем пусть каждый поток интегрирует функцию через свой подинтервал. Наконец, сделайте глобальное уменьшение суммы по всем интегралам, вычисленным на предыдущем этапе. Вы можете поэкспериментировать с тем, сколько потоков использовать на каждом этапе.

2 голосов
/ 16 января 2012

Вот что я придумал.Мне не удалось провести сквозное тестирование этого ядра.Я сделаю обновление, когда у меня будет немного больше времени.

comp_trap - базовый метод «разделяй и властвуй», основанный на коде, который вы использовали выше.comp_trap_multi повышает точность, заставляя каждый рабочий элемент делить его подсекцию

Вам нужно всего лишь выделить массив двойников на хосте, чтобы у каждой рабочей группы был один дубль, чтобы вернуть свой результат.Это должно помочь сократить выделение вектора, которого вы хотели избежать.

Пожалуйста, дайте мне знать, если есть какие-либо проблемы с этим.

Обновлено:

1) изменил вседвойные ссылки на float, потому что double не обязателен в opencl

2) жестко закодировал размер рабочей группы до 64. Это значение является оптимальным в моей системе и должно быть определено экспериментально.Я предпочитаю жестко программировать это значение, а не передавать его в локальный массив с плавающей точкой, потому что в конечном итоге хост-программа будет использовать только оптимальное значение в целевой системе.

3) исправлено неверное вычисление (a1 былнеправильно, должно быть лучше сейчас) 

/*
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];
    }
}
...