Невозможно успешно реализовать параллельное сокращение для графического процессора muti-CUDA

Во-первых, хорошей практикой является использование правильной проверки ошибок CUDA.

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

Если мы начнем с этого, то blocksPerGrid будет «умножено» (т. е. увеличено) на количество графических процессоров, чтобы охватить весь размер вашей проблемы. Давайте пройдемся по вашему коду и «согласуем» расчеты. Например, для двух графических процессоров, размера вектора 1048576 и 512 потоков на блок мы ожидаем, что blocksPerGrid будет равно 1024, потому что 2x1024x512 = 1048576. Это согласуется с вашими расчетами самого blocksPerGrid и вызова ядра.

1. Это неверно:

long sb = blocksPerGrid*sizeof(float);
...
h_C = (float*)malloc(sb);

Хранилище хоста для результата должно соответствовать (по крайней мере) размеру проблемы. Это должен быть один float элемент на блок, умноженный на количество графических процессоров. Но sb — это размер хранилища на один графический процессор. Нам нужно умножить его на количество графических процессоров при расчете необходимого размера для h_C.

2. Это неверно:

  cudaMalloc((void**)&d_C, sb/NGPU);

sb — это уже размер хранилища на один графический процессор согласно вашим расчетам blocksPerGrid. Не следует снова делить его на NGPU. Когда вы это сделаете, у вас теперь будут блоки потоков в каждом графическом процессоре, которые пытаются записать результаты в несуществующую память, и ваше ядро ​​будет выполнять недопустимое поведение. Учитывая достаточно большую проблему и/или использование compute-sanitizer, вы наверняка засвидетельствуете это при правильной проверке ошибок CUDA, о которой я упоминал.

3. Это неверно:

  cudaMemcpy(h_C+blocksPerGrid/NGPU*cpu_thread_id, d_C, sb/NGPU, cudaMemcpyDeviceToHost);

по причинам, которые мы уже рассмотрели. h_C необходимо покрыть весь размер задачи, а размер проблемы на каждый графический процессор уже покрыт blocksPerGrid. Его не следует далее делить на NGPU, а sb — это уже масштабирование на GPU, его не следует дальше делить на NGPU.

4. Это неверно:

  for (int i = 0; i < blocksPerGrid; i++) {
      h_G += (double) h_C[i];
  }

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

5. У вас неправильное расположение расчета h_G. Мы требуем, чтобы все потоки OMP завершили свою работу, прежде чем h_G будет вычислен результат. Следовательно, этот расчет должен выполняться после закрытия параллельной области OMP, чтобы гарантировать, что все потоки обновили свою часть h_C.

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

$ cat t3.cu
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>          // header for OpenMP
#include <cuda_runtime.h>

// Variables
float* h_A;   // host vectors
float* h_B;
float* h_C;
float* h_D;

// Functions
void RandomInit(float*, int);

// Device code
__global__ void VecAdd(const float* A, const float* B, float* C, int N)
{
    extern __shared__ float cache[];

    int i = blockDim.x * blockIdx.x + threadIdx.x;
    int cacheIndex = threadIdx.x;

    float temp = 0.0;  // register for each thread
    while (i < N) {
        temp += A[i]*B[i];
        i += blockDim.x*gridDim.x;
    }

    cache[cacheIndex] = temp;   // set the cache value

    __syncthreads();

    // perform parallel reduction, threadsPerBlock must be 2^m

    int ib = blockDim.x/2;
    while (ib != 0) {
      if (cacheIndex < ib)
        cache[cacheIndex] += cache[cacheIndex + ib];

      __syncthreads();

      ib /=2;
    }

    if (cacheIndex == 0)
      C[blockIdx.x] = cache[0];
}

// Host code

