Блог

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

         num_chunks = 1000;
        int threads_assigned = 0;
        int num_threads = 10;// omp_get_num_threads();
        sum_arr = (double *) malloc(sizeof(double) * num_threads);
        memset(sum_arr, 0, num_threads);
        int sums_each_thread = num_chunks / num_threads;
        printf("num of threads = %dn", num_threads);
        sum = 0.0;
        omp_set_num_threads(num_threads);
        chunk_size = 1.0/(double) num_chunks;
        double start_time = omp_get_wtime();
        #pragma omp parallel
        {
                int idx = omp_get_thread_num();
                if (idx == 0) threads_assigned = omp_get_num_threads();
                int sums_count;
                int loop_size = idx * sums_each_thread;
                for (sums_count = loop_size; sums_count < loop_size   sums_each_thread; sums_count  ) {
                        printf("starting thread no %dn", idx);
                        // double x = (idx * chunk_size)   (chunk_size / 2.0); simplified version in the next line
                        double x = chunk_size * (sums_count   0.5);
                        double y = 4.0 / (1   x*x);
                        sum_arr[idx - 1]  = (y * chunk_size);
                        sum  = y;
                }
                printf("ending thread no %dn", idx);
        }
        double end_time = omp_get_wtime();
        printf("sum = %f, in time %f(s)n", sum * chunk_size, end_time - start_time);
        int new_idx = 0;
        sum = 0.0;
        for(new_idx = 0; new_idx < num_threads; new_idx  ) {
                sum  = sum_arr[new_idx];
        }
        printf("arr sum = %f, with total assigned threads = %dn", sum, threads_assigned);
 

Я рассчитал площадь по:

1 ) вычисление суммы каждой области прямоугольника.

2 ) вычисление суммы всех высот (скалярной суммы), а затем кратной размеру фрагмента.

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

Кто-нибудь может объяснить причину?

В вашем коде есть две ошибки:

Один из них-доступ к недопустимой записи в arr_sum том, где к элементу -1 обращается главный поток. Заменять

 sum_arr[idx - 1]  = (y * chunk_size);
 

с

 sum_arr[idx]  = (y * chunk_size);
 

Другой проблемой является модификация sum без синхронизации.
Извлечение предыдущего sum и сохранение обновленного значения может чередоваться с аналогичной операцией в другом потоке. Это испортит результат. Этот эффект должен быть заметен после удаления printf("starting thread no %dn", idx); из внутреннего цикла. Чтобы исправить это, просто используйте атомарные обновления.

 #pragma omp atomic
sum  = y;
 

     #pragma omp parallel for reduction( :sum)
    for (int i = 0; i < num_chunks;   i) {
        double x = chunk_size * (i   0.5);
        double y = 4.0 / (1   x*x);
        sum  = y;
    }