Оптимизировать байтовые операции CUDA

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

Я получаю только 5 ГБ / с, что не намного лучше, чем быстрый процессор, и даже не приближается к теоретической пропускной способности моего DDR5 GT640 около 80 ГБ / с.

мой набор данных состоит из вектора запроса 1024 байта, базы данных 100 000 x 1024 байта

У меня 100 000 блоков из 128 потоков, если каждый блок обращается к одному и тому же 1024-байтному query_vector, это приведет к снижению производительности? Так как каждый блок имеет доступ к одной и той же ячейке памяти.

blockSize и общая память установлены на 128 и 128 * sizeof (int), 128 — # define’d как THREADS_PER_BLOCK

template<UINT blockSize> __global__ void reduction_sum_abs( BYTE* query_vector, BYTE* db_vector, uint32_t* result )
{
extern __shared__ UINT sum[];
UINT db_linear_index = (blockIdx.y*gridDim.x) + blockIdx.x ;
UINT i = threadIdx.x;

sum[threadIdx.x] = 0;

int* p_q_int = reinterpret_cast<int*>(query_vector);
int* p_db_int = reinterpret_cast<int*>(db_vector);

while( i < VECTOR_SIZE/4 ) {

/* memory transaction */
int q_int = p_q_int[i];
int db_int = p_db_int[db_linear_index*VECTOR_SIZE/4 + i];

uchar4 a0 = *reinterpret_cast<uchar4*>(&q_int);
uchar4 b0 = *reinterpret_cast<uchar4*>(&db_int);

/* sum of absolute difference */
sum[threadIdx.x] += abs( (int)a0.x - b0.x );
sum[threadIdx.x] += abs( (int)a0.y - b0.y );
sum[threadIdx.x] += abs( (int)a0.z - b0.z );
sum[threadIdx.x] += abs( (int)a0.w - b0.w );

i += THREADS_PER_BLOCK;

}

__syncthreads();

if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}

/* reduce the final warp */
if ( threadIdx.x < 32 ) {
if ( blockSize >= 64 ) { sum[threadIdx.x] += sum[threadIdx.x + 32]; } __syncthreads();

if ( blockSize >= 32 ) { sum[threadIdx.x] += sum[threadIdx.x + 16]; } __syncthreads();

if ( blockSize >= 16 ) { sum[threadIdx.x] += sum[threadIdx.x + 8 ]; } __syncthreads();

if ( blockSize >= 8  ) { sum[threadIdx.x] += sum[threadIdx.x + 4 ]; } __syncthreads();

if ( blockSize >= 4  ) { sum[threadIdx.x] += sum[threadIdx.x + 2 ]; } __syncthreads();

if ( blockSize >= 2  ) { sum[threadIdx.x] += sum[threadIdx.x + 1 ]; } __syncthreads();

}/* copy the sum back to global */
if ( threadIdx.x == 0 ) {
result[db_linear_index] = sum[0];
}
}

Я могу получить примерно 4-кратное увеличение пропускной способности, если запустить ядро ​​с четырьмя строками кода, которые закомментируют фактическое вычисление абсолютной разности, очевидно, это приводит к неправильному ответу, но я считаю, что по крайней мере значительная часть времени провел там.

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

Является ли мое использование reinterpret_cast правильный?

Есть ли лучший метод для выполнения 8-битных беззнаковых вычислений?

Какие еще (я бы сказал, многие, так как я начинающий) оптимизации можно сделать?

Спасибо

РЕДАКТИРОВАТЬ:

Мои технические характеристики машины следующие:

Windows XP 2002 SP3

Intel 6600 2,40 ГГц

2 ГБ оперативной памяти

GT640 GDDR5 1 ГБ

Visual C ++ 2010 Express

1

Решение

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

Исправление ошибок:

В вашем коде есть как минимум 2 ошибки, одна из которых @Jez уже указала. После этого шага «частичного сокращения»:

if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}

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