int main(void)
{
    printf("\n");
    printf("Vector Dot Product with multiple GPUs \n");
    int N, NGPU, cpu_thread_id=0;
    int *Dev;
    long mem = 1024*1024*1024;     // 4 Giga for float data type.

    printf("Enter the number of GPUs: ");
    //scanf("%d", &NGPU);
    NGPU = 2;
    printf("%d\n", NGPU);
    Dev = (int *)malloc(sizeof(int)*NGPU);

    int numDev = 0;
    printf("GPU device number: ");
    for(int i = 0; i < NGPU; i++) {
      //scanf("%d", &Dev[i]);
      Dev[i] = i;
      printf("%d ",Dev[i]);
      numDev++;
//      if (getchar() == '\n') break;
    }
    printf("\n");
    if (numDev != NGPU) {
      fprintf(stderr,"Should input %d GPU device numbers\n", NGPU);
      exit(1);
    }

    printf("Enter the size of the vectors: ");
    //scanf("%d", &N);
    N = 1048576;
    printf("%d\n", N);
    if (3*N > mem) {
        printf("The size of these 3 vectors cannot be fitted into 4 Gbyte\n");
        exit(1);
    }
    long size = N*sizeof(float);

    // Set the sizes of threads and blocks
    int threadsPerBlock;
    printf("Enter the number of threads per block: ");
    //scanf("%d", &threadsPerBlock);
    threadsPerBlock = 512;
    printf("%d\n", threadsPerBlock);
    if (threadsPerBlock > 1024) {
      printf("The number of threads per block must be less than 1024 ! \n");
      exit(1);
    }
    int blocksPerGrid = (N + threadsPerBlock*NGPU - 1) / (threadsPerBlock*NGPU);
    printf("The number of blocks is %d\n", blocksPerGrid);
    if (blocksPerGrid > 2147483647) {
      printf("The number of blocks must be less than 2147483647 ! \n");
      exit(1);
    }
    long sb = blocksPerGrid*sizeof(float);
    long sm = threadsPerBlock*sizeof(float);
    // Allocate input vectors h_A and h_B in host memory
    h_A = (float*)malloc(size);
    h_B = (float*)malloc(size);
    h_C = (float*)malloc(sb*NGPU);
    if (! h_A || ! h_B || ! h_C) {
      printf("!!! Not enough memory.\n");
      exit(1);
    }

    // Initialize input vectors

    RandomInit(h_A, N);
    RandomInit(h_B, N);

    // declare cuda event for timer
    cudaEvent_t start, stop;
//    cudaEventCreate(&start);    // events must be created after devices are set
//    cudaEventCreate(&stop);

    float Intime,gputime,Outime;
    double h_G = 0.0;
    omp_set_num_threads(NGPU);

    #pragma omp parallel private(cpu_thread_id)
    {
      float *d_A, *d_B, *d_C;
      cpu_thread_id = omp_get_thread_num();
      cudaSetDevice(Dev[cpu_thread_id]);
    //  cudaSetDevice(cpu_thread_id);

      // start the timer
      if (cpu_thread_id == 0) {
        cudaEventCreate(&start);
        cudaEventCreate(&stop);
        cudaEventRecord(start,0);
      }

      // Allocate vectors in device memory
      cudaMalloc((void**)&d_A, size/NGPU);
      cudaMalloc((void**)&d_B, size/NGPU);
      cudaMalloc((void**)&d_C, sb);

      // Copy vectors from host memory to device memory
      cudaMemcpy(d_A, h_A+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
      cudaMemcpy(d_B, h_B+N/NGPU*cpu_thread_id, size/NGPU, cudaMemcpyHostToDevice);
      #pragma omp barrier

        // stop the timer
      if (cpu_thread_id == 0) {
              cudaEventRecord(stop,0);
              cudaEventSynchronize(stop);
              cudaEventElapsedTime( &Intime, start, stop);
              printf("Data input time for GPU: %f (ms) \n",Intime);
      }

        // start the timer
      if (cpu_thread_id == 0) cudaEventRecord(start,0);

        VecAdd<<<blocksPerGrid, threadsPerBlock, sm>>>(d_A, d_B, d_C, N/NGPU);

        cudaDeviceSynchronize();

        // stop the timer

      if (cpu_thread_id == 0) {
              cudaEventRecord(stop,0);
              cudaEventSynchronize(stop);
              cudaEventElapsedTime( &gputime, start, stop);
              printf("Processing time for GPU: %f (ms) \n",gputime);
              printf("GPU Gflops: %f\n",3*N/(1000000.0*gputime));
      }

        // Copy result from device memory to host memory
        // h_C contains the result in host memory

        // start the timer
      if (cpu_thread_id == 0) cudaEventRecord(start,0);

      cudaMemcpy(h_C+blocksPerGrid*cpu_thread_id, d_C, sb, cudaMemcpyDeviceToHost);

      cudaFree(d_A);
      cudaFree(d_B);
      cudaFree(d_C);
      // stop the timer

      if (cpu_thread_id == 0) {
              cudaEventRecord(stop,0);
              cudaEventSynchronize(stop);
              cudaEventElapsedTime( &Outime, start, stop);
              printf("Data output time for GPU: %f (ms) \n",Outime);
      }
    }
    //compute the solution
    for (int i = 0; i < blocksPerGrid*NGPU; i++) {
          h_G += (double) h_C[i];
    }

    float gputime_tot;
    gputime_tot = Intime + gputime + Outime;
    printf("Total time for GPU: %f (ms) \n",gputime_tot);

    // start the timer
    cudaEventRecord(start,0);

    double h_D = 0.0;     // compute the reference solution
    for (int i = 0; i < N; ++i)
        h_D += (double) h_A[i]*h_B[i];

//        h_D[i] = 1.0/cos(h_A[i]) + 1.0/cos(h_B[i]);

    // stop the timer
    cudaEventRecord(stop,0);
    cudaEventSynchronize(stop);

    float cputime;
    cudaEventElapsedTime( &cputime, start, stop);
    printf("Processing time for CPU: %f (ms) \n",cputime);
    printf("CPU Gflops: %f\n",3*N/(1000000.0*cputime));
    printf("Speed up of GPU = %f\n", cputime/gputime_tot);

    // Destroy timer
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    // check result
    printf("Check result:\n");
    // for (int i = 0; i < N; ++i) {
    //     diff = abs(h_D[i] - h_C[i]);
    //     sum += diff*diff;
    // }
    double diff = abs( (h_D - h_G)/h_D );
    printf("|(h_G - h_D)/h_D|=%20.15e\n",diff);
    printf("h_G =%20.15e\n",h_G);
    printf("h_D =%20.15e\n",h_D);

    for (int i=0; i < NGPU; i++) {
        cudaSetDevice(i);
        cudaDeviceReset();
    }

    return 0;
}


// Allocates an array with random float entries.
void RandomInit(float* data, int n)
{
    for (int i = 0; i < n; ++i)
        data[i] = 1.0f; //rand() / (float)RAND_MAX;
}
$ nvcc -o t3 t3.cu -Xcompiler -fopenmp
$ compute-sanitizer ./t3
========= COMPUTE-SANITIZER

Vector Dot Product with multiple GPUs
Enter the number of GPUs: 2
GPU device number: 0 1
Enter the size of the vectors: 1048576
Enter the number of threads per block: 512
The number of blocks is 1024
Data input time for GPU: 2.405280 (ms)
Processing time for GPU: 8.202272 (ms)
GPU Gflops: 0.383519
Data output time for GPU: 0.429728 (ms)
Total time for GPU: 11.037280 (ms)
Processing time for CPU: 2.361696 (ms)
CPU Gflops: 1.331978
Speed up of GPU = 0.213974
Check result:
|(h_G - h_D)/h_D|=0.000000000000000e+00
h_G =1.048576000000000e+06
h_D =1.048576000000000e+06
========= ERROR SUMMARY: 0 errors
$

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