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

ОБЯЗАТЕЛЬНО

Разработчику Инженеру

Пример - параллельный обход с синхронизацией

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

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

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

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

---

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

Для матриц n × n:

for i = 0 .. n-1
  for j = 0 .. n-1
    C[i,j] = 0
    for k = 0 .. n-1
      C[i,j] += A[i,k] * B[k,j]

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

Порядок циклов влияет на локальность: вариант i-k-j лучше для cache (строка A, столбец B).

---

Элементный параллелизм (embarrassingly parallel output)

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

#pragma omp parallel for collapse(2)
for (int i = 0; i < n; ++i)
  for (int j = 0; j < n; ++j) {
    double sum = 0;
    for (int k = 0; k < n; ++k)
      sum += A[i,k] * B[k,j];
    C[i,j] = sum;
  }

Плюсы: простота, до n² параллельных задач.

Минусы: много обращений к B по столбцам (плохая locality); на distributed — B нужен всем.

Когда хватит: GPU (тысячи потоков), OpenMP на средних n.

---

Блочное умножение (cache-aware)

Разбить на блоки b × b (под L1/L2):

for i_block
  for j_block
    for k_block
      C_block += A_block * B_block   // маленькое matmul в cache

Параллелизм: независимые (i_block, j_block) при фиксированном проходе k — или parallel по k_block с reduction в C.

На одном узле это стандарт GotoBLAS/OpenBLAS.

---

1D row decomposition (MPI)

P процессов: строки A распределены по rank; полная B на каждом (replicate) или broadcast по мере надобности.

rank r владеет строками A[r·n/p .. (r+1)·n/p)
локально считает свои строки C

Коммуникация: 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 блоков по кольцу.

Каждый шаг: все процессоры делают GEMM на своих блоках; коммуникация — nearest neighbor.

Сложность коммуникаций: O(n² / √p) на шаг × √p шагов — лучше, чем naive broadcast.

SUMMA (Scalable Universal Matrix Multiply Algorithm)

Проще Cannon: q итераций; на шаге k broadcast строки блоков A и столбцов B по измерениям решётки; локальный update C.

---

---

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

Шаг 0 (initial skew): сдвиг A влево на coord col, B вверх на coord row.

Шаг 1: локально C += A_local * B_local; shift A left, B up.

Шаг 2: снова multiply-accumulate + shift.

После 2 шагов (для 2×2) каждый блок C накопил нужные произведения. Коммуникация — только с 4 соседями по кольцу, не broadcast всей B.

SUMMA — итерация k

На шаге k:

1. Broadcast блоков A[:, k] по строкам процессорной решётки.
2. Broadcast блоков B[k, :] по столбцам.
3. Локально: C_block += A_slice * B_slice.

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

---

На GPU

CUDA cuBLAS sgemm: пользователь передаёт матрицы в device memory; kernel — tiled matmul с shared memory.

Типичный учебный kernel:

__global__ void matmul_tiled(float* C, const float* A, const float* B, int N) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];
  // load tile → sync → multiply-accumulate → sync
}

Warp выполняет SIMD внутри tile — см. SIMD/MIMD.

---

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

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

---

Верификация

1. Сравнить с последовательным dgemm на малых n.
2. Norm-wise error: ‖C_par − C_seq‖ / ‖C_seq‖ < ε.
3. Масштабирование: plot T_p, S(p), E(p) — законы.

---

Strassen (упоминание)

O(n^2.807) — асymptotically быстрее, но большие константы; на практике для n < 10000 редко выигрывает у блочного O(n³) с BLAS. В параллельном виде — сложная коммуникация.

---

Что дальше

• Инженерия алгоритмов
• Анализ эффективности — Big O и практика
• Нейросети и GPU

---