矩阵乘法性能优化课件

CUDA程序设计主要内容GPGPU及CUDA介绍CUDA编程模型多线程及存储器硬件GPGPU及CUDA介绍多核时代多个适当复杂度、相对低功耗内核并行工作配置并行硬件资源提高处理能力核心时钟频率基本不变nVidiaGT200Quad-coreOpteronIBMCellBroadbandEngineGPU与CPU硬件架构的对比CPU：更多资源用于缓存及流控制GPU：更多资源用于数据计算适合具备可预测、针对数组的计算模式CacheALUControlALUALUALUDRAMCPUDRAMGPU应用范围CPU:controlprocessor不规则数据结构不可预测存取模式递归算法分支密集型算法单线程程序GPU:dataprocessor规则数据结构可预测存取模式油气勘探、金融分析、医疗成像、有限元、基因分析、地理信息系统、…GPGPU(GeneralPurposeComputingonGPU)GPGPU核心思想用图形语言描述通用计算问题把数据映射到vertex或者fragment处理器缺点硬件资源使用不充分存储器访问方式严重受限难以调试和查错高度图形处理和编程技巧CUDA(ComputeUnifiedDeviceArchitecture)CUDA有效结合CPU+GPU编程串行部分在CPU上运行并行部分在GPU上运行CPUSerialCodeGrid0......GPUParallelKernelKernelA<<<nBlk,nTid>>>(args);Grid1CPUSerialCodeGPUParallelKernelKernelB<<<nBlk,nTid>>>(args);CUDA极大提高了现有应用的效果Gridding1FFTCartesianScanData(a)(b)(b)IterativeReconstruction(c)SpiralScanDataSpiralscandata+Gridding+FFTReconstructionrequireslittlecomputationBasedonFig1ofLustigetal,FastSpiralFourierTransformforIterativeMRImageReconstruction,IEEEInt’lSymp.onBiomedicalImaging,2004MRIReconstructionAdvancedMRIReconstructionFFTCartesianScanData(a)SpiralScanDataIterativeReconstruction(c)Gridding(b)(b)Spiralscandata+IterativereconReconstructionrequiresalotofcomputationAdvancedMRIReconstructionComputeQAcquireDataComputeFHdFindρMorethan99.5%oftimeQ只和扫描参数有关FHd是数据相关的使用线性求解器计算ρHaldar,etal,“Anatomically-constrainedreconstructionfromnoisydata,”MRinMedicine.Codefor(p=0;p<numP;p++){for(d=0;d<numD;d++){exp=2*PI*(kx[d]*x[p]+

ky[d]*y[p]+

kz[d]*z[p]);

cArg=cos(exp);

sArg=sin(exp);

rFhD[p]+=rRho[d]*cArg–

iRho[d]*sArg;

iFhD[p]+=iRho[d]*cArg+

rRho[d]*sArg;}}__global__voidcmpFhD(float*gx,gy,gz,grFhD,giFhD){

intp=blockIdx.x*THREADS_PB+threadIdx.x;

//registerallocateimage-spaceinputs&outputsx=gx[p];y=gy[p];z=gz[p];

rFhD=grFhD[p];iFhD=giFhD[p];

for(intd=0;d<SCAN_PTS_PER_TILE;d++){//s(scandata)isheldinconstantmemoryfloatexp=2*PI*(s[d].kx*x+

s[d].ky*y+

s[d].kz*z);

cArg=cos(exp);sArg=sin(exp);rFhD+=s[d].rRho*cArg–s[d].iRho*sArg;

iFhD+=s[d].iRho*cArg+s[d].rRho*sArg;}

grFhD[p]=rFhD;giFhD[p]=iFhD;}CPUGPU性能提升情况S.S.Stone,etal,“AcceleratingAdvancedMRIReconstructionusingGPUs,”ACMComputingFrontierConference2008,Italy,May2008.计算结果对比CUDA成功案例广泛应用于生命科学、机械、石油、金融、数学、天文和通信等行业MRI(磁共振成像)

GRAPPA自动校准

加速网格化快速重建ComputedTomography(CT)

GE

Digisens

SnapCTStone,UIUCBatenburg,Sijbersetal

医疗成像量子化学KYasuda,NagoyaU,Japan双电子积分RI-MP2correlationenergyinQ-Chem3.1LeslieVogt,Harvard

现有的分子动力学软件

NAMD/VMD(alpharelease)GROMACS(alpharelease)HOOMD

OpenMM:分子建模

https:///home/openmm分子动力学

MonteCalo模拟投资组合优化期权及衍生品定价对冲基金风险分析CUDA中的随机数发生器SciFinance的MonteCalo定价模型SciCompCo.金融序列对比蛋白质对接生物系统的随机仿真(SSA)

