Параллельное умножение матриц с помощью OpenMP

Доброго времени суток всем. Недавно я опубликовал коротенькую инструкцию по настройке OpenMP в CLion, а сейчас пришло время для реализаций параллельных алгоритмов. Первым на очереди стоит простенький алгоритм перемножения матриц, он очень хорошо подходит для распараллеливания, потому что состоит из трех вложенных циклов.

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

Реализация алгоритма умножения матриц на C/C++

Пусть матрица хранится в двумерном массиве int **matrix, и доступ к элементам осуществляется двойным индексом matrix[i][j]. Для начала произведем простенькую проверку на то, что матрицы согласованы, после этого можно выделить память и выполнить умножение по формуле.

//Производит умножение матрицы размером n1 x m1
//на матрицу размером n2 x m2
int** matrixMulti(int **matrix1, int n1, int m1, int **matrix2, int n2, int m2) {
    //Если матрицы не согласованы, то не выполняем умножение
    if(m1 != n2) {
        cout << "Error! m1 != n2" << endl;
        return NULL;
    }

    //Выделяем память под результат умножения
    int **result;
    result = (int**)malloc(sizeof(int*)*n1);
    for(int i = 0; i < n1; i++) {
        result[i] = (int*)malloc(sizeof(int)*m2);
    }

    //Умножение по формуле
    for(int i = 0; i < n1; i++) {
        for(int j = 0; j < m2; j++) {
            result[i][j] = 0;
            for(int k = 0; k < m1; k++) {
                result[i][j] += (matrix1[i][k] * matrix2[k][j]);
            }
        }
    }

    return result;
}

Реализация параллельного умножения матриц на C/C++

Нам не придется писать большое количество кода благодаря тому, что мы будем оперировать потоками с общей памятью, а не процессами, как в MPI. Вся суть распараллеливания состоит в том, чтобы дописать директиву #pragma omp parallel for перед внешним циклом. После этого нужное нам количество потоков автоматически разберет вычислительные задачи на себя.

Прилагаю исходный код с функцией случайной генерации матрицы и непосредственно самим параллельным алгоритмом их умножения.

#include <iostream>
#include <cstdlib>
#include <omp.h>

using namespace std;

void randomiseMatrix(int **matrix, int n, int m) {
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < m; j++) {
            matrix[i][j] = rand() % 11;
        }
    }

    return;
}

int main(int argc, char** argv) {
    srand(time(NULL));
    int n1 = 1000;
    int m1 = 500;
    int n2 = 500;
    int m2 = 1200;

    //Матрица n1 x m1
    int **matrix1;
    //Матрица n2 x m2
    int **matrix2;

    matrix1 = (int**)malloc(sizeof(int*)*n1);
    for(int i = 0; i < n1; i++) {
        matrix1[i] = (int*)malloc(sizeof(int)*m1);
    }
    matrix2 = (int**)malloc(sizeof(int*)*n2);
    for(int i = 0; i < n2; i++) {
        matrix2[i] = (int*)malloc(sizeof(int)*m2);
    }

    //Генерируем случайные матрицы для умножения
    randomiseMatrix(matrix1, n1, m1);
    randomiseMatrix(matrix2, n2, m2);

    int **result = (int**)malloc(sizeof(int*)*n1);;
    for(int i = 0; i < n1; i++) {
        result[i] = (int*)malloc(sizeof(int)*m2);
    }

    //Устанавливаем число потоков
    int threadsNum = 2;
    omp_set_num_threads(threadsNum);
    int i, j, k;
#pragma omp parallel for shared(matrix1, matrix2, result) private(i, j, k)
    for (i = 0; i < n1; i++) {
        for (j = 0; j < m2; j++) {
            result[i][j] = 0;
            for (k = 0; k < m1; k++) {
                result[i][j] += (matrix1[i][k] * matrix2[k][j]);
            }
        }
    }

    return 0;
}

Заключение

Кажется, что ничего и не сделали, но если внимательно проследить за диспетчером(например htop) задач во время исполнения, становится видно, как нагрузка распределяется по ядрам процессора. В скором времени опубликую еще реализации параллельных алгоритмов с помощью OpenMP, подписывайтесь на обновления по электронной почте. А на сегодня у меня все, спасибо за внимание!