if ( threadIdx.x < 32 ) {

это нелегальный иметь __syncthreads() оператор внутри блока условного кода:

  if ( blockSize >= 64 ) { sum[threadIdx.x] += sum[threadIdx.x + 32]; } __syncthreads();

(и аналогично для последующих строк, которые делают то же самое). Поэтому рекомендуется это исправить. Есть несколько способов решить эту проблему, один из которых — перейти на использование volatile набранный указатель для ссылки на общие данные. Так как мы сейчас работаем в варпе, volatile квалификатор заставляет компилятор делать то, что мы хотим:

volatile UINT *vsum = sum;
if ( threadIdx.x < 32 ) {
if ( blockSize >= 64 ) vsum[threadIdx.x] += vsum[threadIdx.x + 32];
if ( blockSize >= 32 ) vsum[threadIdx.x] += vsum[threadIdx.x + 16];
if ( blockSize >= 16 ) vsum[threadIdx.x] += vsum[threadIdx.x + 8 ];
if ( blockSize >= 8  ) vsum[threadIdx.x] += vsum[threadIdx.x + 4 ];
if ( blockSize >= 4  ) vsum[threadIdx.x] += vsum[threadIdx.x + 2 ];
if ( blockSize >= 2  ) vsum[threadIdx.x] += vsum[threadIdx.x + 1 ];
}

CUDA пример кода с параллельным сокращением а также связанный PDF может быть хорошим обзором для вас.

ВРЕМЯ / ИДЕАЛЬНЫЙ АНАЛИЗ:

У меня случается GT 640, устройство cc3.5. Когда я бегу bandwidthTest на нем я наблюдаю около 32 ГБ / с для передачи между устройствами. Это число представляет разумную приблизительную верхнюю границу достижимой полосы пропускания, когда ядра устройства обращаются к памяти устройства. Кроме того, когда я добавляю cudaEvent на основе синхронизации и построения примера кода на основе того, что вы показали, с имитацией данных, я наблюдаю пропускную способность около 16 ГБ / с, а не 5 ГБ / с. Таким образом, ваша фактическая методика измерения была бы полезной информацией здесь (фактически, полный код, вероятно, необходим для анализа различий между моей синхронизацией вашего ядра и вашей синхронизацией).

Остается вопрос: можно ли его улучшить? (при условии, что ~ 32 ГБ / с — приблизительная верхняя граница).

ВАШИ ВОПРОСЫ:

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

Поскольку ваше ядро ​​фактически эффективно загружает байты в виде 32-битной величины (uchar4), и каждый поток загружает смежное, последовательное 32-битное количество, я не верю, что есть какие-либо проблемы с конфликтом банков с вашим ядром.

Правильно ли я использую reinterpret_cast?

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

Есть ли лучший метод для выполнения 8-битных беззнаковых вычислений?

Существует, и в этом случае, как указывает @njuffa, SIMD присущи может справиться с этим, как оказалось, с помощью одной инструкции (__vsadu4()см. пример кода ниже).

Какие еще (я бы сказал, многие, так как я начинающий) оптимизации можно сделать?

  1. Используйте cc3.0 метод уменьшения деформации основы, предложенный @MichalHosala

  2. Воспользуйтесь внутренним SIMD __vsadu4() упростить и улучшить обработку количества байтов, предложенного @njuffa.

  3. Реорганизовать векторные данные вашей базы данных, чтобы они находились в главном хранилище столбцов. Это позволяет нам обойтись без обычного метода параллельного сокращения (даже такого, как упомянуто в пункте 1) и переключиться на ядро ​​с прямым циклом чтения, один поток вычисляет полное сравнение векторов. Это позволяет нашему ядру приблизительно достигать пропускной способности памяти устройства в этом случае (cc3.5 GT640).

Вот код и пример выполнения, показывающий 3 реализации: ваша оригинальная реализация (плюс вышеперечисленные «исправления» для получения правильных результатов), ядро ​​opt1, которое модифицирует ваше, чтобы включить пункты 1 и 2 из списка выше, и ядро ​​opt2, которое заменяет ваше на подход, использующий 2 и 3 из списка выше. Согласно моим измерениям, ваше ядро ​​достигает 16 ГБ / с, около половины пропускной способности GT640, ядро ​​opt1 работает со скоростью около 24 ГБ / с (увеличение примерно в равных долях по сравнению с пунктами 1 и 2 выше) и ядро ​​opt2 с реорганизацией данных работает примерно на полной пропускной способности (36 ГБ / с).

$ cat t574.cu
#include <stdio.h>
#include <stdlib.h>
#define THREADS_PER_BLOCK 128
#define VECTOR_SIZE 1024
#define NUM_DB_VEC 100000

typedef unsigned char BYTE;
typedef unsigned int UINT;
typedef unsigned int uint32_t;template<UINT blockSize> __global__ void reduction_sum_abs( BYTE* query_vector, BYTE* db_vector, uint32_t* result )
{
extern __shared__ UINT sum[];
UINT db_linear_index = (blockIdx.y*gridDim.x) + blockIdx.x ;
UINT i = threadIdx.x;

sum[threadIdx.x] = 0;

int* p_q_int = reinterpret_cast<int*>(query_vector);
int* p_db_int = reinterpret_cast<int*>(db_vector);

while( i < VECTOR_SIZE/4 ) {

/* memory transaction */
int q_int = p_q_int[i];
int db_int = p_db_int[db_linear_index*VECTOR_SIZE/4 + i];

uchar4 a0 = *reinterpret_cast<uchar4*>(&q_int);
uchar4 b0 = *reinterpret_cast<uchar4*>(&db_int);

/* sum of absolute difference */
sum[threadIdx.x] += abs( (int)a0.x - b0.x );
sum[threadIdx.x] += abs( (int)a0.y - b0.y );
sum[threadIdx.x] += abs( (int)a0.z - b0.z );
sum[threadIdx.x] += abs( (int)a0.w - b0.w );

i += THREADS_PER_BLOCK;

}

__syncthreads();

if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}
__syncthreads(); // **
/* reduce the final warp */
if ( threadIdx.x < 32 ) {
if ( blockSize >= 64 ) { sum[threadIdx.x] += sum[threadIdx.x + 32]; } __syncthreads();

if ( blockSize >= 32 ) { sum[threadIdx.x] += sum[threadIdx.x + 16]; } __syncthreads();

if ( blockSize >= 16 ) { sum[threadIdx.x] += sum[threadIdx.x + 8 ]; } __syncthreads();

if ( blockSize >= 8  ) { sum[threadIdx.x] += sum[threadIdx.x + 4 ]; } __syncthreads();

if ( blockSize >= 4  ) { sum[threadIdx.x] += sum[threadIdx.x + 2 ]; } __syncthreads();

if ( blockSize >= 2  ) { sum[threadIdx.x] += sum[threadIdx.x + 1 ]; } __syncthreads();

}/* copy the sum back to global */
if ( threadIdx.x == 0 ) {
result[db_linear_index] = sum[0];
}
}

