Распараллеливание алгоритма Гаусса-Жордана с MPI

      Комментарии к записи Распараллеливание алгоритма Гаусса-Жордана с MPI отключены

Главная Форумы Программирование Параллельное программирование Использование MPI Распараллеливание алгоритма Гаусса-Жордана с MPI

В этой теме 0 ответов, 1 участник, последнее обновление  Васильев Владимир Сергеевич 4 мес. назад.

  • Автор
    Сообщения
  • #3090

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

    В алгоритме Гаусса-Жордана удобно распараллеливать обработку элементарных преобразований над матрицей (домножение строки на константу и ее сложение с другой строкой матрицы).

    Приведем исходный код, а затем рассмотрим наиболее интересные части решения:

    #include "mpi.h"
    #include <cstdlib>
    #include <fstream>
    #include <vector>
    
    using namespace std;
    
    int main(int argc, char* argv[])
    {
        double t_start, t_end; // time measure
        int n = 1000; // size of matrix
        double* A = new double[n*(n + 1)]; // matrix
        double* X = new double[n]; // variables
        double buf;
        double* rbuf = new double[n*(n + 1)];
        double* row_buf = new double[n + 1];
        int i, j, k;
        ofstream fout;
    
        int *pSendNum;    // Количество элементов, посылаемых процессу
        int *pSendInd;    // Индекс первого элемента данных, посылаемого процессу
    
        int procs_rank, procs_count;
    
        MPI_Init(&argc, &argv);
    
        MPI_Comm_size(MPI_COMM_WORLD, &procs_count);
        MPI_Comm_rank(MPI_COMM_WORLD, &procs_rank);
        
    
        // Init
        if (procs_rank == 0)
        {
            fout.open("output.txt");
            t_start = MPI_Wtime();
    
            const double zeroEps = 0.0000000001;
            // generate matrix
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n + 1; j++)
                    A[i*(n + 1) + j] =  -500 + rand() % 1001; //i*(n + 1) + j; //
                if (abs(A[i*(n + 1) + i]) < zeroEps)
                    A[i*(n + 1) + i] = 1;
            }
        }
    
        // Gauss
        for (k = 0; k < n; k++)
        {
            if (procs_rank == 0)
                for (j = 0; j < n + 1; j++)
                    row_buf[j] = A[k*(n + 1)+j];
            MPI_Bcast(row_buf, n+1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
            pSendInd = new int[procs_count];
            pSendNum = new int[procs_count];
            int RowsLeft = n - k - 1;
            int Rows = (n - k - 1) / procs_count;
            pSendNum[0] = Rows*(n + 1);
            pSendInd[0] = 0;
            for (j = 1; j < procs_count; j++)
            {
                RowsLeft -= Rows;
                Rows = RowsLeft / (procs_count - j);
                pSendNum[j] = Rows*(n+1);
                pSendInd[j] = pSendInd[j - 1] + pSendNum[j - 1];
            }
    
            MPI_Scatterv(&A[(k + 1)*(n + 1)], pSendNum, pSendInd, MPI_DOUBLE, rbuf,
                pSendNum[procs_rank], MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
            for (i = 0; i < pSendNum[procs_rank]/(n+1); i++)
            {
                buf = rbuf[i*(n + 1) + k] / row_buf[k];
                for (j = k; j < n + 1; j++)
                    rbuf[j] = rbuf[j] - buf*row_buf[j];
            }
    
            MPI_Gatherv(rbuf, pSendNum[procs_rank], MPI_DOUBLE, 
                &A[(k + 1)*(n + 1)], pSendNum, pSendInd, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
            delete[] pSendInd, pSendNum;
        }
        
        if (procs_rank == 0)
        {
            X[n - 1] = A[(n - 1)*(n + 1) + n] / A[(n - 1)*(n + 1) + n - 1];
            for (i = n - 2; i >= 0; i--)
            {
                buf = 0;
                for (j = i + 1; j < n; j++)
                    buf = buf + A[i*(n + 1) + j] * X[j];
                X[i] = 1.0 / A[i*(n + 1) + i] * (A[i*(n + 1) + n] - buf);
            }
            t_end = MPI_Wtime();
        }
    
        // Final
        if (procs_rank == 0)
        {
            fout << "Time elapsed: " << t_end - t_start << '\n';
            for (i = 0; i < n; i++)
                fout << "X" << i + 1 << ": " << X[i] << '\n';
            fout.close();
            delete[] A, X;
        }
    
        MPI_Finalize();
    
        return 0;
    }

    Все процессы в MPI выполняют один и тот же код, для того, чтобы разделить между ними работу нужно использовать ранг (порядковый номер) процесса как показано в этом фрагменте:

    if (procs_rank == 0)
      for (j = 0; j < n + 1; j++)
        row_buf[j] = A[k*(n + 1)+j];
    MPI_Bcast(row_buf, n+1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    Здесь происходит копирование нулевым процессом k-ой строки в вектор row_buf и его рассылка другим процессам. При выполнении MPI_Bcast рассылка будет выполнена нулевым процессом (номер отправителя задается в аргументах), а остальные процессы выполнят прием данных.
    int MPI_Bcast – функция широковещательной рассылки, где процесс с номером root рассылает из своего буфера передачи всем процессам области связи коммуникатора comm.
    Параметры функции:

    1. void* buffer – адрес начала расположения в памяти рассылаемых данных;
    2. int count – число посылаемых элементов;
    3. MPI_Datatype datatype – тип посылаемых элементов;
    4. int root – номер процесса-отправителя;
    5. MPI_Comm comm – коммуникатор.

    pSendInd = new int[procs_count];
    pSendNum = new int[procs_count];
    int RowsLeft = n - k - 1;
    int Rows = (n - k - 1) / procs_count;
    pSendNum[0] = Rows*(n + 1);
    pSendInd[0] = 0;
    for (j = 1; j < procs_count; j++)
    {
      RowsLeft -= Rows;
      Rows = RowsLeft / (procs_count - j);
      pSendNum[j] = Rows*(n+1);
      pSendInd[j] = pSendInd[j - 1] + pSendNum[j - 1];
    }
    
    MPI_Scatterv(&A[(k + 1)*(n + 1)], pSendNum, pSendInd, MPI_DOUBLE, rbuf,
      pSendNum[procs_rank], MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
    for (i = 0; i < pSendNum[procs_rank]/(n+1); i++)
    {
      buf = rbuf[i*(n + 1) + k] / row_buf[k];
      for (j = k; j < n + 1; j++)
        rbuf[j] = rbuf[j] - buf*row_buf[j];
    }
    
    MPI_Gatherv(rbuf, pSendNum[procs_rank], MPI_DOUBLE, 
      &A[(k + 1)*(n + 1)], pSendNum, pSendInd, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
    delete[] pSendInd, pSendNum;

    В данном участке программы происходит разбиение матрицы на строки и их отправка процессам.
    MPI_Scatterv и MPI_Gatherv используются для передачи строк матрицы другим процессам и получения от них новых строк. Для того, чтобы распределить элементы матрицы между процессами, нам нужно для каждого набора строк знать индексы первого элемента в матрице, а также суммарное количество элементов в каждом наборе.
    Массив pSendInd содержит индексы элементов относительно первого элемента k+1 – строки; массив pSendNum – число элементов i-наборе.
    Необходимость использования векторного варианта Scatter и Gather продиктована изменяющимся числом строк в каждой итерации. Если при разделении на строки окажется, что число строк не кратно числу процессов, то разделение произойдет некорректно: число передаваемых элементов будет не кратно числу столбцов.
    В цикле for происходит расчет коэффициента buf, необходимого для получения нуля в k-ом столбце i-ой строки, а также вычитание k-ой строки, умноженной на buf из i-ой строки.

Для ответа в этой теме необходимо авторизоваться.