如何实现一维模板代码的推力版本?

How to implement a Thrust version of a 1D stencil code?

基本上,是否可以使用纯 Thrust 实现如下所示的 1D 模板内核?我希望这个实现尽可能高效,这意味着 Thrust 应该以某种方式知道对相同元素的多重访问,并且需要使用共享内存访问。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#define BLOCK_SIZE 128
#define RADIUS 8
#define SIZE 1024*1024*8
const dim3 DimBlock(BLOCK_SIZE);
const dim3 DimGrid(SIZE/BLOCK_SIZE);

__global__ void stencil_1d(const int * in, int *out) {
    __shared__ int temp[BLOCK_SIZE + 2 * RADIUS];
    int gindex = threadIdx.x + blockIdx.x * blockDim.x;
    int lindex = threadIdx.x + RADIUS;

    // Read input elements into shared memory
    if( gindex < SIZE ) temp[lindex] = in[gindex]; else temp[lindex] = 0;
    if (threadIdx.x < RADIUS) {
    if(gindex - RADIUS>=0 )temp[lindex - RADIUS] = in[gindex - RADIUS]; else  temp[lindex - RADIUS] = 0;
    if(gindex + BLOCK_SIZE < SIZE )  temp[lindex + BLOCK_SIZE] = in[gindex + BLOCK_SIZE]; else  temp[lindex + BLOCK_SIZE] = 0;
    }

    // Synchronize (ensure all the data is available)
    __syncthreads();

    // Apply the stencil
    int result = 0;
    for (int offset = -RADIUS; offset <= RADIUS; offset++)
    if( gindex < SIZE ) result += temp[lindex + offset];

    // Store the result
    if( gindex < SIZE ) out[gindex] = result;
}

int main()
{
    cudaError_t cudaStat;
    thrust::device_vector<int> dev_vec_inp(SIZE,1);
    thrust::device_vector<int> dev_vec_out(SIZE);
    try 
    {
        stencil_1d<<< DimGrid, DimBlock >>>(thrust::raw_pointer_cast(dev_vec_inp.data()) , thrust::raw_pointer_cast(dev_vec_out.data()));
        cudaStat =  cudaGetLastError(); 
        if (cudaStat != cudaSuccess)
        throw cudaGetErrorString(cudaStat);
        else std::cout<<"1D stencil has been executed successfully" << std:: endl;
    }   
    catch(const char* e)
    {
        std::cout<< e;
    }
}

对于这种特殊类型的模板操作 (+),您可以使用前缀求和方法。您将首先执行前缀和,然后使用模板半径从模板右端 window.

减去模板左端 window

这是一个粗略的草图:

$ cat t1777.cu
#include <thrust/device_vector.h>
#include <thrust/scan.h>
#include <thrust/copy.h>
#include <thrust/transform.h>
#include <thrust/host_vector.h>
#include <iostream>

const int ds = 20;
const int stencil_radius = 7;
using namespace thrust::placeholders;
typedef int mt;

int main(){

  thrust::device_vector<mt> data(ds, 1);
  thrust::device_vector<mt> result(ds-(2*stencil_radius));
  thrust::inclusive_scan(data.begin(), data.end(), data.begin());
  thrust::transform(data.begin(), data.end()-(2*stencil_radius),data.begin()+(2*stencil_radius), result.begin(), _2-_1);
  thrust::host_vector<mt> h_result = result;
  thrust::copy(h_result.begin(), h_result.end(), std::ostream_iterator<mt>(std::cout, ","));
  std::cout << std::endl;
}
$ nvcc -o t1777 t1777.cu
$ ./t1777
14,14,14,14,14,14,
$

如果您 运行 在分析器下执行上述代码,您可以确定其中一个为前缀 sum op 推力启动的内核确实使用了共享内存。

我真的没有尝试创建一个合适的模板 window 由左边的半径、右边的半径和中心的模板位置组成,但这应该只需要稍微修改一下我上面显示的内容。

我并没有声称这比“纯 CUDA”方法更好或更快,但它似乎是一种使用“纯推力”方法执行它的方法,推力利用共享内存。

我并不是说此代码 defect-free 或适合任何特定用途。提供它只是为了展示一个想法。使用它需要您自担风险。