Перейти к основному содержимому

Параллельное умножение матриц

ОБЯЗАТЕЛЬНОДЛЯ НОВИЧКОВ
Разработчику Инженеру

Пример — от формулы к параллельной схеме​

C = A × B — умножение матриц. Элемент C[i,j] получается скалярным произведением строки i матрицы A и столбца j матрицы B. Эта операция лежит в основе линейной алгебры, нейросетей и многих физических расчётов.

Алгоритм иллюстрирует параллелизм по данным, разбиение по узлам и цену обмена данными между процессами. Готовые библиотеки (BLAS, LAPACK, cuBLAS, MKL) — результат многолетней оптимизации; ниже — идеи, которые в них заложены, в виде псевдокода и с пояснением эталонных записей на C++/CUDA.

Интуиция размера​

Для квадратных матриц n×n тройной цикл выполняет порядка n³ умножений и хранит n² элементов результата. При n = 10 000 это порядка 10¹² операций — на одном ядре часы и дни. Кластер или GPU делают задачу реалистичной, если данные разрезаны так, чтобы не пересылать всю матрицу B на каждом шаге.

Перед matmul полезно разобрать более простые построения — они же входят в prefix sum на PRAM и в BLAS.

Сумма и префиксные суммы​

Сумма n чисел​

Последовательно — O(n). На p процессорах — дерево редукции за O(log p) раундов обмена плюс локальные куски по n/p элементов.

АЛГОРИТМ СУММА_ПАРАЛЛЕЛЬНО(n, a, p, rank)
  локальная := сумма элементов a на своём отрезке индексов
  частичная[rank] := локальная
  для уровень от 1 до ⌈log2 p⌉
    если rank кратен 2^уровень
      частичная[rank] := частичная[rank] + частичная[rank + 2^(уровень−1)]
    синхронизация
  конец для
  итог := частичная[0]   // у rank 0
КОНЕЦ

Все частичные суммы (scan)​

Задача: s[i] = a[0] + … + a[i]. Нужна в спarse-алгебре, сортировках, построении дерева. На PRAM — O(log n) раундов; в MPI — MPI_Exscan / MPI_Scan.

Умножение матрицы на вектор​

y = A·x — проще matmul: меньше зависимостей, удобный учебный шаг к декомпозиции.

СхемаИдеяОбмены
По строкамПроцесс r считает y[i] для своих строк iВектор x — broadcast или реплика
По столбцамКаждый владеет столбцами A, локально накапливает вклад в y[i]Редукция по y
БлочноПодматрица A и кусок x на узлеHalo по границе блока
АЛГОРИТМ MATVEC_ПО_СТРОКАМ(n, A, x, y, p, rank)
  (i_нач, i_кон) := свой_диапазон_строк(rank, p, n)
  x_лок := получить_копию_вектора(x)    // broadcast
  для i от i_нач до i_кон
    y[i] := 0
    для j от 0 до n − 1
      y[i] := y[i] + A[i,j] * x_лок[j]
    конец для
  конец для
КОНЕЦ

На кластере строки — типичный первый выбор: один MPI_Bcast для x, затем независимые циклы.

Последовательный алгоритм​

Псевдокод (тройной цикл)​

АЛГОРИТМ УМНОЖИТЬ_МАТРИЦЫ(n, A, B, C)
  для i от 0 до n − 1
    для j от 0 до n − 1
      C[i,j] := 0
      для k от 0 до n − 1
        C[i,j] := C[i,j] + A[i,k] * B[k,j]
      конец для
    конец для
  конец для
КОНЕЦ

Построчный разбор​

ЦиклРоль
i, jВыбираем позицию в результате C
C[i,j] := 0Обнуляем накопитель перед суммированием
kПеребираем слагаемые произведения A[i,k]·B[k,j] — "скалярное произведение" строки и столбца

Сложность — O(n³) операций, O(n²) памяти для матриц.

