Программа для умножения матриц

Программа для умножения матриц

В данной статье вы расскажете о программе для умножения матриц.

Введение

Умножение матриц является одним из основных алгоритмов, который широко используется в разных численных методах, включая алгоритмы машинного обучения. Много алгоритмов прямого и обратного перемещения сигнала в сверточных слоях нейронной сети основаны именно на данной операции. Достаточно часто примерно девяносто процентов всего времени, которое затрачивается на машинное обучение, занимает именно эта операция.

Это происходит из-за достаточно эффективной реализации данного алгоритма для процессора, графического ускорителя, а иногда и для специализированных ускорителей матричного умножения. Операция матричного умножения является одним из немногочисленных алгоритмов, позволяющих эффективно использовать всю совокупность вычислительных ресурсов передовых процессоров и графических ускорителей. По этой причине не вызывает удивления, что большинство алгоритмов сводится к операции матричного умножения. При этом возникающие дополнительные расходы, которые связаны с подготовкой данных, обычно с лихвой окупаются за счет общего ускорения алгоритмов.

Программа для умножения матриц

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

Рассмотрим пошагово процесс формирования программы для умножения матриц, и начать следует с постановки задачи. В общем варианте функция матричного умножения может быть представлена в следующем виде:

C[í,j] = a⁕ C$[$í,j$]$ + b⁕ Sum(A$[$í,k$]$⁕ B$[$k,ј$]$);

Здесь матрица A обладает размером M х K. Матрица B имеет размер K х N, а матрица C, соответственно, обладает размером M х N.

Практически без ущерба для изложения, можно принять, что a = 0 и b = 1, тогда:

C$[$í,j$]$ = Sum(A$[$í,k$]$⁕ B$[$k,ј$]$);

Реализация данного выражения на языке программирования С++ просто, как говорится, «в лоб» может быть оформлена следующим образом:

voíd gemm_v0(ínt M, ínt N, ínt K, const float ⁕ A, const float ⁕ B, float ⁕ C)

{

for (ínt í = 0; í ∠ M; ++í)

{

for (ínt ј = 0; ј ∠ N; ++ј)

{

C$[$í⁕ N + ј$]$ = 0;

for (ínt k = 0; k ∠ K; ++k)

C$[$í⁕ N + ј$]$ += A$[$í⁕ K + k$]$ ⁕ B$[$k⁕ N + ј$]$;

}

}

}

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

2⁕ M⁕ N⁕ K = 2⁕ 10^9.

Иначе говоря, производительность составила 1.6 GFLOPS, что можно расценить как очень далекий от теоретического предела результат однопоточной производительности для имеющегося процессора, для которого это значение равно примерно120 GFLOPS.

Сначала следует устранить наиболее очевидные недоработки алгоритма:

1. Определение адресов компонентов массивов следует упростить путем вынесения постоянной части из внутреннего цикла.
2. В оригинальном варианте доступ к компонентам массива B осуществляется не последовательно. Его следует упорядочить путем замены порядка вычисления таким образом, чтобы внутренним циклом был последовательный обход по строкам для каждой из трех матриц.

В результате получаем следующую реализацию:

voíd gemm_v1(ínt M, ínt N, ínt K, const float ⁕ A, const float ⁕ B, float ⁕ C)

{

for (ínt í = 0; í ∠ M; ++í)

{

float ⁕ c = C + í ⁕ N;

for (ínt ј = 0; ј ∠ N; ++ј)

c$[$ј$]$ = 0;

for (ínt k = 0; k ∠ K; ++k)

{ const float ⁕ b = B + k ⁕ N;

float a = A$[$í⁕ K + k$]$;

for (ínt ј = 0; ј ∠ N; ++ј)

c$[$ј$]$ += a ⁕ b$[$ј$]$;

}

}

}

Согласно результатов тестовых замеров время исполнения программы составило 250 мс, что соответствует 11.4 GFLOPS. То есть, эти достаточно небольшие правки позволили получить ускорение в восемь раз!

Следующим шагом будет векторизация внутреннего цикла. Если детально рассмотреть внутренний цикл (по переменной j), то можно увидеть, что процесс вычислений может проводиться блоками (векторами). Фактически любой современный процессор предоставляет возможность проведения вычисления над такими векторами. Например, совокупность инструкций AVX способна оперировать с векторами размерностью 256 бит и это позволяет исполнить восемь операций для вещественных чисел с одинарной точностью за такт. AVX2/FMA сделал еще один шаг вперед, то есть, он способен позволить исполнить слитную операцию умножения и сложения (d = a⁕ b + c) над вектором. Необходимо отметить, что инструкции AVX2/FMA достаточно просто могут быть задействованы напрямую из языка программирования С++ с помощью встроенных функций (intrinsics). Для AVX2/FMA их можно объявить в заголовочном файле ∠ímmintrin.h>:

voíd gemm_v2(ínt M, ínt N, ínt K, const float ⁕ A, const float ⁕ B, float ⁕ C)

{

for (ínt í = 0; í ∠ M; ++í)

{

float ⁕ c = C + í ⁕ N;

for (ínt ј = 0; ј ∠ N; ј += 8)

_mm256_storeu_ps(c + ј + 0, _mm256_setzero_ps());

for (ínt k = 0; k ∠ K; ++k)

{

const float ⁕ b = B + k ⁕ N;

__m256 a = _mm256_set1_ps(A$[$í⁕ K + k$]$);

for (ínt ј = 0; ј ∠ N; ј += 16)

{

_mm256_storeu_ps(c + ј + 0, _mm256_fmadd_ps(a,

_mm256_loadu_ps(b + ј + 0), _mm256_loadu_ps(c + ј + 0)));

_mm256_storeu_ps(c + ј + 8, _mm256_fmadd_ps(a,

_mm256_loadu_ps(b + ј + 8), _mm256_loadu_ps(c + ј + 8)));

}

}

}

}

Проверочные тесты показали время 217 мс или 13.1 GFLOPS, что означает ускорение всего лишь на пятнадцать процентов. Такой результат объясняется следующими факторами:

1. Программы компиляции сегодня достаточно продвинутые, и они способны сами справляться с задачей автоматического осуществления векторизации простых циклов. То есть, уже в первом случае компилятор фактически задействовал инструкции AVX2/FMA, потому ручная оптимизация не принесла практически ощутимых результатов.
2. Скорость выполнения расчетных операций в данном случае ограничивается не вычислительными возможностями процессора, а скоростью загрузки и выгрузки данных. В рассматриваемом случае процессору для задействования 2 256-bít FMA блоков требуется загрузить четыре и выгрузить 2 256-bít вектора за такт.