看起来你的问题有两个部分：
- 为什么我的第二个内核不工作？
- 如何使我的代码运行更快？
为什么我的第二个内核不工作？
你有几个问题：
- 在
i、j 的初始计算以及存储 C 值的索引中存在索引问题。
_syncthreads() 在条件块内的使用违规使用。
第一项是让代码正常工作的关键因素。
如何使我的代码运行更快？
这需要更多的工作。首先，你尝试“增加每个线程的工作量”并没有做到这一点，它只是增加了每个块的线程数（从128个到8 * 128个）。每个线程大约执行相同数量的工作。此外，在尝试使用2D线程块时，我认为发生了一些问题：
- 各种合并和共享内存冲突的读写模式被破坏了。
- 由于每个块所需的共享内存数量增加，有效占用率下降了。
第二个kernel的净效应是将执行时间大约延长了一倍。这不是我们想要的。
然而，增加每个线程的工作量可能是一个好主意，同时使用共享内存，并尝试保持良好的（全局、共享）内存访问模式，以及允许增加占用率。
以下是一项正在进行中的工作。下面的代码修复了您的第二个内核，以及时间基础设施，以及完整的数据验证，以及2个新内核。第一个新内核（＃3）是我所谓的“天真”内核。它只为每个输出点分配一个线程，并且每个线程循环遍历必要的向量，计算其个别结果。没有使用共享内存，甚至没有太多关注协同或任何其他优化。然而，通过对线程块配置进行微调（16,16） ->（8,32）线程，我从@talonmies的答案（现已删除）中观察到，这个内核比您的“快速”内核快得多（3倍）。经过进一步思考（8,32）的观察，我得出结论，下一次优化尝试应该集中在以下方面：
- 消除使用并行归约来计算向量距离的用法（即允许相邻的线程使用直接for循环来循环遍历向量）
- 最大化缓存的收益
- 有效使用共享内存
- 坚持完美的全局协同/完美地使用共享内存进行所有读写
项目4引发了评论中的问题：“我可以转置矩阵吗？”有了这个许可，就可以重新组织数据以便于上面的项目4。我的“快速”内核(#4)通过将B向量加载到共享内存中来解决上述第2项问题，同时允许缓存主要集中在缓存A向量，希望减少缓存抖动（A是两个向量数组中较小的一个，约为2MB - fermi L2为768K，Kepler L2为1.5MB）。通过以转置形式提供A，并有效地从共享内存“转置”B，可以使用直接的for循环计算向量距离，同时允许相邻线程具有完全协同的读写，以及对共享内存的“有效”使用（即非冲突负载和广播读取）。
对于我的特定时序（Quadro5000 cc2.0 GPU，CUDA 6，RHEL 5.5），我看到您的“快速”内核需要约2秒钟，我的“naive”内核需要约0.7秒钟，而我的“快速”内核需要约0.2秒钟，尽管使用了转置的（A，C）数据。
编辑：我进行了一项额外的优化，即每个块一次计算多个（CHKSIZE）B向量。您可以将CHKSIZE设置为1以查看先前的结果（~0.2秒）。我发现CHKSIZE为4时有很好的改进。这是试图利用A的数据重用的攻击。通过在CHKSIZE为4的情况下进行此额外的优化，内核4的内核时间降至约0.1秒。
以下是代码和示例运行：
$ cat t460.cu
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#define SIZE 128
#define N 4000
#define M 20000
#define CHKSIZE 4
__global__ void EuclideanDistances1( float *A, float *B , float *C , int n , int m)
{
__shared__ float accumResult[SIZE];
float sA;
float sB;
int bx = blockIdx.x;
int by = blockIdx.y;
int ty = threadIdx.y;
sA = A [bx * SIZE + ty];
sB = B [by * SIZE + ty];
__syncthreads();
accumResult[ty] = (sA - sB) * (sA - sB);
__syncthreads();
for (int stride = SIZE/2 ; stride > 0 ; stride >>= 1){
if (ty < stride)
{
accumResult[ty] += accumResult [stride + ty];
}
__syncthreads();
}
if ((ty == 0))
C [bx * m + by] = accumResult[ty];
__syncthreads();
}
__global__ void EuclideanDistances2( float *A, float *B , float *C, int n , int m)
{
__shared__ float accumResult[SIZE][8];
__shared__ float sA[SIZE][8];
__shared__ float sB[SIZE][8];
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int i = ((bx*8) + tx) * SIZE + ty;
int j = by * SIZE + ty;
sA[ty][tx] = A[i];
sB[ty][tx] = B[j];
__syncthreads();
accumResult[ty][tx] = (sA[ty][tx] - sB[ty][tx]) * (sA[ty][tx] - sB[ty][tx]);
__syncthreads();
for (int stride = SIZE/2 ; stride > 0 ; stride>>=1){
if (ty < stride)
{
accumResult[ty][tx] += accumResult [stride + ty][tx];
}
__syncthreads();
}
if (ty == 0)
C[((bx*8)+tx) * m + by] = accumResult[0][tx];
}
__global__ void EuclideanDistances3( float *A, float *B , float *C, int n , int m){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
int idy = threadIdx.y+blockDim.y*blockIdx.y;
float result = 0.0f;
if ((idx < n) && (idy < m)){
for (int i = 0; i < SIZE; i++){
float temp = A[(idx*SIZE)+i] - B[(idy*SIZE)+i];
result += temp * temp;}
C[(idx*m) + idy] = result;
}
}
__global__ void EuclideanDistances4( const float *A, const float *B , float *C, const int n , const int m){
__shared__ float my_sB[CHKSIZE*SIZE];
int bx = blockIdx.x;
while ((bx*CHKSIZE) < m){
int tx = threadIdx.x;
for (int i = 0; i < CHKSIZE; i++)
my_sB[(i*SIZE)+tx] = B[(((bx*CHKSIZE)+i)*SIZE)+tx];
__syncthreads();
while (tx < n){
float result[CHKSIZE];
for (int i = 0; i < CHKSIZE; i++)
result[i] = 0.0f;
for (int i = 0; i < SIZE; i++){
float Atemp = A[(n*i)+tx];
for (int j = 0; j < CHKSIZE; j++){
float temp = Atemp - my_sB[i + (j*SIZE)];
result[j] += temp * temp;}}
for (int i = 0; i < CHKSIZE; i++)
C[((i+(bx*CHKSIZE))*n)+ tx] = result[i];
tx += blockDim.x; }
__syncthreads();
bx += gridDim.x;}
}
float comp_euclid_sq(const float *rA, const float *rB, const int size){
float result = 0.0f;
float temp;
for (int i = 0; i < size; i++){
temp = (rA[i] - rB[i]);
result += temp * temp;}
return result;
}
int main()
{
float et1=0.0f, et2=0.0f, et3=0.0f, et4=0.0f;
cudaEvent_t start1, start2, start3,start4, stop1, stop2, stop3, stop4;
cudaEventCreate(&start1);
cudaEventCreate(&start2);
cudaEventCreate(&start3);
cudaEventCreate(&start4);
cudaEventCreate(&stop1);
cudaEventCreate(&stop2);
cudaEventCreate(&stop3);
cudaEventCreate(&stop4);
int n = N;
int m = M;
srand((unsigned)time(0));
float *matrixA = (float *) malloc (n * SIZE * sizeof(float));
for(int i=0; i < n * SIZE; i++)
matrixA[i] = (float) (rand()%100)+1;
float *matrixB = (float *) malloc (m * SIZE * sizeof(float));
for(int i=0; i < m * SIZE; i++)
matrixB[i] = (float) (rand()%100)+1;
float *results_kernel = (float *) malloc (n * m * sizeof(float));
float *cpu_results_kernel = (float *) malloc (n * m * sizeof(float));
for (int i = 0; i< n*m; i++)
cpu_results_kernel[i] = comp_euclid_sq(matrixA + ((i/m)*SIZE), matrixB + (i%m)*SIZE, SIZE);
float *d_matrixA;
float *d_matrixB;
cudaMalloc((void **)&d_matrixA, n * SIZE * sizeof(float));
cudaMalloc((void **)&d_matrixB, m * SIZE * sizeof(float));
cudaMemcpy(d_matrixA , matrixA , n * SIZE * sizeof(float) , cudaMemcpyHostToDevice);
cudaMemcpy(d_matrixB , matrixB , m * SIZE * sizeof(float) , cudaMemcpyHostToDevice);
float *d_results_kernel;
cudaMalloc((void **)&d_results_kernel , n * m * sizeof(float));
dim3 threads1 (1 , SIZE);
dim3 blocks1 (n , m);
cudaEventRecord(start1);
EuclideanDistances1 <<<blocks1 , threads1>>> (d_matrixA , d_matrixB , d_results_kernel , n , m);
cudaEventRecord(stop1);
cudaMemcpy(results_kernel , d_results_kernel , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
for (int i = 0; i< n*m; i++) {
if (results_kernel[i] != cpu_results_kernel[i]) {printf("cpu/kernel1 mismatch at %d, cpu: %f, kernel1: %f\n", i, cpu_results_kernel[i], results_kernel[i]); return 1;}}
cudaMemset(d_results_kernel, 0, n*m*sizeof(float));
cudaEventSynchronize(stop1);
cudaEventElapsedTime(&et1, start1, stop1);
dim3 threads2 (8 , SIZE);
dim3 blocks2 (n/8 , m);
cudaEventRecord(start2);
EuclideanDistances2 <<<blocks2 , threads2>>> (d_matrixA , d_matrixB , d_results_kernel , n , m);
cudaEventRecord(stop2);
cudaMemcpy(results_kernel , d_results_kernel , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
for (int i = 0; i< n*m; i++) {
if (results_kernel[i] != cpu_results_kernel[i]) {printf("cpu/kernel2 mismatch at %d, cpu: %f, kernel1: %f\n", i, cpu_results_kernel[i], results_kernel[i]); return 1;}}
cudaMemset(d_results_kernel, 0, n*m*sizeof(float));
cudaEventSynchronize(stop2);
cudaEventElapsedTime(&et2, start2, stop2);
cudaFuncSetCacheConfig(EuclideanDistances3, cudaFuncCachePreferL1);
dim3 threads3 (8, 32);
dim3 blocks3 (n/threads3.x , m/threads3.y);
cudaEventRecord(start3);
EuclideanDistances3 <<<blocks3 , threads3>>> (d_matrixA , d_matrixB , d_results_kernel , n , m);
cudaEventRecord(stop3);
cudaMemcpy(results_kernel , d_results_kernel , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
for (int i = 0; i< n*m; i++) {
if (results_kernel[i] != cpu_results_kernel[i]) {printf("cpu/kernel3 mismatch at %d, cpu: %f, kernel3: %f\n", i, cpu_results_kernel[i], results_kernel[i]); return 1;}}
cudaMemset(d_results_kernel, 0, n*m*sizeof(float));
cudaEventSynchronize(stop3);
cudaEventElapsedTime(&et3, start3, stop3);
float *matrixA_T = (float *) malloc (n * SIZE * sizeof(float));
for (int i = 0; i < n; i++)
for (int j = 0; j < SIZE; j++)
matrixA_T[(j*n)+i] = matrixA[(i*SIZE)+j];
cudaMemcpy(d_matrixA , matrixA_T , n * SIZE * sizeof(float) , cudaMemcpyHostToDevice);
cudaFuncSetCacheConfig(EuclideanDistances4, cudaFuncCachePreferL1);
dim3 threads4(SIZE);
dim3 blocks4(m/CHKSIZE);
cudaEventRecord(start4);
EuclideanDistances4 <<<blocks4 , threads4>>> (d_matrixA , d_matrixB , d_results_kernel , n , m);
cudaEventRecord(stop4);
cudaMemcpy(results_kernel , d_results_kernel , n * m *sizeof(float) , cudaMemcpyDeviceToHost);
for (int i = 0; i< n; i++)
for (int j = 0; j < m; j++)
if (results_kernel[(j*n)+i] != cpu_results_kernel[(i*m)+j]) {printf("cpu/kernel4 mismatch at %d,%d, cpu: %f, kernel4: %f\n", i,j, cpu_results_kernel[(i*m)+j], results_kernel[(j*n)+i]); return 1;}
cudaEventSynchronize(stop4);
cudaEventElapsedTime(&et4, start4, stop4);
cudaFree(d_results_kernel);
printf("Success!\n");
printf("kernel1 : %.fms, kernel2 : %.fms, kernel3 : %.fms, kernel4 : %.fms\n", et1, et2, et3, et4);
free(matrixA);
free(matrixB);
free(results_kernel);
return 0;
}
$ nvcc -O3 -arch=sm_20 -o t460 t460.cu
$ ./t460
Success!
kernel1 : 2213ms, kernel2 : 4660ms, kernel3 : 691ms, kernel4 : 99ms
$
希望这些内容能够给您更多的工作灵感。当然，您在使用cc3.0设备时可能会得到不同的时间。
是否有进一步的优化可能？可能是。我首先会瞄准如何利用向量A中的数据重用机会。 （向量B的数据重用已经在内核4中通过将其加载到共享内存中处理。可能有办法使用一些共享内存来存储A的部分，以使代码运行得更快。）
我想我也应该提到，遵循您提供的代码的方法，此代码正在计算
欧几里得距离的平方。对内核进行微小修改即可使其计算实际的欧几里得距离（
C [ ... ] = sqrtf(...);）。但是，我包含的验证假定结果在
float中完美存储为整数数量的“范围内”。您的测试案例符合此要求，但否则验证代码需要进行修改（如果使用
sqrtf）。
for循环中的if看起来是不一致的，因此它不应该包含__syncthreads();。 - Maciej Piechotka