Кратко: исходный алгоритм имеет сложность $O(nm)$. Главное — минимизировать промахи кеша и позволить компилятору/процессору векторизовать внутренние вычисления. Ниже — рекомендации с примерами.
1) Анализ доступа к памяти
- Строки: если массивы в C-хранении (row-major), выражение $A[i][j]$ и $C[i][j]$ при внутреннем цикле по $j$ даёт последовательный доступ (хорошая пространственная локальность). $B[j]$ также читается последовательно.
- Переупорядочивание циклов (j вне, i внутри) делает $B[j]$ константой в внутреннем цикле (в регистре), но доступ к $A[i][j]$/$C[i][j]$ станет страйдным (шаг $m$) — плохая локальность, если $m$ велика.
2) Алгоритмические изменения
- Логика операций независима для каждой пары $(i,j)$; асимптотика не меняется ($O(nm)$), но можно изменить порядок/формат хранения:
- Если часто масштабируете столбцы, рассмотреть хранение матриц в column-major или транспонирование $A$/$C$ заранее, тогда цикл по столбцам будет эффективен.
- Альтернатива: выполнять операцию «по блокам столбцов» (tiling) для кеш-локальности.
3) Блокирование (tiling) по $j$ - Выбор блока $T$ так, чтобы блокы $A[i][j..j+T-1]$, $C[i][j..j+T-1]$ и $B[j..j+T-1]$ помещались в L1/L2.
- Структура:
for (j0 = 0; j0 < m; j0 += T)
for (i = 0; i < n; ++i)
for (j = j0; j < min(j0+T,m); ++j)
C[i][j] += A[i][j] * B[j];
- Блокирование сохраняет хорошую строковую проходность и увеличивает вероятность того, что блок B и кусок A/C останутся в кеше при повторном использовании.
4) Векторизация
- Внутренний цикл по $j$ при row-major очень векторизуем. Требования:
- Убедиться в отсутствии алиасинга (используйте ключевое слово restrict в C и/или компилять с -fno-strict-aliasing/соответствующими гарантиями).
- Выравнивание данных (aligned malloc) улучшает производительность.
- Дать подсказки компилятору: #pragma omp simd или #pragma GCC ivdep.
- Пример (простая версия, допускающая векторизацию):
for (i = 0; i < n; ++i) {
float *Ci = C[i];
float *Ai = A[i];
#pragma omp simd
for (j = 0; j < m; ++j)
Ci[j] += Ai[j] * B[j];
}
- Для максимума — писать SIMD-интринсики (AVX2/AVX-512). Для float + AVX2 шаг = $w=8$ элементов, для double — $w=4$. Используйте FMA-инструкции (_mm256_fmadd_ps) если есть.
5) Параллелизм по потокам
- Ряды независимы → параллелить по $i$ (OpenMP): #pragma omp parallel for — хорошо сочетается с внутренней векторизацией.
- Если параллелить по $j$, нужно учитывать гонки при записи в C (каждый j пишет в разные столбцы — безопасно, но доступы к A/C станут страйдными).
6) Выбор порядка циклов — резюме
- Row-major (C): лучший порядок — for i outer, for j inner (текущий код) для пространственной локальности и векторизации по j.
- Column-major или если хотите держать $B[j]$ в регистре и A/ C устроены по столбцам — наоборот: for j outer, for i inner.
- При больших размерах используйте блокирование, чтобы сохранить преимущества row-major и при этом повысить кеш-локальность B.
7) Микрооптимизации и настройки компилятора
- Компилятор: -O3 -march=native -ffast-math (осторожно с точностью) или -Ofast, включить vectorize.
- Выравнивание: _mm_malloc / posix_memalign, использовать __attribute__((aligned(64))).
- Разворачивание циклов вручную в критичных местах, префетч (использовать sparingly).
- Проверять производительность профайлером и командой компилятора, смотреть vectorization reports (например -Rpass=loop-vectorize в clang/gcc).
8) Пример оптимизированной комбинации (блокирование + параллелизм + pragma simd)
#define TILE 128
#pragma omp parallel for
for (int j0 = 0; j0 < m; j0 += TILE) {
int j1 = (j0 + TILE < m) ? j0 + TILE : m;
for (int i = 0; i < n; ++i) {
float *Ai = A[i];
float *Ci = C[i];
#pragma omp simd
for (int j = j0; j < j1; ++j)
Ci[j] += Ai[j] * B[j];
}
}
Выбор TILE подберите экспериментально (значения типично в диапазоне от $32$ до $1024$ в зависимости от L1/L2 и типа).
9) Замечания по точности
- Векторизация и FMA/перестановки могут изменить ассоциативность и дать немного другие плавающие результаты — учитывайте, если важна детерминированная ошибка.
Вывод: для стандартного row‑major C поддерживать текущий порядок циклов, добавить блокирование по $j$ для больших $m$, включить векторизацию (restrict, align, pragma), распараллелить по $i$ и при необходимости применять SIMD‑интринсики. Эти меры обычно дают наибольший выигрыш без изменения алгоритмической сложности $O(nm)$.