[2016/06/25] CUDA矩阵乘法 简单实现

0 随便说两句

嗯，因为项目需要，所以学CUDA。感觉对之后的计算机图形学会稍微有点帮助？如果毕业要换电脑，我一定要记得买一个有Nvidia显卡的电脑，否则这CUDA有些算是白学了哈哈。

刚接触，感觉还是挺难的，主要是一些概念没有理清。这个矩阵乘法，我也算是迷迷糊糊弄懂，特别是最后计算index的部分，实在是。。。线性代数、数学、空间想象能力真的不管在哪个领域都重要到不行啊！以后如果有孩子了一定要好好教育（呵呵。

真的是懵逼了半个星期，今晚终于看懂了点。分享下我觉得入门CUDA（真的是入门！入门）的一些步骤和概念吧。

0  先下好各种驱动程序，安装包，然后在VS里面跑Sample

1 知道CPU，GPU是什么东西，还有一些关键词 host, device, global, shared 等。

2 知道三个基本过程， host到device，device跑核心函数，device到host，已经这三个过程所对应的函数、参数的意义、用法

3 仔细参透grid block 和 thread之间的关系，熟记并理解blockDim.x/y, blockIdx.x/y, threadIdx.x/y各代表什么含义，注意它跟数组的下标是有区别的，我一开始真的想了很久，因为我把这个索引当成是数组下标，就怎么都不明白它为什么是这样表达的。（入门就先别想三维了吧。。。

4 第一个例子：实现向量加法，并且修改device核心函数中的num_blocks, num_threads，以作为理解(3)的参考。

5 第二个例子：实现矩阵乘法，对就是我现在这步。本来跑完了向量加法的你可能会觉得哈哈，挺简单。然后这个矩阵的优化算法（我是指优化算法！不是指每个线程算一个矩阵元素这种简单的。。。不过你可以先实现简单的，再使用优化的。。）的各种Index就把我弄晕了。这不是编程的问题，这就是数学。。。数学不好是硬伤啊！

6 第三个例子：接下来我要实现的FFT。嗯既然还没开始做我就先不写了，调试好了我会整理上来。

1 main函数

首先你要实现的，只是简单的给A矩阵、B矩阵赋值。

为了避免一上来就过大给大脑造成负担：） 让A矩阵的大小为16*16(ha*wa)（16是BLOCK_SIZE，所以A的维度至少也得是16的倍数），B(wa*wb)也是，乘出来的C(ha*wb)显然也是。

为了方便看出结果到底对不对，就不随机赋值矩阵了，每个值都赋1好了。

最后调用一下

MatrixMulCPU(c,a,b,wa,ha,wb);

这个函数，这个函数是在CPU上运行的，所以不需要那个格式<<<>>>。

