▶ 矩阵乘法，按照书里的内容进行了几方面的优化，包括局部内存，矢量数据类型，寄存器，流水线等。
● 最直接的乘法。调用时 main.c 中使用 size_t globalSize[2] = { rowA, colB }, localSize[2] = { 16, 16 }; 。rowA 蕴含在 get_global_id(0) 中了，不再出现在函数中，后面的几种方法也如此。
1 // multiply.cl 2 __kernel void multiply01(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 3 { 4 const int row = get_global_id(0), col = get_global_id(1); 5 int k; 6 float sum; 7 for (k = 0, sum = 0.0f; k < colA; k++) 8 sum += inputA[row * colA + k] * inputB[k * colB + col]; 9 outputC[row * colB + col] = sum; 10 return; 11 }
1 // main.c 2 #include <stdio.h> 3 #include <math.h> 4 #include <stdlib.h> 5 #include <time.h> 6 #include <cl.h> 7 8 const int rowA = 4096, colA = 1024, colB = 2048; 9 //const int rowA = 128, colA = 128, colB = 128; // 测试用，刚够 multiply05 的 1 组 10 const char *sourceText = "D:\Code\OpenCL\multiply.cl"; 11 12 bool floatEq(const float a, const float b)// 相等返回 1 13 { 14 if (b == 0) 15 return fabs(a) < 0.001; 16 return fabs(a / b - 1) < 0.001; 17 } 18 19 int readText(const char* kernelPath, char **pcode)// 读取文本文件放入 pcode，返回字符串长度 20 { 21 FILE *fp; 22 int size; 23 //printf("<readText> File: %s ", kernelPath); 24 fopen_s(&fp, kernelPath, "rb"); 25 if (!fp) 26 { 27 printf("Open kernel file failed "); 28 getchar(); 29 exit(-1); 30 } 31 if (fseek(fp, 0, SEEK_END) != 0) 32 { 33 printf("Seek end of file failed "); 34 getchar(); 35 exit(-1); 36 } 37 if ((size = ftell(fp)) < 0) 38 { 39 printf("Get file position failed "); 40 getchar(); 41 exit(-1); 42 } 43 rewind(fp); 44 if ((*pcode = (char *)malloc(size + 1)) == NULL) 45 { 46 printf("Allocate space failed "); 47 getchar(); 48 exit(-1); 49 } 50 fread(*pcode, 1, size, fp); 51 (*pcode)[size] = ' '; 52 fclose(fp); 53 return size + 1; 54 } 55 56 int main() 57 { 58 int i, j, k, correct; 59 float *A, *B, *C, tempSum; 60 //char info[1024] = { 0 }; 61 clock_t time; 62 63 A = (float*)malloc(sizeof(float) * rowA * colA); 64 B = (float*)malloc(sizeof(float) * colA * colB); 65 C = (float*)malloc(sizeof(float) * rowA * colB); 66 srand(2); 67 for (i = 0; i < rowA * colA; A[i] = rand() & 0xF, i++); 68 for (i = 0; i < colA * colB; B[i] = rand() & 0xF, i++); 69 for (i = 0; i < rowA * colB; C[i] = 0, i++); 70 71 // 初始化平台到创建命令队列 72 cl_int status; 73 cl_uint nPlatform; 74 clGetPlatformIDs(0, NULL, &nPlatform); 75 cl_platform_id *listPlatform = (cl_platform_id*)malloc(nPlatform * sizeof(cl_platform_id)); 76 clGetPlatformIDs(nPlatform, listPlatform, NULL); 77 cl_uint nDevice; 78 clGetDeviceIDs(listPlatform[0], CL_DEVICE_TYPE_ALL, 0, NULL, &nDevice); 79 cl_device_id *listDevice = (cl_device_id*)malloc(nDevice * sizeof(cl_device_id)); 80 clGetDeviceIDs(listPlatform[0], CL_DEVICE_TYPE_ALL, nDevice, listDevice, NULL); 81 cl_context context = clCreateContext(NULL, nDevice, listDevice, NULL, NULL, &status); 82 cl_command_queue queue = clCreateCommandQueue(context, listDevice[0], NULL, &status); 83 84 // 缓冲区 85 cl_mem bufferA = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * rowA * colA, A, &status); 86 cl_mem bufferB = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA * colB, B, &status); 87 cl_mem bufferC = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(float) * rowA * colB, NULL, &status); 88 89 // 程序与内核参数 90 char *code; 91 size_t codeLength = readText(sourceText, &code); 92 cl_program program = clCreateProgramWithSource(context, 1, (const char**)&code, &codeLength, &status); 93 status = clBuildProgram(program, nDevice, listDevice, NULL, NULL, NULL); 94 //clGetProgramBuildInfo(program, listDevice[0], CL_PROGRAM_BUILD_LOG, 1024, info, NULL); 95 //printf(" %s ", info); 96 cl_kernel kernel = clCreateKernel(program, "multiply01", &status); 97 clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufferA); 98 clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufferB); 99 clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufferC); 100 clSetKernelArg(kernel, 3, sizeof(int), &colA); 101 clSetKernelArg(kernel, 4, sizeof(int), &colB); 102 size_t globalSize[2] = { rowA, colB }, localSize[2] = { 16, 16 };// 注意不同函数需要调整工作组网格参数 103 104 // 执行内核 105 time = clock(); 106 status = clEnqueueNDRangeKernel(queue, kernel, 2, NULL, globalSize, localSize, 0, NULL, NULL); 107 clFinish(queue); 108 printf(" Time kernel : %d ms ", clock() - time); 109 110 // 返回并检查结果 111 clEnqueueReadBuffer(queue, bufferC, CL_TRUE, 0, sizeof(float) * rowA * colB, C, 0, NULL, NULL); 112 for (i = 0, correct = 1; i < rowA && correct; i++) 113 { 114 for (j = 0; j < colB && correct; j++) 115 { 116 for (k = 0, tempSum = 0.0f; k < colA; tempSum += A[i * colA + k] * B[k * colB + j], k++); 117 if (!floatEq(tempSum, C[i * colB + j])) 118 { 119 printf("Error at [%d, %d], calculation: %f, reference: %f ", i, j, C[i*colA + j], tempSum); 120 correct = 0; 121 } 122 } 123 } 124 if (correct) 125 printf("Result correct. "); 126 printf("Time total: %d ms ", clock() - time); 127 128 // 释放资源 129 free(A); 130 free(B); 131 free(C); 132 free(listPlatform); 133 free(listDevice); 134 clReleaseContext(context); 135 clReleaseCommandQueue(queue); 136 clReleaseProgram(program); 137 clReleaseKernel(kernel); 138 clReleaseMemObject(bufferA); 139 clReleaseMemObject(bufferB); 140 clReleaseMemObject(bufferC); 141 getchar(); 142 return 0; 143 }
● 输出结果。纯 CPU 串行计算的版本花了 145730 ms，开 O3 后优化为 9157 ms（真是 9 秒）
Time kernel : 228 ms Result correct. Time total: 112593 ms// 时间全花在检查结果上了
● 使用局部内存优化，一个工作组负责输出 outputC 中大小为 TILE_DIM * TILE_DIM 的矩阵。每次从 inputA 和 inputB 中分别取出一个该大小的子矩阵放入局部内存中（Asub 和 Bsub），完成该部分的乘法运算，再抓取下一对子矩阵（inputA 中向右挪动，inputB 中向下挪）。调用时在 mian.c 中修改核函数名字就行
1 // multiply.cl 2 #define TILE_DIM 16 3 4 __kernel void multiply02(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 5 { 6 __local float Asub[TILE_DIM][TILE_DIM], Bsub[TILE_DIM][TILE_DIM]; 7 const int localRow = get_local_id(0), localCol = get_local_id(1); 8 const int groupRow = get_group_id(0) * get_local_size(0), groupCol = get_group_id(1) * get_local_size(1);// 当前工作项的首元素位置 9 const int nTile = colA / TILE_DIM; // 行程长除以块长，即需要计算的方块数 10 int t, k; 11 float acc; 12 for (t = 0, acc = 0.0f; t < nTile; t++) 13 { 14 Asub[localRow][localCol] = inputA[(groupRow + localRow) * colA + TILE_DIM * t + localCol]; // 读取块A，注意列要用 t 来算 15 Bsub[localCol][localRow] = inputB[(TILE_DIM * t + localRow) * colB + groupCol + localCol]; // 读取块B，注意行要用 t 来算，注意交换了 Bsub 的行列，优化访问 16 //Bsub[localRow][localCol] = inputB[(t * TILE_DIM + localRow) * colB + groupCol + localCol];// 不交换 Bsub 行列的情形，与上一行相互替换 17 barrier(CLK_LOCAL_MEM_FENCE); 18 19 for (k = 0; k < TILE_DIM; k++) // 在局部内存上计算累进，acc 为各工作项私有，不会有冲突 20 acc += Asub[localRow][k] * Bsub[localCol][k]; 21 //acc += Asub[localRow][k] * Bsub[k][localCol]; // 不交换 Bsub 行列的情形，与上一行相互替换 22 barrier(CLK_LOCAL_MEM_FENCE); // 保证整个方块算完，才能进下一个方块 23 } 24 outputC[(groupRow + localRow) * colB + groupCol + localCol] = acc; 25 }
● 输出结果，有明显提高。不交换 Bsub 行列的情形时间稍长，约 91 ms。
Time kernel : 81 ms Result correct. Time total: 111484 ms
● 在 multiply02 的基础上，将 Asub 扩展为 float4 类型，四个元素分布在同一列隔 TILE_DIM 的位置上，如 Asub[0][0] == (float4)(A[0][0], A[16][0], A[32][0], A[48][0])，Bsub保持不变（为了方便后面 multiply04 的理解，这里先变A不变B），计算乘法时，4个 Asub 中的元素同时与 Bsub 相乘，最后输出到 outputC 时再分开，相当于一个工作项负责计算 outputC 中同一列中 4 个分散位置上的元素。调用时 main.c 中使用 size_t globalSize[2] = { rowA / 4, colB }, localSize[2] = { 16, 16 }; 。
1 // multiply.cl 2 #define TILE_DIM 16 3 4 __kernel void multiply03(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 5 { 6 __local float4 Asub[TILE_DIM][TILE_DIM]; // Asub 扩展为 float4，相当于原来的 4 倍大小 7 __local float Bsub[TILE_DIM][TILE_DIM]; // Bsub 保持不变 8 const int localRow = get_local_id(0), localCol = get_local_id(1); 9 const int groupRow = get_group_id(0) * get_local_size(0) * 4, globalCol = get_global_id(1); // 注意不同的工作组在行方向上相差为 get_local_size(0) * 4 10 const int nTile = colA / TILE_DIM; 11 int t, k; 12 float4 acc = (float4)(0.f), temp; 13 for (t = 0; t < nTile; t++) 14 { 15 Asub[localRow][localCol] = (float4)(inputA[(groupRow + TILE_DIM * 0 + localRow) * colA + TILE_DIM * t + localCol], // 从 A 中离散的 4 个位置上取元素放入 Asub 16 inputA[(groupRow + TILE_DIM * 1 + localRow) * colA + TILE_DIM * t + localCol], 17 inputA[(groupRow + TILE_DIM * 2 + localRow) * colA + TILE_DIM * t + localCol], 18 inputA[(groupRow + TILE_DIM * 3 + localRow) * colA + TILE_DIM * t + localCol]); 19 Bsub[localCol][localRow] = inputB[(TILE_DIM * t + localRow) * colB + globalCol]; 20 barrier(CLK_LOCAL_MEM_FENCE); 21 22 for (k = 0; k < TILE_DIM; k++) // 在局部内存上计算累进，每次计算四个位置乘法 23 acc += Asub[localRow][k] * Bsub[localCol][k]; 24 barrier(CLK_LOCAL_MEM_FENCE); 25 } 26 outputC[(groupRow + TILE_DIM * 0 + localRow) * colB + globalCol] = acc.x; // 输出时分散到各位置 27 outputC[(groupRow + TILE_DIM * 1 + localRow) * colB + globalCol] = acc.y; 28 outputC[(groupRow + TILE_DIM * 2 + localRow) * colB + globalCol] = acc.z; 29 outputC[(groupRow + TILE_DIM * 3 + localRow) * colB + globalCol] = acc.w; 30 }
● 输出结果，不如 multiply02 的正方形局部内存表现好，主要原因是全局内存的读取和写入时位置太分散。
Time kernel : 103 ms Result correct. Time total: 34299 ms
● 在 multiply03 的基础上，将 Bsub 也扩展为 float4 类型，四个元素分布在同一行连续位置上，与 Asub 不同，如 Bsub[0][0] == (float4)(B[0][0], B[0][1], B[0][2], B[0][3])，计算乘法时，Asub的各元素分别与 Bsub 的各元素相乘（可以使用矢量乘法加快速度），最后输出到 outputC 时再分开。由于 Bsub 的四个元素是横向相邻的，所以可以使用矢量读取和写入函数 vload4 和 vstore4（注意这两个函数的偏移量参数以 (float4 *) 为单位），所以在找准了目标数组中的位置（需要读取或写入的数组的一维索引）后，把该索引除以 4，即为偏移量。相当于一个工作项负责计算 outputC 中连续 4 列的 4 个分散行上共 16 个的元素。调用时 main.c 中使用 size_t globalSize[2] = { rowA / 4, colB / 4 }, localSize[2] = { 16, 16 }; 。
1 // multiply.cl 2 #define TILE_DIM 16 3 4 __kernel void multiply04(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 5 { 6 __local float4 Asub[TILE_DIM][TILE_DIM], Bsub[TILE_DIM][TILE_DIM]; // Asub 是离散的，Bsub是横向连续的 4 个元素 7 const int localRow = get_local_id(0), localCol = get_local_id(1); 8 const int groupRow = get_group_id(0) * get_local_size(0) * 4, groupCol = get_group_id(1) * get_local_size(1) * 4; 9 const int nTile = colA / TILE_DIM; 10 int t, k; 11 float4 acc[4] = { (float4)(0.f) }, tempA, tempB; 12 for (t = 0; t < nTile; t++) 13 { 14 Asub[localRow][localCol] = (float4)(inputA[(groupRow + TILE_DIM * 0 + localRow) * colA + TILE_DIM * t + localCol], 15 inputA[(groupRow + TILE_DIM * 1 + localRow) * colA + TILE_DIM * t + localCol], 16 inputA[(groupRow + TILE_DIM * 2 + localRow) * colA + TILE_DIM * t + localCol], 17 inputA[(groupRow + TILE_DIM * 3 + localRow) * colA + TILE_DIM * t + localCol]); 18 Bsub[localCol][localRow] = vload4(((TILE_DIM * t + localRow) * colB + groupCol + localCol * 4) / 4, inputB);// 向量读取 19 barrier(CLK_LOCAL_MEM_FENCE); 20 21 for (k = 0; k < TILE_DIM; k++) 22 { 23 tempA = Asub[localRow][k], tempB = Bsub[localCol][k]; 24 acc[0] += tempA.x * tempB; // 注意这里是一个标量乘以一个向量，得到一个向量 25 acc[1] += tempA.y * tempB; 26 acc[2] += tempA.z * tempB; 27 acc[3] += tempA.w * tempB; 28 } 29 barrier(CLK_LOCAL_MEM_FENCE); 30 } 31 for (k = 0; k < 4; k++) 32 vstore4(acc[k], ((groupRow + TILE_DIM * k + localRow) * colB + groupCol + localCol * 4) / 4, outputC); // 向量写入 33 }
● 输出结果，有显著进步
Time kernel : 40 ms Result correct. Time total: 32563 ms
● 作为 multiply03 和 multiply04 的一个补充，我一开始看到 multiply04 中 Asub 这么跳着取值的时候是很懵的，直到看了后面的 multiply06 才知道这么做是方便对该方法进行扩展。在写 multiply06 之前，我按自己的理解，让 Asub 不跳着取，而是类似 Bsub 那样，在竖方向上连续取 4 个元素，看看计算效率如何。调用时 main.c 使用 size_t globalSize[2] = { rowA / 4, colB / 4 }, localSize[2] = { 16, 16 }; 。
1 // multiply.cl 2 #define TILE_DIM 16 3 __kernel void multiply05(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 4 { 5 __local float4 Asub[TILE_DIM][TILE_DIM], Bsub[TILE_DIM][TILE_DIM]; // Asub，Bsub是连续的 4 个元素 6 const int localRow = get_local_id(0), localCol = get_local_id(1); 7 const int groupRow = get_group_id(0) * get_local_size(0) * 4, groupCol = get_group_id(1) * get_local_size(1) * 4; 8 const int nTile = colA / TILE_DIM; 9 int t, k; 10 float4 acc[4] = { (float4)(0.f) }, tempA, tempB; 11 for (t = 0; t < nTile; t++) 12 { 13 Asub[localRow][localCol] = (float4)(inputA[(groupRow + localRow * 4 + 0) * colA + TILE_DIM * t + localCol], // 改了读取位置 14 inputA[(groupRow + localRow * 4 + 1) * colA + TILE_DIM * t + localCol], 15 inputA[(groupRow + localRow * 4 + 2) * colA + TILE_DIM * t + localCol], 16 inputA[(groupRow + localRow * 4 + 3) * colA + TILE_DIM * t + localCol]); 17 Bsub[localCol][localRow] = vload4(((TILE_DIM * t + localRow) * colB + groupCol + localCol * 4) / 4, inputB); 18 barrier(CLK_LOCAL_MEM_FENCE); 19 20 for (k = 0; k < TILE_DIM; k++) // 乘法计算上没有变化 21 { 22 tempA = Asub[localRow][k], tempB = Bsub[localCol][k]; 23 acc[0] += tempA.x * tempB; 24 acc[1] += tempA.y * tempB; 25 acc[2] += tempA.z * tempB; 26 acc[3] += tempA.w * tempB; 27 } 28 barrier(CLK_LOCAL_MEM_FENCE); 29 } 30 for (k = 0; k < 4; k++) 31 vstore4(acc[k], ((groupRow + localRow * 4 + k) * colB + groupCol + localCol * 4) / 4, outputC); // 改了写入位置 32 }
● 输出结果，跟 multiply04 差不多
Time kernel : 38 ms Result correct. Time total: 33787 ms
● 高潮部分来了， multiply04 虽然使用了局部内存，利用 float4 类型的存储和内建函数，但认为它没有充分利用工作项的寄存器数量，所以最简单的改善方法就是要求每个工作项去计算 outputC 中更多的元素（multiply04 中是 4 * 4 = 16 个）。考虑将 Asub 在竖方向上扩展为原来的 nAsub 倍，将 Bsub 在横方向上扩展为原来的 nBsub 倍（为了优化访存，右边扩出来的部分以此放到原来 Bsub 的下方，看起来就像 Asub 那样在竖方向上扩展了，其实不是），一个工作项负责计算 outputC 中 4 * 4 * nAsub * nBsub 个元素，适当调整扩展倍数，可以提高计算效率。注意这里当 nAsub == nBsub == 1 时该算法退化为 multiply04。调用时main.c 中使用 size_t globalSize[2] = { rowA / 4 / 2, colB / 4 / 2}, localSize[2] = { 16, 16 }; ，其中的两个 2 分别等于 nAsub 和 nBsub。
1 // multiply.cl 2 #define TILE_DIM 16 3 #define nAsub 2 // Asub 扩展倍数 4 #define nBsub 2 // Bsub 扩展倍数 5 6 __kernel void multiply06(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 7 { 8 __local float4 Asub[TILE_DIM * nAsub][TILE_DIM], Bsub[TILE_DIM * nBsub][TILE_DIM];// 两个局部内存矩阵都在竖方向上扩展，但是扩展的定义不同 9 const int localRow = get_local_id(0), localCol = get_local_id(1); 10 const int groupRow = get_group_id(0) * get_local_size(0) * 4 * nAsub, groupCol = get_group_id(1) * get_local_size(1) * 4 * nBsub;// groupRow 的跨度需要调整 11 const int nTile = colA / TILE_DIM; 12 int t, k, subi, subj; 13 float4 acc[4 * nAsub][nBsub], tempA, tempB; // acc 扩展为二维 float4 数组 14 for (k = 0; k < 4 * nAsub; k++) // acc 清零 15 { 16 for (subj = 0; subj < nBsub; subj++) 17 acc[k][subj] = (float4)(0.f); 18 } 19 for (t = 0; t < nTile; t++) 20 { 21 for (subi = 0; subi < nAsub; subi++) // 对 Asub 添加一层循环，每次迭代填入 TILE_DIM * TILE_DIM 个元素 22 { 23 Asub[TILE_DIM * subi + localRow][localCol] = (float4)(inputA[(groupRow + TILE_DIM * (4 * subi + 0) + localRow) * colA + TILE_DIM * t + localCol], 24 inputA[(groupRow + TILE_DIM * (4 * subi + 1) + localRow) * colA + TILE_DIM * t + localCol], 25 inputA[(groupRow + TILE_DIM * (4 * subi + 2) + localRow) * colA + TILE_DIM * t + localCol], 26 inputA[(groupRow + TILE_DIM * (4 * subi + 3) + localRow) * colA + TILE_DIM * t + localCol]); 27 } 28 for (subj = 0; subj<nBsub; subj++) // 对 Bsub 类似 29 Bsub[TILE_DIM * subj + localRow][localCol] = vload4(((TILE_DIM * t + localRow) * colB + groupCol + TILE_DIM * 4 * subj + localCol * 4) / 4, inputB); 30 barrier(CLK_LOCAL_MEM_FENCE); 31 32 for (k = 0; k < TILE_DIM; k++) 33 { 34 for (subi = 0; subi < nAsub; subi++)// 计算时添加两层循环，分别完成 Asub 和 Bsub 中各分块的元素的乘法（比较纠结 k 和 subi/subj 哪个循环放里边比较好） 35 { 36 for (subj = 0; subj < nBsub; subj++) 37 { 38 tempA = Asub[TILE_DIM * subi + localRow][k]; 39 tempB = Bsub[TILE_DIM * subj + k][localCol]; 40 acc[4 * subi + 0][subj] += tempA.x * tempB; 41 acc[4 * subi + 1][subj] += tempA.y * tempB; 42 acc[4 * subi + 2][subj] += tempA.z * tempB; 43 acc[4 * subi + 3][subj] += tempA.w * tempB; 44 } 45 } 46 } 47 barrier(CLK_LOCAL_MEM_FENCE); 48 } 49 for (k = 0; k < 4 * nAsub; k++) // 写入时添加一层循环（实际上是添加了两层循环，但是 k 和 subi 的循环合并了，就像对 acc 清零的时候一样） 50 { 51 for (subj = 0; subj < nBsub; subj++) 52 vstore4(acc[k][subj], ((groupRow + TILE_DIM * k + localRow) * colB + groupCol + TILE_DIM * 4 * subj + localCol * 4) / 4, outputC); 53 } 54 }
● 输出结果，大有进步。当扩展倍数为（2，1）时没有明显提高（40 ms）；倍数为（1，2）或（2，2）时提高最大，如下所示；倍数为（2，4），（4，2），（4，4）时与下面类似（26 ~ 28 ms）；倍数为（8，4）时超出局部内存大小，clCreateKernel 返回 -45（CL_INVALID_PROGRAM_EXECUTABLE）。
Time kernel : 26 ms Result correct. Time total: 34380 ms
● 优化流水线。multiply04，multiply05，multiply06 的思想是：取第 0 块子矩阵放入局部内存 → 计算局部内存中第 0 块子矩阵的积累 acc → 取第 1 块子矩阵放入局部内存 → 计算局部内存中第 1 块子矩阵的积累 acc → ……。每次迭代都需要两个栅栏，即等待所有数据填入局部内存以及等待所有工作项完成计算。这里先考虑添加一个预读步骤，变成：预读第 0 块子矩阵到私有存储器（寄存器）→ 将预读的第 0 块子矩阵存入局部内存 → 计算局部内存中第 0 块子矩阵长度积累 acc → 预读第 1 块子矩阵到私有存储器（寄存器）→ 将预读的第 1 块子矩阵存入局部内存 → 计算局部内存积累 acc → ……，看起来暂时没有什么改善，但是注意到，在局部内存中进行计算的时候，虽然局部内存被占用了不能修改，但是寄存器可以空出来去预读下一块子矩阵，一旦计算完成，且寄存器预读完成，就可以用寄存器中预读好的数据去覆盖局部内存，相当于部分掩盖了从全局内存中读取数据的延迟，注意这时栅栏存在于 “用寄存器数据去覆盖局部内存数据” 的前后，分别用于保证 “上一个子矩阵的计算完成”，以及 “可以开始本次子矩阵的计算”。根据上面这种预读的思想，调整流水线为：预读第 0 块子矩阵到私有存储器 → 将预读的第 0 块子矩阵存入局部内存 → 预读第 1 块子矩阵到私有存储器 → 计算局部内存中第 0 块子矩阵长度积累 acc → 将预读的第 1 块子矩阵存入局部内存 → 预读第 2 块子矩阵到私有存储器 → 计算局部内存中第 1 块子矩阵长度积累 acc → ……。调用时参数同 multiply06。
1 // multiply.cl 2 #define TILE_DIM 16 3 #define nAsub 2 4 #define nBsub 2 5 6 #define PRE_LOAD(startA, startB) 7 { 8 for (subi = 0; subi < nAsub; subi++) 9 { 10 fetchA[subi].x = inputA[startA + (4 * subi + 0) * TILE_DIM * colA + localRow * colA + localCol]; 11 fetchA[subi].y = inputA[startA + (4 * subi + 1) * TILE_DIM * colA + localRow * colA + localCol]; 12 fetchA[subi].z = inputA[startA + (4 * subi + 2) * TILE_DIM * colA + localRow * colA + localCol]; 13 fetchA[subi].w = inputA[startA + (4 * subi + 3) * TILE_DIM * colA + localRow * colA + localCol]; 14 } 15 for (subj = 0; subj < nBsub; subj++) 16 fetchB[subj] = vload4((startB + localRow * colB + subj * TILE_DIM * 4 + localCol * 4) / 4, inputB); 17 } 18 19 #define STORE 20 { 21 barrier(CLK_LOCAL_MEM_FENCE); 22 for (int subi = 0; subi < nAsub; subi++) 23 Asub[subi * TILE_DIM + localRow][localCol] = fetchA[subi]; 24 for (int subj = 0; subj < nBsub; subj++) 25 Bsub[subj * TILE_DIM + localRow][localCol] = fetchB[subj]; 26 barrier(CLK_LOCAL_MEM_FENCE); 27 }; 28 29 #define COMPUTE 30 { 31 for(int k = 0; k < TILE_DIM; k++) 32 { 33 for (int subi = 0; subi < nAsub; subi++) 34 { 35 for (int subj = 0; subj < nBsub; subj++) 36 { 37 tempA = Asub[subi * TILE_DIM + localRow][k]; 38 tempB = Bsub[subj * TILE_DIM + k][localCol]; 39 acc[4 * subi + 0][subj] += tempA.x * tempB; 40 acc[4 * subi + 1][subj] += tempA.y * tempB; 41 acc[4 * subi + 2][subj] += tempA.z * tempB; 42 acc[4 * subi + 3][subj] += tempA.w * tempB; 43 } 44 } 45 } 46 } 47 48 __kernel void multiply07(__global float *inputA, __global float *inputB, __global float *outputC, int colA, int colB) 49 { 50 __local float4 Asub[TILE_DIM * nAsub][TILE_DIM], Bsub[TILE_DIM * nBsub][TILE_DIM]; 51 const int localRow = get_local_id(0), localCol = get_local_id(1); 52 const int groupRow = get_group_id(0) * get_local_size(0) * 4 * nAsub, groupCol = get_group_id(1) * get_local_size(1) * 4 * nBsub; 53 const int nTile = colA / TILE_DIM; 54 int k, startA, startB, subi, subj; 55 float4 acc[4 * nAsub][nBsub], fetchA[nAsub], fetchB[nBsub], tempA, tempB; // fetchA 和 fetchB 为从 inputA 和 inputB 中预读的值 56 for (k = 0; k < 4 * nAsub; k++) 57 { 58 for (subj = 0; subj < nBsub; subj++) 59 acc[k][subj] = (float4)(0.f); 60 } 61 62 startA = groupRow * colA, startB = groupCol; 63 PRE_LOAD(startA, startB) // 从 inputA 和 inputB 中预读第一个子矩阵到 tempA 和 tempB中 64 for (; startA < groupRow + colA - TILE_DIM; startA += TILE_DIM, startB += TILE_DIM * colB)// 推进 TILE，最后一次不需要预读，循环减少一次 65 { 66 STORE // 将预读数据写入 Asub 和 Bsub 67 PRE_LOAD(startA + TILE_DIM, startB + TILE_DIM * colB) // 预读下一个子矩阵，inputA 中向右挪动，inputB 中向下挪 68 COMPUTE // 利用当前 Asub 和 Bsub 中的值计算积累 acc 69 } 70 STORE // 存储最后一个子矩阵 71 COMPUTE // 计算最后一个子矩阵的积累 acc 72 73 for (k = 0; k < 4 * nAsub; k++) // 数据写入 outputC，与以前相同 74 { 75 for (subj = 0; subj < nBsub; subj++) 76 vstore4(acc[k][subj], ((groupRow + TILE_DIM * k + localRow) * colB + groupCol + TILE_DIM * 4 * subj + localCol * 4) / 4, outputC); 77 } 78 }
● 输出结果，仍然有进步
Time kernel : 13 ms Result correct. Time total: 34598 ms
● 笔记
■ 可惜向量数据类型没有内建的点积之类的运算，不然可以把 Asub 取成横向连续 4 个元素，Bsub 取成竖向连续 4 个元素
■ 上面所有结果没有考虑 warming up 的问题。试了一下，multiply04 首次跑需要 40 ms 左右，但如果在其之后再重复跑 100 次，总共只要 2198 ms
■ 书上代码非常难看，变量名各种 ab，ae，ii，jj （我知道是 A_begin，A_end，还有用于区别 i，j），还有 bx，by，tx，ty（移植自 CUDA 的 blockIdx，threadIdx？）；输入 B 和输出 C 的形参莫名其妙就变成 float4 类型了，用着还挺直率
1 tb[ty][tx] = b[j + ty * N_float4 + tx];// 这种代码看着真恶心，关键还是向量操作（索引已经从 float* 转化为 float4*）
■ 中括号记法，这是我在这次矩阵乘法编程中使用的重要方法，记录下来方便以后使用和分析，并以 multiply06 中复杂的索引计算为例进行说明。写到另一篇（http://www.cnblogs.com/cuancuancuanhao/p/9063167.html）里面去了。