人体视觉皮层的自组织计算模型分析基因表达的DNA微阵列工具Schatzetal,UMaryland生物信息学和生命科学3DLattice-Boltzman解算器基于Lattice-Boltzman的PDE解算器用于照明的LatticeBoltzman

Navier-Stokes解算器等离子体湍流建模ThibaultandSenocakTolkeandKrafczy流体动力学

GPMAD:离子束动力学模拟

FDTD法进行的光散射模拟

Acceleware的解算器FDTD加速Accelerware电磁学和电磁力学天气研究与预测模型(WRF)

25%~30%的性能提升海啸模拟天气,大气,海洋科学与空间建模加密编码模式匹配CUDA编程模型CUDA设备与线程计算设备（device）作为CPU（host）的协处理器有独立的存储设备（devicememory）同时启动大量线程计算密集部分使用大量线程并行的kernelGPU与CPU线程的区别GPU的线程非常轻量，线程切换~1cycle，而CPU需要~1000cycleGPU上的线程数>1000时才能有效利用GPU的计算能力StreamingProcessor（SP）Afullypipelined,single-issue,inorder

microprocessor2ALUsandaFPURegisterfile32-bitscalarprocessingNoinstructionfetchand

schedulingNocacheStreamingMultiprocessor（SM）AnarrayofSPs8streamingprocessor2SpecialFunctionUnits(SFU)A16KBread/writesharedmemoryNotacacheButasoftware-manageddatastoreMultithreadingissuingunitInstructionandconstantcacheCUDA程序基本结构串行部分在CPU上运行（host）并行部分在GPU上运行（device）CPUSerialCode(host)Grid0......GPUParallelKernel(device)KernelA<<<nBlk,nTid>>>(args);Grid1CPUSerialCode(host)GPUParallelKernel(device)

KernelB<<<nBlk,nTid>>>(args);C扩展Declspecsglobal,device,shared,local,constantKeywordsthreadIdx,blockIdxIntrinsics__syncthreadsRuntimeAPIMemory,symbol,executionmanagementFunctionlaunch__device__floatfilter[N];__global__voidconvolve(float*image){__shared__floatregion[M];...

region[threadIdx]=image[i];__syncthreads()...

image[j]=result;}//AllocateGPUmemoryvoid*myimage=cudaMalloc(bytes)//100blocks,10threadsperblockconvolve<<<100,10>>>(myimage);CUDA程序的编译使用nvcc编译工具

nvcc<filename>.cu[-oexcutable]调试选项：-g（debug）、-deviceemu（CPU模拟GPU）并行线程组织并行性的维度一维 y=a+b二维

P=MN三维

CTorMRI并行线程组织结构Thread:并行的基本单位Threadblock:互相合作的线程组CooperativeThreadArray(CTA)允许彼此同步通过快速共享内存交换数据以1维、2维或3维组织最多包含512个线程Grid:一组threadblock以1维、2维或3维组织共享全局内存Kernel:在GPU上执行的核心程序Onekernel↔onegrid线程层次BlockandThreadIDsBlocks和Threads具有IDsthreadIdx,blockIdxBlockID:1Dor2DThreadID:1D,2Dor3D由此决定相应处理数据CUDA线程组织CUDAkernel函数由一系列线程组成单指令多数据流（SPMD）通过IDs确定处理的数据线程可划分为不同的Block在同一个block中，可以通过sharememory、atomicoperation和barriersynchronization进行协同…floatx=input[threadID];floaty=func(x);output[threadID]=y;…threadIDThreadBlock0……floatx=input[threadID];floaty=func(x);output[threadID]=y;…ThreadBlock1…floatx=input[threadID];floaty=func(x);output[threadID]=y;…ThreadBlockN-1012345670123456701234567一个简单的例子——IncrementArrayElements//CPUprogramvoidinc_cpu(float*a,floatb,intN){for(intidx=0;idx<N;idx++)a[idx]=a[idx]+b;}voidmain(){…

inc_cpu(a,b,N);}//CUDAprogram__global__voidinc_gpu(float*a,floatb,intN){

intidx=blockIdx.x*blockDim.x+threadIdx.x;

if(idx<N)a[idx]=a[idx]+b;}voidmain(){

…dim3dimBlock(blocksize);dim3dimGrid(ceil(N/(float)blocksize));

inc_gpu<<<dimGrid,dimBlock>>>(a,b,N);}CUDA线程的同步void__syncthreads();Barriersynchronization同步thread

block之内的所有线程避免访问共享内存时发生RAW/WAR/WAW冒险__shared__floatscratch[256];scratch[threadID]=begin[threadID];__syncthreads();intleft=scratch[threadID-1];在此等待，直至所有线程到达才开始执行下面的代码存储器模型与内存分配R/Wper-threadregisters1-cyclelatencyR/Wper-threadlocalmemorySlow–registerspillingtoglobalmemoryR/Wper-blocksharedmemory1-cyclelatencyButbankconflictsmaydragdownR/Wper-gridglobalmemory~500-cyclelatencyButcoalescingaccessingcouldhidelatencyReadonlyper-gridconstantandtexture

