如何写出正确CUDA程序(持续更新中)

本文不谈CUDA优化，只谈如何正确写出CUDA程序，先写对，再说优化。最近写CUDA程序也是断断续续的，经常程序写完了需要很长的时间调试，各种错误，调试的方法也是简单的printf()，目前没找到更好的调试方法，如果有更好的办法的化希望告知。
昨天总结了下为什么写的CUDA程序总是出错（就先别说优化呢？写对再再谈别的吧），总结出一些心得和经验(算是经验)吧，结合例子看吧。

首先谈谈调试方法吧，目前我遇到的主要是两种错误：内存错误和逻辑错误，其他错误。

其他错误

对于其他错误，一般加上错误处理即可，错误处理见CUDA编程的错误处理
这样便于快速定位，不加的话出错误了定位更难受。

内存错误

目前对于内存错误，编译的时候加入下面选项

-ftz=true -Xcompiler -rdynamic -lineinfo

运行的前可以先用cuda-memcheck检查下

cuda-memcheck a.out

如果有内存错误，会输出报告，告知你哪一行访问了非法内存，往往问题是出在kernel里，这样就可以快速定位到行，进行排错了。cuda-memcheck还有其他很多功能，目前我就用到这个，其他的等需要的时候再说吧。

逻辑错误

确定好没有内存错误，没有其他错误后，发现结果不对，这也是最头疼的地方，目前为止一直用的printf打印方法来进行观察，没用过cuda-gdb，也在网上搜索了很多关于CUDA的调试方法，也没有个好的结果，用cuda-gdb一是感觉还得学cuda-gdb的使用，而是使用起来也不是那么方便(没有细看cuda-gdb的用法，有兴趣的可以看看)，这里对于printf方法，也是摸索出一些方法。首先写程序的时候这样写CUDA Pro Tip：Write Flexible Kernels with Grid-Stride Loops,这样是为了用printf调试方便。用printf调试的时候，把block数和thread数都设置为1，就可以很好的用printf进行调试了。用printf的时候也可以一个功能一个功能(一段代码一段代码)的进行打印，这样来进行定位。

自己总结的一些小经验，有更好的方法，希望大家指出。

---

如何写出正确的CUDA程序

每次写完CUDA程序后总是各种错误，然后拍错会花很长的时间，好吧，我承认是自己太笨的原因。总结出一些小方法(大神勿喷)
先看一个例子，矩阵*向量，看看调试方法以及自己的思考
kernel1（错误）

__global__ void MxvBlock(const float *h_a, const float *h_b, float *h_c, const int row_size, const int col_size)
{
    int tid = threadIdx.x;
    extern __shared__ float s_tmp[];
    float tmp = 0.0f;
    for (int i=tid; i< row_size; i+=blockDim.x)
        tmp += h_a[i] * h_b[i]; 
    s_tmp[tid] = tmp;
    __syncthreads();

    //reduction，当然对于规约有更好的办法，这里就不写了
    float sum = 0.0f;
    if (tid == 0)
    {   
        for (int j=0; j<blockDim.x; ++j)
        {   
            sum += s_tmp[j];
        }   
    }   
    h_c[bid] = sum;
}

kernel2(错误2)

__global__ void MxvBlock(const float *h_a, const float *h_b, float *h_c, const int row_size, const int col_size)
{
    int tid = threadIdx.x;
    int bid = blockIdx.x;

    extern __shared__ float sr[];
    float tmpSum = 0.0f;
    for (int i=tid; i<col_size; i+=blockDim.x)
        tmpSum += h_a[bid*col_size + i] * h_b[i];
    sr[tid] = tmpSum;
    __syncthreads();

    //reduction, fixed me 
    float sum = 0.0f;
    if (tid == 0)
    {
        for (int i=0; i<blockDim.x; ++i)
            sum += sr[i];
    }
    h_c[bid] = sum;
}

先解释下，本kernel实现矩阵h_a(维度row_size*col_size)乘以向量h_b（维度col_size*1），结果存在h_c（维度row_size*1）中。
本kernel是一个block处理计算一个h_a的行乘以一个h_b的列得到h_c的一个元素
(注意数据划分很重要，只有知道了数据如何划分才知道如何写kernel)
数据划分可分为以下三种情况
1. 一个thread处理一个h_a的行乘以一个h_b的列得到h_c的一个元素
2. 一个block处理一个h_a的行乘以一个h_b的列得到h_c的一个元素
3. 一个warp处理一个h_a的行乘以一个h_b的列得到h_c的一个元素
当然这三种都是可以的，这里选择的是第二种方法。
思考方法如下
1. 首先既然是以block为单位处理h_a的每一行，那blockIdx.x就代表着h_a的行号，就是处理h_a哪一行(明确bid(blockIdx.x),tid(threadIdx.x)的物理含义很重要)，只有知道了h_a的行号和列号才可以定位到h_a的某一个元素，才能进行计算；可以看看上面的错误kernel例子，明显就错了，都没有blockIdx.x出现。
2. kernel2已经很接近正确版本了(结果我还是一步步printf很久才找出问题所在，哎)。
先说说2的思路(以后这样记录下思路，不然每次写CUDA程序的时候都很乱，不像些CPU程序，思路比较明朗)
1. 既然确定是一个block处理一行乘一列，block里的threads来进行这些乘加运算，那就写上下面两句吧，肯定不会出错，然后明确其物理意义

