Привет, Хабр! На связи Кирилл Колодяжный, разработчик систем хранения данных в YADRO, ML-энтузиаст, автор книги Hands-on Machine Learning with C++.

Сегодня я расскажу про реализацию матричного умножения и особенности разработки для GPU. Познакомлю вас с устройством GPU, объясню, чем отличается программирование от привычного для CPU, какие нюансы нужно учитывать для эффективной реализации операций GEMM. А затем сравним производительность разных подходов к реализации.

Произведение матриц

Для начала вспомним, что такое произведение матриц. 

Если у нас есть матрица M на K, умноженная на матрицу K на N, мы получаем матрицу C размером M на N. 

Результирующие элементы матрицы C — это сумма произведения элементов соответствующей строки матрицы А на столбец матрицы B.

Если реализовывать это обычным способом для CPU, все выходит достаточно просто. Мы получаем три вложенных цикла, в которых перемножаем элементы и сохраняем результирующий элемент матрицы C:

void matmul(int M, int N, int K,
            const float* A, const float* const B, float* C) {
    for (int m = 0; m < M; m++)
        for (int n = 0; n < N; n++)
            for (int k = 0; k < K; k++) 
                C[m * M  + n] += A[m * M + k] * B[k * K + n];
}

Матричное умножение — элементарная операция из линейной алгебры, но используют ее повсеместно: 

• В моделирование процессов, управлении системами и планировании.

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

• В алгоритмах на графах, которые используются в навигации роботов и для построения систем SLAM. 

Обобщенное матричное умножение GEMM (general matrix multiplication)

Особенность GEMM в том, что это взвешенное произведение матриц, которое складывается с аккумулятором, тоже умноженным на какой-то коэффициент. Алгоритмическая сложность такого подхода с тремя циклами — это n³, когда у нас все размерности матриц одинаковые. Есть линейные алгоритмы, при которых можно уменьшить логарифмическую сложность. Она будет примерно n^2,4.

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

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

Далее я буду показывать только важные моменты из кода, но вы можете изучить проект с полноценными бенчмарками. Его можно скомпилировать, посмотреть примеры и проверить у себя на оборудовании.

Полный проект →

Чем отличается CPU от GPU 

И как это будет влиять на нашу дальнейшую разработку?

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

На архитектуре NVIDIA такие ядра сгруппированы еще и в стриминг-мультипроцессоры. В терминологии CUDA или в OpenCL это называется Compute Units.

Каждый такой стриминг-мультипроцессор (streaming multiprocessor) делится на блоки. Для архитектуры NVIDIA A100 стриминг-мультипроцессор состоит из четырех блоков (streaming processor), в которых находится по 16 ядер для операций единичной точности (float32), 8 ядер для операций двойной точности (float64), 16 ядер для целочисленных операций int32 и одно TensorCore для матричных операций. Аналог streaming multiprocessor в OpenCL — ComputeUnit. Однако иногда Compute Units, считают за один из стриминг-процессоров, четвертая часть мультипроцессора. 

В терминологии OpenCL Processing element — ядро, на котором работает Work Item-поток, выполняющий программу.

Модель памяти на GPU похожа на привычную CPU, но есть и существенные отличия.

Чтобы загрузить данные в GPU, мы сначала загружаем их в оперативную память, а затем в глобальную память GPU. Глобальная память GPU — это то, что мы знаем по спецификациям GPU и видеокарт: 32, 80, 100 гигабайт и так далее.

В рамках одного compute unit или стриминг-процессора есть еще локальная память (shared memory в терминологии CUDA). Она разделяется между теми потоками, которые запустились на стриминг-процессоре. Там эти потоки могут обмениваться данными в рамках локальной памяти.

Одному потоку еще доступна private memory или регистры в терминологии CUDA. Также они разделяют пространство памяти с кешем первого уровня: можно регулировать, сколько отойдет на регистры, сколько на кеш первого уровня. Размер файла регистров ограничен для каждого стриминг-процессора. Например для A100 — 16384х32bit. Количество параллельно работающих потоков ограничено использованием этих типов памяти. Соответственно если программа (далее называю ее kernel), которая работает в рамках потока, использует слишком много ресурсов, то количество параллельных потоков будет ограничено.

Ndrange как топология вычислительной задачи