__global__ void reduction_sum_abs_opt1( BYTE* query_vector, BYTE* db_vector, uint32_t* result )
{
__shared__ UINT sum[THREADS_PER_BLOCK];
UINT db_linear_index = (blockIdx.y*gridDim.x) + blockIdx.x ;
UINT i = threadIdx.x;

sum[threadIdx.x] = 0;

UINT* p_q_int = reinterpret_cast<UINT*>(query_vector);
UINT* p_db_int = reinterpret_cast<UINT*>(db_vector);

while( i < VECTOR_SIZE/4 ) {

/* memory transaction */
UINT q_int = p_q_int[i];
UINT db_int = p_db_int[db_linear_index*VECTOR_SIZE/4 + i];
sum[threadIdx.x] += __vsadu4(q_int, db_int);

i += THREADS_PER_BLOCK;

}
__syncthreads();
// this reduction assumes THREADS_PER_BLOCK = 128
if (threadIdx.x < 64) sum[threadIdx.x] += sum[threadIdx.x+64];
__syncthreads();

if ( threadIdx.x < 32 ) {
unsigned localSum = sum[threadIdx.x] + sum[threadIdx.x + 32];
for (int i = 16; i >= 1; i /= 2)
localSum = localSum + __shfl_xor(localSum, i);
if (threadIdx.x == 0) result[db_linear_index] = localSum;
}
}

__global__ void reduction_sum_abs_opt2( BYTE* query_vector, UINT* db_vector_cm, uint32_t* result)
{
__shared__ UINT qv[VECTOR_SIZE/4];
if (threadIdx.x < VECTOR_SIZE/4) qv[threadIdx.x] = *(reinterpret_cast<UINT *>(query_vector) + threadIdx.x);
__syncthreads();
int idx = threadIdx.x + blockDim.x*blockIdx.x;
while (idx < NUM_DB_VEC){
UINT sum = 0;
for (int i = 0; i < VECTOR_SIZE/4; i++)
sum += __vsadu4(qv[i], db_vector_cm[(i*NUM_DB_VEC)+idx]);
result[idx] = sum;
idx += gridDim.x*blockDim.x;}
}

unsigned long compute_host_result(BYTE *qvec, BYTE *db_vec){

unsigned long temp = 0;
for (int i =0; i < NUM_DB_VEC; i++)
for (int j = 0; j < VECTOR_SIZE; j++)
temp += (unsigned long) abs((int)qvec[j] - (int)db_vec[(i*VECTOR_SIZE)+j]);
return temp;
}

