Краткий план оптимизации умножения матриц для CPU с кэшем и объяснения.
1) Цель и метрика
- Цель — максимизировать реальное использование FLOPs при минимизации пропуска в память; измерять пропускную способность памяти и полезную нагрузку (GFLOPS). Теоретическая сложность: $\Theta (n^3)$ FLOPs (для классического алгоритма), но производительность ограничена памятью без кеш‑оптимизаций.
2) Почему блочная (tiled / blocked) реализация нужна
- Блочная реализация повышает локальность данных и повторное использование кеш‑строк (temporal locality). Без блокировки элементы матриц часто перезапрашиваются из памяти → много кеш‑промахов.
- Арифметическая плотность (arithmetic intensity) увеличивается: для блока размеров $(m_b,n_b,k_b)$ количество операций и байт памяти приблизительно
$$\text{FLOPs}\approx 2m_b{n}_b{k}_b,\text{Bytes}\approx s(m_b{k}_b+k_b{n}_b+m_b{n}_b),$$ где $s$ — размер элемента (например $4$ для float, $8$ для double). Большая дробь $\frac{\text{FLOPs}}{\text{Bytes}}$ = выше загрузка CPU и меньше зависимость от памяти.
3) Общая схема оптимизации (пошагово)
- Профилирование: измерить базовую реализацию и понять, что ограничивает (memory bound vs compute bound).
- Выбор порядка циклов и блокирования: использовать трёхуровневое блокирование (register / L1 / L2/L3):
- микроблок (микро‑ядро) размера $(m_r,n_r,k_r)$ — помещается в регистры и кеш L1.
- макроблоки $(M_B,N_B,K_B)$ — заполняют L1/L2/L3.
- Пакование (packing): копировать блоки из A и B в смежные буферы перед вычислением, чтобы получить последовательный доступ и уменьшить TLB/кеш‑конфликты.
- Оптимизация микрокернела: вручную (intrinsics/ASM) использовать векторные инструкции и FMA; раскрытие циклов (unroll) и минимизация загрузок в регистры.
- Векторизация: использовать SIMD (AVX2/AVX‑512) — загрузки, умножения и накопления в векторах.
- Параллелизация: разделить макроблоки между потоками (OpenMP/TBB), избегать false sharing, привязка потоков к ядрам и NUMA‑распределение памяти.
- Выбор размеров блоков по кешам: подобрать $M_B,N_B,K_B$ по размеру L1/L2/L3 (см формулу дальше).
- Претфетчинг / порядок записи: программный prefetch, non‑temporal stores для больших выходов, выравнивание данных и паддинг для обхода remainder.
- Использование оптимизированных BLAS (OpenBLAS, MKL) либо авто‑тюнинг (ATLAS) если готовое решение допустимо.
4) Правила выбора размеров блоков (приближённо)
- Пусть доступно в L1 памяти $C_{L1}$ байт для данных (после учёта прочего). Требуется вместить пакеты A и B и частично C:
$$s(m_r{k}_b+k_b{n}_r+m_r{n}_r)\le C_{L1}.$$ - Обычно выбирают микроблоки так, чтобы $m_r,n_r$ соответствовали числу скалярных регистрах/векторных регистров (пример: для AVX2 double типичные $(m_r,n_r)=(4,4)$, для AVX‑512 — $(8,8)$), а $k_b$ подбирают так, чтобы $k_b$ уместился в L1 с учётом пакетов.
- K‑блок (часть суммирования) обычно выбирают так, чтобы пакованный блок B (или A) помещался в L2/L3: $s(M_B{K}_B+K_B{N}_B)\le C_{L2/L3}$.
5) Детали реализации, которые дают большой выигрыш
- Packing: копировать A‑блоки и B‑блоки в непрерывные буферы в порядке, удобном для микрокернела.
- Микрокернел на инстр.-уровне: использовать FMA и держать накопления в регистрах; минимизировать обращения к памяти.
- Векторные инструкции и выравнивание: выравнивание буферов на границы векторных регистров, использовать aligned loads/stores.
- Разворачивание циклов (loop unrolling) для уменьшения ветвлений.
- Prefetch: программный prefetch для будущих адресов, но осторожно — слишком ранний/поздний prefetch вреден.
- Non‑temporal stores (NT) при записи больших выходных блоков, чтобы не портить кешы.
- Использование больших страниц (huge pages) для уменьшения TLB‑промахов.
- NUMA и affinity: выделять память на локальной ноде, фиксировать потоки на ядрах.
- Избегать false sharing: каждая нить работает с собственным регионом C.
6) Алгоритмические и аппаратные приёмы дополнительно
- Алгоритмические:
- Strassen/Cannon и т. п.: снижает сложность (<$\Theta (n^3)$), но усложняет реализацию, ухудшает численную устойчивость и увеличивает память; полезно для очень больших размеров.
- Cache‑oblivious алгоритмы: рекурсивное деление, автоматически даёт блокирование на всех уровнях кеша.
- Аппаратные/системные:
- Использовать специализированные библиотеки (MKL, OpenBLAS) — они содержат тщательно оптимизированные микрокернелы.
- Аппаратная поддержка: AVX‑512, FMA, большой набор регистров — ускоряет микрокернел.
- Использование GPU или ускорителей для матрицы большого размера, если доступно.
- Настройка OS (huge pages), управление частотой CPU (P‑states), отключение турбо если неконсистентность.
7) Практические рекомендации и порядок работ
- Начните с реализации блочного GEMM с packing и базовым микрокернелом.
- Добавьте SIMD‑интринсики и регистровое блокирование.
- Подберите параметры блокирования экспериментально (auto‑tuning) согласно L1/L2/L3.
- Параллелизуйте по блокам, проверьте NUMA и affinity.
- Сравните с оптимизированной BLAS; если требуется — взять её реализацию или вдохновиться.
8) Короткие эмпирические подсказки
- Элемент размера $s$: float $s=4$, double $s=8$.
- Типично: микрокернел для AVX2 double: $m_r=4$, $n_r=4$; AVX‑512: $m_r=8$, $n_r=8$.
- Подбирать $K_B$ так, чтобы $K_B\approx \frac{C_{L1}}{s(m_r+n_r)}$.
Заключение: блочная реализация — фундамент: она повышает локальность и arithmetic intensity; дальше — микрокернел с SIMD/FMA, packing, правильный выбор размеров по кешам, параллелизация и системные настройки дают основной выигрыш. Использование готовых оптимизированных BLAS/авто‑тюнинга часто быстрее, чем самостоятельная оптимизация с нуля.