В API, таких как CUDA, OpenCL, Vulkan, мы формулируем задачу для GPU в рамках какого-то N-мерного грида (сетки). Это может быть одномерная, двумерная и трехмерная сетка. Сетка задается количеством ворк-групп или блоков на верхнем уровне.

Такой ворк-групп будет выполняться на одном compute unit или стриминг-мультипроцессоре. В рамках него работают потоки — ворк-айтемы, количество и размеры которых тоже нужно задавать самостоятельно.

Напишем реализацию матричного умножения на GPU 

Теперь попробуем написать простую реализацию матричного умножения с учетом того, что мы узнали о GPU. И познакомимся с описанием вспомогательного кода, необходимой для OpenCL. OpenCL предоставляет два API — стандартный C-API, и С++-API для конфигурирования устройств, запуска программ и так далее.

Инициализируем буферы девайсов

Сначала мы должны инициализировать объект типа девайс. Он инициализируется индексом доступного устройства в системе. Для OpenCL такое устройство может быть как для GPU, так и интерфейс OpenCL может предоставлять CPU-реализацию, допустим, для Intel есть реализация OpenCL под центральные процессоры. Можно строить гетерогенные вычисления.

На основе девайса мы строим объект-контекст, который дальше используем для инициализации CommandQueue. 

CommandQueue — это очередь команд, в которой мы будем задавать наши команды выполнения на GPU: работа с памятью, запуск kernel и так далее.

И таких CommandQueue мы можем сделать несколько. У нас получаются еще более высокоуровневые потоки. Это своего рода аналог CUDA stream:

auto device    = cl::Device(…);
auto context   = cl::Context(device);
Auto cmd_queue = cl::CommandQueue(context, device);

buffer_a = cl::Buffer(context, CL_MEM_READ_ONLY,  matrix_data.mem_size_a());
buffer_b = cl::Buffer(context, CL_MEM_READ_ONLY,  matrix_data.mem_size_b());
buffer_c = cl::Buffer(context, CL_MEM_READ_WRITE, matrix_data.mem_size_c());

cmd_queue.enqueueWriteBuffer(buffer_a, CL_TRUE, 0, mem_size_a(), matrix_a.data());
cmd_queue.enqueueWriteBuffer(buffer_b, CL_TRUE, 0, mem_size_b(), matrix_b.data());
cmd_queue.enqueueWriteBuffer(buffer_c, CL_TRUE, 0, mem_size_c(), matrix_c.data());

Для переноса данных из оперативной памяти (хоста) системы на GPU мы должны описать объекты типа cl::Buffer, для которого надо задать его размер и тип работы с памятью: 

• read-only, 

• write-only, 

• read-write. 

Далее мы ставим команду в очередь через enqueueWriteBuffer и передаем туда ссылку на объект буфера, размер данных и указатель на данные, расположенные на хост-системе. Параметр CL_TRUE означает, что мы используем блокирующую операцию. А matrix_a — это просто вектор с данными.

Загружаем и компилируем kernel 

Далее нужна сама программа, которая будет вызываться на ядрах. 

Для этого в OpenCL есть два подхода: мы можем предварительно скомпилировать ее внешним компилятором или сделать это динамически прямо в рантайме. Чтобы скомпилировать программу в рантайме, используем объект типа cl::Program::Sources — на самом деле просто вектор строк:

cl::Program::Sources sources;
sources.push_back(load_kernel_source(kernel_file_name));

auto program = cl::Program(context, sources);
program.build(device, "-Werror");

Туда загружаются текстовые файлы с исходным кодом, после чего строится объект cl::Program на основе контекста и вызывается метод build для конкретного устройства, которое вызывает компилятор, а он предоставляет драйвер данного устройства. Вторым параметром передаются почти все стандартные для GCC флаги компиляции. В данном случае — Werror.

Задаем аргументы для запуска kernel 

После того, как скомпилировали набор наших исходных файлов, нужна точка входа программы (конкретная функция), которая будет подаваться в ядро.

auto gemm_kernel = cl::Kernel(program, kernel_name.c_str());
gemm_kernel.setArg(0, matrix_data.m);
gemm_kernel.setArg(1, matrix_data.n);
gemm_kernel.setArg(2, matrix_data.k);
gemm_kernel.setArg(3, buffer_a);
gemm_kernel.setArg(4, buffer_b);
gemm_kernel.setArg(5, buffer_c);

