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

Question

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

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

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

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

blockSize и совместно используемая память установлены на 128 и 128*sizeof(int), 128 определяется как 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

Answer 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 или выше, как я настроил тестовый пример. Это потому, что я создаю 100000 блоков в одной одномерной сетке, поэтому мы не можем запустить это как есть, например, на устройстве cc2.0.
2. Возможна небольшая дополнительная настройка, особенно если она выполняется на разных устройствах для ядра opt2, путем изменения параметров сетки и блока. У меня они установлены на 64 и 512, и эти значения не должны быть критическими (хотя блок должен быть VECTOR_SIZE/4 или более), так как алгоритм использует цикл шага сетки, чтобы покрыть весь векторный набор. GT640 имеет только 2 SM, поэтому в этом случае 64 потоковых блоков достаточно, чтобы держать устройство занятым (вероятно, даже 32 будет в порядке). Вы можете изменить их, чтобы получить максимальную производительность на больших устройствах.

Answer 2

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

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 );