Опубликован: 30.05.2014 | Уровень: для всех | Доступ: платный | ВУЗ: Нижегородский государственный университет им. Н.И.Лобачевского

Самостоятельная работа 5: Оптимизация вычислений в задаче матричного умножения. Оптимизация работы с памятью

Использование блоков прямоугольной формы

Разработаем еще одну реализацию блочного умножения матриц, в которой в качестве блока матрицы B берется не квадратная область, а горизонтальная полоса. Блок матрицы A остается квадратным. Также сохраним предположение о том, что порядок матрицы делится на размер блока без остатка.

//singleBlock_t.cpp
#include "mult.h"
#include "assert.h"
#include <omp.h>

void mult(ELEMENT_TYPE * A, ELEMENT_TYPE * B, 
          ELEMENT_TYPE * C, int n, int bSize)
{
  ELEMENT_TYPE s, err;
  
  assert(n % bSize == 0);
  
  for(int j = 0; j < n; j++ )
    for(int i = 0; i < n; i++ )  
      C[j * n + i] = 0;
  
  for(int jk = 0; jk < n / bSize; jk++ )
    for(int ik = 0; ik < n / bSize; ik++ )
      for(int j = jk * bSize; j < jk * bSize + bSize; j++)
        for(int k = ik * bSize; k < ik * bSize + bSize;k++)
#pragma simd
          for(int i = 0; i < n; i++ )
            C[j * n + i] += A[j * n + k] * B[k * n + i];
}

Данный вариант обладает рядом преимуществ. Во-первых, элементы блока матрицы A повторно используются несколько раз, как и при реализации с квадратными блоками. Как следствие, часть данных читается не из памяти, а из кэш-памяти. Во-вторых, если матрица сравнительно небольшого размера, то полоса матрицы B и C может поместиться в кэш-память последнего уровня, тем самым повышая эффективность использования памяти. В-третьих, внутренний цикл содержит много итераций, что, по всей видимости, лучше соответствует эвристикам компилятора. Подробнее о преимуществах данного подхода можно прочитать в работах [12.9,12.10,12.11].

Попробуем скомпилировать и запустить полученный код. На рис. 10.10 представлен результат запуска кода с измененным размером блока матрицы B.

Результат запуска алгоритма с блоком не квадратной формы

Рис. 10.10. Результат запуска алгоритма с блоком не квадратной формы

В таблице 10.9 представлено сравнение времени выполнения последовательного блочного алгоритма с разным размером блока и Intel MKL. Как и в случае с блочным алгоритмом для размера 1024 в начале размер блока задавался через параметр, а затем была использована константа. Для матриц размера 2048 и 3073 размер блока задавался константой.

Таблица 10.9. Сравнение времени работы блочной реализации алгоритма и Intel MKL на Intel Xeon Phi (время в секундах).
Размер блока: 16 32 64 128 jki MKL seq.
N=1024\ Размер блока\параметр 1,29 1,28 1,26 1,51 1,95 0,207
N=1024\ Размер блока\константа 1,12 1,11 1,09 1,32 1,95 0,207
N=2048\ Размер блока\ константа 8,53 8,34 15,05 15,05 15,4 1,363
N=3072\ Размер блока\ константа 28,97 28,07 49,4 50,6 50,8 4,411

Время работы алгоритма в зависимости от размера блока (N = 1024) на Intel Xeon Phi.

Рис. 10.11. Время работы алгоритма в зависимости от размера блока (N = 1024) на Intel Xeon Phi.

Новая блочная реализация алгоритма позволила получить лучшие результаты, чем версия с квадратными блоками.

Анализируя таблицу можно заметить следующий факт – при вычислениях с матрицами большого размера происходить резкое замедление при использовании размера блока, равного 64. По всей видимости, при переходе с размера блока равного 32 на 64 полоса матрицы B начинает не помещаться в кеш-память. Для решения данной проблемы на матрицах еще большего размера необходимо в качестве блока матрицы B взять только часть полосы, т.е. использовать прямоугольные блоки. В качестве дополнительного задания предлагается реализовать блочное умножение матриц с прямоугольными блоками и подобрать оптимальный размер блока.

Параллельная реализация

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

#include "mult.h"
#include "assert.h"
#include <omp.h>