Чтобы получить конкретную точку входа, мы используем объект cl::Program, который содержит все скомпилированные файлы. Вторым параметром указываем строковое название функции, которая будет выполняться. После чего для объекта kernel надо установить, откуда он будет брать аргументы, и определить их количество. Устанавливаются аргументы по индексу, методом setArg.

В нашем случае первые три аргумента — это размерности матриц, целочисленные значения. Дальше по коду ссылки на объекты типы cl::Buffer, которые мы определили перед этим.

Задаем топологию задачи для kernel 

Далее, чтобы запустить kernel, нам нужно описать топологию нашей задачи или NDRange в терминологии OpenCL. Обычно задаются глобальный NDRange, глобальный размер сетки и размер блока — local range.

auto local_range  = cl::NDRange{32, 32};
auto global_range = cl::NDRange{matrix_data.m, matrix_data.n};

Здесь я описываю local range как 32 на 32. Использовать его напрямую мы сейчас не будем, но задать все равно нужно — это повысит производительность. Размер local_range надо подбирать так, чтобы с одной стороны загрузить compute unit, а с другой оставаться в рамках ограничений по local и private memory. А в глобальном range сетка будет равно размерности выходной матрицы — m на n. Мы в нашей программе (одном потоке) будем считать один элемент результирующей матрицы.

Если посмотреть, как задается наша сетка, то global range — это один элемент в нашей матрице. Индексацию по Local range мы прямо не используем — только по необходимости.

Запускаем kernel 

Подготовив все нужные нам объекты, мы опять ставим следующую команду в очередь, используя метод enqueueNDRangeKernel. Объект kernel отправляется выполняться на GPU с учетом заданных global и local range.

Чтобы синхронизироваться, мы используем объект cl::Event:

cl::Event event;
cmd_queue.enqueueNDRangeKernel(gemm_kernel,
                               cl::NullRange, 
                               global_range,
                               local_range, 
                               nullptr, 
                               &event);
event.wait();

Использовать его постоянно при запуске kernel необязательно, так как мы можем запускать код и на хосте, и на GPU. Объекты cl::Event нужны, чтобы задать точки синхронизации именно там, где вам нужно.

И вот так выглядит код простого подхода для вычисления одного элемента матрицы C:

__kernel void naive_gemm(int m, int n, int k,
                         const __global float* a, 
                         const __global float* b,  
                        __global float* c) {
  const int global_row = get_global_id(0);
  const int global_col = get_global_id(1);
  float acc = 0.0f;
  for (int i = 0; i < k; ++i) {
    acc += a[i * m + global_row] * b[global_col * i + k];
  }
  c[global_col * m + global_row] = acc;
}

Как мы задавали ранее, первые три аргумента — размеры матриц и следующие три аргумента — это буферы данных. Через встроенные функции OpenCL get_global_id мы получаем индексы потока, в котором сейчас выполняется код. get_global_id — это ID потоков в двумерных координатах внутри global range.

Что делаем далее:

• используя эти индексы, вычисляем сумму произведения элементов по внутренней размерности К — одной строки и одного столбца соответствующих матриц, 

• складываем в переменную «аккумулятор»,

• передаем это в глобальную память C.

Измеряем производительность

Я буду приводить дальше сравнение производительности, которую я делал на своем ноутбуке. Его характеристики:

• Intel Core i9-13900HX: 8 / 16 ядер, 2.2 ГГц - 1.6 ГГц,

• RAM 32 ГБ,

• GeForce RTX 4070 8ГБ:

  • 46 SMs,

  • 184 Tensor Cores,

  • L1 Cache 128 KB (per SM),

  • L2 Cache 36 MB.

Производительность простого подхода представлена на графике ниже. Я сравнивал его с двумя профессиональными библиотеками под CPU — OpenBLAS и Intel MKL — и двумя библиотеками под GPU — clBlast, реализованная на OpenCL, и cuBLAS от NVIDIA.

Вот такой простой подход вполне можно сравнить с реализацией Intel MKL. Он дает существенный прирост в производительности на относительно небольших размерностях матриц. Однако на больших размерностях подход отстает:

Проблема простого подхода и как его можно улучшить

