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

[artey]

Здравствуйте, написал параллельную программу на OpenMP для решения системы разностных уравнений методом переменных направлений. Проблема в следующем - последовательная версия программы выполняется на сервере за 0.03 секунды, а кода запускаю на 2 и более потоках, то программа превышает лимит времени для выполнения на сервере и файл вывода, в котором выводится время выполнения программы оказывается пустой. Поэтому я даже не знаю сколько она выполняется. Но по идее на 2 и более потоках она должна работать быстрее. В чем дело?

Вот код:

Используется синтаксис 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;
}
double max(double a, double b)
{
        if (a > b)
                return a;
        else
                return b;
}
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.99;
        

        double f1[n + 1][n + 1];
        //решение на n-шаге
        double C0[n + 1][n + 1];
        //решение на (n+1/2)-шаге
        double C1[n + 1][n + 1];
        //решение на (n+1)-шаге
        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();
        while (error > 1e-7)
        {

                #pragma omp parallel for private(i,j) shared(C0)
                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];

                
                //прогонка по столбцам
                #pragma omp parallel for private(j) shared(C1)
                for (j = 1; j < n; j++)
                {
                        a = -a1; b = -(1 / tau + e2); c = -c1;
                        i = 1;
                        d[1] = d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
                        i = n-1;
                        d[n - 1] = b1 * C0[i][j-1] +  (1 / tau - e3)*C0[i][j] + f1[i][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+1/2)
                #pragma omp parallel for private(i) shared(C1)
                for (i = 0; i < n + 1; i++)
                        C1[n][i] = C0[n - 1][i];
                #pragma omp parallel for private(i) shared(C1)
                for (i = 0; i < n + 1; i++)
                        C1[i][n] = C0[i][n - 1];
                #pragma omp parallel for private(i) shared(C1)
                for (i = 0; i < n + 1; i++)
                        C1[0][i] = 0.0;
                #pragma omp parallel for private(i) shared(C1)
                for (i = 0; i < n + 1; i++)
                        C1[i][0] = 0.0;


                //прогонка по строкам
                #pragma omp parallel for private(i) shared(C2)
                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];
                        }
                }



                //граничные условия шаг (n+1)
                #pragma omp parallel for private(i) shared(C2)
                for (i = 0; i < n + 1; i++)
                        C2[n][i] = C1[n - 1][i];
                #pragma omp parallel for private(i) shared(C2)
                for (i = 0; i < n + 1; i++)
                        C2[i][n] = C1[i][n - 1];
                #pragma omp parallel for private(i) shared(C2)
                for (i = 0; i < n + 1; i++)
                        C2[0][i] = 0.0;
                #pragma omp parallel for private(i) shared(C2)
                for (i = 0; i < n + 1; i++)
                        C2[i][0] = 0.0;



                error = 0.0;
                #pragma omp parallel for private(i,j) reduction(+:error)
                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);
                

        }
        t=omp_get_wtime()-t;
        FILE *f;
        f=fopen("result.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;
}
 

[mihaild]

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

[artey]

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

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

Нет возможности запускать локально. Я технологию эту недавно начал изучать, поэтому как тут что работает еще много не знаю, к сожалению.

[Dmitriy40]

Я для отладки часто сую во все дыры счётчики, потом в конце (а иногда и в промежуточных точках) вывожу их в файл/экран и анализирую, что куда сколько раз попало и сколько я ожидал. Отлавливается много разных ошибок, в том числе и незаметных по времени исполнения кода.
Насчёт совместной записи в одно место, я думаю в openmp это как раз решено на уровне компилятора, ведь это один из основных принципов использования распареллеливания. Но не уверен.
Попробуйте удалить (закомментировать) все #pragma и возвращать их по одной — найдёте где подвисает. А дальше уже думать почему.

[mihaild]

Dmitriy40 в сообщении #1433862 писал(а):

Насчёт совместной записи в одно место, я думаю в openmp это как раз решено на уровне компилятора

Нет, не решено, в чем легко убедиться на практике. И собственно непонятно, как это можно было бы сделать на уровне компилятора, кроме взятия лока на каждое обращение к shared переменным, но это понятно с какой скоростью работать будет.

