华南理工大学 陈虎 博士 tommychen74@yahoo.com.cn CUDA例子程序——矩阵乘法 华南理工大学 陈虎 博士 tommychen74@yahoo.com.cn

一个简单的矩阵乘法例子 用矩阵乘法说明CUDA编程中内存和线程管理的基本特性： 共享存储器的用法 本地存储器、寄存器的用法 线程ID的用法 主机和设备之间数据传输的API 为了方便，以方形矩阵说明

方形矩阵乘法 矩阵P = M * N 大小为 WIDTH x WIDTH 在没有采用分片优化算法的情况下： 一个线程计算P矩阵中的一个元素

串行版本的矩阵乘法 // 宿主机的双精度矩阵乘法 void MatrixMulOnHost(float* M, float* N, float* P, int Width){ for (int i = 0; i < Width; ++i) for (int j = 0; j < Width; ++j){ double sum = 0; for (int k = 0; k < Width; ++k){ double a = M[i * width + k]; double b = N[k * width + j]; sum += a * b; } P[i * Width + j] = sum; N k j WIDTH M P i WIDTH k WIDTH WIDTH

向GPU传输矩阵数据 void MatrixMulOnDevice(float* M, float* N, float* P, int Width) { int size = Width * Width * sizeof(float); float* Md, Nd, Pd; //设置调用内核函数时的线程数目 dim3 dimBlock(Width, Width); dim3 dimGrid(1, 1); //在设备存储器上给M和N矩阵分配空间，并将数据复制到设备存储器中 cudaMalloc(&Md, size); cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice); cudaMalloc(&Nd, size); cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice); //在设备存储器上给P矩阵分配空间 cudaMalloc(&Pd, size);

计算结果向主机传输 //内核函数调用，将在后续部分说明 //只使用了一个线程块(dimGrid)，此线程块中有Width*Width个线程 MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd，Width); // 从设备中读取P矩阵的数据 cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost); // 释放设备存储器中的空间 cudaFree(Md); cudaFree(Nd); cudaFree (Pd); }

矩阵乘法的内核函数 // 矩阵乘法的内核函数——每个线程都要执行的代码 __global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width) { // 2维的线程ID号 int tx = threadIdx.x; int ty = threadIdx.y; // Pvalue用来保存被每个线程计算完成后的矩阵的元素 float Pvalue = 0;

内核函数 (继上.) //每个线程计算一个元素 for (int k = 0; k < Width; ++k) { 内核函数 (继上.) //每个线程计算一个元素 for (int k = 0; k < Width; ++k) { float Melement = Md[ty * Width + k]; float Nelement = Nd[k * Width + tx]; Pvalue += Melement * Nelement; } // 将计算结果写入设备存储器中 Pd[ty * Width + tx] = Pvalue; Nd k WIDTH tx Md Pd This should be emphasized! Maybe another slide on “This is the first thing” ty ty WIDTH tx k WIDTH WIDTH

只使用了一个线程块 一个线程块中的每个线程计算Pd中的一个元素 每个线程 缺点： 载入矩阵Md中的一行 载入矩阵Nd中的一列 Grid 1 一个线程块中的每个线程计算Pd中的一个元素 每个线程 载入矩阵Md中的一行 载入矩阵Nd中的一列 为每对Md和Nd元素执行了一次乘法和加法 缺点： 计算和片外存储器存访问比例接近1：1，受存储器延迟影响很大； 矩阵的大小受到线程块所能容纳最大线程数（512个线程）的限制 Block 1 Thread (2, 2) 48 WIDTH Md Pd

处理任意大小的方形矩阵 让每个线程块计算结果矩阵中的一个大小为(TILE_WIDTH)2的子矩阵 总共有(WIDTH/TILE_WIDTH)2 个线程块 Nd WIDTH Pd Md 需要注意的是：当WIDTH/TILE_WIDTH大于最大的网格数量（64K）时，需要在内核函数附近设置一个循环！（见最后的总结） by TILE_WIDTH ty WIDTH bx tx WIDTH WIDTH

G80显卡的存储器瓶颈 所有的线程都要访问全局存储器获取输入矩阵元素 G80显卡的峰值速度为346.5GFlops 效率仅为6% 每一次的单精度浮点乘法和加法需要两次的内存访问 (8 bytes) 全局存储器的访问带宽为86.4 GB/s 每秒钟可以读取21.6G个浮点数 每秒钟最多可以完成21.6GFlops G80显卡的峰值速度为346.5GFlops 效率仅为6% 全局存储器成为计算瓶颈 要充分使用高带宽的片上局部存储器 Grid Block (0, 0) Block (1, 0) Shared Memory Shared Memory Registers Registers Registers Registers Thread (0, 0) Thread (1, 0) Thread (0, 0) Thread (1, 0) Host Global Memory Constant Memory

使用共享存储器以便重用全局存储器中的数据 每个输入元素都需要被WIDTH个线程读取 将每个元素都装载到共享存储器中，让很多线程都使用本地数据以便减少存储带宽 使用分片算法 N WIDTH P M ty WIDTH tx WIDTH WIDTH

分片矩阵乘法 每个线程块计算一个大小为TILE_WIDTH的方形子矩阵Pdsub 每个线程计算Pdsub子矩阵中的一个元素 Md Nd Pd Pdsub TILE_WIDTH WIDTH bx tx 1 TILE_WIDTH-1 2 by ty TILE_WIDTHE 每个线程块计算一个大小为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 全局存储器不再是性能瓶颈！

内核函数线程数配置 //每个线程块有TILE_WIDTH2个线程 dim3 dimBlock(TILE_WIDTH, TILE_WIDTH); //有(Width/TILE_WIDTH)2个线程块 dim3 dimGrid(Width/TILE_WIDTH, Width/TILE_WIDTH); //调用内核函数 MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd，Width);

内核函数 //获得线程块号 int bx = blockIdx.x; int by = blockIdx.y; //获得块内的线程号 int tx = threadIdx.x; int ty = threadIdx.y; //Pvalue：线程计算完成后的子矩阵元素——自动变量 float Pvalue = 0; //循环，遍历M和N的所有子矩阵 for (int m = 0; m < Width/TILE_WIDTH; ++m) { //此处代码在下面 };

将数据装入共享存储器 // 获取指向当前矩阵M子矩阵的指针Msub Float* Mdsub = GetSubMatrix(Md, m, by, Width); //获取指向当前矩阵N的子矩阵的指针Nsub Float* Ndsub = GetSubMatrix(Nd, bx, m, Width); //共享存储器空间声明 __shared__float Mds[TILE_WIDTH][TILE_WIDTH]; __shared__float Nds[TILE_WIDTH][TILE_WIDTH]; // 每个线程载入M的子矩阵的一个元素 Mds[ty][tx] = GetMatrixElement(Mdsub, tx, ty); //每个线程载入N的子矩阵的一个元素 Nds[ty][tx] = GetMatrixElement(Ndsub, tx, ty);

从局部存储器中计算结果 //同步，在计算之前，确保子矩阵所有的元素都已载入共享存储器中__syncthreads(); //每个线程计算线程块内子矩阵中的一个元素 for (int k = 0; k < TILE_WIDTH; ++k) Pvalue += Mds[ty][k] * Nds[k][tx]; //同步，确保重新载入新的M和N子矩阵数据前，上述计算操作已全部完成 __syncthreads();

一些其它代码 GetSubMatrix(Md, x, y, Width) x*TILE_WIDTH GetSubMatrix(Md, x, y, Width) 获取第(x, y)号子矩阵的起始地址 Md + y*TILE_WIDTH*Width + x*TILE_WIDTH); GetMatrixElement(Mdsub,tx,ty,Width) 获取子矩阵中某个元素的地址 *(Mdsub+ty*Width+tx)); Md y*TILE_WIDTH MSub TILE_WIDTH ty TILE_WIDTH tx Width

CUDA 代码 – 保存结果 上面的代码在G80上运行的速度是45 GFLOP！ // 获取指向矩阵P的子矩阵的指针 Matrix Psub = GetSubMatrix(P, bx, by); //向全局存储器写入线程块计算后的结果子矩阵 //每个线程写入一个元素 SetMatrixElement(Psub, tx, ty, Pvalue); 上面的代码在G80上运行的速度是45 GFLOP！