Оптимизация выполнения ядра CUDA для расчета треугольной матрицы
Я занимаюсь разработкой своего первого приложения Cuda, и у меня есть ядро с «ниже ожидаемой пропускной способностью», что на данный момент является самым узким местом.
Задача ядра — вычислить матрицу размером N на N (DD) содержит квадраты расстояний между всеми элементами в матрице данных. Матрица данных (Y) имеет размер N на D (для поддержки многомерных данных) и хранится как мажорная строка.
Источник:
__global__ void computeSquaredEuclideanDistance(const float * __restrict__ Y, float * __restrict__ DD, const int N, const int D) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for (int i = index; i < N * N; i += stride) {
const int m = i / N;
const int n = i % N;
float tmp = 0;
for (int d = 0; d < D; ++d) {
const float Ynd = Y[d + D * n];
const float Ymd = Y[d + D * m];
const float Ydiff = Ynd - Ymd;
tmp += Ydiff * Ydiff;
}
DD[n + N * m] = tmp;
}
}
Это называется с size_t blockSize = 256 а также size_t numBlocks = (N*N + blockSize - 1)/blockSize,
Как я могу оптимизировать это ядро? Сначала я думал, что трудоемкая часть — это чтение данных без использования какой-либо разделяемой памяти, но может ли кто-нибудь подсказать мне, как к этому подойти?
Замечания от nvvc инструмент профилирования:
- Анализ латентности:
- Использование вычислений на уровне около 40%
- Использование памяти (кэш-память второго уровня) около 35%
- Занятость не проблема
- Активные деформации на 57,59 от теоретического 64
- Заполняемость на уровне 90% от теоретического 100
Для моего приложения типичные значения:
- 5k <
N< 30k D2 или 3
Решение
Я обычно игнорирую эти типы вопросов оптимизации, потому что они находятся на грани не по теме, по моему мнению. Что еще хуже, вы не предоставляете MCVE поэтому любой, кто пытается ответить, должен написать весь собственный код поддержки для компиляции и тестирования вашего ядра. И такая работа требует бенчмаркинга и анализа кода. Но поскольку ваша проблема — это в основном проблема линейной алгебры (а мне нравится линейная алгебра), я ответил на нее, а не закрыл голосование за слишком широкую ……
С этим с моей груди. Есть несколько вещей, которые сразу же появляются в коде, которые можно улучшить и которые, вероятно, окажут существенное влияние на время выполнения.
Во-первых, счетчик срабатываний внутреннего цикла известен априори. В любое время, когда у вас возникнет такая ситуация, сообщите об этом компилятору. Развертывание цикла и переупорядочение кода — очень мощная оптимизация компилятора, и компилятор NVIDIA в этом очень хорош. Если вы переместите D в параметр шаблона, вы можете сделать что-то вроде этого:
template<int D>
__device__ float esum(const float *x, const float *y)
{
float val = 0.f;
#pragma unroll
for(int i=0; i<D; i++) {
float diff = x[i] - y[i];
val += diff * diff;
}
return val;
}
template<int D>
__global__
void vdistance0(const float * __restrict__ Y, float * __restrict__ DD, const int N)
{
int index = blockIdx.x * blockDim.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for (int i = index; i < N * N; i += stride) {
const int m = i / N;
const int n = i % N;
DD[n + N * m] = esum<D>(Y + D * n, Y + D * m);
}
}
template __global__ void vdistance0<2>(const float *, float *, const int);
template __global__ void vdistance0<3>(const float *, float *, const int);
Компилятор будет встроен esum и разверните внутренний цикл, и он может затем использовать свою переупорядоченную эвристику для лучшего чередования нагрузок и флопов, чтобы улучшить пропускную способность. Полученный код также имеет более низкую площадь регистра. Когда я запускаю это для N = 10000 и D = 2, я получаю ускорение примерно на 35% (7,1 мс против 4,5 мс на GTX 970 с CUDA 9.1).
Но есть еще более очевидная оптимизация, чем эта. Выполняемый вами расчет даст симметричную выходную матрицу. Вам нужно только сделать (N*N)/2 операции для вычисления полной матрицы, а не N*N вы делаете в своем коде [технически N(N/2 -1) потому что диагональные элементы равны нулю, но давайте забудем диагональ для целей этого обсуждения].
Поэтому, используя другой подход и используя один блок для вычисления каждой строки верхней треугольной выходной матрицы, вы можете сделать что-то вроде этого:
struct udiag
{
float *p;
int m;
__device__ __host__ udiag(float *_p, int _m) : p(_p), m(_m) {};
__device__ __host__ float* get_row(int i) { return p + (i * (i + 1)) / 2; };
};template<int D>
__global__
void vdistance2(const float * __restrict__ Y, float * __restrict__ DD, const int N)
{
int rowid = blockIdx.x;
int colid = threadIdx.x;
udiag m(DD, N);
for(; rowid < N; rowid += gridDim.x) {
float* p = m.get_row(rowid);
const float* y = Y + D * rowid;
for(int i=colid; i < (N-rowid); i += blockDim.x) {
p[i] = esum<D>(y, y + D * i);
}
}
}
template __global__ void vdistance2<2>(const float *, float *, const int);
template __global__ void vdistance2<3>(const float *, float *, const int);
При этом используется небольшой вспомогательный класс для инкапсуляции номеров треугольников, необходимых для схемы адресации для верхней треугольной выходной матрицы. Это позволяет сэкономить огромное количество памяти и пропускную способность памяти, а также уменьшить общее количество FLOP для расчета. Если после этого вам нужно сделать что-то другое, BLAS (и CUBLAS) поддерживает вычисления на верхних или нижних треугольных матрицах. Используй их. Когда я запускаю это, я получаю ускорение примерно на 75% (7,1 мс против 1,6 мс на том же GTX 970).
Огромный отказ от ответственности: Весь код, который вы видите здесь, был написан во время 45-минутного перерыва на обед и как было очень слегка проверено. Я абсолютно не утверждаю, что что-то в этом ответе действительно правильно. Я подтвердил, что он компилируется и не выдает ошибку времени выполнения при запуске для получения данных профилирования. Вот и все. Кавает Эмптор и все такое.
Другие решения
Других решений пока нет …