Ядро CUDA не возвращает значения

Question

Я работаю с сервером, на котором установлено несколько графических процессоров.Я использую openMP для запуска ядра через несколько графических процессоров одновременно.Проблема, которую я вижу, состоит в том, что ядро, которое я использую, похоже, не обновляет значения в векторах упорных устройств, которым оно передается.Приведенный ниже код должен вывести значение 1 для всех элементов в векторах устройства, но вместо этого вывести значение 0. Код компилируется и запускается и показывает мне, что ядро ​​выполняется успешно.

Я не понимаю, почему этокод работает не так, как ожидалось.

#include <iostream>
#include <cmath>
#include <omp.h>
#include <vector>
#include <thrust/host_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/device_malloc.h>
#include <thrust/device_free.h>
#include <thrust/device_vector.h>



using namespace::std;


const long N_R1 = 100;
const long N_R2 = 100;


__global__ void kernel(long* ND, long* NR1, 
                       float* a, float* b, float* c, float* d)

{
    // Calculate Global index (Generic 3D block, 3D thread)
    long idx = ( blockIdx.x + blockIdx.y * gridDim.x * gridDim.y * blockIdx.z )
              * ( threadIdx.z * ( blockDim.x*blockDim.y ) ) + threadIdx.y 
              * blockDim.x + threadIdx.x;

    //Values correspond to 2D array limits
    long idxR1 = idx / ND[0];
    long idxR2 = idx % ND[0];

    if(idxR1 >= NR1[0] || idxR2 >= ND[0])
    {
        return;
    }

        a[idx] =1.0;
        b[idx] =1.0;
        c[idx] =1.0;
        d[idx] =1.0;

}


void kernel_wrapper()
{
    // GPU Count
    int num_gpus = 0;
    cudaGetDeviceCount(&num_gpus);
    omp_set_num_threads(num_gpus);

    //Calculate Dimensioning
    long D_total = N_R1 * N_R2;
    //Region 1 coordinates are loaded on to each GPU
    //Region 2 coordinates are divided up onto GPUs
    long R2_stride = ceil(float(N_R2)/float(num_gpus));

    //Distance arrays need to be split longo whole sections of region 1. 
    //(Distances size = N_R1 * N_R2) subset of distance size needs to be N_R1
    long D_stride = R2_stride * N_R1;


#pragma omp parallel
    {

        // Get CPU thread number
        long cpu_thread_id = omp_get_thread_num();

        cudaSetDevice(cpu_thread_id);

        // Set up Local Arrays for distance and potential
        // Step 1: Calculate rough Array Limits
        // If array spaces divide evenly between threads then beginnings and endings can be calculated below
        long R2_begin = cpu_thread_id * R2_stride;
        long D_begin  = cpu_thread_id * D_stride;

        long R2_end = R2_begin + R2_stride;
        long D_end  = D_begin + D_stride;

        // Step 2: Check Ends are not out of bounds
        //         The last thread in the calculation is likely to have array sizings that are out of bounds
        //         if this is the case then the ends need to be clipped:
        if(R2_end >= N_R2)
        {
            R2_end = N_R2;
        }
        if(D_end >= D_total)
        {
            D_end = D_total;
        }

        // Local aray sizes are (end - begin)
        long l_R2 = R2_end - R2_begin;
        long l_D     = D_end - D_begin;

        float zero = 0.0;
        // Create Region 2 potential components
        thrust::host_vector<float > a(l_D,zero);
        thrust::host_vector<float > b(l_D,zero);
        thrust::host_vector<float > c(l_D,zero);
        thrust::host_vector<float > d(l_D,zero);

        long* p_NR1;
        long nr1 = N_R1;
        cudaMalloc( (void**)&p_NR1, sizeof(long) );
        cudaMemcpy( p_NR1, &nr1, sizeof(long), cudaMemcpyHostToDevice);

        long* p_NR2;
        cudaMalloc( (void**)&p_NR2, sizeof(long) );
        cudaMemcpy( p_NR2, &l_D, sizeof(long), cudaMemcpyHostToDevice);

        //Generate Device Side Data for region 2 potential components
        thrust::device_vector< float > d_a = a;
        thrust::device_vector< float > d_b = b;
        thrust::device_vector< float > d_c = c;
        thrust::device_vector< float > d_d = d;
        // Generate pointers to Device Side Data for region 2 potential components
        float* p_a = thrust::raw_pointer_cast(d_a.data());
        float* p_b = thrust::raw_pointer_cast(d_b.data());
        float* p_c = thrust::raw_pointer_cast(d_c.data());
        float* p_d = thrust::raw_pointer_cast(d_d.data());

        dim3 blocks = N_R1;
        dim3 threads = l_R2;
        kernel<<<blocks,threads>>>(p_NR2, p_NR1,
                                   p_a, p_b, p_c, p_d);
        cudaDeviceSynchronize();
        if(cudaGetLastError() == cudaSuccess)
        {
            cout << "Kernel Successful!" << cudaGetErrorString(cudaGetLastError()) << endl;
            cin.ignore(1);
        }

        a = d_a;
        b = d_b;
        c = d_c;
        d = d_d;

        for(long j = 0; j != a.size(); j++)
        {
            cout << "a[" << j << "] = " << a[j] << endl;
        }
        for(long j = 0; j != b.size(); j++)
        {
            cout << "b[" << j << "] = " << b[j] << endl;
        }
        for(long j = 0; j != c.size(); j++)
        {
            cout << "c[" << j << "] = " << c[j] << endl;
        }
        for(long j = 0; j != c.size(); j++)
        {
            cout << "c[" << j << "] = " << c[j] << endl;
        }
}
        cin.ignore(1);
}

