Параллельная программа на OpenMP

[artey]

mihaild в сообщении #1433988 писал(а):

artey, у вас по-прежнему гонка: разные потоки пишут в одно и то же место, и в итоге получается непонятно что.

Может тогда как то вычисления по-другому оформить можно, чтобы можно было распараллелить удобней? Что посоветуете?

[Mihaylo]

Ну все же вроде просто:
1. Инициализируются различные переменные и массивы C0, C1, C2, f1.
2. Потенциальный поток 1: Присваиваются значения a, b, c, пересчитываются вектора d, P, Q, вычисляется массив C1
3. Потенциальный поток 2: Присваиваются значения a, b, c, пересчитываются вектора d, P, Q, вычисляется массив C2
4. Вычисление error.

Проблема в том, что поток 2 пишет в те же переменные a, b, c, d, P, Q, что и поток 1. Когда потоки выполняются последовательно, то проблема не возникает, так как после выполнения потока 1 данные в этих переменных не нужны, а поток 2 никак не использует. Но когда параллельный поток еще выполняется, то они мешают друг другу. Выход: использовать разные переменные a, b, c, d, P, Q. Красиво это делается путем вынесения повторного кода в отдельную функцию.

[artey]

Mihaylo в сообщении #1434156 писал(а):

Ну все же вроде просто:
1. Инициализируются различные переменные и массивы C0, C1, C2, f1.
2. Потенциальный поток 1: Присваиваются значения a, b, c, пересчитываются вектора d, P, Q, вычисляется массив C1
3. Потенциальный поток 2: Присваиваются значения a, b, c, пересчитываются вектора d, P, Q, вычисляется массив C2
4. Вычисление error.

Проблема в том, что поток 2 пишет в те же переменные a, b, c, d, P, Q, что и поток 1. Когда потоки выполняются последовательно, то проблема не возникает, так как после выполнения потока 1 данные в этих переменных не нужны, а поток 2 никак не использует. Но когда параллельный поток еще выполняется, то они мешают друг другу. Выход: использовать разные переменные a, b, c, d, P, Q. Красиво это делается путем вынесения повторного кода в отдельную функцию.

А если потоков больше двух, например 3 или 4, то как тогда?

[mihaild]

artey в сообщении #1434138 писал(а):

Может тогда как то вычисления по-другому оформить можно, чтобы можно было распараллелить удобней?

Делать вспомогательные массивы приватными для потоков. В случае omp, кажется, самый простой способ - объявить parallel секцию, внутри которой их создавать, и внутри которой же будет for. Либо (лучше) заранее выделить память для всех потоков (переменные превращаются в массивы, одномерные массивы в двумерные и т.д.), и использовать get_thread_num в качестве индекса.

[artey]

Внес изменения в код:

Используется синтаксис C++

