Матричный метод решения СЛАУ (распараллеливание с MPI)

Тема в разделе "Общие вопросы по С и С++", создана пользователем rrrFer, 4 окт 2016.

  1. rrrFer

    rrrFer Well-Known Member
    Команда форума C\C++ Team

    Регистрация:
    6 сен 2011
    Сообщения:
    1.324
    Симпатии:
    36
    Алгоритм матричного метода решения СЛАУ подробно описан в теме: Матричный метод решения СЛАУ (распараллеливание с openMP).
    Распараллелим этот алгоритм с использованием библиотеки MPI (в отлчии от OpenMP она реализует параллелизм процессов, а не потоков и ориентирована на архитектуру ЭВМ с распределенной памятью (в т.ч. кластерного типа).

    Программа тестировалась на СЛАУ из 200 уравнений.

    В данной лабораторные использовались коллективные функции, такие как MPI_Scatter и MPI_Gather.

    Было принято решение распараллелить циклы:


    1) Деления элементов строки на элемент, лежащий на главной диагонали.


    Исходный код данного распараллеленного цикла представлен в листинге 1.1:
    Код (Text):
      MPI_Scatter(matrix[k], num_iter, MPI_DOUBLE, mini_segmentM, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E[k], num_iter, MPI_DOUBLE, mini_segmentE, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&div, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            for (int j = 0; j < num_iter; j++)
            {
                mini_segmentM[j] /= div;
                mini_segmentE[j] /= div;
            }

            MPI_Gather(mini_segmentM, num_iter, MPI_DOUBLE, matrix[k], num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(mini_segmentE, num_iter, MPI_DOUBLE, E[k], num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    Здесь, в первом MPI_Scatter каждому процессу отправляется часть строки матрицы A равной num_iter, в которой необходимо поделить элементы.

    Во втором MPI_Scatter каждому процессу отправляется часть строки единичной матрицы E равной num_iter, в которой необходимо поделить элементы.

    Сам элемент, который является делителем отправляется всем процессам с помощью MPI_Bcast.

    Далее с помощью функций MPI_Gather происходит “сбор” строки обратно в процессе 0.

    2) Зануление элементов ниже главной диагонали.

    Исходный код данного распараллеленного цикла представлен в листинге 1.2:

    Код (Text):
    if (rank == 0)
            {
                for (int i = 0; i < size; i++)
                {
                    segmentMK[i] = matrix[k][i];
                    segmentEK[i] = E[k][i];
                }
            }

            MPI_Bcast(segmentMK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(segmentEK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            MPI_Scatter(matrix, num_iter * size, MPI_DOUBLE, segmentM, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            for (int i = 0; i < num_iter; i++)
            {
                if ((rank * num_iter) + i <= k)
                    continue;

                multi = segmentM[i][k];
                for (int j = 0; j < size; j++)
                {
                    segmentM[i][j] -= multi * segmentMK[j];
                    segmentE[i][j] -= multi * segmentEK[j];
                }
            }
           
            MPI_Gather(segmentM, size * num_iter, MPI_DOUBLE, matrix, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(segmentE, size * num_iter, MPI_DOUBLE, E, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    Здесь в процессе с номером 0 в массивах segmentMK и segmentEK запоминаются строки, которые будут отниматься от строк, лежащих ниже данной (k-ой), матриц A и E соответственно.

    Таким образом, данные массивы отправляются всем процессам с помощью функций MPI_Bcast.

    С помощью функций MPI_Scatter каждый процесс получает свою часть матриц A и Е.

    В цикле с параметром i происходит обход строк полученных матриц. Если строка лежит выше и равна той, которая была отправлена функцией MPI_Bcast, то данная итерация цикла пропускается. Иначе от очередной строки отнимается строка, которая была отправлена функцией MPI_Bcast умноженная на элемент очередной строки, который находится ниже элемента, лежащего на главной диагонали. Таким образом, зануляются элементы, находящиеся ниже элемента, лежащего на главной диагонали, в этом же столбце.

    Далее с помощью функций MPI_Gather, в процессе 0 матрицы A и E собираются из частей, которые были изменены в процессах.

    3) Зануление элементов выше главной диагонали.


    Исходный код данного распараллеленного цикла представлен в листинге 1.3:

    Код (Text):
    for (int k = size - 1; k > 0; k--)
        {
            if (rank == 0)
            {
                for (int i = 0; i < size; i++)
                {
                    segmentMK[i] = matrix[k][i];
                    segmentEK[i] = E[k][i];
                }
            }

            MPI_Bcast(segmentMK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(segmentEK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(matrix, num_iter * size, MPI_DOUBLE, segmentM, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            for (int i = num_iter - 1; i > -1; i--)
            {
                if ((rank * num_iter) + i >= k)
                    continue;

                multi = segmentM[i][k];
                for (int j = 0; j < size; j++)
                {
                    segmentM[i][j] -= multi * segmentMK[j];
                    segmentE[i][j] -= multi * segmentEK[j];
                }
            }
           
            MPI_Gather(segmentM, size * num_iter, MPI_DOUBLE, matrix, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(segmentE, size * num_iter, MPI_DOUBLE, E, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        }

    p { text-indent: 1.25cm; margin-bottom: 0.25cm; direction: ltr; color: rgb(0, 0, 0); line-height: 120%; text-align: justify; orphans: 2; widows: 2; }p.western { font-family: "Times New Roman",serif; font-size: 14pt; }p.cjk { font-family: "Times New Roman"; font-size: 14pt; }p.ctl { font-family: "Times New Roman"; font-size: 14pt; }

    Здесь происходит “обратный ход”. То есть зануляются элементы выше главной диагонали. В результате прямого хода была получена верхнетреугольная матрица из матрицы A (элементы на главной диагонали равны 1, элементы ниже главной диагонали равны 0).

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

    Таким образом, функциями MPI_Bcast отправляются строки, которые будет отниматься от очередных строк. Функциями MPI_Scatter каждому процессу отправляется своя часть матриц A и E.

    После вычислений новая матрица собирается в процессе с номером 0 из частей, с помощью функции MPI_Gather.

    4) Вычисление неизвестных переменных X.

    Исходный код вычисления неизвестных переменных X представлен в листинге 1.4:
    Код (Text):
    MPI_Bcast(B, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        // Вычисление X
        for (int i = 0; i < num_iter; i++)
        {
            segmentX[i] = 0;
            for (int j = 0; j < size; j++)
                segmentX[i] += segmentE[i][j] * B[j];
        }
    MPI_Gather(segmentX, num_iter, MPI_DOUBLE, X, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    Здесь, с помощью функции MPI_Bcast каждому процессу отправляется матрица-столбец свободных членов B.

    Так как в результате элементарных преобразований матрица A трансформировалась в единичную, а матрица E трансформировалась в обратную матрицу A, то с помощью функции MPI_Scatter каждому процессу отправляется своя часть обратной матрицы.

    После того как каждый процесс вычислил свою часть неизвестных X, происходит “сбор” частей в одну часть в процессе 0, с помощью функции MPI_Gather.

    Тестирование программы:

    Для начала программа тестировалась на 5 процессах:

    upload_2016-10-4_13-52-45.png

    Рисунок 1 – результат выполнения программы на 5 процессах

    Позже программа была протестирована на 8 процессах (рис.2), 4х (рис.3) и 2х (рис.4):

    upload_2016-10-4_13-52-45.png

    Рисунок 2 – результат выполнения программы на 8 процессах

    upload_2016-10-4_13-52-45.png

    Рисунок 3 – результат выполнения программы на 4 процессах

    upload_2016-10-4_13-52-45.png

    Рисунок 4 – результат выполнения программы на 2 процессах

    Как видно из рисунков 1-4, быстрее всего программа работала на 2х процессах (скорее всего потому, что использовалась на двуядерном процессоре).

    Исходный код программы целиком:
    Код (Text):
    #include <iostream>
    #include <time.h>
    #include <vector>
    #include <mpi.h>

    using namespace std;

    double random(const int min, const int max)
    {
        if (min == max)
            return min;
        return min + rand() % (max - min);
    }

    int main(int argc, char* argv[])
    {
        setlocale(LC_ALL, "RUS");
        const int size = 200;

        int size_proc, rank;
        MPI_Status status;
        MPI_Request request;

        MPI_Init(&argc, &argv);
        MPI_Comm_size(MPI_COMM_WORLD, &size_proc);
        MPI_Comm_rank(MPI_COMM_WORLD, &rank);

        const int num_iter = 100;

        double matrix[size][size];
        double B[size];
        double E[size][size];

        double segmentM[num_iter][size];
        double segmentE[num_iter][size];

        double segmentMK[size];
        double segmentEK[size];

        double mini_segmentM[num_iter];
        double mini_segmentE[num_iter];

        //Заполнение матрицы A, B и E
        if (rank == 0)
        {
            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                {
                    matrix[i][j] = random(0, 100);
                    if (i == j) E[i][j] = 1.0;
                    else E[i][j] = 0.0;
                }
                B[i] = random(0, 100);
            }
        }
       
        // Прямой ход
        double t = clock();

        double div, multi;
        for (int k = 0; k < size; k++)
        {
            if (rank == 0)
            {
                if (matrix[k][k] == 0.0)
                {
                    bool changed = false;
                    for (int i = k + 1; i < size; i++)
                    {
                        if (matrix[i][k] != 0)
                        {
                            swap(matrix[k], matrix[i]);
                            swap(E[k], E[i]);
                            changed = true;
                            break;
                        }
                    }
                    if (!changed)
                    {
                        cout << endl << "Error: матрица не может быть найдена" << endl;
                        return -1;
                    }
                }
                div = matrix[k][k];
            }

            MPI_Scatter(matrix[k], num_iter, MPI_DOUBLE, mini_segmentM, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E[k], num_iter, MPI_DOUBLE, mini_segmentE, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&div, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            for (int j = 0; j < num_iter; j++)
            {
                mini_segmentM[j] /= div;
                mini_segmentE[j] /= div;
            }

            MPI_Gather(mini_segmentM, num_iter, MPI_DOUBLE, matrix[k], num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(mini_segmentE, num_iter, MPI_DOUBLE, E[k], num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            if (rank == 0)
            {
                for (int i = 0; i < size; i++)
                {
                    segmentMK[i] = matrix[k][i];
                    segmentEK[i] = E[k][i];
                }
            }

            MPI_Bcast(segmentMK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(segmentEK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            MPI_Scatter(matrix, num_iter * size, MPI_DOUBLE, segmentM, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            for (int i = 0; i < num_iter; i++)
            {
                if ((rank * num_iter) + i <= k)
                    continue;

                multi = segmentM[i][k];
                for (int j = 0; j < size; j++)
                {
                    segmentM[i][j] -= multi * segmentMK[j];
                    segmentE[i][j] -= multi * segmentEK[j];
                }
            }
           
            MPI_Gather(segmentM, size * num_iter, MPI_DOUBLE, matrix, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(segmentE, size * num_iter, MPI_DOUBLE, E, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        }

        //Обратный ход
        for (int k = size - 1; k > 0; k--)
        {
            if (rank == 0)
            {
                for (int i = 0; i < size; i++)
                {
                    segmentMK[i] = matrix[k][i];
                    segmentEK[i] = E[k][i];
                }
            }

            MPI_Bcast(segmentMK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(segmentEK, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(matrix, num_iter * size, MPI_DOUBLE, segmentM, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            for (int i = num_iter - 1; i > -1; i--)
            {
                if ((rank * num_iter) + i >= k)
                    continue;

                multi = segmentM[i][k];
                for (int j = 0; j < size; j++)
                {
                    segmentM[i][j] -= multi * segmentMK[j];
                    segmentE[i][j] -= multi * segmentEK[j];
                }
            }
           
            MPI_Gather(segmentM, size * num_iter, MPI_DOUBLE, matrix, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Gather(segmentE, size * num_iter, MPI_DOUBLE, E, size * num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        }
        /*if (rank == 0)
        {
            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                    cout << E[i][j] << ' ';

                cout << '\n';
            }
        }*/

        double X[size];
        double segmentX[num_iter];

        MPI_Bcast(B, size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        MPI_Scatter(E, num_iter * size, MPI_DOUBLE, segmentE, num_iter * size, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        // Вычисление X
        for (int i = 0; i < num_iter; i++)
        {
            segmentX[i] = 0;
            for (int j = 0; j < size; j++)
                segmentX[i] += segmentE[i][j] * B[j];
        }

        MPI_Gather(segmentX, num_iter, MPI_DOUBLE, X, num_iter, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        /*if (rank == 0)
        {
            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                    cout << matrix[i][j] << " ";
                cout << endl;
            }
        }*/

        if (rank == 0)
        {
            cout << "\nThe system of equations:";
            for (int i = 0; i < size; i++)
                cout << "\nx" << i + 1 << " = " << X[i];

            t = (clock() - t) / 1000;
            cout << "\n\nThe time spent on computation: " << t << "s.";
        }
        MPI_Finalize();
        return 0;
    }
     
    WebWare Team нравится это.
Загрузка...
Похожие Темы - Матричный метод решения
  1. rrrFer
    Ответов:
    0
    Просмотров:
    307
  2. DeuS7
    Ответов:
    0
    Просмотров:
    92
  3. Sander
    Ответов:
    1
    Просмотров:
    510
  4. kuzduk
    Ответов:
    0
    Просмотров:
    350
  5. kir2700

    Заплачу Метод крылова c#

    kir2700, 20 дек 2015, в разделе: C/C++/C#
    Ответов:
    0
    Просмотров:
    793

Поделиться этой страницей