Как я сказал в начале, операция матричного умножения должна быть арифметически более интенсивна, чем обращение к памяти. Но у нас не совсем получилось этого достичь, потому что каждый поток достаточно интенсивно загружает элементы из глобальной памяти, а потом туда еще и пишет. К тому же, операция умножения и сложения компилятором объединится в оптимизированную команду под GPU. Скоре всего, это будет команда FMA (fused multiply add).

Поэтому хотелось бы снизить количество обращений к памяти. Дело в том, что у GPU есть разные виды памяти: глобальная, локальная и регистры (private memory). Чем ближе память к регистрам, тем быстрее к ней обращаться.

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

Сделаем реализацию GPU лучше

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

Как реализовать обработку тайлов

Здесь мы начнем использовать индексы локального NDRange, которых раньше не использовали. Он тоже задан размером 32 на 32 элемента — это квадратный блок (тайл). Внутренний цикл по общей размерности матриц К мы разделим  на две части: сначала загрузим этот тайл в локальную память, а потом сделаем умножение. Мы будем обрабатывать тайлы по строкам и столбцам: каждую строку тайла будем умножать на соответствующий столбец.

Нам потребуются немножко другие, более детальные индексы внутри каждого потока. То есть когда мы разрабатываем какой-то кернел для GPU, достаточно большая часть работы — это правильно организовать индексацию в потоке для доступа к данным. В данном случае я использую функцию get_local_id:

int row = get_local_id(0); // Local row ID (max: TS)
int col = get_local_id(1); // Local col ID (max: TS)
int global_row = TS * get_group_id(0) + row; // Row ID of C (0..M) 
int global_lol = TS * get_group_id(1) + col; // Col ID of C (0..N)

Это получение индекса потока внутри одного из блоков (внутри work-group), в котором он работает. get_group_id возвращает индекс самого блока, в котором работает текущий поток. Таким образом, мы получаем два индекса — глобальный и локальный.

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

__local float a_sub[TS][TS]; // Local Memory для A tile
__local float b_sub[TS][TS]; // Local Memory для B tile

float acc = 0.0f; // Аккумулятор для результата вычислений

И изменяем наш цикл, который был достаточно простым:

// Цикл по tiles
const int num_tiles = K / TS;
for (int t = 0; t < num_tiles; t++) {
    // Загружаем tile для A и B в local memory
    …
    // Выполняем вычисления для одного tile
    …
}
// Сохраняем результат в C
C[global_col * M + global_row] = acc;

Как это работает:

• посчитаем количество тайлов в зависимости от размера тайла — это будет верхнеуровневый цикл,

• загрузим все в локальную память текущего тайла,

• посчитаем результат, накопим его и сохраним результат.

Как организована загрузка в локальную память

Каждый поток загружает один элемент тайла в зависимости от его индекса. Задача не представляется последовательно, как в императивном стиле, а разбивается на обработку отдельных элементов. Тут мы должны немножко поменять восприятие написания программы и учесть, что kernel стартует с индекса — он и будет использоваться для обращения к данным.

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

Загружаем тайлы в локальную память

Перенос данных из глобальной памяти в локальную — это просто копирование, оператор присваивания из элементов одного массива в другой с соответствующими индексами. А далее мы используем функцию barrier с параметром CLK_LOCAL_MEM_FENCE. Каждый поток во время выполнения будет дожидаться завершения всех операций c локальной памятью:

const int tiled_row = TS * t + row;
const int tiled_col = TS * t + col;

a_sub[col][row] = A[tiled_col * M + global_row];
b_sub[col][row] = B[global_col * K + tiled_row];
 
// синхронизируем что быть уверенным, что tile загружен
barrier(CLK_LOCAL_MEM_FENCE);

Такой barrier можно делать не только на локальную память, но и на глобальную и память изображений.

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

for (int k=0; k < TS; k++) {
    acc += a_sub[k][row] * b_sub[col][k];
}
 
// синхронизируем перед следующим tile
barrier(CLK_LOCAL_MEM_FENCE);

Почему тайлинг работает?

Дело в том, что мы сократили количество доступа к глобальной памяти — у нас существенно уменьшилась задержка (latency) работы с памятью. К тому же, загружая данные из потоков поэлементно, мы и оптимизировали доступ к глобальной памяти. Он у нас стал сгруппированным (coalesced).

Я дальше покажу, в чем суть. 

Оптимизация работы с глобальной памятью в рамках локальной

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

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

Как мы оптимизировали еще и доступ к глобальной памяти