Порядок циклов влияет на локальность в кэше. Вариант внешние i и k, внутренний j часто лучше для доступа к строке A и строке B (в памяти матрицы хранятся по строкам) — это тема оптимизации на этапе языка, здесь достаточно знать, что один и тот же математический алгоритм можно переставить для скорости.

Элементный параллелизм (независимые элементы C)​

Каждый C[i,j] независим при фиксированных A и B — разные пары (i,j) можно считать параллельно (чтения A и B общие, записи в разные ячейки C).

Псевдокод​

АЛГОРИТМ УМНОЖИТЬ_ПАРАЛЛЕЛЬНО_ПО_ВЫХОДУ(n, A, B, C)
  параллельно для i от 0 до n − 1
    параллельно для j от 0 до n − 1
      сумма := 0
      для k от 0 до n − 1
        сумма := сумма + A[i,k] * B[k,j]
      конец для
      C[i,j] := сумма
    конец параллельно
  конец параллельно
КОНЕЦ
ИдеяСледствие
Два уровня параллельно дляДо n² независимых задач
Внутренний kОстаётся последовательным внутри одной пары (i,j)

Плюсы — простота, огромный параллелизм на GPU.

Минусы — обход B по столбцам плох для кэша; в MPI матрица B часто должна быть доступна на каждом узле (replicate или broadcast).

Справочно на C++ (OpenMP)​

Код ITЗагрузка примера кода…

ЭлементСмысл
collapse(2)Объединить два цикла i и j в одно пространство итераций для распределения между потоками
double sumЛокальная переменная на итерацию (i,j) — у каждой задачи своя sum
A[i,k]В C++ индексация с нуля; запятая в учебниках иногда пишут как A[i][k]

Когда хватает — GPU, OpenMP на средних n, учебные задачи.

Блочное умножение (учёт кэша)​

Разбиваем матрицы на блоки размера b×b, чтобы маленькое умножение блоков помещалось в кэш L1/L2.

Псевдокод​

АЛГОРИТМ УМНОЖИТЬ_БЛОКАМИ(n, b, A, B, C)
  для i_блок от 0 до n − 1 шаг b
    для j_блок от 0 до n − 1 шаг b
      обнулить подматрицу C_блок[i_блок, j_блок]
      для k_блок от 0 до n − 1 шаг b
        C_блок += УМНОЖИТЬ_МАЛЕНЬКИЕ_БЛОКИ(
          A[i_блок..i_блок+b), k_блок..),
          B[k_блок.., j_блок..j_блок+b)
        )
      конец для
    конец для
  конец для
КОНЕЦ

Параллелизм — независимые пары (i_блок, j_блок) при фиксированном проходе k_блок, либо параллель по k_блок с аккуратным накоплением в C.

На одном узле так работают GotoBLAS/OpenBLAS.

1D row decomposition (MPI)​

P процессов. Процесс rank владеет полосой строк A и соответствующими строками C; матрица B реплицируется или рассылается broadcast.

Псевдокод​

АЛГОРИТМ УМНОЖИТЬ_MPI_ПО_СТРОКАМ(n, A, B, C, rank, size)
  (строка_начало, строка_конец) := РАЗБИТЬ_1D(n, size, rank)
  если rank = 0 то
    РАЗОСЛАТЬ_ВСЕМ(B)
  иначе
    ПОЛУЧИТЬ(B)
  конец если

  для i от строка_начало до строка_конец
    для j от 0 до n − 1
      C[i,j] := СКАЛЯРНОЕ(A_строка_i, B_столбец_j)
    конец для
  конец для
КОНЕЦ

Коммуникация — broadcast B стоит O(n²) данных; приемлемо при малом p.

2D block decomposition (Cannon, SUMMA)​

Процессоры образуют решётку q×q, p = q². Блоки A, B, C — подматрицы на каждом rank.

Cannon (квадратная решётка)​

  1. Initial skew — сдвиг блоков A влево, B вверх по координатам процесса.
  2. q шагов — локальное умножение-накопление + циклический shift соседям.