int main(){

float et;
cudaEvent_t start, stop;
BYTE *h_qvec, *d_qvec, *h_db_vec, *d_db_vec;
uint32_t *h_res, *d_res;
h_qvec =   (BYTE *)malloc(VECTOR_SIZE*sizeof(BYTE));
h_db_vec = (BYTE *)malloc(VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE));
h_res = (uint32_t *)malloc(NUM_DB_VEC*sizeof(uint32_t));
for (int i = 0; i < VECTOR_SIZE; i++){
h_qvec[i] = rand()%256;
for (int j = 0; j < NUM_DB_VEC; j++) h_db_vec[(j*VECTOR_SIZE)+i] = rand()%256;}
cudaMalloc(&d_qvec, VECTOR_SIZE*sizeof(BYTE));
cudaMalloc(&d_db_vec, VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE));
cudaMalloc(&d_res, NUM_DB_VEC*sizeof(uint32_t));
cudaMemcpy(d_qvec, h_qvec, VECTOR_SIZE*sizeof(BYTE), cudaMemcpyHostToDevice);
cudaMemcpy(d_db_vec, h_db_vec, VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE), cudaMemcpyHostToDevice);
cudaEventCreate(&start); cudaEventCreate(&stop);

// initial run

cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs<THREADS_PER_BLOCK><<<NUM_DB_VEC, THREADS_PER_BLOCK, THREADS_PER_BLOCK*sizeof(int)>>>(d_qvec, d_db_vec, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
unsigned long h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("1: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if (h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("1: mismatch!\n");

// optimized kernel 1
cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs_opt1<<<NUM_DB_VEC, THREADS_PER_BLOCK>>>(d_qvec, d_db_vec, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("2: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if(h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("2: mismatch!\n");

// convert db_vec to column-major storage for optimized kernel 2

UINT *h_db_vec_cm, *d_db_vec_cm;
h_db_vec_cm = (UINT *)malloc(NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT));
cudaMalloc(&d_db_vec_cm, NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT));
for (int i = 0; i < NUM_DB_VEC; i++)
for (int j = 0; j < VECTOR_SIZE/4; j++)
h_db_vec_cm[(j*NUM_DB_VEC)+i] = *(reinterpret_cast<UINT *>(h_db_vec + (i*VECTOR_SIZE))+j);
cudaMemcpy(d_db_vec_cm, h_db_vec_cm, NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT), cudaMemcpyHostToDevice);
cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs_opt2<<<64, 512>>>(d_qvec, d_db_vec_cm, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("3: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if(h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("3: mismatch!\n");

return 0;
}

$ nvcc -O3 -arch=sm_35 -o t574 t574.cu
$ ./run35 t574
1: et: 6.34ms, bw: 16.14GB/s
Success!
2: et: 4.16ms, bw: 24.61GB/s
Success!
3: et: 2.83ms, bw: 36.19GB/s
Success!
$

Несколько заметок:

  1. Приведенный выше код, в частности ваше ядро, должен быть скомпилирован для cc3.0 или выше, как я настроил тестовый пример. Это потому, что я создаю 100 000 блоков в одной сетке 1D, поэтому мы не можем запустить это как есть, например, на устройстве cc2.0.
  2. Возможна небольшая дополнительная настройка, особенно если она выполняется на разных устройствах для ядра opt2, путем изменения параметров сетки и блока. У меня они установлены на 64 и 512, и эти значения не должны быть критическими (хотя блок должен быть VECTOR_SIZE / 4 потока или больше), так как алгоритм использует цикл шага сетки, чтобы покрыть весь векторный набор. GT640 имеет только 2 SM, поэтому в этом случае 64 потоковых блоков достаточно, чтобы держать устройство занятым (вероятно, даже 32 будет в порядке). Вы можете изменить их, чтобы получить максимальную производительность на больших устройствах.
7

Другие решения

Одна вещь сразу привлекла мое внимание:

if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}

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

if ( threadIdx.x < 64 ) {
if ( blockSize >= 128 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}

Это позволило бы всем варпам, кроме первых двух, завершить их выполнение раньше.

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

/* reduce the final warp */
if ( threadIdx.x < 32 ) {
auto localSum = sum[threadIdx.x] + sum[threadIdx.x + 32]);
for (auto i = 16; i >= 1; i /= 2)
{
localSum = localSum + __shfl_xor(localSum, i);
}

if (threadIdx.x == 0) result[db_linear_index] = localSum;
}

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

Редактировать:
Также кажется, что вы без необходимости пишете в общую память четыре раза:

/* sum of absolute difference */
sum[threadIdx.x] += abs( (int)a0.x - b0.x );
sum[threadIdx.x] += abs( (int)a0.y - b0.y );
sum[threadIdx.x] += abs( (int)a0.z - b0.z );
sum[threadIdx.x] += abs( (int)a0.w - b0.w );

Почему бы просто не сделать следующее?

    /* sum of absolute difference */
sum[threadIdx.x] += abs( (int)a0.x - b0.x )
+ abs( (int)a0.y - b0.y )
+ abs( (int)a0.z - b0.z );
+ abs( (int)a0.w - b0.w );
1

По вопросам рекламы [email protected]