На GPU операции с памятью делаются группами — транзакциями. Каждая транзакция должна быть определенного размера и правильно выровнена. Например, для доступа к 32-байтовым элементам одна транзакция будет размером 128 байт. Если обращения к памяти идут по одному элементу подряд (несгруппированно), то транзакций будет больше, и пропускная способность снизится. А если обращаться к элементам последовательно, то транзакций будет меньше, и работа с памятью станет быстрее.

Как это визуализировать: есть трехмерная матрица и три потока, которые работают с ее строками. В памяти она расположена линейно, и каждый поток при доступе к первому элементу соответствующей строки будет это делать через смещение (stride). Если размер этого смещения не соответствует требованиям размера транзакции и выравнивания. то у нас будет хуже производительность.

Одним из подходов к решению проблемы в данном случае может быть обработка не строк, а столбцов. 

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

Результаты замера производительности

Naïve GEMM (прошлая реализация) — это серая линия. Более темным серым выделена текущая реализация Tiling GEMM. Видно что она уже в несколько раз превзошла наше предыдущее решение и решения на CPU, приблизившись к решениям на GPU. Значит, подход с учетом архитектурных особенностей уже сыграл существенную роль.

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

Улучшаем решение и соблюдаем баланс в работе с памятью

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

Что мы можем сделать? Переместить данные в регистры. И нужно нагрузить поток, чтобы он больше выполнял арифметических действий. В каждом тайле мы обрабатывали строку и столбец. А теперь сделаем так, что каждый поток GPU будет обрабатывать двумерный тайл, то есть все строки и столбцы конкретного тайла, а данные перенесем в регистры.

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

Что нам для этого нужно? Заведем переменные локальные, которые попадут в регистры и в которые мы перенесем данные из локальной памяти:

#define WPTM 8           // фактор уменьшения тайла по M
#define WPTN 8           // фактор уменьшения тайла по N
…
float a_reg;
float b_reg[WPTN];
float acc[WPTM][WPTN];   // 2D аккумуляторы
…

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

Визуализировать алгоритм можно так:

Берем один элемент тайла матрицы А, загружаем  строку матрицы В уножаем и складываем  с результатом. Меняя один из элементов матрицы А, мы не меняем строку, — она остается в private памяти, и накапливаем результат в аккумуляторе.

Реализуется это так:

// цикл по тайлам
for (int k = 0; k < TSK; k++) {
  for (int wn = 0; wn < WPTN; wn++) {
    int col = local_row + wn * RTSN;
    b_reg[wn] = b_sub[col][k];
  }
…

Остается тот же самый цикл по тайлам, однако добавляются циклы, которые загружают промежуточные значения (строку матрицы В) из локальной памяти в массивы в регистрах:

…
// циклы по элементам тайла
for (int wm = 0; wm < WPTM; wm++) {
    int row = local_row + wm * RTSM;
    a_reg = a_sub[k][row];
…

Далее мы идем по элементам тайла, где выбираем соответствующий элемент тайла А:

for (int k = 0; k < TSK; k++) {
 … b_reg = …
 for (int wm = 0; wm < WPTM; wm++) {
    … a_reg = …
    for (int wn = 0; wn < WPTN; wn++) {
      acc[wm][wn] += a_reg * b_reg[wn];
    }
  }
} 

И уже с помощью уже предзагруженных элементов для тайла В умножаем на соответствующий элемент тайла А и сохраняем наш двумерный аккумулятор.

Как сделать еще лучше

Решение выше снова можно улучшить.

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

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

Однако, скорее всего, такой подход приведет к конфликту банков памяти в локальной памяти. Так как у нас элементы в тайле А и в тайле В будут расположены с одинаковой индексацией.

Почему? Дело в том, что локальная память GPU поделена на банки памяти, и последовательные адреса попадают в следующие банки памяти.

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

Если же у нас все потоки, которые выполняются в рамках Compute Unit, обращаются к одному банку, такие обращения к памяти будут последовательными — существенно менее эффективными.

Визуализировать конфликт банков памяти можно следующим образом:

Допустим, у нас есть четыре банка памяти, на которые ложится матрица 3 на 4, и опять три потока, которые работают со строками. И на первом же такте мы обратимся к одному банку памяти. Такой доступ будет последовательный. 

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

Программно реализуется это достаточно просто. Мы добавляем какое-то количество элементов к размерности локальной памяти:

__local float b_sub[TSN][TSK + 2];

Размерность padding зависит от типа данных, алгоритма, индексации и так далее. Стоит продумывать, сколько нужно добавлять. 

Результаты замера производительности

Черным цветом отмечены результаты optim GEMM для решения со всеми улучшениями, про которые мы говорили. Мы получили производительность, сравнимую с хорошими библиотеками для GPU.

На больших размерах результаты уже вполне сравнимы с clBlast, производительность существенно выросла.

На что еще нужно обращать внимание при разработке эффективных ядер для GPU-программ

Давайте познакомимся с понятием occupancy. Occupancy в вычислениях для GPU— сущность, которая показывает, насколько эффективно потоковые мультипроцессоры (SMs, Compute Units) на GPU используются параллельными потоками. Occupancy эквивалентно отношению количества активных Warp (wave fronts) к максимальному количеству Warp, которое теоретически может поддерживать SM. Факторы, которые влияют на GPU occupancy:

• Размер блока потоков: количество потоков в блоке влияет на то, сколько блоков может одновременно находиться в SM.

• Использование регистров: каждый поток использует регистры, а увеличение использования регистров в потоке уменьшает количество одновременных потоков.

• Использование общей памяти: общая память выделяется на каждый блок, поэтому при более высоком уровне использования количество активных блоков ограничено.

• Архитектура графического процессора: у разных графических процессоров NVIDIA — разные ограничения на количество потоков, блоков и регистров на SM.

При этом стоит учитывать, что эти ресурсы делятся с кешем первого уровня, что также влияет на эффективность использования, а потоки в рамках work-group организуются в Wavefront (Warp). 

И compute unit оперирует именно волнами потоков, а не единичными потоками. Когда формируются транзакции доступа к памяти они как раз собираются из операций волны. 

Хорошая GPU Occupancy влияет на:

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

• Использование ресурсов: обеспечивает эффективное использование ресурсов графического процессора, таких как вычислительные блоки и пропускная способность памяти.

Скрытие задержки доступа к памяти можно проиллюстрировать следующим рисунком:

У планировщика всегда должно быть достаточное количество Wavefront, пока одни работают с памятью другие занимаются вычислениями.

Еще стоит учитывать, что у потоков в Wavefront одна программа — Program Address. В современных GPU у каждого потока свой Program Counter, это позволяет более эффективно выполнять потоки исполняющие разные ветки ветвления. Но ветвления стоит избегать, так как часть потоков все равно не будет работать.

Как было показано ранее, в системе могут быть доступны специальные вычислительные ядра, предназначенные для выполнения определенных операций более оптимально. Например, Tensor Cores на NVIDIA GPU выполняют умножение небольших матриц и позволяют значительно ускорить соответствующие вычисления. Аналогичные инструкции существуют и у Intel GPU — их тоже следует учитывать при оптимизации.

Важно помнить о нюансах точности таких специальных ядер: они могут выполнять вычисления в половинной точности (FP16), что подходит не для всех задач. Поэтому необходимо учитывать требования к точности ваших алгоритмов и решаемых задач.

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

Использование intrinsics в OpenCL позволяет напрямую обращаться к возможностям аппаратного обеспечения, например, для загрузки данных из глобальной памяти с помощью текстурных объектов, которые могут выступать в роли дополнительного кэша. Это может существенно повысить пропускную способность и общую производительность.

Для эффективной реализации алгоритмов на GPU необходимо хорошо разбираться в архитектуре устройства: как работают механизмы памяти (глобальная, локальная, регистры), как взаимодействуют различные компоненты системы и какие подходы наиболее подходят под конкретную архитектуру. Однако даже без глубокого знания архитектуры можно добиться высокой эффективности — GPU все равно обеспечивают значительно лучшие показатели по сравнению с CPU при выполнении задач требующего большого объёма относительно однотипных вычислений.

Если ваша задача предполагает интенсивные вычисления и использование GPU оправдано — стоит попробовать реализовать ее именно там, чтобы максимально использовать преимущества параллельной обработки.

Что еще изучить по теме:

• OpenCL SGEMM implementaiton

• Beating cuBLAS in Single-Precision GEMM

• Optimized FP32 matrix multiplication on AMD RDNA3 GPU