使用Cuda并行实现计算阵列中连续子序列的总和

Question

可以考虑以下数组： tab = [80,12,14,5,70,9,26,30,8,12,16 ，15 我要计算使用CUDA大小为4的所有可能序列的总和： 例如：

S1=80+12+14+5=111 
S2=12+14+5+70 =101 
S3=14+5+70+9 =98 
.... 

您有一个高效的想法parallise使用CUDA此任务。前面的表只是一个例子，在我的情况下，我将使用巨大的一个。

Answer 1

我们可以使用推力在单个操作（thrust::transform）中执行此操作。在CUDA中，这可以被认为是相当简单的一维模板操作。

可以在幻灯片49-58上找到here的1-D模板操作的很好描述。

这实际上是一个简化的情况，因为模板宽度是4并且它只位于中心点的一侧。

这里的一个工作实例比较2点的方法：

$ cat t88.cu 
#include <thrust/device_vector.h> 
#include <thrust/transform.h> 
#include <thrust/iterator/zip_iterator.h> 
#include <thrust/copy.h> 
#include <iostream> 

const int nTPB=256; 
typedef float mytype; 
const int ds = 1048576*32; 

struct sum4 
{ 
    template <typename T> 
    __host__ __device__ 
    mytype operator()(const T t){ 
    return thrust::get<0>(t) + thrust::get<1>(t) + thrust::get<2>(t) + thrust::get<3>(t); 
    } 
}; 

template <typename T> 
__global__ void sum4kernel(const T * __restrict__ in, T * __restrict__ out, const unsigned dsize) 
{ 

    __shared__ T sdata[nTPB+3]; 
    unsigned idx = threadIdx.x+blockDim.x*blockIdx.x; 
    if (idx < dsize) sdata[threadIdx.x] = in[idx]; 
    if ((threadIdx.x < 3) && ((idx+blockDim.x) < dsize)) sdata[threadIdx.x + blockDim.x] = in[idx + blockDim.x]; 
    __syncthreads(); 
    T temp = sdata[threadIdx.x]; 
    temp += sdata[threadIdx.x+1]; 
    temp += sdata[threadIdx.x+2]; 
    temp += sdata[threadIdx.x+3]; 
    if (idx < dsize - 4) out[idx] = temp; 
} 

int main(){ 

    mytype hdata1[] = {80,12,14,5,70,9,26,30,8,12,16,15}; 
    unsigned ds1 = sizeof(hdata1)/sizeof(hdata1[0]); 
    mytype hres1[ds1-4]; 
    thrust::device_vector<mytype> ddata1(hdata1, hdata1+ds1); 
    thrust::device_vector<mytype> dres1(ds1-4); 
    thrust::transform(thrust::make_zip_iterator(thrust::make_tuple(ddata1.begin(), ddata1.begin()+1, ddata1.begin()+2, ddata1.begin()+3)), thrust::make_zip_iterator(thrust::make_tuple(ddata1.end()-3, ddata1.end()-2, ddata1.end()-1, ddata1.end())), dres1.begin(), sum4()); 
    thrust::copy(dres1.begin(), dres1.end(), std::ostream_iterator<mytype>(std::cout, ",")); 
    std::cout << std::endl; 
    sum4kernel<<<(ds1+nTPB-1)/nTPB, nTPB>>>(thrust::raw_pointer_cast(ddata1.data()), thrust::raw_pointer_cast(dres1.data()), ds1); 
    cudaMemcpy(hres1, thrust::raw_pointer_cast(dres1.data()), (ds1-4)*sizeof(mytype), cudaMemcpyDeviceToHost); 
    for (int i = 0; i < ds1-4; i++) 
    std::cout << hres1[i] << ","; 
    std::cout << std::endl; 

    thrust::device_vector<mytype> ddata2(ds, 1); 
    thrust::device_vector<mytype> dres2(ds-4); 

    cudaEvent_t start, stop; 
    cudaEventCreate(&start); cudaEventCreate(&stop); 

    cudaEventRecord(start); 
    thrust::transform(thrust::make_zip_iterator(thrust::make_tuple(ddata2.begin(), ddata2.begin()+1, ddata2.begin()+2, ddata2.begin()+3)), thrust::make_zip_iterator(thrust::make_tuple(ddata2.end()-3, ddata2.end()-2, ddata2.end()-1, ddata2.end())), dres2.begin(), sum4()); 
    cudaEventRecord(stop); 
    thrust::host_vector<mytype> hres2 = dres2; 
    float et; 
    cudaEventElapsedTime(&et, start, stop); 
    std::cout << "thrust time: " << et << "ms" << std::endl; 
// validate 
    for (int i = 0; i < ds-4; i++) if (hres2[i] != 4) {std::cout << "thrust validation failure: " << i << "," << hres2[i] << std::endl; return 1;} 
    cudaEventRecord(start); 
    sum4kernel<<<(ds+nTPB-1)/nTPB, nTPB>>>(thrust::raw_pointer_cast(ddata2.data()), thrust::raw_pointer_cast(dres2.data()), ds); 
    cudaEventRecord(stop); 
    cudaMemcpy(&(hres2[0]), thrust::raw_pointer_cast(dres2.data()), (ds-4)*sizeof(mytype), cudaMemcpyDeviceToHost); 
    cudaEventElapsedTime(&et, start, stop); 
    std::cout << "cuda time: " << et << "ms" << std::endl; 
    for (int i = 0; i < ds-4; i++) if (hres2[i] != 4) {std::cout << "cuda validation failure: " << i << "," << hres2[i] << std::endl; return 1;} 
} 


$ nvcc -arch=sm_61 -o t88 t88.cu 
$ ./t88 
111,101,98,110,135,73,76,66, 
111,101,98,110,135,73,76,66, 
thrust time: 0.902464ms 
cuda time: 0.76288ms 
$ 

对于这个特定的GPU（泰坦X帕斯卡）没有推力时间之间太大的差别（约15％）用于32M元件数据集和CUDA时间。我们希望这个算法是内存限制的。

对于这个pascal titan x，bandwidthTest报告了关于可测量内存带宽的345 GB/s。

的CUDA实现必须加载整个数据集大小和存储整个数据集大小（约）=每元件2点的操作，因此所获得的带宽计算为这CUDA代码是：

(32*1048576 elements * 2 ops/element * 4 bytes/op)/0.00076288 s = ~350GB/s 

所以它看起来CUDA实现实现了大约最大可用带宽。