CUDA平行十字产品

Question

免责声明:我对CUDA和并行编程都很陌生 - 所以如果你不想回答我的问题,请忽略这一点,或者至少指出我正确的资源,这样我就可以自己找到答案了.

这是我想要使用并行编程解决的特殊问题.我具有存储3D矢量以这种格式一些一维数组- > [v0x, v0y, v0z, ... vnx, vny, vnz]其中n是所述载体,和x,y,z是各自的组分.

假设我想找到[v0, v1, ... vn]一个数组中向量[v0, v1, ... vn]与另一个数组中相应向量之间的叉积.

没有并行化,计算非常简单:

result[x] = vec1[y]*vec2[z] - vec1[z]*vec2[y];

result[y] = vec1[z]*vec2[x] - vec1[x]*vec2[z];

result[z] = vec1[x]*vec2[y] - vec1[y]*vec2[x];

我遇到的问题是了解如何为我目前拥有的数组实现CUDA并行化.由于结果向量中的每个值都是单独的计算,因此我可以有效地并行运行每个向量的上述计算.由于得到的叉积的每个分量都是单独的计算,因此它们也可以并行运行.我将如何设置块和线程/考虑为这样的问题设置线程？

Answer 1

任何CUDA程序员的前2个优化优先级都是有效地使用内存,并暴露足够的并行性来隐藏延迟.我们将使用它们来指导我们的算法选择.

一个非常简单的线程策略(线程策略回答问题,"每个线程将做什么或负责什么？")在任何转换(而不是简化)类型问题中是让每个线程负责1个输出值.您的问题符合转换描述- 输出数据集大小是输入数据集大小的顺序.

我假设您打算有两个相等长度的矢量包含您的3D矢量,并且您想要获取每个中的第一个3D矢量和每个中的第二个3D矢量的叉积,依此类推.

如果我们选择每个线程1个输出点的线程策略(即result[x]或者,result[y]或者result[z]一起将是3个输出点),那么我们将需要3个线程来计算每个向量交叉积的输出.如果我们有足够的向量来增加,那么我们将有足够的线程来保持我们的机器"忙"并且很好地隐藏延迟.根据经验,如果线程数为10000或更多,您的问题将开始在GPU上变得有趣,所以这意味着我们希望您的1D向量包含大约3000个3D向量或更多.我们假设是这样的.

为了解决内存效率目标,我们的首要任务是从全局内存中加载矢量数据.理想情况下,我们希望将其合并,这大致意味着相邻线程访问内存中的相邻元素.我们也希望输出存储也能够合并,并且我们为每个线程选择一个输出点/一个向量组件的线程策略将很好地支持它.

为了有效使用内存,我们希望理想情况下只从全局内存加载一次项目.您的算法自然会涉及少量数据重用.数据重用是显而易见的,因为计算result[y]取决于vec2[z]并且计算result[x]也取决于vec2[z]仅选择一个示例.因此,存在数据重用时的典型策略是首先将数据加载到CUDA共享内存中,然后允许线程基于共享内存中的数据执行其计算.正如我们将要看到的,这使得我们可以相当容易/方便地安排来自全局存储器的合并负载,因为全局数据加载安排不再紧密地耦合到线程或用于计算的数据的使用.

最后一个挑战是找出一个索引模式,以便每个线程从共享内存中选择适当的元素来相乘.如果我们查看您在问题中描述的计算模式,我们会看到第一个加载来自vec1计算结果的索引的偏移模式+1(模3).所以x- > y,y- > z和z- > x.同样地,我们看到下一个负载的+2(模3),下一个负载的vec2另一个+2(模3)模式vec1和最终负载的另一个+1(模3)模式vec2.

如果我们将所有这些想法结合起来,我们就可以编写一个应具有通常有效特性的内核:

$ cat t1003.cu
#include <stdio.h>

#define TV1 1
#define TV2 2
const size_t N = 4096;    // number of 3D vectors
const int blksize = 192;  // choose as multiple of 3 and 32, and less than 1024
typedef float mytype;
//pairwise vector cross product
template <typename T>
__global__ void vcp(const T * __restrict__ vec1, const T * __restrict__ vec2, T * __restrict__ res, const size_t n){

  __shared__ T sv1[blksize];
  __shared__ T sv2[blksize];
  size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
  while (idx < 3*n){ // grid-stride loop
    // load shared memory using coalesced pattern to global memory
    sv1[threadIdx.x] = vec1[idx];
    sv2[threadIdx.x] = vec2[idx];
    // compute modulo/offset indexing for thread loads of shared data from vec1, vec2
    int my_mod = threadIdx.x%3;   // costly, but possibly hidden by global load latency
    int off1 = my_mod+1;
    if (off1 > 2) off1 -= 3;
    int off2 = my_mod+2;
    if (off2 > 2) off2 -= 3;
    __syncthreads();
    // each thread loads its computation elements from shared memory
    T t1 = sv1[threadIdx.x-my_mod+off1];
    T t2 = sv2[threadIdx.x-my_mod+off2];
    T t3 = sv1[threadIdx.x-my_mod+off2];
    T t4 = sv2[threadIdx.x-my_mod+off1];
    // compute result, and store using coalesced pattern, to global memory
    res[idx] = t1*t2-t3*t4;
    idx += gridDim.x*blockDim.x;}  // for grid-stride loop
}