int main(){    const int wa = 16;const int ha = 16;const int wb = 16;    float a[wa * ha], b[wa * wb], c[ha * wb];int i=0;printf("A:\n");for(;i<wa*ha;i++){if((i+1)%wa==0)printf("\n[%d]",(i+1)/wa);//a[i]=rand()*1.0 / RAND_MAX;a[i]=1;printf("%f ",a[i]);}printf("B:\n");for(i=0;i<wa*wb; i++){if((i+1)%wb==0)printf("\n[%d]",(i+1)/wb);//b[i] = rand()*1.0 / RAND_MAX;b[i]=1;printf("%f ",b[i]);}MatrixMulCPU(c,a,b,wa,ha,wb);getchar();    return 0;}

2 CPU函数

我已经写了比较详细的注释了，基本上跟着以下几个步骤，在明白向量加法之后，这个函数是没有什么问题的了。

void MatrixMulCPU(float* C, float* A, float* B, int wa, int ha, int wb){int size,i=0;cudaError_t err = cudaSuccess;//1 将A和B加载到内存中float* Ad;size = ha * wa * sizeof(float);//A矩阵所占的内存空间err = cudaMalloc((void **)&Ad, size);//在GPU上分配A矩阵的内存空间if(err!=cudaSuccess)printf("A分配内存失败");err = cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);//将CPU的A传给GPU的Adif(err!=cudaSuccess)printf("A to device 失败");float* Bd;size = wa * wb * sizeof(float);//B矩阵所占的内存空间err = cudaMalloc((void**)&Bd, size);//在GPU上分配B矩阵的内存空间if(err!=cudaSuccess)printf("B分配内存失败");err = cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);//将CPU的B传给GPU的Bdif(err!=cudaSuccess)printf("B to device 失败");//2 为GPU上的C分配内存float* Cd;size = ha * wb * sizeof(float);//C矩阵所占的内存空间err = cudaMalloc((void**)&Cd, size);//在GPU上分配C矩阵的内存空间if(err!=cudaSuccess)printf("C分配内存失败");int hc = ha, wc = wb;//3 计算矩阵运算所对应的线程块维度和线程布局dim3 Block(wc / BLOCK_SIZE, hc / BLOCK_SIZE);//根据hc和wc计算在grid中的线程块安排dim3 Thread(BLOCK_SIZE, BLOCK_SIZE);//16X16的线程安排(一个线程块中 thread呈16*16分布)//4 运行GPU计算matrixMulGPU<<<Block,Thread>>>(Ad, Bd, Cd, wa, ha, wb);//5 从GPU中的C读取回CPUerr = cudaMemcpy(C, Cd, ha * wb * sizeof(float), cudaMemcpyDeviceToHost);if(err!=cudaSuccess)printf("C to device 失败");//6 释放Ad,Bd,Cd在GPU上的内存cudaFree(Ad);cudaFree(Bd);cudaFree(Cd);//7 结果输出size = size/sizeof(float);printf("C:\n");for(;i<size;i++){if((i+1)%wb==0)printf("\n[%d]",(i+1)/wb);printf("%f ",C[i]);}}

3 GPU函数

最关键最核心的当然是在于这个函数。。我画了一张图，但是懒得传上来。根据图来计算相应的索引，理解起来会相对容易一些

//每个线程块 计算一个Csub(16*16)Block_size*Block_size//每个线程 计算Csub中的一个元素__global__ void matrixMulGPU(const float* A,const float* B,float *C,int wa,int ha,int wb){<span style="white-space:pre"></span>int bx = blockIdx.x;//block的x索引    int by = blockIdx.y;//block的y索引    int tx = threadIdx.x;//block内的thread的x索引    int ty = threadIdx.y;//block内的thread的y索引//该block要处理的Aint aBegin = wa * (by * BLOCK_SIZE);//A(0,by)    int aEnd = aBegin + wa - 1;    int aStep = BLOCK_SIZE;//offsetA    int bBegin = BLOCK_SIZE * bx;//B(bx,0)    int bStep = BLOCK_SIZE * wb;//offsetB        float Csub = 0;    for (int a = aBegin,b = bBegin; a <= aEnd; a += aStep,b += bStep)    {//As 和 Bs 是共享的、被用来作块内计算的A、B子矩阵        __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];        __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];        //每个线程负责一个元素拷贝        As[ty][tx] = A[a + wa * ty + tx];        Bs[ty][tx] = B[b + wb * ty + tx];        __syncthreads();                //每个线程负责计算一个子块i 和 子块j的子乘积        for (int k = 0; k < BLOCK_SIZE; ++k)        {            Csub += (float)(As[ty][k] * Bs[k][tx]);        }        __syncthreads();    }//全局地址，向全局寄存器写回去    //一个线程负责一个元素，一个block负责一个子块    //int c = wb * BLOCK_SIZE * by + BLOCK_SIZE * bx;    int c= (by * BLOCK_SIZE + ty)*wb + (bx * BLOCK_SIZE + tx);C[c]=Csub;//C[c + wb * ty + tx] = Csub;}

我本来以为，写这个矩阵乘法的程序是我的任务，结果我发现，我要写的是加固程序。。而且还是要在3个级别上实现的。。。。。。我很崩溃，因为我感觉我连入门都很艰难。。。。。。。。。。。。。。