在 CUDA 中查找矩阵的最大值

Question

我刚开始使用 CUDA。现在我有一个问题。我有 N*N 矩阵，window 比例是 8x8。我想将这个矩阵细分为多个子矩阵并找到它的最大值。例如，如果我有 64*64 矩阵，那么我将有 8 个 8*8 比例的小矩阵，并找出 8 个最大值。最后我将所有最大值保存到新数组中，但它的顺序总是改变。我想找到让它们保持正确顺序的解决方案

__global__ void calculate_emax_kernel(float emap[],float emax[], int img_height, int img_width,int windows_size)
{
    int x_index = blockIdx.x*blockDim.x+threadIdx.x;
    int y_index = blockIdx.y*blockDim.y+threadIdx.y;

    int num_row_block = img_height/windows_size;
    int num_col_block = img_width/windows_size;
    __shared__ float window_elements[256];
    __shared__ int counter;
    __shared__ int emax_count;

    if (threadIdx.x == 0) emax_count = 0;
    __syncthreads();
    int index;
    int emax_idx = 0;


    if(y_index >= img_height|| x_index >= img_width) return;
    for(int i = 0; i < num_row_block; i++)
    {
        for(int j = 0; j < num_col_block; j++)
        {
            counter = 0;
            if(y_index >= i*windows_size && y_index < (i+1)*windows_size
                    && x_index >= j*windows_size && x_index < (j+1)*windows_size)
            {
                int idx = y_index*img_height + x_index;
                index = atomicAdd(&counter, 1);

                window_elements[index] = emap[idx];
                __syncthreads();


                // reduction
                unsigned int k = (windows_size*windows_size)/2;
                while(k != 0)
                {
                    if(index < k)
                    {
                        window_elements[index] = fmaxf(window_elements[index], window_elements[index+k]);

                    }
                    k /= 2;
                }
                if(index == 0)
                {
                    emax[i*num_row_block+j] = window_elements[index];
                }
            }
            __syncthreads();
        }
        __syncthreads();
    }
    __syncthreads();
}

这是我的配置

void construct_emax(float *input,float *output, int img_height, int img_width)
{
    int windows_size = 4;
    float * d_input, * d_output;
    cudaMalloc(&d_input, img_width*img_height*sizeof(float));
    cudaMalloc(&d_output, img_width*img_height*sizeof(float));

    cudaMemcpy(d_input, input, img_width*img_height*sizeof(float), cudaMemcpyHostToDevice);
    dim3 blocksize(16,16);
    dim3 gridsize;

    gridsize.x=(img_width+blocksize.x-1)/blocksize.x;
    gridsize.y=(img_height+blocksize.y-1)/blocksize.y;

    calculate_emax_kernel<<<gridsize,blocksize>>>(d_input,d_output,img_height,img_width,windows_size);

}

Answer 1

使用 CUDA，parallel reduction is tricky; segmented parallel reduction 更棘手。现在你在二维中做，你的 segment/window 比线程块小。

对于大 window 尺寸，我认为这不是问题。您可以使用一个线程块来减少一个 window。例如，如果您有一个 16x16 window，您可以简单地使用 16x16 线程块。如果您有更大的 window 大小，例如 64x64，您仍然可以使用 16x16 线程块。首先在数据加载期间将 64x64 window 减少到 16x16 元素，然后在线程块内减少到 1 个标量。

对于小于块大小的 window 大小，您必须减少每个线程块的多个 windows 以获得更高的性能。您可以使用当前的 block/grid 配置，其中每个 256 线程块 (16x16) 负责 16 个 4x4 windows。但这不是最佳方案，因为每个 32 线程环绕都分为两部分 (2x16)。这对 coalesced global memory access 不利，并且很难将 2x16 扭曲映射到一个或多个 4x4 windows 以实现有效的并行缩减。

或者，我建议您使用具有 256 个线程的一维线程块。每 m 个线程减少一个 mxm window。然后你可以使用二维网格覆盖整个图像。

const int m = window_size;
dim3 blocksize(256);
dim3 gridsize((img_width+255)/256, (img_height+m-1)/m);

在内核函数中，你可以

在全局数据加载期间将每个 mxm window 减少为 1xm 向量；
使用树减少方法将 1xm 向量减少为标量。

下面的代码是一个概念性演示，它在 m 是 2 和 m <= 32 的幂时起作用。您可以进一步修改它以实现任意 m 和更好的边界检查。

#include <assert.h>
#include <cuda.h>
#include <thrust/device_vector.h>

__global__ void calculate_emax_kernel(const float* input, float* output,
                                      int height, int width, int win_size,
                                      int out_width) {
  const int tid = threadIdx.x;
  const int i = blockIdx.y * win_size;
  const int j = blockIdx.x * 256 + tid;
  const int win_id = j % win_size;

  __shared__ float smax[256];

  float tmax = -1e20;
  if (j < width) {
    for (int tile = 0; tile < win_size; tile++) {
      if (i + tile < height) {
        tmax = max(tmax, input[(i + tile) * width + j]);
      }
    }
  }
  smax[tid] = tmax;
  for (int shift = win_size / 2; shift > 0; shift /= 2) {
    if (win_id < shift) {
      smax[tid] = max(smax[tid], smax[tid + shift]);
    }
  }
  if (win_id == 0 && j < width) {
    output[blockIdx.y * out_width + (j / win_size)] = smax[tid];
  }
}

int main() {
  const int height = 1024;
  const int width = 1024;
  const int m = 4;
  thrust::device_vector<float> in(height * width);
  thrust::device_vector<float> out(
      ((height + m - 1) / m) * ((width + m - 1) / m));

  dim3 blocksize(256);
  dim3 gridsize((width + 255) / 256, (height + m - 1) / m);

  assert(m == 2 || m == 4 || m == 8 || m == 16 || m == 32);
  calculate_emax_kernel<<<gridsize, blocksize>>>(
      thrust::raw_pointer_cast(in.data()),
      thrust::raw_pointer_cast(out.data()),
      height, width, m, (width + m - 1) / m);

  return 0;
}

Answer 2

如果您愿意使用图书馆，请提供几点建议：

使用 NPP，原语集（来自 nvidia） https://docs.nvidia.com/cuda/npp/group__image__filter__max.html
一个较低级别的库，用于其他 reduce 操作和更细化的硬件使用方式（来自 nvidia / nvlabs） http://nvlabs.github.io/cub/

在 CUDA 中查找矩阵的最大值

Find max of matrix in CUDA

cuda

reduction