对于任何 CUDA 程序员来说，最重要的 2 个优化优先级是有效地使用内存，并提供足够的并行性来隐藏延迟。我们将使用这些来指导我们的算法选择。
一个非常简单的线程策略（线程策略回答了问题，“每个线程将做什么或负责什么？”）在任何转换（相对于reduction) 类型问题是让每个线程负责 1 个输出值。您的问题符合 transformation 的描述 - 输出数据集大小与输入数据集大小相同。
我假设您打算拥有两个包含 3D 向量的长度相等的向量，并且您想要取每个中的第一个 3D 向量和每个中的第二个 3D 向量的叉积，依此类推。
如果我们选择每个线程1个输出点的线程策略（即result[x]或result[y]或result[z]，加起来就是3个输出点），那么我们将需要3个线程来计算每个线程的输出矢量叉积。如果我们有足够的向量来相乘，那么我们将有足够的线程来保持我们的机器“忙碌”并很好地隐藏延迟。根据经验，如果线程数为 10000 或更多，您的问题将开始在 GPU 上变得有趣，因此这意味着我们希望您的 1D 向量包含大约 3000 个或更多的 3D 向量。让我们假设是这种情况。
为了解决内存效率目标，我们的首要任务是从全局内存中加载您的矢量数据。理想情况下，我们希望它是 coalesced，这大致意味着相邻线程访问内存中的相邻元素。我们还希望合并输出存储，并且我们为每个线程选择一个输出点/一个向量分量的线程策略将很好地支持这一点。
为了有效地使用内存，我们希望从全局内存中只加载每个项目一次。您的算法自然会涉及少量数据重用。数据重用是显而易见的，因为result[y] 的计算依赖于vec2[z]，而result[x] 的计算也依赖于vec2[z]，仅举一个例子。因此，当存在数据重用时，典型的策略是先将数据加载到 CUDA 共享内存中，然后让线程根据共享内存中的数据执行计算。正如我们将看到的，这使得我们从全局内存中安排合并加载变得相当容易/方便，因为全局数据加载安排不再与线程或计算数据的使用紧密耦合。
最后一个挑战是找出一种索引模式，以便每个线程从共享内存中选择适当的元素来相乘。如果我们查看您在问题中描述的计算模式，我们会看到来自vec1 的第一个负载遵循计算结果的索引的 +1（模 3）偏移模式。所以x->y，y->z，和z -> x。同样，我们看到来自vec2 的下一次加载的+2（模3），来自vec1 的下一次加载的另一个+2（模3）模式和来自@ 的最终加载的另一个+1（模3）模式987654338@。
如果我们将所有这些想法结合起来，我们就可以编写一个应该具有普遍高效特征的内核：
$ 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 loop 编写内核。这对于本次讨论并不是非常重要，并且与此问题无关，因为我选择了一个等于问题大小 (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%。