Коммуникация — только ближайшие соседи, не полный broadcast B.

SUMMA​

На шаге k —

  1. Разослать блоки столбца k матрицы A по строкам решётки.
  2. Разослать блоки строки k матрицы B по столбцам.
  3. Локально C_block += A_slice * B_slice.

Проще Cannon в реализации; q итераций для решётки q×q.

Cannon — пошагово на решётке 2×2​

Матрицы 4×4, блоки 2×2 на 4 процессах:

P(0,0)  P(0,1)        A00 A01    B00 B01
P(1,0)  P(1,1)        A10 A11    B10 B11
ШагДействие
0Initial skew — сдвиг A влево на номер столбца процесса, B вверх на номер строки
1Локально C += A_local * B_local; shift A влево, B вверх
2Снова multiply-accumulate и shift

После q шагов каждый блок C накопил нужные слагаемые. Обмен — с соседями по кольцу, не пересылка всей B.

SUMMA — итерация k (псевдокод)​

для k от 0 до q − 1
  РАЗОСЛАТЬ_ПО_СТРОКАМ(блок_столбца_A[:,k])
  РАЗОСЛАТЬ_ПО_СТОЛБЦАМ(блок_строки_B[k,:])
  C_локально += УМНОЖИТЬ(полученный_кусок_A, полученный_кусок_B)
конец для

На GPU (идея без обязательного знания CUDA)​

cuBLAS sgemm — готовая реализация. Учебный tiled kernel делит матрицу на плитки, загружает плитки в быструю shared memory на чипе, синхронизирует потоки блока, считает частичные произведения.

Псевдокод уровня "плитка"​

АЛГОРИТМ ЯДРО_GPU_ПЛИТКА(блок_потоков, плитка_T)
  загрузить плитку A_T в быструю_память_блока
  загрузить плитку B_T в быструю_память_блока
  синхронизация_внутри_блока
  для каждого элемента плитки, закреплённого за потоком
    накопить произведение в регистре
  синхронизация_внутри_блока
  записать результат в C
КОНЕЦ

Справочно на CUDA (скелет)​

__global__ void matmul_tiled(float* C, const float* A, const float* B, int N) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];
  // загрузка плитки → __syncthreads → MAC → __syncthreads
}
ЭлементСмысл
__global__Функция выполняется на GPU, вызывается с CPU
__shared__Память, общая для потоков одного блока на GPU
TILEРазмер плитки (константа компиляции)

Warp из 32 потоков выполняет одну инструкцию над разными данными — см. SIMD и SIMT. Полный курс CUDA — после базового C/C++.

Сравнение схем​

СхемаПамять на узелКоммуникацияСложность реализации
Replicate BO(n²)Broadcast BНизкая
1D rowO(n²/p)Bcast / haloСредняя
2D Cannon/SUMMAO(n²/p)O(n²/√p) на фазуВысокая
Strassen (редко MPI)Меньше opsСложный обменИсследовательская

Верификация​

  1. Сравнить с последовательным эталоном на малых n (псевдокод или однопоточная программа).
  2. Относительная ошибка — ‖C_par − C_seq‖ / ‖C_seq‖ < ε для floating point.
  3. Графики T_p, speedup, efficiency — законы.

Модифицированный блочный метод​

В учебной литературе по HPC различают классическое блочное умножение и модифицированную схему — блоки A, B, C закрепляют за узлами так, чтобы внутренний тройной цикл по i, j, k выполнялся локально, а обмены шли только по границам блоков и по расписанию фаз (аналог SUMMA).

Идея для инженера — согласовать размер блока с L3 и топологией сети: блок 64×64 double ≈ 32 KiB на подматрицу, плюс halo при необходимости.

Strassen​

Асимптотика порядка O(n^2.807) — меньше операций при огромных n, но большие константы. На практике при n < 10 000 чаще выигрывает блочный O(n³) с BLAS. Параллельный Strassen — сложная коммуникация.

Что дальше​

Содержание