CUDA实例-矩阵乘法.

本课程得到Intel®大学合作计划支持课程网站：例子程序——矩阵乘法华南理工大学陈虎博士[e-mail]一个简单的矩阵乘法例子用矩阵乘法说明CUDA编程中内存和线程管理的基本特性：共享存储器的用法本地存储器、寄存器的用法线程ID的用法主机和设备之间数据传输的API为了方便，以方形矩阵说明M方形矩阵乘法矩阵P=M*N大小为WIDTHxWIDTH在没有采用分片优化算法的情况下：一个线程计算P矩阵中的一个元素M和N需要从全局存储器载入WIDTH次NPWIDTHWIDTHWIDTHWIDTH串行版本的矩阵乘法MNPWIDTHWIDTHWIDTHWIDTH//宿主机的双精度矩阵乘法voidMatrixMulOnHost(float*M,float*N,float*P,intWidth){for(inti=0;iWidth;++i)for(intj=0;jWidth;++j){doublesum=0;for(intk=0;kWidth;++k){doublea=M[i*width+k];doubleb=N[k*width+j];sum+=a*b;}P[i*Width+j]=sum;}}ikkjvoidMatrixMulOnDevice(float*M,float*N,float*P,intWidth){intsize=Width*Width*sizeof(float);float*Md,Nd,Pd;//设置调用内核函数时的线程数目dim3dimBlock(Width,Width);dim3dimGrid(1,1);//在设备存储器上给M和N矩阵分配空间，并将数据复制到设备存储器中cudaMalloc(&Md,size);cudaMemcpy(Md,M,size,cudaMemcpyHostToDevice);cudaMalloc(&Nd,size);cudaMemcpy(Nd,N,size,cudaMemcpyHostToDevice);//在设备存储器上给P矩阵分配空间cudaMalloc(&Pd,size);向GPU传输矩阵数据计算结果向主机传输//内核函数调用，将在后续部分说明//只使用了一个线程块(dimGrid)，此线程块中有Width*Width个线程MatrixMulKerneldimGrid,dimBlock(Md,Nd,Pd，Width);//从设备中读取P矩阵的数据cudaMemcpy(P,Pd,size,cudaMemcpyDeviceToHost);//释放设备存储器中的空间cudaFree(Md);cudaFree(Nd);cudaFree(Pd);}矩阵乘法的内核函数//矩阵乘法的内核函数——每个线程都要执行的代码__global__voidMatrixMulKernel(float*Md,float*Nd,float*Pd,intWidth){//2维的线程ID号inttx=threadIdx.x;intty=threadIdx.y;//Pvalue用来保存被每个线程计算完成后的矩阵的元素floatPvalue=0;NdMdPdWIDTHWIDTHWIDTHWIDTH内核函数(继上.)//每个线程计算一个元素for(intk=0;kWidth;++k){floatMelement=Md[ty*Width+k];floatNelement=Nd[k*Width+tx];Pvalue+=Melement*Nelement;}//将计算结果写入设备存储器中Pd[ty*Width+tx]=Pvalue;}tytxtytxkk只使用了一个线程块一个线程块中的每个线程计算Pd中的一个元素每个线程载入矩阵Md中的一行载入矩阵Nd中的一列为每对Md和Nd元素执行了一次乘法和加法缺点：计算和片外存储器存访问比例接近1：1，受存储器延迟影响很大；矩阵的大小受到线程块所能容纳最大线程数（512个线程）的限制Grid1Block13254242648Thread(2,2)WIDTHMdPdNd处理任意大小的方形矩阵让每个线程块计算结果矩阵中的一个大小为(TILE_WIDTH)2的子矩阵每个线程块中有(TILE_WIDTH)2个线程总共有(WIDTH/TILE_WIDTH)2个线程块MdNdPdWIDTHWIDTHWIDTHWIDTHtytxbybx需要注意的是：当WIDTH/TILE_WIDTH大于最大的网格数量（64K）时，需要在内核函数附近设置一个循环！（见最后的总结）TILE_WIDTH11GridGlobalMemoryBlock(0,0)SharedMemoryThread(0,0)RegistersThread(1,0)RegistersBlock(1,0)SharedMemoryThread(0,0)RegistersThread(1,0)RegistersHostConstantMemoryG80显卡的存储器瓶颈所有的线程都要访问全局存储器获取输入矩阵元素每一次的单精度浮点乘法和加法需要两次的内存访问(8bytes)全局存储器的访问带宽为86.4GB/s每秒钟可以读取21.6G个浮点数每秒钟最多可以完成21.6GFlopsG80显卡的峰值速度为346.5GFlops效率仅为6%全局存储器成为计算瓶颈要充分使用高带宽的片上局部存储器M12使用共享存储器以便重用全局存储器中的数据每个输入元素都需要被WIDTH个线程读取将每个元素都装载到共享存储器中，让很多线程都使用本地数据以便减少存储带宽使用分片算法NPWIDTHWIDTHWIDTHWIDTHtytx13MdNdPdPdsubTILE_WIDTHWIDTHWIDTHTILE_WIDTHTILE_WIDTHbxtx01TILE_WIDTH-12012byty210TILE_WIDTH-1210TILE_WIDTHTILE_WIDTHTILE_WIDTHEWIDTHWIDTH分片矩阵乘法每个线程块计算一个大小为TILE_WIDTH的方形子矩阵Pdsub每个线程计算Pdsub子矩阵中的一个元素假设Md和Nd的大小都是TILE_WIDTH的倍数G80中首先需要考虑的事项每个线程块内应该有较多的线程TILE_WIDTH=16时有16*16=256个线程应该有分解为若干个线程块一个1024*1024大小的Pd矩阵有64*64=4096个线程块每个线程块从全局存储器将矩阵M和N的一小块读入到局部存储器中，然后完成计算从全局存储器中读出2*256=512个单精度浮点数；完成256*(2*16)=8,192次浮点计算操作；浮点操作：全局存储器读出操作＝16:1全局存储器不再是性能瓶颈！15内核函数线程数配置//每个线程块有TILE_WIDTH2个线程dim3dimBlock(TILE_WIDTH,TILE_WIDTH);//有(Width/TILE_WIDTH)2个线程块dim3dimGrid(Width/TILE_WIDTH,Width/TILE_WIDTH);//调用内核函数MatrixMulKerneldimGrid,dimBlock(Md,Nd,Pd，Width);内核函数//获得线程块号intbx=blockIdx.x;intby=blockIdx.y;//获得块内的线程号inttx=threadIdx.x;intty=threadIdx.y;//Pvalue：线程计算完成后的子矩阵元素——自动变量floatPvalue=0;//循环，遍历M和N的所有子矩阵for(intm=0;mWidth/TILE_WIDTH;++m){//此处代码在下面};将数据装入共享存储器//获取指向当前矩阵M子矩阵的指针MsubFloat*Mdsub=GetSubMatrix(Md,m,by,Width);//获取指向当前矩阵N的子矩阵的指针NsubFloat*Ndsub=GetSubMatrix(Nd,bx,m,Width);//共享存储器空间声明__shared__floatMds[TILE_WIDTH][TILE_WIDTH];__shared__floatNds[TILE_WIDTH][TILE_WIDTH];//每个线程载入M的子矩阵的一个元素Mds[ty][tx]=GetMatrixElement(Mdsub,tx,ty);//每个线程载入N的子矩阵的一个元素Nds[ty][tx]=GetMatrixElement(Ndsub,tx,ty);18从局部存储器中计算结果//同步，在计算之前，确保子矩阵所有的元素都已载入共享存储器中__syncthreads();//每个线程计算线程块内子矩阵中的一个元素for(intk=0;kTILE_WIDTH;++k)Pvalue+=Mds[ty][k]*Nds[k][tx];//同步，确保重新载入新的M和N子矩阵数据前，上述计算操作已全部完成__syncthreads();19一些其它代码GetSubMatrix(Md,x,y,Width)获取第(x,y)号子矩阵的起始地址Md+y*TILE_WIDTH*Width+x*TILE_WIDTH);GetMatrixElement(Mdsub,tx,ty,Width)获取子矩阵中某个元素的地址*(Mdsub+ty*Width+tx));MdTILE_WIDTHx*TILE_WIDTHy*TILE_WIDTHTILE_WIDTHWidthMSubtytx20CUDA代码–保存结果//获取指向矩阵P的子矩阵的指针MatrixPsub=GetSubMatrix(P,bx,by);//向全局存储器写入线程块计算后的结果子矩阵//每个线程写入一个元素SetMatrixElement(Psub,tx,ty,Pvalue);上面的代码在G80上运行的速度是45GFLOP！