CUDA求任意长度向量和debug实录

症状：点击运行，程序一直跑，却没有输出

问题代码：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
 
#include <stdio.h>
 
#define N 256 //向量长度，自行设置
#define BLOCK 128//线程块数目，可根据硬件限制随意设置
#define BLOCKDIM 128 //线程块内线程数目，可根据硬件限制随意设置
 
__global__ void addKernel(int *c, const int *a, const int *b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    while (i<N)
    {
        c[i] = a[i] + b[i];
        i += blockDim.x * blockIdx.x;
    }
}
 
int main()
{
    int* a = new int[N];
    int* b = new int[N];
    int* c = new int[N];
    int *dev_a = 0;
    int *dev_b = 0;
    int *dev_c = 0;
 
    for (int i = 1; i <= N; ++i) {
        a[i-1] = i;
        b[i-1] = i & 1;
    }
    cudaMalloc((void**)&dev_c, N * sizeof(int));
    cudaMalloc((void**)&dev_a, N * sizeof(int));
    cudaMalloc((void**)&dev_b, N * sizeof(int));
    cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice);
 
    addKernel<<<BLOCK, BLOCKDIM>>>(dev_c, dev_a, dev_b);
    cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost);
    
    for (int i = 0; i < N; ++i) {
        printf("第%d维:%d+%d=%d\n", i+1, a[i], b[i], c[i]);
    }
 
    delete[] a;
    delete[] b;
    delete[] c;
    cudaFree(dev_c);
    cudaFree(dev_a);
    cudaFree(dev_b);
 
    return 0;
}

第一步：在main函数中各个可疑程序段后添加printf语句，运行查看是哪段代码的问题，修改后的main函数

int main()
{
    int* a = new int[N];
    int* b = new int[N];
    int* c = new int[N];
    int* dev_a = 0;
    int* dev_b = 0;
    int* dev_c = 0;

    for (int i = 1; i <= N; ++i) {
        a[i - 1] = i;
        b[i - 1] = i & 1;
    }
    cudaMalloc((void**)&dev_c, N * sizeof(int));
    cudaMalloc((void**)&dev_a, N * sizeof(int));
    cudaMalloc((void**)&dev_b, N * sizeof(int));
    cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice);
    printf("111111111111111\n");

    addKernel << <BLOCK, BLOCKDIM >> > (dev_c, dev_a, dev_b);
    printf("222222222222222\n");
    cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost);
    printf("333333333333333\n");

    for (int i = 0; i < N; ++i) {
        printf("第%d维:%d+%d=%d\n", i + 1, a[i], b[i], c[i]);
    }

    delete[] a;
    delete[] b;
    delete[] c;
    cudaFree(dev_c);
    cudaFree(dev_a);
    cudaFree(dev_b);

    return 0;
}

运行修改main函数后的代码，程序输出

111111111111111
222222222222222

说明问题出在语句cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost);   上网查询cudaMemcpyDeviceToHost模式运行失败的原因后，我尝试了重启电脑，没有用，说明是程序代码有问题。这个模式运行失败可能的原因还有并行处理数据的时候有部分数据超过了内存的边界，会导致处理结果无法拷贝出来。于是我仔细地检查了add函数

__global__ void addKernel(int* c, const int* a, const int* b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    while (i < N)
    {
        c[i] = a[i] + b[i];
        i += blockDim.x * blockIdx.x;
    }
}

add函数中我已经对边界做了判断（i<N），因而不可能是越界的问题。但是既然无法把并行处理结果拷贝出来，肯定是并行处理的过程有问题。

第二步：在add中增加printf语句，运行程序查看各个线程块的运行情况，修改后的add函数如下

__global__ void addKernel(int* c, const int* a, const int* b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    while (i < N)
    {
        c[i] = a[i] + b[i];
        i += blockDim.x * blockIdx.x;
        printf("%d号线程块", blockIdx.x);
    }
}

运行修改add函数后的代码，程序疯狂输出“0号线程块”，而且没有要停止的迹象，因此可以判断，0号线程块出现了死循环。为什么会这样呢？仔细看语句i += blockDim.x * blockIdx.x;   0号block的blockIdx.x是0，所以在0号block中的每个thread的循环变量i每次增量都为0，永远不会破坏i<N这个条件，while循环就一直进行。

这个程序的意思是重复利用固定数量的线程运算矢量的若干维度求和。实际上，GPU一次同时运行的线程数目，是一个grid里面横着的一排，所以一次同时计算的矢量维度序号范围是grid里一横排包括的序号范围，所以对于每个重复利用的线程，下次计算的维度序号和这次相差一横排的宽度。
thread、block、grid之间到底是什么关系？详情可以看这篇博客的图画。总结来说，就是一个grid里包含若干block，一个block里包含若干thread，blockDim是一个线程格在横排上的宽度，即横排上分布多少个thread，gridDim是一个grid在横排上的宽度，即横排上分布多少个block。

经过以上分析，add代码应修正为

__global__ void addKernel(int* c, const int* a, const int* b)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    while (i < N)
    {
        c[i] = a[i] + b[i];
        i += blockDim.x * gridDim.x;
    }
}

运行修正后代码，程序正常运行。