memories~500-cyclelatencyButcachedGPUGlobalMemory分配cudaMalloc()分配显存中的globalmemory两个参数对象数组指针数组尺寸cudaFree()释放显存中的globalmemory一个参数对象数组指针代码实例分配6464单精度浮点数组数组指针Md建议用“d”表示GPU显存数据结构GPUGlobalMemory分配intBLOCK_SIZE=64;float*Md;intsize=BLOCK_SIZE*BLOCK_SIZE*sizeof(float);cudaMalloc((void**)&Md,size);cudaFree(Md);Host-Device数据交换cudaMemcpy()在存储器直接传输数据四个参数目的对象数组指针源对象数组指针数组尺寸传输方向Host到HostHost到DeviceDevice到HostDevice到Device代码实例M.elements:CPU主存Md:GPU显存符号常数:cudaMemcpyHostToDevice和cudaMemcpyDeviceToHostHost-Device数据交换cudaMemcpy(Md,M.elements,size,cudaMemcpyHostToDevice);cudaMemcpy(M.elements,Md,size,cudaMemcpyDeviceToHost);CUDA变量与函数CUDA引入的变量修饰词__device__储存于GPU上的globalmemory空间和应用程序具有相同的生命期(lifetime)可被grid中所有线程存取,CPU代码通过runtime函数存取__constant__储存于GPU上的constantmemory空间和应用程序具有相同的生命期(lifetime)可被grid中所有线程存取,CPU代码通过runtime函数存取__shared__储存于GPU上threadblock内的共享存储器和threadblock具有相同的生命期(lifetime)只能被threadblock内的线程存取无修饰（Local变量）储存于SM内的寄存器和localmemory和具有相同的生命期(lifetime)Thread私有Built-indim3Type定义grid和thread

block的组织dim3dimGrid(2,2);dim3dimBlock(4,2,2);kernelFunction<<<dimGrid,dimBlock>>>(…);CUDA函数定义Executedon

the:Onlycallablefrom

the:__device__floatDeviceFunc()devicedevice__global__voidKernelFunc()devicehost__host__floatHostFunc()hosthost__global__定义kernel函数必须返回void__device__和__host__可以组合使用则被定义的函数在CPU和GPU上都被编译__device__函数不能用&运算符取地址限制不支持递归调用不支持静态变量(staticvariable)不支持可变长度参数函数调用typeva_list(stdarg.h)doubleaverage(intcount,...)CUDA函数定义Kernel函数调用调用时必须给出线程配置方式__global__voidKernelFunc(...);dim3DimGrid(100,50);//5000threadblocksdim3DimBlock(4,8,8);//256threadsperblocksize_t

SharedMemBytes=64;//64bytesofsharedmemoryKernelFunc<<<DimGrid,DimBlock,SharedMemBytes>>>(...);CUDA数学函数pow,sqrt,cbrt,hypot,exp,exp2,expm1,log,log2,log10,

log1p,sin,cos,tan,asin,acos,atan,atan2,sinh,cosh,tanh,asinh,acosh,atanh,ceil,floor,trunc,round,etc.只支持标量运算许多函数有一个快速、较不精确的对应版本以“__”为前缀，如__sin()编译开关-use_fast_math强制生成该版本的目标码每个多处理器包含两个超越函数计算单元CUDA程序设计实例——方阵相乘P=M*N（长宽均为WIDTH）计算策略每个线程计算矩阵P中的一个元素MNPWIDTHWIDTHWIDTHWIDTH第一步：CPU实现MNPWIDTHWIDTHWIDTHWIDTHikkj//Matrixmultiplicationonthe(CPU)hostindoubleprecisionvoidMatrixMulOnHost(float*M,float*N,float*P,intWidth){for(inti=0;i<Width;++i)for(intj=0;j<Width;++j){doublesum=0;for(intk=0;k<Width;++k){doublea=M[i*width+k];doubleb=N[k*width+j];sum+=a*b;}

P[i*Width+j]=sum;}}第二步：将矩阵数据传给显存voidMatrixMulOnDevice(float*M,float*N,float*P,intWidth){

intsize=Width*Width*sizeof(float);

float*Md,Nd,Pd;

…1.//AllocateandLoadM,Ntodevicememory

cudaMalloc(&Md,size);

cudaMemcpy(Md,M,size,cudaMemcpyHostToDevice);

cudaMalloc(&Nd,size);

cudaMemcpy(Nd,N,size,cudaMemcpyHostToDevice);//AllocatePonthedevice

cudaMalloc(&Pd,size);第三步：将计算结果传回内存2.//Kernelinvocationcode–tobeshownlater…3.//ReadPfromthedevice

cudaMemcpy(P,Pd,size,cudaMemcpyDeviceToHost);//Freedevicematrices

cudaFree(Md);cudaFree(Nd);cudaFree(Pd);}第四步：kernel函数//Matrixmultiplicationkernel–perthreadcode__global__voidMatrixMulKernel(float*Md,float*Nd,float*Pd,intWidth){//2DThreadID

