Суммирование с OpenMP с использованием C

Question

Я пытался распараллелить этот кусок кода около двух дней и продолжаю иметь логические ошибки. Программа заключается в том, чтобы найти площадь интеграла, используя сумму очень малых значений dx, и рассчитать каждое дискретное значение интеграла. Я пытаюсь реализовать это с помощью openmp, но на самом деле у меня нет опыта работы с openmp. Я хотел бы вашей помощи, пожалуйста. Фактическая цель - распараллелить переменную suma в потоках, чтобы каждый поток вычислял меньше значений интеграла. Программа успешно компилируется, но когда я ее выполняю, она возвращает неверные результаты.

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char *argv[]){
    float down = 1, up = 100, dx, suma = 0, j;
    int steps, i, nthreads, tid;
    long starttime, finishtime, runtime; 

    starttime = omp_get_wtime();
    steps = atoi(argv[1]);
    dx = (up - down) / steps;

    nthreads = omp_get_num_threads();
    tid = omp_get_thread_num();
    #pragma omp parallel for private(i, j, tid) reduction(+:suma)
    for(i = 0; i < steps; i++){
        for(j = (steps / nthreads) * tid; j < (steps / nthreads) * (tid + 1); j += dx){
            suma += ((j * j * j) + ((j + dx) * (j + dx) * (j + dx))) / 2 * dx;
        }
    }
    printf("For %d steps the area of the integral  3 * x^2 + 1 from %f to %f is: %f\n", steps, down, up, suma);
    finishtime = omp_get_wtime();
    runtime = finishtime - starttime;
    printf("Runtime: %ld\n", runtime);
    return (0);
}

Answer 1

Проблема лежит в вашем цикле. Если вы используете for-pragma, OpenMP сделает разделение цикла за вас:

#pragma omp parallel for private(i) reduction(+:suma)
for(i = 0; i < steps; i++) {
    // recover the x-position of the i-th step
    float x = down + i * dx;
    // evaluate the function at x
    float y = (3.0f * x * x + 1)
    // add the sum of the rectangle to the overall integral
    suma += y * dx
}

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

Вам также следует рассмотреть возможность переключения на удвоение для повышения точности.

Answer 2

Давайте просто рассмотрим случай с нитями = 1.Это:

#pragma omp parallel for private(i, j, tid) reduction(+:suma)
for(i = 0; i < steps; i++){
    for(j = (steps / nthreads) * tid; j < (steps / nthreads) * (tid + 1); j += dx){
        suma += ((j * j * j) + ((j + dx) * (j + dx) * (j + dx))) / 2 * dx;
    }
}

превращается в это:

for(i = 0; i < steps; i++){
    for(j = 0; j < steps; j += dx){
        suma += ((j * j * j) + ((j + dx) * (j + dx) * (j + dx))) / 2 * dx;
    }
}

, и вы можете начать видеть проблему;вы в основном зацикливаетесь на шаги ² .

Кроме того, ваш второй цикл не имеет никакого смысла, так как вы увеличиваете на dx.Та же самая путаница между признаками (i, j) и местоположениями в физической области (i * dx) проявляется в вашем приращении.j+dx не имеет никакого смысла.Предположительно, вы хотите увеличить suma на (f (x) + f (x ')) * dx / 2 (например, правило трапеции);это должно быть

        float x = down + i*dx;
        suma += dx * ((3 * x * x + 1) + (3 * (x + dx) * (x + dx) + 1)) / 2;

Как указывает эбо, вы хотите суммировать integrand , а не его антидериватив.

Теперь, если мы включим проверку ответа:

printf("For %d steps the area of the integral  3 * x^2 + 1 from %f to %f is: %f (expected: %f)\n",
            steps, down, up, suma, up*up*up-down*down*down + up - down);

и мы запускаем его последовательно, мы начинаем получать правильный ответ:

$ ./foo 10
For 10 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1004949.375000 (expected: 1000098.000000)
Runtime: 0
$ ./foo 100
For 100 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000146.562500 (expected: 1000098.000000)
Runtime: 0
$ ./foo 1000
For 1000 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000098.437500 (expected: 1000098.000000)
Runtime: 0

Нет никакого смысла беспокоиться о случае OpenMP, пока не сработает случай с последовательным интерфейсом.

Как только ebo указывает на OpenMP, самое простое, что нужно сделать, это просто позволить OpenMP выполнить декомпозицию вашего цикла: например,

#pragma omp parallel for reduction(+:suma)
    for(i = 0; i < steps; i++){
        float x = down + i*dx;
        suma += dx * ((3 * x * x + 1) + (3 * (x + dx) * (x + dx) + 1)) / 2;
    }

Запуск этого,каждый получает

$ setenv OMP_NUM_THREADS 1
$ ./foo 1000
For 1000 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000098.437500 (expected: 1000098.000000)
Runtime: 0
$ setenv OMP_NUM_THREADS 2
$ ./foo 1000
For 1000 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000098.437500 (expected: 1000098.000000)
Runtime: 0
$ setenv OMP_NUM_THREADS 4
$ ./foo 1000
For 1000 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000098.625000 (expected: 1000098.000000)
Runtime: 0
$ setenv OMP_NUM_THREADS 8
$ ./foo 1000
For 1000 steps the area of the integral  3 * x^2 + 1 from 1.000000 to 100.000000 is: 1000098.500000 (expected: 1000098.000000)

Можно явно блокировать в OpenMP, если вы действительно этого хотите, но у вас должна быть причина для этого.