Возвращение минимального и максимального элементов массива в CUDA

Question

Я выполняю некоторые манипуляции/вычисления массивов в CUDA (через Cudafy.NET library, хотя меня одинаково интересуют методы CUDA/C++), и вам нужно вычислить минимальные и максимальные значения, которые находятся в массиве. Одно из ядер выглядит следующим образом:

[Cudafy] 
    public static void UpdateEz(GThread thread, float time, float ca, float cb, float[,] hx, float[,] hy, float[,] ez) 
    { 
     var i = thread.blockIdx.x; 
     var j = thread.blockIdx.y; 

     if (i > 0 && i < ez.GetLength(0) - 1 && j > 0 && j < ez.GetLength(1) - 1) 
      ez[i, j] = 
       ca * ez[i, j] 
       + cb * (hx[i, j] - hx[i - 1, j]) 
       + cb * (hy[i, j - 1] - hy[i, j]) 
       ; 
    } 

Я хотел бы сделать что-то вроде этого:

[Cudafy] 
    public static void UpdateEz(GThread thread, float time, float ca, float cb, float[,] hx, float[,] hy, float[,] ez, out float min, out float max) 
    { 
     var i = thread.blockIdx.x; 
     var j = thread.blockIdx.y; 

     min = float.MaxValue; 
     max = float.MinValue; 

     if (i > 0 && i < ez.GetLength(0) - 1 && j > 0 && j < ez.GetLength(1) - 1) 
     { 
      ez[i, j] = 
       ca * ez[i, j] 
       + cb * (hx[i, j] - hx[i - 1, j]) 
       + cb * (hy[i, j - 1] - hy[i, j]) 
       ; 

      min = Math.Min(ez[i, j], min); 
      max = Math.Max(ez[i, j], max); 

     } 
    } 

кто знает удобный способ вернуть минимальное и максимальное значения (для всего массива , а не только на поток или блок)?

Answer 1

Основываясь на вашем комментарии к вашему вопросу, вы пытались найти значения max и min при их вычислении; в то время как это возможно, это не самый эффективный. Если вы настроитесь на это, тогда вы можете провести атомарное сравнение с каким-то глобальным минимумом и глобальным максимумом, а недостатком будет то, что каждый поток будет сериализован, что, вероятно, станет значительным узким местом.

Для более канонического подхода к нахождению максимума или минимума в массиве с помощью сокращения, вы можете сделать что-то вдоль линий:

#define MAX_NEG ... //some small number 

template <typename T, int BLKSZ> __global__ 
void cu_max_reduce(const T* d_data, const int d_len, T* max_val) 
{ 
    volatile __shared__ T smem[BLKSZ]; 

    const int tid = threadIdx.x; 
    const int bid = blockIdx.x; 
     //starting index for each block to begin loading the input data into shared memory 
    const int bid_sidx = bid*BLKSZ; 

    //load the input data to smem, with padding if needed. each thread handles 2 elements 
    #pragma unroll 
    for (int i = 0; i < 2; i++) 
    { 
       //get the index for the thread to load into shared memory 
     const int tid_idx = 2*tid + i; 
     const int ld_idx = bid_sidx + tid_idx; 
     if(ld_idx < (bid+1)*BLKSZ && ld_idx < d_len) 
      smem[tid_idx] = d_data[ld_idx]; 
     else 
      smem[tid_idx] = MAX_NEG; 

     __syncthreads(); 
    } 

    //run the reduction per-block 
    for (unsigned int stride = BLKSZ/2; stride > 0; stride >>= 1) 
    { 
     if(tid < stride) 
     { 
      smem[tid] = ((smem[tid] > smem[tid + stride]) ? smem[tid]:smem[tid + stride]); 
     } 
     __syncthreads(); 
    } 

    //write the per-block result out from shared memory to global memory 
    max_val[bid] = smem[0]; 
} 


//assume we have d_data as a device pointer with our data, of length data_len 
template <typename T> __host__ 
T cu_find_max(const T* d_data, const int data_len) 
{ 
    //in your host code, invoke the kernel with something along the lines of: 
    const int thread_per_block = 16; 
    const int elem_per_thread = 2; 
    const int BLKSZ = elem_per_thread*thread_per_block; //number of elements to process per block 
    const int blocks_per_grid = ceil((float)data_len/(BLKSZ)); 

    dim3 block_dim(thread_per_block, 1, 1); 
    dim3 grid_dim(blocks_per_grid, 1, 1); 

    T *d_max; 
    cudaMalloc((void **)&d_max, sizeof(T)*blocks_per_grid); 

    cu_max_reduce <T, BLKSZ> <<<grid_dim, block_dim>>> (d_data, data_len, d_max); 

    //etc.... 
} 

Найдет поблочно максимальное значение. Вы можете запустить его снова на своем выходе (например, с d_max в качестве входных данных и с обновленными параметрами запуска) на 1 блоке, чтобы найти глобальный максимум - запуск его в многопроходном режиме, как это необходимо, если ваш набор данных слишком велик (в этом случае, выше 2 * 4096 элементов, так как у нас есть каждый поток процесса 2 элемента, хотя вы можете просто обрабатывать больше элементов на поток, чтобы увеличить это).

Следует отметить, что это не особенно эффективно (вы хотите использовать более интеллектуальный шаг при загрузке разделяемой памяти, чтобы избежать конфликтов с банками), и я не уверен на 100%, что это правильно (это сработало на нескольких небольших тестовых ящиках, которые я пробовал), но я попытался написать его для максимальной ясности. Кроме того, не забудьте добавить код проверки ошибок, чтобы убедиться, что ваши вызовы CUDA успешно завершены, я оставил их здесь, чтобы сохранить его коротким (er).

Я также должен направить вас к более углубленной документации; вы можете взглянуть на сокращение выборки CUDA на http://docs.nvidia.com/cuda/cuda-samples/index.html, хотя это не делает вычисление min/max, это та же самая общая идея (и более эффективная). Кроме того, если вы ищете для простоты, вы можете просто хотите использовать функции Thrust в thrust::max_element и thrust::min_element и документацию по адресу: thrust.github.com/doc/group__extrema.html

Answer 2

Вы можете разработать собственный алгоритм min/max с использованием метода divide and conquer.

Если у вас есть возможность использовать npp, эта функция может быть полезна: nppsMinMax_32f.

Answer 3

Если вы пишете электромагнитную волну симулятор и не хотите изобретать колесо, вы можете использовать thrust::minmax_element. Ниже я представляю простой пример того, как его использовать. Добавьте свою собственную проверку CUDA.

#include <stdio.h> 

#include <cuda_runtime_api.h> 

#include <thrust\pair.h> 
#include <thrust\device_vector.h> 
#include <thrust\extrema.h> 

int main() 
{ 
    const int N = 5; 

    const float h_a[N] = { 3., 21., -2., 4., 5. }; 

    float *d_a;  cudaMalloc(&d_a, N * sizeof(float)); 
    cudaMemcpy(d_a, h_a, N * sizeof(float), cudaMemcpyHostToDevice); 

    float minel, maxel; 
    thrust::pair<thrust::device_ptr<float>, thrust::device_ptr<float>> tuple; 
    tuple = thrust::minmax_element(thrust::device_pointer_cast(d_a), thrust::device_pointer_cast(d_a) + N); 
    minel = tuple.first[0]; 
    maxel = tuple.second[0]; 

    printf("minelement %f - maxelement %f\n", minel, maxel); 

    return 0; 
}