int main()
{

    kernel_wrapper();

    return 0;
}

Любая помощь будет принята с благодарностью.

Answer 1

Некоторые выходные значения устанавливаются на 1, некоторые нет.Проблема заключается в следующем:

// Calculate Global index (Generic 3D block, 3D thread)
long idx = ( blockIdx.x + blockIdx.y * gridDim.x * gridDim.y * blockIdx.z )
          * ( threadIdx.z * ( blockDim.x*blockDim.y ) ) + threadIdx.y 
          * blockDim.x + threadIdx.x;

Это не то, что я бы назвал правильным общим преобразованием трехмерной сетки / блока в глобально уникальный одномерный индекс, который, как я полагаю, является вашим намерением.Давайте просто выберем один пример, чтобы доказать, что он сломан.Предположим, вы запускаете 1D сетку из 1D блоков (что вы и делаете).Тогда все переменные (block, thread) Idx.y и .z будут равны нулю.Только blockIdx.x и threadIdx.x могут принимать ненулевые значения в этой конфигурации запуска.

В этом случае ваше выражение уменьшается до:

// Calculate Global index (Generic 3D block, 3D thread)
long idx = ( blockIdx.x + 0 * gridDim.x * gridDim.y * 0 )
          * ( 0 * ( blockDim.x*blockDim.y ) ) + 0 
          * blockDim.x + threadIdx.x;

, то есть оно уменьшается до:

long idx = threadIdx.x;

Итак, первые (размером с блок) элементы ваших массивов (a, b, c, d) установлены правильно, остальные - нет.Поскольку threadIdx.x не является уникальным от одного блока к другому, это не правильный глобально уникальный идентификатор потока, и поэтому каждый блок записывает одинаковые выходные местоположения, а не каждый заботится о отдельномчасть массива.

Итак, что такое возможное (правильное) общее преобразование индекса 3D-в-1D?

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

Когда я заменяю ваши вычисления в ядре idx с этим кодом ваше ядро ​​заполняет все записи массива 1.0 согласно моему тестированию.