[photon]

С телефона неудобно смотреть код - посмотрю потом внимательнее, но насколько я увидел, там есть циклы, где идет переписывание из одного массива в другой, причем подряд по памяти. Подозреваю, что распараллеливание только замедлит работу. Быстрее просто копировать память стандартными средствами

[Dmitriy40]

(Оффтоп)

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

Dmitriy40 в сообщении #1433862 писал(а):

Насчёт совместной записи в одно место, я думаю в openmp это как раз решено на уровне компилятора

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

Спасибо за уточнение. Просто это один из базовых примеров использования openmp что я везде видел (заполнение большого массива), вот и подумал что предусмотрели деление цикла на N одинаковых больших частей, заполняемых потоками параллельно. Т.е. не вперемешку, а с большим смещением, а значит и конфликтов будет ничтожно мало.
Впрочем я с openmp никогда не разбирался.

[mihaild]

(Оффтоп)

Dmitriy40 в сообщении #1433889 писал(а):

предусмотрели деление цикла на N одинаковых больших частей, заполняемых потоками параллельно. Т.е. не вперемешку, а с большим смещением

Я не про то. Итерации бьются на чанки, и чанки как-то распределяются по потокам (это может происходить разными способами). Но в приведенном коде треды пишут не в соседние места, а ровно в одно и то же без синхронизации, что является UB.

[Dmitriy40]

(Оффтоп)

mihaild
А, Вы видимо про массивы P,Q,d и даже некоторые переменные, это да, я не обратил внимания. Тогда полностью согласен, это должен решать программист. Оказывается программа гораздо хуже чем мне казалось.

[photon]

Да, посмотрел теперь внимательнее. Результат с превышением времени "немного предсказуем". В параллельных потоках массивы P, Q, C2, d портят друг друга, так что вместо корректных результатов будет какой-то хлам, соответственно, ожидать, что выполнится условие малости ошибки для выхода из цикла while не стоит. Ускорить работу этой программы, уверен, можно, но, думаю, не путём применения omp.

[artey]

photon, к сожалению, надо именно с помощью openmp. Я попытался внести изменения в программу, но они не дали результатов:

Используется синтаксис 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;
}
double max(double a, double b)
{
        if (a > b)
                return a;
        else
                return b;
}
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.9;
        

        double f1[n + 1][n + 1];
        //решение на n-шаге
        double C0[n + 1][n + 1];
        //решение на (n+1/2)-шаге
        double C1[n + 1][n + 1];
        //решение на (n+1)-шаге
        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();
        while (error > 1e-7)
        {

                #pragma omp parallel for private(i,j) shared(C0,C2)
                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)
                #pragma omp parallel for private(i) shared(C0,C1)
                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 for private(j) shared(a,b,c,i,d,C0,P,Q,C1)
                for (j = 1; j < n; j++)
                {
                        a = -a1; b = -(1 / tau + e2); c = -c1;
                        i = 1;
                        d[1] = d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
                        i = n-1;
                        d[n - 1] = b1 * C0[i][j-1] +  (1 / tau - e3)*C0[i][j] + f1[i][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+2)
                #pragma omp parallel for private(i) shared(C1,C2)
                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 for private(i) shared(a,b,c,j,d,C1,P,Q,C2)
                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;
                #pragma omp parallel for private(i,j) reduction(+:error)
                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);
                

        }
        t=omp_get_wtime()-t;
        FILE *f;
        f=fopen("result.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;
}
 

[photon]

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

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

к сожалению, надо именно с помощью openmp

Вам надо шашечки или ехать? - в смысле, это учебная задача на использование OpenMP или вам надо добиться максимального ускорения?

[mihaild]

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

а чтобы избежать совместной записи их надо делать private

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

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

[photon]

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

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

Ну да, "заприватить" надо не только указатель, наверное этого можно добиться отказавшись от массивов, а перейдя к векторам и уже их объявлять как private

[artey]

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

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

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

к сожалению, надо именно с помощью openmp

Вам надо шашечки или ехать? - в смысле, это учебная задача на использование OpenMP или вам надо добиться максимального ускорения?

это учебная задача