使用 CUDA 提取矩阵列？
Extracting matrix columns with CUDA?
使用nvprof，我发现以下内核是我的CUDA应用程序的瓶颈
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
它打算将 indices 中列出的矩阵 src 的列提取到矩阵 tgt 中。请注意，矩阵 src 和 tgt 都有 numRows 行，并且以列主维度存储。此外，len = length(indices)*numRows 是矩阵 tgt.
的条目总数
我的问题：有没有更有效的方法来做到这一点？也感谢参考较旧的问题。我很惊讶我找不到之前问过的这个问题，因为它是 MATLAB 中使用的非常常见的操作tgt = src(:,indices(:));。
非常感谢！
作为一个复制内核，最佳性能大约受可用内存带宽的限制。可以通过 运行 宁 bandwidthTest cuda sample code 并参考设备到设备传输报告的数量（将因 GPU 而异）获得对此的粗略估计。
您的内核已经编写得相当好，加载和存储操作应该很好地合并。通过检查代码可以明显看出这一点，但是您可以通过 运行ning nvprof 和 --metrics gld_efficiency 以及 运行 和 --metrics gst_efficiency 来向自己证明这一点。这两个数字都应该接近 100%。 （根据我的测试。）
所以在我的例子中，当我 运行 你的内核在 Quadro5000 GPU 上时，将传输大小除以内核执行时间，我得到的数字大约是可用带宽的 60%。您的内核中没有太多其他事情发生，所以我们只关注循环中的这两行：
int colId = j / numRows;
int rowId = j % numRows;
事实证明，整数除法和模运算在 GPU 上都相当昂贵；它们是由编译器生成的指令序列创建的——没有本地除法或模机器指令。因此，如果我们能想出一种在主循环中摆脱这些的方法，我们就可以更接近 bandwidthTest（设备到设备）报告的 100% 带宽的目标。
由于你在循环中的增量是固定的（stride），我相信我们可以预先计算（大部分）需要添加到 colId 和 rowId 的增量对于循环的每次迭代，并在循环内使用加法和减法，而不是除法。修改后的内核如下所示：
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int div = stride/numRows;
int rem = stride%numRows;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
因此，通过预先计算循环每次迭代的增量，我们可以避免主循环中的 "expensive" 除法操作。那么性能呢？该内核更接近 100% 带宽目标。这是完整的测试代码：
#include <stdio.h>
#define NUMR 1000
#define NUMC 20000
#define DSIZE (NUMR*NUMC)
#define EXTC 10000
#define EXSZ (NUMR*EXTC)
#define nTPB 256
#define nBLK 64
typedef float real_t;
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int div = stride/numRows;
int rem = stride%numRows;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
__global__ void copy_kernel(real_t *tgt, real_t *src, int len){
int tid = threadIdx.x+blockDim.x*blockIdx.x;
while (tid < len){
tgt[tid] = src[tid];
tid+=blockDim.x*gridDim.x;
}
}
int main(){
real_t *h_a, *d_a, *h_b, *d_b, *h_bi;
h_a = (real_t *) malloc(DSIZE*sizeof(real_t));
cudaMalloc(&d_a, DSIZE*sizeof(real_t));
h_b = (real_t *) malloc(EXSZ*sizeof(real_t));
cudaMalloc(&d_b, EXSZ*sizeof(real_t));
h_bi = (real_t *) malloc(EXSZ*sizeof(real_t));
int *h_ind, *d_ind;
h_ind = (int *) malloc(EXTC*sizeof(int));
cudaMalloc(&d_ind, EXTC*sizeof(int));
for (int i = 0; i < EXTC; i++) h_ind[i] = i;
for (int i = 0; i < DSIZE; i++) h_a[i] = i;
cudaMemcpy(d_a, h_a, DSIZE*sizeof(real_t), cudaMemcpyHostToDevice);
cudaMemcpy(d_ind, h_ind, EXTC*sizeof(int), cudaMemcpyHostToDevice);
extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_b, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
copy_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR*EXTC);
cudaDeviceSynchronize();
my_extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_bi, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
for (int i = 0; i < EXSZ; i++)
if (h_bi[i] != h_b[i]) {printf("mismatch at %d, was: %f, should be: %f\n", i, h_bi[i], h_b[i]); return 1;}
printf("Success!\n");
return 0;
}
我已经包含了对你的内核的测试，我的内核，以及对 "copy kernel" 的中间测试，它只是对相同数量的数据进行纯复制。这有助于我们确认可用带宽上限的想法（见下文）。
现在是性能数据。 bandwidthTest 在这个 GPU 上告诉我们：
$ /usr/local/cuda/samples/bin/x86_64/linux/release/bandwidthTest
[CUDA Bandwidth Test] - Starting...
Running on...
Device 0: Quadro 5000
Quick Mode
Host to Device Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 5855.8
Device to Host Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 6334.8
Device to Device Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 101535.4
Result = PASS
NOTE: The CUDA Samples are not meant for performance measurements. Results may vary when GPU Boost is enabled.
$
所以我们有大约 100GB/s 的可用带宽。现在 运行ning nvprof --print-gpu-trace 我们看到：
$ nvprof --print-gpu-trace ./t822
==17985== NVPROF is profiling process 17985, command: ./t822
Success!
==17985== Profiling application: ./t822
==17985== Profiling result:
Start Duration Grid Size Block Size Regs* SSMem* DSMem* Size Throughput Device Context Stream Name
781.98ms 29.400ms - - - - - 80.000MB 2.7211GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
811.40ms 9.0560us - - - - - 40.000KB 4.4170GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
811.44ms 1.3377ms (64 1 1) (256 1 1) 15 0B 0B - - Quadro 5000 (0) 1 7 extractColumn_kernel(float*, float*, int*, int, int) [188]
812.78ms 21.614ms - - - - - 40.000MB 1.8507GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
834.94ms 816.10us (64 1 1) (256 1 1) 9 0B 0B - - Quadro 5000 (0) 1 7 copy_kernel(float*, float*, int) [194]
835.77ms 911.39us (64 1 1) (256 1 1) 18 0B 0B - - Quadro 5000 (0) 1 7 my_extractColumn_kernel(float*, float*, int*, int, int) [202]
836.69ms 20.661ms - - - - - 40.000MB 1.9360GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
$
此处的传输大小为 1000 行 * 10000 列 * 4 bytes/element * 每个元素 2 次传输（一次读取，一次写入）= 80,000,000 字节。您的原始内核在 1.34 毫秒内传输此数据，平均带宽约为 60GB/s。 "pure copy" 内核在 0.816 毫秒内传输相同的数据，平均带宽为 98GB/s——非常接近我们 100GB/s 的目标。我修改后的 "column-copy" 内核需要 0.911 毫秒，因此它提供了大约 88GB/秒的速度。
想要更高的性能？由于我的计算变量 rem 和 div 对于每个线程都是相同的，您可以在主机代码中预先计算这些量，并将它们作为内核参数传递。我不确定它会有多大的不同，但你可以试一试。您现在有了评估性能效果的路线图（如果有的话）。
备注：
请注意，我相信我在修改后的内核中更新索引的逻辑是合理的，但我还没有对它进行详尽的测试。它通过了我在这里介绍的简单测试用例。
我省略了proper cuda error checking。但是，如果您遇到问题，您应该使用它，and/or 运行 您的代码 cuda-memcheck.
正如我提到的，从合并的角度来看，你的内核已经写得很好，所以它已经达到了 "best case" 的大约 60%。因此，如果您在这里寻找 2 倍、5 倍或 10 倍的加速，您将找不到，并且期望它是不合理的。
编辑：进一步改进
在这种情况下我们没有更接近纯复制内核的一个可能原因是由于这种间接：
tgt[j] = src[indices[colId]*numRows + rowId];
^^^^^^^^^^^^^^
从 indices（一个全局变量）读取代表了一个 "extra" 数据访问，我们的纯拷贝内核不需要这样做。也可能有一些聪明的方法来优化这种访问的处理。由于它将是重复的（与 src 的读取和 tgt 的写入不同），这表明可能使用缓存或专用内存可能会有所帮助。
如果我们仔细检查访问的性质，我们可以观察到（对于相当大的矩阵），在大多数情况下，对 indices 的访问可能是 uniform 在 warp 中跨线程。这意味着通常，在给定的 warp 中，所有线程都将具有相同的 colId 值，因此将从 indices 请求相同的元素。这种类型的模式表明可以使用 CUDA __constant__ 内存进行优化。这里所需的更改并不广泛；我们基本上需要将 indices 数据移动到 __constant__ 数组。
这是修改后的代码：
$ cat t822.cu
#include <stdio.h>
#define NUMR 1000
#define NUMC 20000
#define DSIZE (NUMR*NUMC)
#define EXTC 10000
#define EXSZ (NUMR*EXTC)
#define nTPB 256
#define nBLK 64
typedef float real_t;
__constant__ int my_indices[EXTC];
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int numRows, int len, int div, int rem) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[my_indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
__global__ void copy_kernel(real_t *tgt, real_t *src, int len){
int tid = threadIdx.x+blockDim.x*blockIdx.x;
while (tid < len){
tgt[tid] = src[tid];
tid+=blockDim.x*gridDim.x;
}
}
int main(){
real_t *h_a, *d_a, *h_b, *d_b, *h_bi;
h_a = (real_t *) malloc(DSIZE*sizeof(real_t));
cudaMalloc(&d_a, DSIZE*sizeof(real_t));
h_b = (real_t *) malloc(EXSZ*sizeof(real_t));
cudaMalloc(&d_b, EXSZ*sizeof(real_t));
h_bi = (real_t *) malloc(EXSZ*sizeof(real_t));
int *h_ind, *d_ind;
h_ind = (int *) malloc(EXTC*sizeof(int));
cudaMalloc(&d_ind, EXTC*sizeof(int));
for (int i = 0; i < EXTC; i++) h_ind[i] = i;
for (int i = 0; i < DSIZE; i++) h_a[i] = i;
cudaMemcpy(d_a, h_a, DSIZE*sizeof(real_t), cudaMemcpyHostToDevice);
cudaMemcpy(d_ind, h_ind, EXTC*sizeof(int), cudaMemcpyHostToDevice);
extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_b, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
copy_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR*EXTC);
cudaDeviceSynchronize();
cudaMemcpyToSymbol(my_indices, h_ind, EXTC*sizeof(int));
int mydiv = (nBLK*nTPB)/NUMR;
int myrem = (nBLK*nTPB)%NUMR;
my_extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR, NUMR*EXTC, mydiv, myrem);
cudaMemcpy(h_bi, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
for (int i = 0; i < EXSZ; i++)
if (h_bi[i] != h_b[i]) {printf("mismatch at %d, was: %f, should be: %f\n", i, h_bi[i], h_b[i]); return 1;}
printf("Success!\n");
return 0;
}
$
从性能结果来看：
$ nvprof --print-gpu-trace ./t822
==18998== NVPROF is profiling process 18998, command: ./t822
Success!
==18998== Profiling application: ./t822
==18998== Profiling result:
Start Duration Grid Size Block Size Regs* SSMem* DSMem* Size Throughput Device Context Stream Name
773.01ms 28.300ms - - - - - 80.000MB 2.8269GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
801.33ms 9.0240us - - - - - 40.000KB 4.4326GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
801.38ms 1.3001ms (64 1 1) (256 1 1) 15 0B 0B - - Quadro 5000 (0) 1 7 extractColumn_kernel(float*, float*, int*, int, int) [188]
802.68ms 20.773ms - - - - - 40.000MB 1.9256GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
823.98ms 811.75us (64 1 1) (256 1 1) 9 0B 0B - - Quadro 5000 (0) 1 7 copy_kernel(float*, float*, int) [194]
824.82ms 8.9920us - - - - - 40.000KB 4.4484GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
824.83ms 824.65us (64 1 1) (256 1 1) 13 0B 0B - - Quadro 5000 (0) 1 7 my_extractColumn_kernel(float*, float*, int, int, int, int) [204]
825.66ms 21.023ms - - - - - 40.000MB 1.9027GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
我们看到我们现在非常接近 100% 使用可用带宽的目标（824 us 对比复制内核时间 811 us）。 __constant__ 内存总量限制为 64KB，因此这意味着只有当索引（实际上需要复制的列数）大约小于 16,000 时，它才能以这种方式使用。
使用nvprof，我发现以下内核是我的CUDA应用程序的瓶颈
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
它打算将 indices 中列出的矩阵 src 的列提取到矩阵 tgt 中。请注意，矩阵 src 和 tgt 都有 numRows 行，并且以列主维度存储。此外，len = length(indices)*numRows 是矩阵 tgt.
我的问题：有没有更有效的方法来做到这一点？也感谢参考较旧的问题。我很惊讶我找不到之前问过的这个问题，因为它是 MATLAB 中使用的非常常见的操作tgt = src(:,indices(:));。
非常感谢！
作为一个复制内核，最佳性能大约受可用内存带宽的限制。可以通过 运行 宁 bandwidthTest cuda sample code 并参考设备到设备传输报告的数量（将因 GPU 而异）获得对此的粗略估计。
您的内核已经编写得相当好，加载和存储操作应该很好地合并。通过检查代码可以明显看出这一点，但是您可以通过 运行ning nvprof 和 --metrics gld_efficiency 以及 运行 和 --metrics gst_efficiency 来向自己证明这一点。这两个数字都应该接近 100%。 （根据我的测试。）
所以在我的例子中，当我 运行 你的内核在 Quadro5000 GPU 上时，将传输大小除以内核执行时间，我得到的数字大约是可用带宽的 60%。您的内核中没有太多其他事情发生，所以我们只关注循环中的这两行：
int colId = j / numRows;
int rowId = j % numRows;
事实证明，整数除法和模运算在 GPU 上都相当昂贵；它们是由编译器生成的指令序列创建的——没有本地除法或模机器指令。因此，如果我们能想出一种在主循环中摆脱这些的方法，我们就可以更接近 bandwidthTest（设备到设备）报告的 100% 带宽的目标。
由于你在循环中的增量是固定的（stride），我相信我们可以预先计算（大部分）需要添加到 colId 和 rowId 的增量对于循环的每次迭代，并在循环内使用加法和减法，而不是除法。修改后的内核如下所示：
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int div = stride/numRows;
int rem = stride%numRows;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
因此，通过预先计算循环每次迭代的增量，我们可以避免主循环中的 "expensive" 除法操作。那么性能呢？该内核更接近 100% 带宽目标。这是完整的测试代码：
#include <stdio.h>
#define NUMR 1000
#define NUMC 20000
#define DSIZE (NUMR*NUMC)
#define EXTC 10000
#define EXSZ (NUMR*EXTC)
#define nTPB 256
#define nBLK 64
typedef float real_t;
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int div = stride/numRows;
int rem = stride%numRows;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
__global__ void copy_kernel(real_t *tgt, real_t *src, int len){
int tid = threadIdx.x+blockDim.x*blockIdx.x;
while (tid < len){
tgt[tid] = src[tid];
tid+=blockDim.x*gridDim.x;
}
}
int main(){
real_t *h_a, *d_a, *h_b, *d_b, *h_bi;
h_a = (real_t *) malloc(DSIZE*sizeof(real_t));
cudaMalloc(&d_a, DSIZE*sizeof(real_t));
h_b = (real_t *) malloc(EXSZ*sizeof(real_t));
cudaMalloc(&d_b, EXSZ*sizeof(real_t));
h_bi = (real_t *) malloc(EXSZ*sizeof(real_t));
int *h_ind, *d_ind;
h_ind = (int *) malloc(EXTC*sizeof(int));
cudaMalloc(&d_ind, EXTC*sizeof(int));
for (int i = 0; i < EXTC; i++) h_ind[i] = i;
for (int i = 0; i < DSIZE; i++) h_a[i] = i;
cudaMemcpy(d_a, h_a, DSIZE*sizeof(real_t), cudaMemcpyHostToDevice);
cudaMemcpy(d_ind, h_ind, EXTC*sizeof(int), cudaMemcpyHostToDevice);
extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_b, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
copy_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR*EXTC);
cudaDeviceSynchronize();
my_extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_bi, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
for (int i = 0; i < EXSZ; i++)
if (h_bi[i] != h_b[i]) {printf("mismatch at %d, was: %f, should be: %f\n", i, h_bi[i], h_b[i]); return 1;}
printf("Success!\n");
return 0;
}
我已经包含了对你的内核的测试，我的内核，以及对 "copy kernel" 的中间测试，它只是对相同数量的数据进行纯复制。这有助于我们确认可用带宽上限的想法（见下文）。
现在是性能数据。 bandwidthTest 在这个 GPU 上告诉我们：
$ /usr/local/cuda/samples/bin/x86_64/linux/release/bandwidthTest
[CUDA Bandwidth Test] - Starting...
Running on...
Device 0: Quadro 5000
Quick Mode
Host to Device Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 5855.8
Device to Host Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 6334.8
Device to Device Bandwidth, 1 Device(s)
PINNED Memory Transfers
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 101535.4
Result = PASS
NOTE: The CUDA Samples are not meant for performance measurements. Results may vary when GPU Boost is enabled.
$
所以我们有大约 100GB/s 的可用带宽。现在 运行ning nvprof --print-gpu-trace 我们看到：
$ nvprof --print-gpu-trace ./t822
==17985== NVPROF is profiling process 17985, command: ./t822
Success!
==17985== Profiling application: ./t822
==17985== Profiling result:
Start Duration Grid Size Block Size Regs* SSMem* DSMem* Size Throughput Device Context Stream Name
781.98ms 29.400ms - - - - - 80.000MB 2.7211GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
811.40ms 9.0560us - - - - - 40.000KB 4.4170GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
811.44ms 1.3377ms (64 1 1) (256 1 1) 15 0B 0B - - Quadro 5000 (0) 1 7 extractColumn_kernel(float*, float*, int*, int, int) [188]
812.78ms 21.614ms - - - - - 40.000MB 1.8507GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
834.94ms 816.10us (64 1 1) (256 1 1) 9 0B 0B - - Quadro 5000 (0) 1 7 copy_kernel(float*, float*, int) [194]
835.77ms 911.39us (64 1 1) (256 1 1) 18 0B 0B - - Quadro 5000 (0) 1 7 my_extractColumn_kernel(float*, float*, int*, int, int) [202]
836.69ms 20.661ms - - - - - 40.000MB 1.9360GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
$
此处的传输大小为 1000 行 * 10000 列 * 4 bytes/element * 每个元素 2 次传输（一次读取，一次写入）= 80,000,000 字节。您的原始内核在 1.34 毫秒内传输此数据，平均带宽约为 60GB/s。 "pure copy" 内核在 0.816 毫秒内传输相同的数据，平均带宽为 98GB/s——非常接近我们 100GB/s 的目标。我修改后的 "column-copy" 内核需要 0.911 毫秒，因此它提供了大约 88GB/秒的速度。
想要更高的性能？由于我的计算变量 rem 和 div 对于每个线程都是相同的，您可以在主机代码中预先计算这些量，并将它们作为内核参数传递。我不确定它会有多大的不同，但你可以试一试。您现在有了评估性能效果的路线图（如果有的话）。
备注：
请注意，我相信我在修改后的内核中更新索引的逻辑是合理的，但我还没有对它进行详尽的测试。它通过了我在这里介绍的简单测试用例。
我省略了proper cuda error checking。但是，如果您遇到问题，您应该使用它，and/or 运行 您的代码
cuda-memcheck.正如我提到的，从合并的角度来看，你的内核已经写得很好，所以它已经达到了 "best case" 的大约 60%。因此，如果您在这里寻找 2 倍、5 倍或 10 倍的加速，您将找不到，并且期望它是不合理的。
编辑：进一步改进
在这种情况下我们没有更接近纯复制内核的一个可能原因是由于这种间接：
tgt[j] = src[indices[colId]*numRows + rowId];
^^^^^^^^^^^^^^
从 indices（一个全局变量）读取代表了一个 "extra" 数据访问，我们的纯拷贝内核不需要这样做。也可能有一些聪明的方法来优化这种访问的处理。由于它将是重复的（与 src 的读取和 tgt 的写入不同），这表明可能使用缓存或专用内存可能会有所帮助。
如果我们仔细检查访问的性质，我们可以观察到（对于相当大的矩阵），在大多数情况下，对 indices 的访问可能是 uniform 在 warp 中跨线程。这意味着通常，在给定的 warp 中，所有线程都将具有相同的 colId 值，因此将从 indices 请求相同的元素。这种类型的模式表明可以使用 CUDA __constant__ 内存进行优化。这里所需的更改并不广泛；我们基本上需要将 indices 数据移动到 __constant__ 数组。
这是修改后的代码：
$ cat t822.cu
#include <stdio.h>
#define NUMR 1000
#define NUMC 20000
#define DSIZE (NUMR*NUMC)
#define EXTC 10000
#define EXSZ (NUMR*EXTC)
#define nTPB 256
#define nBLK 64
typedef float real_t;
__constant__ int my_indices[EXTC];
__global__ void extractColumn_kernel(real_t *tgt, real_t *src, int *indices, int numRows, int len) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
for (int j = tid; j < len; j += stride) {
int colId = j / numRows;
int rowId = j % numRows;
tgt[j] = src[indices[colId]*numRows + rowId];
}
}
__global__ void my_extractColumn_kernel(real_t *tgt, real_t *src, int numRows, int len, int div, int rem) {
int stride = gridDim.x * blockDim.x;
int tid = blockDim.x * blockIdx.x + threadIdx.x;
int colId = tid / numRows;
int rowId = tid % numRows;
for (int j = tid; j < len; j += stride) {
tgt[j] = src[my_indices[colId]*numRows + rowId];
colId += div;
rowId += rem;
if (rowId >= numRows) {rowId-=numRows; colId++;}
}
}
__global__ void copy_kernel(real_t *tgt, real_t *src, int len){
int tid = threadIdx.x+blockDim.x*blockIdx.x;
while (tid < len){
tgt[tid] = src[tid];
tid+=blockDim.x*gridDim.x;
}
}
int main(){
real_t *h_a, *d_a, *h_b, *d_b, *h_bi;
h_a = (real_t *) malloc(DSIZE*sizeof(real_t));
cudaMalloc(&d_a, DSIZE*sizeof(real_t));
h_b = (real_t *) malloc(EXSZ*sizeof(real_t));
cudaMalloc(&d_b, EXSZ*sizeof(real_t));
h_bi = (real_t *) malloc(EXSZ*sizeof(real_t));
int *h_ind, *d_ind;
h_ind = (int *) malloc(EXTC*sizeof(int));
cudaMalloc(&d_ind, EXTC*sizeof(int));
for (int i = 0; i < EXTC; i++) h_ind[i] = i;
for (int i = 0; i < DSIZE; i++) h_a[i] = i;
cudaMemcpy(d_a, h_a, DSIZE*sizeof(real_t), cudaMemcpyHostToDevice);
cudaMemcpy(d_ind, h_ind, EXTC*sizeof(int), cudaMemcpyHostToDevice);
extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, d_ind, NUMR, NUMR*EXTC);
cudaMemcpy(h_b, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
copy_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR*EXTC);
cudaDeviceSynchronize();
cudaMemcpyToSymbol(my_indices, h_ind, EXTC*sizeof(int));
int mydiv = (nBLK*nTPB)/NUMR;
int myrem = (nBLK*nTPB)%NUMR;
my_extractColumn_kernel<<<nBLK, nTPB>>>(d_b, d_a, NUMR, NUMR*EXTC, mydiv, myrem);
cudaMemcpy(h_bi, d_b, EXSZ*sizeof(real_t), cudaMemcpyDeviceToHost);
for (int i = 0; i < EXSZ; i++)
if (h_bi[i] != h_b[i]) {printf("mismatch at %d, was: %f, should be: %f\n", i, h_bi[i], h_b[i]); return 1;}
printf("Success!\n");
return 0;
}
$
从性能结果来看：
$ nvprof --print-gpu-trace ./t822
==18998== NVPROF is profiling process 18998, command: ./t822
Success!
==18998== Profiling application: ./t822
==18998== Profiling result:
Start Duration Grid Size Block Size Regs* SSMem* DSMem* Size Throughput Device Context Stream Name
773.01ms 28.300ms - - - - - 80.000MB 2.8269GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
801.33ms 9.0240us - - - - - 40.000KB 4.4326GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
801.38ms 1.3001ms (64 1 1) (256 1 1) 15 0B 0B - - Quadro 5000 (0) 1 7 extractColumn_kernel(float*, float*, int*, int, int) [188]
802.68ms 20.773ms - - - - - 40.000MB 1.9256GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
823.98ms 811.75us (64 1 1) (256 1 1) 9 0B 0B - - Quadro 5000 (0) 1 7 copy_kernel(float*, float*, int) [194]
824.82ms 8.9920us - - - - - 40.000KB 4.4484GB/s Quadro 5000 (0) 1 7 [CUDA memcpy HtoD]
824.83ms 824.65us (64 1 1) (256 1 1) 13 0B 0B - - Quadro 5000 (0) 1 7 my_extractColumn_kernel(float*, float*, int, int, int, int) [204]
825.66ms 21.023ms - - - - - 40.000MB 1.9027GB/s Quadro 5000 (0) 1 7 [CUDA memcpy DtoH]
Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
我们看到我们现在非常接近 100% 使用可用带宽的目标（824 us 对比复制内核时间 811 us）。 __constant__ 内存总量限制为 64KB，因此这意味着只有当索引（实际上需要复制的列数）大约小于 16,000 时，它才能以这种方式使用。