int tid = threadIdx.x;//tid代表h_a列号
int bid = blockIdx.x;//bid代表h_a的行号

1. 计算h_a[i][j] *h_b[j]，这也是kernel的核心，定位h_a[i][j]和h_b[j]

h_a[i][j] = h_a[i*列长+j]
h_b[j] = h_b[j]

既然是以一个block为单位处理h_a的一行乘以h_b的一列，以block的threads来处理乘加运算，那么for循环变量就应该以tid来表示，如下

    extern __shared__ float sr[];
    float tmpSum = 0.0f;
    //以block为单位运行，则stride为blockDim.x
    for (int i=tid; i<col_size; i+=blockDim.x)
        tmpSum += h_a[bid*col_size + i] * h_b[i];
    sr[tid] = tmpSum;
    //用到共享内存一定一定一定别忘了要同步呀
    __syncthreads();

注意for循环是以CUDA Pro Tip：Write Flexible Kernels with Grid-Stride Loops形式运行的,如下还有别的形式

//以thread为单位处理
int tid = threadIdx.x + blockIdx.x * blockDim.x;
for (int i = tid; tid < SIZE; tid+=blockDim.x * gridDim.x)
...

然后注意到tmpSum是每个线程都有的，每个线程分别计算他们自己的tmpSum，之后我们需要对这些线程的tmpSum进行规约，因此需要一个数组先保存这些tmpSum,需要多大呢，有多少个线程就有多少个tmpSum,比如说如果每个block里有256个threads，则需要arr[256]来进行保存。
如上代码所示，用了shared mem，则一定需要使用__syncthreads()来进行同步（PS:有一次调试一个程序，结果变来变去，我就纳闷了好久。。。一步步检查，发现没有对共享变量进行同步，哎只有当自己去写程序的时候才会发现各种问题，道理谁都懂，做出来就各种乱七八糟的错误了）
3. 到这一步了，可以进行规约了，注意shared mem是每个block都有的，如下代码

    float sum = 0.0f;
    if (tid == 0)
    {
        for (int i=0; i<blockDim.x; ++i)
            sum += sr[i];
    }

1. 之后将sum赋值给h_c，如下

    float sum = 0.0f;
    if (tid == 0)
    {
        for (int i=0; i<blockDim.x; ++i)
            sum += sr[i];
    }
    h_c[bid] = sum;

之后各种调试，怎么不对，h_c[]的值都是0呀，一步步printf，看到sum的值是正确的呀，可是为什么到h_c[]这里就出错呢？终于找到原因了，下面这行代码，是每个thread都要运行的，而sum是每个thread都有的变量，比如说thread1其sum就为0，只有thread0的sum为正确的值，这样h_c的值就被其他线程掩盖了

h_c[bid] = sum

因此需要明确每个thread的执行路径，变量是不是每个thread都有，或着说每个block有(shared mem)

kernel3(正确)

__global__ void MxvBlock(const float *h_a, const float *h_b, float *h_c, const int row_size, const int col_size)
{
    int tid = threadIdx.x;                                                                                                                                                    
    int bid = blockIdx.x;

    extern __shared__ float sr[];
    float tmpSum = 0.0f;
    for (int i=tid; i<col_size; i+=blockDim.x)
        tmpSum += h_a[bid*col_size + i] * h_b[i];
    sr[tid] = tmpSum;
    __syncthreads();

    //reduction, fixed me 
    float sum = 0.0f;
    if (tid == 0)
    {   
        for (int i=0; i<blockDim.x; ++i)
            sum += sr[i];
        h_c[bid] = sum;
    }   
}

总结

1. 首先明确数据划分方式，数据分割处理方式，之后明确threadIdx.x，blockIdx.x的物理意义；
2. 明确变量在每个thread或block里的
3. 程序出现逻辑错误时，可以试着在纸上上运行下thread0的路径、thread1的路径、block0的路径、block1的路径，来发现错误

好了，以后再更吧。。。。。
都不知道下次什么时候会写CUDA程序了