Блог

Я новичок в OpenMP и отчаянно пытаюсь учиться. Я попытался написать пример кода на C в Visual Studio 2012 для реализации умножения матриц. Я надеялся, что кто-нибудь с опытом работы в OpenMP может взглянуть на этот код и помочь мне получить максимальную скорость / распараллеливание для этого:

 #include <iostream>
#include <stdlib.h>
#include <omp.h>
#include <random>
using namespace std;

#define NUM_THREADS 4

// Program Variables
double**        A;
double**        B;
double**        C;
double          t_Start;
double          t_Stop;
int             Am;
int             An;
int             Bm;
int             Bn;

// Program Functions
void            Get_Matrix();
void            Mat_Mult_Serial();
void            Mat_Mult_Parallel();
void            Delete_Matrix();


int main()
{
    printf("Matrix Multiplication Programnn");
    cout << "Enter Size of Matrix A: ";
    cin >> Am >> An;
    cout << "Enter Size of Matrix B: ";
    cin >> Bm >> Bn;

    Get_Matrix();
    Mat_Mult_Serial();
    Mat_Mult_Parallel();


    system("pause");
    return 0;

}


void Get_Matrix()
{
    A = new double*[Am];
    B = new double*[Bm];
    C = new double*[Am];
    for ( int i=0; i<Am; i   ){A[i] = new double[An];}
    for ( int i=0; i<Bm; i   ){B[i] = new double[Bn];}
    for ( int i=0; i<Am; i   ){C[i] = new double[Bn]; }

    for ( int i=0; i<Am; i   )
    {
         for ( int j=0; j<An; j   )
         {
             A[i][j]= rand() % 10   1;
         }
    }

    for ( int i=0; i<Bm; i   )
    {
        for ( int j=0; j<Bn; j   )
        {
            B[i][j]= rand() % 10   1;
        }
    }
    printf("Matrix Create Complete.n");
}


void Mat_Mult_Serial()
{
    t_Start = omp_get_wtime();
    for ( int i=0; i<Am; i   )
    {
        for ( int j=0; j<Bn; j   )
        {
            double temp = 0;
            for ( int k=0; k<An; k   )
            {
                temp  = A[i][k]*B[k][j];
            }
        }
    }
    t_Stop = omp_get_wtime() - t_Start;
    cout << "Serial Multiplication Time: " << t_Stop << " seconds" << endl;
    }


void Mat_Mult_Parallel()
{
    int i,j,k;
    t_Start = omp_get_wtime();

    omp_set_num_threads(NUM_THREADS);
    #pragma omp parallel for private(i,j,k) schedule(dynamic)
    for ( i=0; i<Am; i   )
    {
        for ( j=0; j<Bn; j   )
        {
            //double temp = 0;
            for ( k=0; k<An; k   )
            {
                C[i][j]  = A[i][k]*B[k][j];
            }
        }
    }

    t_Stop = omp_get_wtime() - t_Start;
    cout << "Parallel Multiplication Time: " << t_Stop << " seconds." << endl;
}


void Delete_Matrix()
{
    for ( int i=0; i<Am; i   ){ delete [] A[i]; }
    for ( int i=0; i<Bm; i   ){ delete [] B[i]; }
    for ( int i=0; i<Am; i   ){ delete [] C[i]; }

    delete [] A;
    delete [] B;
    delete [] B;
}
  

Мои примеры основаны на классе matrix, который я создал для параллельного обучения. Если вы заинтересованы, не стесняйтесь обращаться ко мне. Существует несколько способов ускорить умножение матриц :

Хранение

Используйте одномерный массив в порядке следования строк для более быстрого доступа к элементу.
Вы можете получить доступ к A(i,j) с помощью [i * An j]

Используйте циклическую инвариантную оптимизацию

 for (int i = 0; i < m; i   )
    for (int j = 0; j < p; j   )
    {
        Scalar sigma = C(i, j);
        for (int k = 0; k < n; k   )
            sigma  = (*this)(i, k) * B(k, j);
        C(i, j) = sigma;
    }
  

Это предотвращает повторное вычисление C (i, j) несколько раз в самом внутреннем цикле.

Измените порядок цикла «для k <-> для i»

 for (int i = 0; i < m; i   )
    for (int k = 0; k < n; k   )
    {
        Aik = (*this)(i, k);
        for (int j = 0; j < p; j   )
            C(i, j)  = Aik * B(k, j);
    }
  

Это позволяет играть с локальностью пространственных данных

Используйте блокировку / разбиение цикла

 for(int ii = 0; ii < m; ii  = block_size)
    for(int jj = 0; jj < p; jj  = block_size)
        for(int kk = 0; kk < n; kk  = block_size)
            #pragma omp parallel for // I think this is the best place for this case
            for(int i = ii; i < ii   block_size; i   )
                for(int k = kk; k < kk   block_size; k   )
                {
                    Scalar Aik = (*this)(i, k);
                    for(int j = jj; j < jj   block_size; j   )
                        C(i, j)  =  Aik * B(k, j);
                }
  

Это может использовать лучшую временную локальность данных. Оптимальный размер блока зависит от вашей архитектуры и размера матрицы.

Затем распараллеливайте!

Как правило, параллельный #pragma omp для должен выполняться в самом внешнем цикле. Возможно, использование двух параллельных циклов в двух первых внешних циклах может дать лучшие результаты. Тогда это зависит от используемой вами архитектуры, размера матрицы… Вы должны протестировать! Поскольку умножение матриц имеет статическую рабочую нагрузку, я бы использовал статическое расписание.

Оптимизация Moar !

Вы можете выполнить оптимизацию вложенности цикла. Вы можете векторизовать свой код. Вы можете посмотреть, как это делают BLAS.

Я очень новичок в OpenMP, и этот код очень поучителен. Однако я обнаружил ошибку в последовательной версии, которая дает ей несправедливое преимущество в скорости по сравнению с параллельной версией.

Вместо того, чтобы писать C[i][j] = A[i][k]*B[k][j]; , как вы делаете в параллельной версии, вы написали temp = A[i][k]*B[k][j]; в последовательной версии. Это намного быстрее (но не помогает вам вычислить матрицу C). Таким образом, вы не сравниваете яблоки с яблоками, что делает параллельный код более медленным по сравнению. Когда я исправил эту строку и запустил ее на своем ноутбуке (что позволяет использовать 2 потока), параллельная версия была почти в два раза быстрее. Неплохо!