406130914138刘志伟GPU上矩阵乘法的设计与实现

刘志伟：GPU上矩阵乘法的设计与实现1GPU上矩阵乘法的设计与实现刘志伟（南昌大学信息与工程学院计算机科学与技术406130914138）1概述矩阵乘法是科学计算中的最基本的操作，在许多领域中有广泛的应用。对于矩阵乘法的研究有几个方向。一个是研究矩阵乘法的计算复杂度，研究矩阵乘法的时间复杂度的下界，这方面的工作有strassen算法等。另外一个方向是根据不同的处理器体系结构，将经典的矩阵乘法高效的实现出来，这方面的结果体现在许多高效的BLAS库。许多高效的BLAS库都根据体系结构的特点高效的实现了矩阵乘法，比如GotoBLAS，ATLAS等。众所周知，矩阵乘法是一种大计算量的算法，也是很耗时的运算。CPU提高单个核心性能的主要手段比如提高处理器工作频率及增加指令级并行都遇到了瓶颈，当遇到运算量大的计算，CPU进行大矩阵的乘法就变得相当耗时，运算效率很低下。随着多核CPU和众核GPU的快速发展，计算行业正在从只使用CPU的“中央处理”向CPU与GPU并用的“协同处理”发展，并行系统已成为主流处理器芯片。传统的GPU架构受其硬件架构的影响不能有效利用其资源进行通用计算，NVIDIA(英伟达)公司推出的统一计算设备架构CUDA(ComputeUnifiedDeviceArchitectures)，使得GPU具备更强的可编程性，更精确和更高的性能，应用领域也更加广泛。2CUDA简介2.1CUDA编程模型NVIDIA的CUDA架构通过对硬件的重新组织把GPU带到了更加一般的应用领域。它试图通过提供一般的高层和低层的API来访问GPU的并行元素以减轻问题映射到GPU上的不方便。现在的GPU特别适合计算密集型、高度并行化的计算。CUDA提供了对显卡的抽象，把显卡作为能够同时执行成千上万轻量级线程的设备，如图1所示。这些线程组织成块，块中的每个线程能够访问它们自身的寄存器和块的共享存储器，每个线程也能够和它相邻的线程进行同步。一个内核函数的代码由一个或多个这样的线程执行。由于具有CUDA能力的设备允许DRAM的类属存取(聚集和分散)，所以每个线程都能访问GPU板卡上的显存和纹理存储器。刘志伟：GPU上矩阵乘法的设计与实现2图1CUDA编程模型在实际应用中，首先对问题进行分析，哪些部分可以在GPU上并行实现，哪些在CPU上执行。一旦确定程序中的并行部分，就可以考虑把这部分的计算工作交给GPU。在GPU上运行的CUDA并行计算函数称为kernel(内核函数)。一个kernel函数并不是一个完整的程序，而是整个CUDA程序中的一个可以被并行执行的步骤。如图1所示，一个完整的CUDA程序是由GPU中一系列的kernel函数并行步骤和CPU端串行处理步骤共同组成的，这些处理步骤会按照程序中相应语句的顺序依次执行，满足顺序一致性。2.2CUDA存储器模型除了编程模型外，CUDA还规定了存储器模型，如图2所示。线程在执行时将会访问到处于不同存储空间中的数据。刘志伟：GPU上矩阵乘法的设计与实现3图2GPU的存储器层次结构每个线程都拥有自己私有的存储器、寄存器和局部存储器；每个block拥有一块共享存储器（SharedMemory)；grid中所有的线程都可以访问同一块全局存储器（GlobalMemory）。除此之外，还有两种可以被所有线程访问的只读存储器：常数存储器（ConstantMemory）和纹理存储器（TextureMemory）。这几种存储器的访问速度不同，大小不同，在程序中使用不同的存储器对速度的影响较大。寄存器速度最快，但容量最小，较少使用到；共享存储器容量较大，共享存储器的特点是速度较快且对于同一个block中的所有线程能够共享，这对程序中数据的分配十分有利。常数存储器和纹理存储器分别用于存放常数和图像。3矩阵乘法在GPU上的实现本文主要探讨矩阵乘法如何在GPU上实现，故设计如下3个矩阵，如表1所示。表1实验所用示例矩阵矩阵数据类型维度matrixAfloat10241024matrixBfloat10241024matrixCfloat10241024设有矩阵ijAann，ijBbnn，ijCcnn则矩阵乘法的计算式为：1nijikkjkcab。图3给出了此矩阵乘法在CPU上实现的代码，从中可以看出：矩阵乘刘志伟：GPU上矩阵乘法的设计与实现4法的计算量复杂度为O(n3)。而矩阵乘法的访存复杂度为O(n2)，所以矩阵乘法的计算访存比复杂度为O(n)，是一个典型的计算密集型应用。fori=1toi=ndoforj=1toj=ndoC(i,j)=0fork=1tok=ndoC(i,j)=A(i,k)*B(k,j)+C(i,j)ENDdoENDdoENDdo图3矩阵乘法在CPU上的实现方阵matrixA和matrixB相乘，并把相应的结果传给matrixC方阵。在CPU上实现的源码如图4所示，其中matrixA和matrixB方阵的元素均为随机生成。voidmultiMatrix(float*matrixA,float*matrixB,float*matrixC,intn){inti,j,k;for(i=0;in;i++){for(j=0;jn;j++){floattemp=0;for(k=0;kn;k++){temp+=matrixA[i*n+k]*matrixB[k*n+j];}matrixC[i*n+j]=temp;}}}图4在CPU上运行的multiMatrix函数下面将对multiMatrix函数进行改写，使其能够在GPU上运行。改写之后的代码如图5所示。#includecuda_runtime.h#includedevice_launch_parameters.h#defineN16//内核函数__global__voidmultiMatrix(float*matrixA,float*matrixB,float*matrixC,intn){刘志伟：GPU上矩阵乘法的设计与实现5intRow=blockIdx.y*N+threadIdx.y;intCol=blockIdx.x*N+threadIdx.x;floatsum=0.0;for(intk=0;kn;k++){sum+=matrixA[Row*N+k]*matrixB[k*N+Col];}//每个线程负责计算matrixC中的一个元素matrixC[Row*N+Col]=sum;}//multiMatrix函数（CUDA中）//在外部调用，使用externexternCvoidmultiMatrix_CUDA(float*matrixA,float*matrixB,float*matrixC,intn){cudaSetDevice(0);//设置目标GPUfloat*matrixA_cuda,*matrixB_cuda,*matrixC_cuda;intsize=n*n*sizeof(float);cudaMalloc((void**)&matrixA_cuda,size);cudaMalloc((void**)&matrixB_cuda,size);cudaMalloc((void**)&matrixC_cuda,size);cudaMemcpy(matrixA_cuda,matrixA,size,cudaMemcpyHostToDevice);cudaMemcpy(matrixB_cuda,matrixB,size,cudaMemcpyHostToDevice);dim3dimGrid(n/N,n/N);//网格的维度dim3dimBlock(N,N);//块的维度multiMatrixKerneldimGrid,dimBlock(matrixA_cuda,matrixB_cuda,matrixC_cuda,n);cudaMemcpy(matrixC,matrixC_cuda,size,cudaMemcpyHostToDevice);//释放设备上的矩阵cudaFree(matrixA_cuda);cudaFree(matrixB_cuda);cudaFree(matrixC_cuda);}图5在GPU上运行的multiMatrix函数4实验结果与分析刘志伟：GPU上矩阵乘法的设计与实现6实验环境如表2所示。表2实验环境参数表CPUIntel(R)Core(TM)i3-3240GPUNVIDIAQuadro410内存6G操作系统Windows7(ServicePack1)开发软件VisualStudio2010专业版CUDA版本CUDADriverVersion6.0在相同的测试环境下，本文给出了multiMatrix函数在CPU和GPU上的对比，如表3所示。表3相同测试环境下，multiMatrix函数在CPU和GPU上的对比N值GPU时间(ms)16074622062345538814351676432795在相同测试条件下，CPU的时间维持在18135ms注：N值只影响GPU里面的线程分布，并不影响CPU的性能。故六次测量中，CPU的时间基本一致。根据表3的数据，建立了不同N值，代码在GPU上运行时间如下图。刘志伟：GPU上矩阵乘法的设计与实现7图6不同N值，代码在GPU上运行的时间不同N值，所对应加速比情况，如表4所示。表4不同N值，所对应加速比情况N值加速比10.29920.87943.275812.6381623.7373222.811根据表4的数据，建立了不同N值，所对应加速比情况，如图7所示。图7不同N值，加速比情况可以看出，随着N值的增大，代码在GPU上的运行时间逐渐减少，当增加到16时，运行时间趋于稳定，变化不大。同样的，随着N值增大，加速比逐渐增大，但N到16时，加速比趋于稳定，不再增加。5总结对于计算密集型且数据耦合度低的应用，使用CPU-GPU异构计算平台即采用CUDA架构进行优化，效果会非常明显。矩阵乘法就是这样的一种应用，可以参考表4中的数据，加速比可高达23.737。CPU-GPU异构计算平台中，由CPU负责执行事务管理和复杂逻辑处理等并行性较弱的计算，由GPU负责进行计算密集型的大规模数据并行计算。想要发挥异构平台的性能，必须针对具体的应用，具体分析，并不是每一个应用都适合CPU-GPU异构并行化的。