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

Это хорошая практика для таких вопросов, чтобы обеспечить полный код, который кто-то может скомпилировать и запустить без добавления чего-либо или изменения чего-либо. Вообще говоря, 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 будет в порядке). Вы можете изменить их, чтобы получить максимальную производительность на больших устройствах.