int main(){

  mytype *h_v1, *h_v2, *d_v1, *d_v2, *h_res, *d_res;
  h_v1  = (mytype *)malloc(N*3*sizeof(mytype));
  h_v2  = (mytype *)malloc(N*3*sizeof(mytype));
  h_res = (mytype *)malloc(N*3*sizeof(mytype));
  cudaMalloc(&d_v1,  N*3*sizeof(mytype));
  cudaMalloc(&d_v2,  N*3*sizeof(mytype));
  cudaMalloc(&d_res, N*3*sizeof(mytype));
  for (int i = 0; i<N; i++){
    h_v1[3*i]    = TV1;
    h_v1[3*i+1]  = 0;
    h_v1[3*i+2]  = 0;
    h_v2[3*i]    = 0;
    h_v2[3*i+1]  = TV2;
    h_v2[3*i+2]  = 0;
    h_res[3*i]   = 0;
    h_res[3*i+1] = 0;
    h_res[3*i+2] = 0;}
  cudaMemcpy(d_v1, h_v1, N*3*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_v2, h_v2, N*3*sizeof(mytype), cudaMemcpyHostToDevice);
  vcp<<<(N*3+blksize-1)/blksize, blksize>>>(d_v1, d_v2, d_res, N);
  cudaMemcpy(h_res, d_res, N*3*sizeof(mytype), cudaMemcpyDeviceToHost);
  // verification
  for (int i = 0; i < N; i++) if ((h_res[3*i] != 0) || (h_res[3*i+1] != 0) || (h_res[3*i+2] != TV1*TV2)) { printf("mismatch at %d, was: %f, %f, %f, should be: %f, %f, %f\n", i, h_res[3*i], h_res[3*i+1], h_res[3*i+2], (float)0, (float)0, (float)(TV1*TV2)); return -1;}
  printf("%s\n", cudaGetErrorString(cudaGetLastError()));
  return 0;
}


$ nvcc t1003.cu -o t1003
$ cuda-memcheck ./t1003
========= CUDA-MEMCHECK
no error
========= ERROR SUMMARY: 0 errors
$

请注意,我选择使用grid-stride循环编写内核.这对于这个讨论并不是非常重要,而且与此问题无关,因为我选择的网格大小等于问题大小(4096*3).但是,对于更大的问题大小,您可以选择比整体问题大小更小的网格大小,以获得一些可能的小效率增益.

对于这样一个简单的问题,定义"最优性"相当容易.然而,最佳方案是加载输入数据(仅一次)并写入输出数据所需的时间.如果我们考虑上面的测试代码的更大版本,更改N为40960(并且不进行其他更改),则读取和写入的总数据将为40960*3*4*3字节.如果我们分析该代码,然后将其bandwidthTest作为峰值可实现内存带宽的代理进行比较,我们观察到:

$ CUDA_VISIBLE_DEVICES="1" nvprof ./t1003
==27861== NVPROF is profiling process 27861, command: ./t1003
no error
==27861== Profiling application: ./t1003
==27861== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:   65.97%  162.22us         2  81.109us  77.733us  84.485us  [CUDA memcpy HtoD]
                   30.04%  73.860us         1  73.860us  73.860us  73.860us  [CUDA memcpy DtoH]
                    4.00%  9.8240us         1  9.8240us  9.8240us  9.8240us  void vcp<float>(float const *, float const *, float*, unsigned long)
      API calls:   99.10%  249.79ms         3  83.263ms  6.8890us  249.52ms  cudaMalloc
                    0.46%  1.1518ms        96  11.998us     374ns  454.09us  cuDeviceGetAttribute
                    0.25%  640.18us         3  213.39us  186.99us  229.86us  cudaMemcpy
                    0.10%  255.00us         1  255.00us  255.00us  255.00us  cuDeviceTotalMem
                    0.05%  133.16us         1  133.16us  133.16us  133.16us  cuDeviceGetName
                    0.03%  71.903us         1  71.903us  71.903us  71.903us  cudaLaunchKernel
                    0.01%  15.156us         1  15.156us  15.156us  15.156us  cuDeviceGetPCIBusId
                    0.00%  7.0920us         3  2.3640us     711ns  4.6520us  cuDeviceGetCount
                    0.00%  2.7780us         2  1.3890us     612ns  2.1660us  cuDeviceGet
                    0.00%  1.9670us         1  1.9670us  1.9670us  1.9670us  cudaGetLastError
                    0.00%     361ns         1     361ns     361ns     361ns  cudaGetErrorString
$ CUDA_VISIBLE_DEVICES="1" /usr/local/cuda/samples/bin/x86_64/linux/release/bandwidthTest
[CUDA Bandwidth Test] - Starting...
Running on...

 Device 0: Tesla K20Xm
 Quick Mode

 Host to Device Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     6375.8

 Device to Host Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     6554.3

 Device to Device Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     171220.3

Result = PASS

NOTE: The CUDA Samples are not meant for performance measurements. Results may vary when GPU Boost is enabled.
$

内核需要9.8240us来执行,并且在那时加载或存储总共40960*3*4*3字节的数据.因此,内核实现的内存带宽为40960*3*4*3/0.000009824或150 GB/s.此GPU上可实现峰值的代理测量值为171 GB/s,因此该内核可实现88%的最佳吞吐量.通过更仔细的基准测试连续两次运行内核,第二次执行只需要执行8.99us.这使得在这种情况下实现的带宽达到峰值可实现吞吐量的96%.