void mult(ELEMENT_TYPE * A, ELEMENT_TYPE * B, 
          ELEMENT_TYPE * C, int n, int bSize)
{
  ELEMENT_TYPE s, err;
  
  assert(n % bSize == 0);
  
  for(int j = 0; j < n; j++ )
    for(int i = 0; i < n; i++ )  
      C[j * n + i] = 0;

#pragma omp parallel for   
  for(int jk = 0; jk < n / bSize; jk++ )
    for(int ik = 0; ik < n / bSize; ik++ )
      for(int j = jk * bSize; j < jk * bSize + bSize; j++)
        for(int k = ik * bSize; k < ik * bSize + bSize;k++)
#pragma simd
          for(int i = 0; i < n; i++ )
            C[j * n + i] += A[j * n + k] * B[k * n + i];
}

Скомпилируем и выполним код. На рис. 10.12 представлен результат вычислений.

Результат выполнения  параллельного алгоритма

увеличить изображение
Рис. 10.12. Результат выполнения параллельного алгоритма

В таблице 10.10 приведено сравнение времени выполнения параллельного алгоритма с Intel MKL.

Таблица 10.10. Сравнение времени работы параллельной блочной реализации алгоритма и Intel MKL на Intel Xeon Phi. (время в секундах)
Размер блока: 32 64 MKL parallel
N=1024 0,040 0,070 0,004
N=2048 0,250 0,370 0,019
N=3072 0,760 1,000 0,062

Из таблицы 10.10 видно, что время вычислений отличается от Intel MKL на порядок.

Можно ли еще ускорить работу реализованного алгоритма?

Умножение матриц очень интенсивно использует память. При большом количестве потоков существенно возрастает нагрузка на шину данных. Для уменьшения нагрузки можно попробовать сократить количество потоков, используя функцию OpenMP – omp_set_num_threads. Сократим количество потоков до 120:

#include "mult.h"
#include "assert.h"
#include <omp.h>

void mult(ELEMENT_TYPE * A, ELEMENT_TYPE * B, 
          ELEMENT_TYPE * C, int n, int bSize)
{
  ELEMENT_TYPE s, err;
  
  assert(n % bSize == 0);
  
  for(int j = 0; j < n; j++ )
    for(int i = 0; i < n; i++ )  
      C[j * n + i] = 0;

omp_set_num_threads(120);

#pragma omp parallel for   
  for(int jk = 0; jk < n / bSize; jk++ )
    for(int ik = 0; ik < n / bSize; ik++ )
      for(int j = jk * bSize; j < jk * bSize + bSize; j++)
        for(int k = ik * bSize; k < ik * bSize + bSize;k++)
#pragma simd
          for(int i = 0; i < n; i++ )
            C[j * n + i] += A[j * n + k] * B[k * n + i];
}

В таблице 10.11 приведено сравнение времени выполнения параллельного алгоритма в зависимости от размера блока с Intel MKL при 120 потоках OpenMP.

Таблица 10.11. Сравнение времени работы параллельной блочной реализации алгоритма и Intel MKL на Intel Xeon Phi. (время в секундах)
Размер блока: 244 потока 120 потоков MKL parallel
32 64 32 64
1024 0,04 0,07 0,04 0,08 0,004
2048 0,25 0,37 0,16 0,35 0,019
3048 0,76 1 0,44 0,94 0,062

Как видно из таблицы 10.11, отставание от Intel MKL на больших матрицах несколько сократилось. В качестве дополнительного задания предлагается подобрать оптимальное количество потоков.

В целом работа над кодом может быть продолжена. Желающие достигнуть результаты, показываемые Intel MKL, могут обратиться к соответствующей литературе (см., например, 12.10). При этом, по всей вероятности, обойтись исключительно использованием языка C не удастся, потребуется программировать на ассемблере (или в интринсиках).

Дополнительные задания

  1. Подберите оптимальные размеры блоков для разных версий блочной реализации алгоритма умножения матриц (с квадратным блоком, с прямоугольным блоком).
  2. Подберите оптимальное количество создаваемых OpenMP потоков для реализованного блочного алгоритма умножения матриц в зависимости от размера матриц и блоков.
  3. Исследуйте влияние прогрева на время работы созданных программ.
Svetlana Svetlana
Svetlana Svetlana

Здравствуйие! Я хочу пройти курс Введение в принципы функционирования и применения современных мультиядерных архитектур (на примере Intel Xeon Phi), в презентации самостоятельной работы №1 указаны логин и пароль для доступ на кластер и выполнения самостоятельных работ, но войти по такой паре логин-пароль не получается. Как предполагается выполнение самосоятельных работ в этом курсе?

Егор Кузьмин
Егор Кузьмин
Россия, г. Москва
Вера Борисова
Вера Борисова
Россия