int

tx=threadIdx.x;

int

ty=threadIdx.y;//Pvalueisusedtostoretheelementofthematrix//thatiscomputedbythethreadfloatPvalue=0;NdMdPdWIDTHWIDTHWIDTHWIDTHtytxtytxkk

for(intk=0;k<Width;++k)

{

floatMelement=Md[ty*Width+k];

floatNelement=Nd[k*Width+tx];

Pvalue+=Melement*Nelement;

}

Pd[ty*Width+tx]=Pvalue;}第四步：kernel函数（续）第五步：调用kernel函数2.//Kernelinvocationcode

//Setuptheexecutionconfigurationdim3dimBlock(Width,Width);dim3dimGrid(1,1);//Launchthedevicecomputationthreads!

MatrixMulKernel<<<dimGrid,dimBlock>>>(Md,Nd,Pd);局限性每个线程都需要读：Md矩阵的一行Nd矩阵的一列计算与访存比约为1:1矩阵规模受限于每个block允许的thread数目

Grid1Block148Thread(2,2)WIDTHMdPdNd参考资料——CUDASDKSDK中包含许多CUDA范例多线程及存储器硬件StreamingMultiprocessor执行ThreadBlocks线程以block为单位分配到SM视资源需求,一个SM分配至多8个blockSMinG80可以接受768个线程256(threads/block)*3blocks128(threads/block)*6blocks,etc线程并发运行SM分配并维护线程IDSM管理并调度线程ThreadBlockSizeConsiderations对于矩阵乘法,哪个thread

block尺寸最好:8X8,16X16或者32X32?8X8:64threads/block.每个SM至多接受768threads,即12

blocks。但是,SM至多接受8blocks,所以实际上仅有512

threads16X16:256threads/block.每个SM至多接受768threads,即3

blocks→只要其它计算资源许可，可以满负荷工作32X32:1024threads/block.SM无法处理线程调度和执行Threadblock内部线程组织为32-threadwarpsAnimplementationdecision-notpartofCUDAWarp是SM调度的基本单位Warp就是一条32路SIMD指令Half-warp是warp的前一半或后一半访问存储器的基本单位Warp调度和执行下条指令中全部操作数就位的warps拥有执行资格Warp中全部线程执行同一指令同一warp内的分支语句可能执行不同的指令路径不同指令路径的线程只能顺序执行每次执行warp中一条可能的路径N条指令路径→1/Nthroughput应尽量避免在同一warp内出现分支SM存储器资源Registerandlocalmemory:per-thread线程私有编译器自行分配e.g.floata;Sharedmemory:per-blockBlock内所有线程共享使数据尽量靠近处理器动态分配到blockse.g.__shared__floatregion[M];ConstantcacheTexturecache寄存器阵列G80中每个SM配置8192个寄存器当前设计选择，不属于CUDA寄存器动态分配到划归SM的blocks中一旦分配到某一block,不能被其它blocks访问同一block内部的线程只能使用分配给该线程的寄存器实例：假设每个block有1616个thread每个thread使用10个寄存器，那么每个block需要使用2560个寄存器，因此每个SM能容纳3个block，也就是768个thread假如每个thread多使用1个寄存器，那么每个block需要使用2816个寄存器，SM就只能容纳2个block，造成并行度下降1/3应该综合考虑并行度与访存开销的影响。假如在上面的情况下多使用1个寄存器能够使访存次数减少一半，那么实际性能反而有所提高存储器模型回顾R/Wper-blocksharedmemory1-cyclelatencyButbankconflictsmaydragdownR/Wper-gridglobalmemory~500-cyclelatencyButcoalescingaccessingcouldhidelatency性能优化思路SharedMemory比GlobalMemory快几百倍线程之间通过SharedMemory合作使用一个或少量线程装载和计算threadblock内全部线程共享的数据利用SharedMemory提高性能每个元素都需要被多个线程重复使用将元素存入sharedmemory供线程共享分块计算实例：矩阵乘法性能优化MNPWIDTHWIDTHWIDTHWIDTHtytx每个block计算一块小方矩阵Pd每个thread读入Pd的一个元素假设M和N的大小是小矩阵大小的整数倍矩阵乘法的分块计算MdNdPdPdsubTILE_WIDTHWIDTHWIDTHTILE_WIDTHTILE_WIDTHbxtx01TILE_WIDTH-12