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

rrrFer · 04.10.2016

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

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

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

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

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

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

Код:

  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:

Код:

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:

Код:

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:

Код:

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 процессах:

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

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

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

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

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

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

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

Код:

#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;
}

sfdsfgadgsdfghsdfh · 05.01.2018

День добрый. Компилирую так "mpicxx -o prog mpi.cpp" запускаю так "mpirun -np 5 ./prog" в итоге выводит это "
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= EXIT CODE: 139
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES"
Как решить проблему?
[doublepost=1515153742,1515152654][/doublepost]Решил проблему указав 6 уравнений (const int size = 6

и (const int num_iter = 3

, но в итоге я не понял, какую роль несёт параметр num_iter, то есть как правильно подбирать под количество уравнений

Все сервисы Codeby

Поиск

Поиск

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

rrrFer

sfdsfgadgsdfghsdfh