#include <iostream>
#include <fstream>
#include <math.h>
#include <omp.h>
using namespace std;
const int n = 100;
int i00 = n / 4;
int j00 = n / 4;
//определение источника
double diraca(int x)
{
        if (x == 0)
                return 1.0;
        else
                return 0.0;
}
int main()
{
        double V, U;
        U = 1.0; V = 1.0;
        double D = 0.01;
        const double h = 100.0 / double(n);
        double a1 = D / (h*h) + U / h;
        double b1 = D / (h*h) + V / h;
        double c1 = D / (h*h);
        double d1 = D / (h*h);
        double Q1 = 10.0;
        int i, j;
        double e1 = 2.0*D / (h*h) + U / h + 2.0*D / (h*h) + V / h;
        double e2 = a1 + c1;
        double e3 = b1 + d1;
        double tau = 0.95;
        

        double f1[n + 1][n + 1];
        //решение на n-1шаге
        double C0[n + 1][n + 1];
        //решение на (n-1/2)-шаге
        double C1[n + 1][n + 1];
        //решение на n-шаге
        double C2[n + 1][n + 1];
        //прогоночные коэфф.-ты
        double P[n + 1], Q[n + 1];
        //начальные условия распределения примеси
        for (i = 0; i < n + 1; i++)
                for (j = 0; j < n + 1; j++) 
                {
                        f1[i][j] = Q1 * diraca(i - i00) * diraca(j - j00);
                        C2[i][j] = 0.0;
                        C1[i][j] = 0.0;
                        C0[i][j] = 0.0;
                }

        double error = 1.0;
        double t=omp_get_wtime();
        int iter=0;
        while (error > 1e-7)
        {

                for (i = 0; i < n; i++)
                    for (j = 0; j < n; j++)
                        C0[i][j] = C2[i][j];
                                
                
                double a, b, c, d[n + 1];
                
                
                //граничные условия шаг (n-1/2)
                for (i = 0; i < n+1; i++)
                { 
                  C1[n][i] = C0[n - 1][i];
                  C1[i][n] = C0[i][n - 1];
                  C1[0][i] = 0.0;
                  C1[i][0] = 0.0;
                }

                
                //прогонка по столбцам
                #pragma omp parallel private(a,b,c,d,P,Q) shared(f1,C0,C1,n)
                {
                #pragma omp parallel for  
                for (j = 1; j < n; j++)
                {
                        a = -a1; b = -(1 / tau + e2); c = -c1;
                        d[1] = d1 * C0[1][j + 1] + (1 / tau - e3)*C0[1][j] + f1[1][j];
                        d[n - 1] = b1 * C0[n-1][j-1] +  (1 / tau - e3)*C0[n-1][j] + f1[n-1][j];
                        
                        for (i = 2; i < n-1; i++)
                                d[i] = b1 * C0[i][j - 1] + d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
                        P[1] = c / b;
                        Q[1] = -d[1] / b;
                        for (i = 2; i < n; i++)
                        {
                                P[i] = c / (-a * P[i - 1] + b);
                                Q[i] = (-d[i] + a * Q[i - 1]) / (-a * P[i - 1] + b);
                        }
                        C1[n - 1][j] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
                        for (i = n - 2; i >= 1; i--)
                        {
                                C1[i][j] = P[i] * C1[i + 1][j] + Q[i];
                        }
                }
                }
                
                //граничные условия шаг (n)
                for (i = 0; i < n+1; i++)
                { C2[n][i] = C1[n - 1][i];
                  C2[i][n] = C1[i][n - 1];
                  C2[0][i] = 0.0;
                  C2[i][0] = 0.0;
                }

                                
                //прогонка по строкам
                #pragma omp parallel private(a,b,c,d,P,Q) shared(f1,C1,C2,n)
                {
                #pragma omp parallel for 
                for (i = 1; i < n; i++)
                {
                        a = -b1; b = -(1 / tau + e3); c = -d1;
                        j = 1;
                        d[1] = C1[i+1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
                        j = n-1;
                        d[n - 1] = a1 * C1[i - 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j]; 
                        for (j = 2; j < n-1; j++)
                                d[j] = a1 * C1[i - 1][j] + c1 * C1[i + 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
                        P[1] = c / b;
                        Q[1] = -d[1] / b;
                        for (j = 2; j < n; j++)
                        {
                                P[j] = c / (-a * P[j - 1] + b);
                                Q[j] = (-d[j] + a * Q[j - 1]) / (-a * P[j - 1] + b);
                        }
                        C2[i][n - 1] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
                        for (j = n - 2; j >= 1; j--)
                        {
                                C2[i][j] = P[j] * C2[i][j + 1] + Q[j];
                        }
                }
                }       
                
                error = 0.0;
                
                for (i = 1; i < n; i++)
                        for (j = 1; j < n; j++)
                                error += pow(fabs(C2[i][j] - C0[i][j]), 2);
                error = pow(error / (n*n), 0.5);
                iter++;
                printf("Iteration=  %i  %lf\n",iter,error);

        }
        t=omp_get_wtime()-t;
        FILE *f;
        f=fopen("result_2.txt","w");
        for(i=1;i<n;i++)
            for(j=1;j<n;j++)
                fprintf(f,"%i %i %lf\n",i,j,C2[i][j]);
        fclose(f);      
        printf("Time=%lf\n",t);
    return 0;
}

 

Результаты:
1 поток - 0.032040 секунды;
2 потока - 0.053495 секунды;
3 потока - 0.074289 секунды;
4 потока 0.084838 секунды;
Нет никакого ускорения. В чем дело?

[artey]

Извиняюсь, опечатка в предыдущем коде - для цикла тоже указал директиву parallel. Исправил:

Используется синтаксис C++

#include <iostream>
#include <fstream>
#include <math.h>
#include <omp.h>
using namespace std;
const int n = 100;
int i00 = n / 4;
int j00 = n / 4;
//определение источника
double diraca(int x)
{
        if (x == 0)
                return 1.0;
        else
                return 0.0;
}
int main()
{
        double V, U;
        U = 1.0; V = 1.0;
        double D = 0.01;
        const double h = 100.0 / double(n);
        double a1 = D / (h*h) + U / h;
        double b1 = D / (h*h) + V / h;
        double c1 = D / (h*h);
        double d1 = D / (h*h);
        double Q1 = 10.0;
        int i, j;
        double e1 = 2.0*D / (h*h) + U / h + 2.0*D / (h*h) + V / h;
        double e2 = a1 + c1;
        double e3 = b1 + d1;
        double tau = 0.95;
        

        double f1[n + 1][n + 1];
        //решение на n-1шаге
        double C0[n + 1][n + 1];
        //решение на (n-1/2)-шаге
        double C1[n + 1][n + 1];
        //решение на n-шаге
        double C2[n + 1][n + 1];
        //прогоночные коэфф.-ты
        double P[n + 1], Q[n + 1];
        //начальные условия распределения примеси
        for (i = 0; i < n + 1; i++)
                for (j = 0; j < n + 1; j++) 
                {
                        f1[i][j] = Q1 * diraca(i - i00) * diraca(j - j00);
                        C2[i][j] = 0.0;
                        C1[i][j] = 0.0;
                        C0[i][j] = 0.0;
                }

        double error = 1.0;
        double t=omp_get_wtime();
        int iter=0;
        int nt;
        while (error > 1e-7)
        {

                for (i = 0; i < n; i++)
                    for (j = 0; j < n; j++)
                        C0[i][j] = C2[i][j];
                                
                
                double a, b, c, d[n + 1];
                
                
                //граничные условия шаг (n-1/2)
                for (i = 0; i < n+1; i++)
                { 
                  C1[n][i] = C0[n - 1][i];
                  C1[i][n] = C0[i][n - 1];
                  C1[0][i] = 0.0;
                  C1[i][0] = 0.0;
                }

                //прогонка по столбцам
                #pragma omp parallel default(shared) 
                { nt = omp_get_num_threads();
                #pragma omp for private(a,b,c,d,P,Q)   
                for (j = 1; j < n; j++)
                {
                        a = -a1; b = -(1 / tau + e2); c = -c1;
                        d[1] = d1 * C0[1][j + 1] + (1 / tau - e3)*C0[1][j] + f1[1][j];
                        d[n - 1] = b1 * C0[n-1][j-1] +  (1 / tau - e3)*C0[n-1][j] + f1[n-1][j];
                        
                        for (i = 2; i < n-1; i++)
                                d[i] = b1 * C0[i][j - 1] + d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
                        P[1] = c / b;
                        Q[1] = -d[1] / b;
                        for (i = 2; i < n; i++)
                        {
                                P[i] = c / (-a * P[i - 1] + b);
                                Q[i] = (-d[i] + a * Q[i - 1]) / (-a * P[i - 1] + b);
                        }
                        C1[n - 1][j] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
                        for (i = n - 2; i >= 1; i--)
                        {
                                C1[i][j] = P[i] * C1[i + 1][j] + Q[i];
                        }
                }
                }
                
                //граничные условия шаг (n)
                for (i = 0; i < n+1; i++)
                { C2[n][i] = C1[n - 1][i];
                  C2[i][n] = C1[i][n - 1];
                  C2[0][i] = 0.0;
                  C2[i][0] = 0.0;
                }

                                
                //прогонка по строкам
                #pragma omp parallel default(shared)
                {
                #pragma omp for private(a,b,c,d,P,Q)  
                for (i = 1; i < n; i++)
                {
                        a = -b1; b = -(1 / tau + e3); c = -d1;
                        j = 1;
                        d[1] = C1[i+1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
                        j = n-1;
                        d[n - 1] = a1 * C1[i - 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j]; 
                        for (j = 2; j < n-1; j++)
                                d[j] = a1 * C1[i - 1][j] + c1 * C1[i + 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
                        P[1] = c / b;
                        Q[1] = -d[1] / b;
                        for (j = 2; j < n; j++)
                        {
                                P[j] = c / (-a * P[j - 1] + b);
                                Q[j] = (-d[j] + a * Q[j - 1]) / (-a * P[j - 1] + b);
                        }
                        C2[i][n - 1] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
                        for (j = n - 2; j >= 1; j--)
                        {
                                C2[i][j] = P[j] * C2[i][j + 1] + Q[j];
                        }
                }
                }       
                
                error = 0.0;
                
                for (i = 1; i < n; i++)
                        for (j = 1; j < n; j++)
                                error += pow(fabs(C2[i][j] - C0[i][j]), 2);
                error = pow(error / (n*n), 0.5);
                iter++;
                printf("Iteration=  %i  %lf\n",iter,error);

        }
        t=omp_get_wtime()-t;
        FILE *f;
        if(nt==1) f=fopen("result1.txt","w");
        if(nt==2) f=fopen("result2.txt","w");
        if(nt==3) f=fopen("result3.txt","w");
        if(nt==4) f=fopen("result4.txt","w");
        for(i=1;i<n;i++)
            for(j=1;j<n;j++)
                fprintf(f,"%i %i %lf\n",i,j,C2[i][j]);
        fclose(f);      
        printf("Time=%lf\n",t);
        printf("num_threads=%d\n",nt);
        
    return 0;
}
 

Результаты:
1 поток - 0.029425 секунды;
2 потока - 0.019436 секунды;
3 потока - 0.023503секунды;
4 потока 0.021606 секунды;

Ещё можно что-то тут улучшить? А то на 2 потоках быстрее вычисления получились, чем на 3 и 4.

[Dmitriy40]

artey
Задам тупой, но необходимый вопрос: а Вам известно с какой погрешностью измеряется выдаваемое вам время? Может тут погрешность 10мс и все эти числа равны друг другу? Я просто не в курсе как измеряется время в openmp на серверах.

[photon]

Если точность неизвестна, то можно оценить по разбросу результатов при нескольких запусках

[artey]

photon в сообщении #1434389 писал(а):

Если точность неизвестна, то можно оценить по разбросу результатов при нескольких запусках

На каждом потоке запускал по 3 раза, брал минимальное значение.

artey в сообщении #1434368 писал(а):

1 поток - 0.029425 секунды;
2 потока - 0.019436 секунды;
3 потока - 0.023503 секунды;
4 потока 0.021606 секунды;

[photon]

По возможности, обработку чего угодно, имеющего вид таблицы (сетка какая-то, изображение и т.п.) следует стараться делать так, чтобы внешние циклы шли по строкам, а внутренние по столбцам, тогда во внутренних циклах будут обращения по памяти подряд, а не вразброс. При распараллеливании внешнего цикла тоже выйдет, что каждый поток работает с кусочком памяти, локализованным в одном месте. Дёргание каждым потоком памяти в разных местах может стать бутылочным горлышком, из-за которого распараллеливание не даст эффекта или вовсе ухудшит результат.

[Mihaylo]

artey в сообщении #1434291 писал(а):

А если потоков больше двух, например 3 или 4, то как тогда?

В Вашей задаче очевидное оптимальное число потоков - 2. Увеличение числа потоков можно сделать после тщательного анализа и перестройки алгоритма, но я думаю, что это будет сложно и не очень эффективно.

[photon]

Mihaylo в сообщении #1434511 писал(а):

В Вашей задаче очевидное оптимальное число потоков - 2

Для меня не очевидно. Расскажите, почему именно два?

[Mihaylo]

Имеется два нагруженных цикла for с похожими вычислениями, но разными выходными данными. Они легче всего распараллеливаются. Поэтому 2 потока очевидны.

[photon]

Очевидно, вы не разбираетесь 1) в теме распараллеливания. В частности распараллеливания при помощи OpenMP: parallel for разбивает цикл на куски, которые выполняются в разных потоках, неважно, сколько раз в тексте программы встречаются циклы - исходя из вашей логики, программа с одним циклом однозначно оптимальна без распараллеливания, что, очевидно, верно далеко не всегда; 2) В оценке нагруженности. Скорее всего эти два цикла выполняются сильно разное время, несмотря на внешнюю похожесть и количество операций.

[Mihaylo]

Странно, что два потока не очевидны.