C++实战之OpenCL矩阵相乘

简单概念理解

在opencl中，有个索引空间NDRange的概念，NDRange是一个N维的索引空间，N可以是1，2，3。NDRange由一个长度为N的整数阵列来定义，他指定了索引空间各个维度的宽度，每个work-item的全局id和局部id，都是N维元组。

有多个work-item构成的叫做work-group，作業組的 ID 跟作業項的全局 ID差不多。一個長度為 N 的陣列定義了每個維度上作業組的數 目。作業項在所隸屬的作業組中有一個局部 ID，此 ID 中各維度的取值範圍為0到作業組在相應維 度上的大小減一。因此，作業組的 ID 加上其中一個局部 ID可以唯一確定一個作業項。有兩種途徑 來標識一個作業項:根據全局索引，或根據作業組索引加一個局部索引。

下面一张简单的图介绍下在代码中的应用：

具体实例介绍

==—–矩阵的相乘。==

串行实现代码：

int a[midle*heightA];int b[widthB*midle];//初始化省略for(int l=0;l<heightA;l++){        for(int n = 0;n<widthB;n++){            for(int q=0;q<midle;q++){                result[l*widthB+n] +=a [l*midle+q]*b[q*widthB+n];            }            //std::cout<<"r = "<<result[l*widthB+n]<<std::endl;        }    }

简单从网上找几张图片理解下矩阵乘法的基本步骤和原理：

可以得知矩阵A其中的第i行和矩阵B中的第j列对应项相乘并相加的结果就是结果矩阵C[i][j]的元素了。

opencl实现矩阵相乘

a*b=c

下面是opencl的运行在gpu上内核具体实现，也就是矩阵相乘代码中在gpu中并行执行的单位代码：

__kernel void hello_kernel(__global const int *a,                           __global const int *b,                           __global int *result_matrix,int  result_matrix_row,                           int  result_matrix_col,int  compute_size){    int row = get_global_id(0);    int col = get_global_id(1);    int sum = 0;    for(int i=0;i<compute_size;i++)    {        sum += a[row*compute_size+i] * b[i*result_matrix_col+col];    }    result_matrix[row*result_matrix_col+col] = sum;}

这边的参数：a（M*P）,b（P*N）是输入矩阵，result_matrix 是结果矩阵，result_matrix_row是结果矩阵的行（M），result_matrix_col是结果矩阵的列（N）。compute_size是a和b两个矩阵相同的行列，也就是P，是矩阵相乘中元素想加的个数。

kernel的核心代码即为上述串行代码中去掉外面两个for循环的代码，row和col分别是对应结果矩阵行和列的id，是变化的。目前这个kernel是也是这个work-item每次计算出结果矩阵中的一个元素，可以根据需要，改成每次计算出结果矩阵的一行数据，或其他块数据。我这边有做实验，计算出一个元素和一行元素在我的macbook pro上时间没有什么区别，网上看别的大神用代码说明改成计算一行的的元素，时间上减少了2/3，这部分有待验证。

如上所说，一次执行这个kernel算出结果矩阵当中的一个元素，那么需要执行M*N个这个kernel，意思就是说会有M *N个这样的（work-item）线程同时执行这个kernel，也就是并行执行。

这边为了加深对opencl 运行时的理解分析，对这个kernel做了打印测试分析。不过opencl spec上说printf是从缓冲区flush 出来的，没有顺序的，所以打印的结果仅供参考。
在cl文件中加入以下两句：

#pragma OPENCL EXTENSION cl_khr_fp64: enable#pragma OPENCL EXTENSION cl_amd_printf : enable

就可以在kernel 中使用printf了，语法跟平时一样。

位了方便测试，把矩阵大小改成 6*6 ，10*10，20*20，对比他们的结果：

kernel中的改动：

__kernel void hello_kernel(__global const int *a,                           __global const int *b,                           __global int *result_matrix,int  result_matrix_row,                           int  result_matrix_col,int  compute_size){    int row = get_global_id(0);    int col = get_global_id(1);    printf("row =%d,col=%d\n",row,col);    int sum = 0;    for(int i=0;i<compute_size;i++)    {        sum += a[row*compute_size+i] * b[i*result_matrix_col+col];    }    printf("result_matrix[%d]=%d\n",row*result_matrix_col+col,sum);    result_matrix[row*result_matrix_col+col] = sum;}

下面是结果：

row =0,col=0row =1,col=0row =2,col=0row =3,col=0row =4,col=0row =5,col=0row =0,col=1row =1,col=1row =2,col=1row =3,col=1row =4,col=1row =5,col=1row =0,col=2row =1,col=2row =2,col=2row =3,col=2row =4,col=2row =5,col=2row =0,col=3row =1,col=3row =2,col=3row =3,col=3row =4,col=3row =5,col=3row =0,col=4row =1,col=4row =2,col=4row =3,col=4row =4,col=4row =5,col=4row =0,col=5row =1,col=5row =2,col=5row =3,col=5row =4,col=5row =5,col=5result_matrix[17]=36result_matrix[23]=36result_matrix[29]=36result_matrix[35]=36result_matrix[0]=36result_matrix[6]=36result_matrix[12]=36result_matrix[18]=36result_matrix[24]=36result_matrix[30]=36result_matrix[1]=36result_matrix[7]=36result_matrix[13]=36result_matrix[19]=36result_matrix[25]=36result_matrix[31]=36result_matrix[2]=36result_matrix[8]=36result_matrix[14]=36result_matrix[20]=36result_matrix[26]=36result_matrix[32]=36result_matrix[3]=36result_matrix[9]=36result_matrix[15]=36result_matrix[21]=36result_matrix[27]=36result_matrix[33]=36result_matrix[4]=36result_matrix[10]=36result_matrix[16]=36result_matrix[22]=36result_matrix[28]=36result_matrix[34]=36result_matrix[5]=36result_matrix[11]=36cpu t = 0.00000300gpu t = 0.00001500 Executed program succesfully.

这个结果先是上面一句printf的全部打印冲出来，然后
再是下面打印矩阵值的printf。不急着下结论，再看一组实验结果10*10：

row =6,col=9row =7,col=9row =8,col=9row =9,col=9row =0,col=0row =1,col=0row =2,col=0row =3,col=0row =4,col=0row =5,col=0row =6,col=0row =7,col=0row =8,col=0row =9,col=0row =0,col=1row =1,col=1row =2,col=1row =3,col=1row =4,col=1row =5,col=1row =6,col=1row =7,col=1row =8,col=1row =9,col=1row =0,col=2row =1,col=2row =2,col=2row =3,col=2row =4,col=2row =5,col=2row =6,col=2row =7,col=2row =8,col=2row =9,col=2row =0,col=3row =1,col=3row =2,col=3row =3,col=3row =4,col=3row =5,col=3row =6,col=3row =7,col=3row =8,col=3row =9,col=3row =0,col=4row =1,col=4row =2,col=4row =3,col=4row =4,col=4row =5,col=4row =6,col=4row =7,col=4row =8,col=4row =9,col=4row =0,col=5row =1,col=5row =2,col=5row =3,col=5row =4,col=5row =5,col=5row =6,col=5row =7,col=5row =8,col=5row =9,col=5row =0,col=6row =1,col=6row =2,col=6row =3,col=6row =4,col=6row =5,col=6row =6,col=6row =7,col=6row =8,col=6row =9,col=6row =0,col=7row =1,col=7row =2,col=7row =3,col=7row =4,col=7row =5,col=7row =6,col=7row =7,col=7row =8,col=7row =9,col=7row =0,col=8row =1,col=8row =2,col=8row =3,col=8row =4,col=8row =5,col=8row =6,col=8row =7,col=8row =8,col=8row =9,col=8row =0,col=9row =1,col=9row =2,col=9row =3,col=9row =4,col=9row =5,col=9result_matrix[0]=60result_matrix[10]=60result_matrix[20]=60result_matrix[30]=60result_matrix[40]=60result_matrix[50]=60result_matrix[60]=60result_matrix[70]=60result_matrix[80]=60result_matrix[90]=60result_matrix[1]=60result_matrix[11]=60result_matrix[21]=60result_matrix[31]=60result_matrix[41]=60result_matrix[51]=60result_matrix[61]=60result_matrix[71]=60result_matrix[81]=60result_matrix[91]=60result_matrix[2]=60result_matrix[12]=60result_matrix[22]=60result_matrix[32]=60result_matrix[42]=60result_matrix[52]=60result_matrix[62]=60result_matrix[72]=60result_matrix[82]=60result_matrix[92]=60result_matrix[3]=60result_matrix[13]=60result_matrix[69]=60result_matrix[79]=60result_matrix[89]=60result_matrix[99]=60result_matrix[46]=60result_matrix[56]=60result_matrix[66]=60result_matrix[76]=60result_matrix[86]=60result_matrix[96]=60result_matrix[7]=60result_matrix[17]=60result_matrix[27]=60result_matrix[37]=60result_matrix[47]=60result_matrix[57]=60result_matrix[67]=60result_matrix[77]=60result_matrix[87]=60result_matrix[97]=60result_matrix[8]=60result_matrix[18]=60result_matrix[28]=60result_matrix[38]=60result_matrix[48]=60result_matrix[58]=60result_matrix[68]=60result_matrix[78]=60result_matrix[88]=60result_matrix[98]=60result_matrix[9]=60result_matrix[19]=60result_matrix[29]=60result_matrix[39]=60result_matrix[49]=60result_matrix[59]=60result_matrix[23]=60result_matrix[33]=60result_matrix[43]=60result_matrix[53]=60result_matrix[63]=60result_matrix[73]=60result_matrix[83]=60result_matrix[93]=60result_matrix[4]=60result_matrix[14]=60result_matrix[24]=60result_matrix[34]=60result_matrix[44]=60result_matrix[54]=60result_matrix[64]=60result_matrix[74]=60result_matrix[84]=60result_matrix[94]=60result_matrix[5]=60result_matrix[15]=60result_matrix[25]=60result_matrix[35]=60result_matrix[45]=60result_matrix[55]=60result_matrix[65]=60result_matrix[75]=60result_matrix[85]=60result_matrix[95]=60result_matrix[6]=60result_matrix[16]=60result_matrix[26]=60result_matrix[36]=60cpu t = 0.00000600gpu t = 0.00001600 Executed program succesfully.

看了上面就会发现，也是乱序的，不过仔细一看，其实还是有一点规律的，这边只是让自己对opencl runtime有个初步的认识。后面在对矩阵相乘的优化中，继续学习work-group和work-item等抽象概念。

下面是在cpu上跑的代码：

////  main.cpp//  OpenCL////  Created by wmy on 2017/9/19.//  Copyright © 2017年 wmy. All rights reserved.//#include <OpenCL/OpenCL.h>#include <iostream>#include <fstream>#include <sstream>#include <unistd.h>#include <sys/time.h>#include<time.h>#include<stdio.h>#include<stdlib.h>#include <mach/mach_time.h>#include <boost/algorithm/string.hpp>using namespace std;//const int ARRAY_SIZE = 100000;//4*3---3*5const int midle = 640;const int heightA = 480;const int widthB = 640;//const int heightB = 3;//一、 选择OpenCL平台并创建一个上下文cl_context CreateContext(){    cl_int errNum;    cl_uint numPlatforms;    cl_platform_id firstPlatformId;    cl_context context = NULL;    //选择可用的平台中的第一个    errNum = clGetPlatformIDs(1, &firstPlatformId, &numPlatforms);    if (errNum != CL_SUCCESS || numPlatforms <= 0)    {        std::cerr << "Failed to find any OpenCL platforms." << std::endl;        return NULL;    }    //创建一个OpenCL上下文环境    cl_context_properties contextProperties[] =    {        CL_CONTEXT_PLATFORM,        (cl_context_properties)firstPlatformId,        0    };    context = clCreateContextFromType(contextProperties, CL_DEVICE_TYPE_GPU,                                      NULL, NULL, &errNum);    return context;}//二、 创建设备并创建命令队列cl_command_queue CreateCommandQueue(cl_context context, cl_device_id *device){    cl_int errNum;    cl_device_id *devices;    cl_command_queue commandQueue = NULL;    size_t deviceBufferSize = -1;    // 获取设备缓冲区大小    errNum = clGetContextInfo(context, CL_CONTEXT_DEVICES, 0, NULL, &deviceBufferSize);    if (deviceBufferSize <= 0)    {        std::cerr << "No devices available.";        return NULL;    }    // 为设备分配缓存空间    devices = new cl_device_id[deviceBufferSize / sizeof(cl_device_id)];    errNum = clGetContextInfo(context, CL_CONTEXT_DEVICES, deviceBufferSize, devices, NULL);    //选取可用设备中的第一个    commandQueue = clCreateCommandQueue(context, devices[0], 0, NULL);    *device = devices[0];    delete[] devices;    return commandQueue;}// 三、创建和构建程序对象cl_program CreateProgram(cl_context context, cl_device_id device, const char* fileName){    cl_int errNum;    cl_program program;    std::ifstream kernelFile(fileName, std::ios::in);    if (!kernelFile.is_open())    {        std::cerr << "Failed to open file for reading: " << fileName << std::endl;        return NULL;    }    std::ostringstream oss;    oss << kernelFile.rdbuf();    std::string srcStdStr = oss.str();    const char *srcStr = srcStdStr.c_str();    program = clCreateProgramWithSource(context, 1,                                        (const char**)&srcStr,                                        NULL, NULL);    errNum = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);    return program;}//创建和构建程序对象bool CreateMemObjects(cl_context context, cl_mem memObjects[3],                      int *a, int *b){    memObjects[0] = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,                                   sizeof(int) * midle*heightA, a, NULL);    memObjects[1] = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,                                   sizeof(int) * widthB*midle, b, NULL);    memObjects[2] = clCreateBuffer(context, CL_MEM_READ_WRITE,                                   sizeof(int) * widthB*heightA, NULL, NULL);    return true;}// 释放OpenCL资源void Cleanup(cl_context context, cl_command_queue commandQueue,             cl_program program, cl_kernel kernel, cl_mem memObjects[3]){    for (int i = 0; i < 3; i++)    {        if (memObjects[i] != 0)            clReleaseMemObject(memObjects[i]);    }    if (commandQueue != 0)        clReleaseCommandQueue(commandQueue);    if (kernel != 0)        clReleaseKernel(kernel);    if (program != 0)        clReleaseProgram(program);    if (context != 0)        clReleaseContext(context);}int main(int argc, char** argv){    cl_context context = 0;    cl_command_queue commandQueue = 0;    cl_program program = 0;    cl_device_id device = 0;    cl_kernel kernel = 0;    cl_mem memObjects[3] = { 0, 0, 0 };    cl_int errNum;   // uint64_t t1,t2,t3;    clock_t t1,t2,t3;    const char* filename = "/Users/wangmingyong/Projects/OpenCL/OpenCL/HelloWorld.cl";    // 一、选择OpenCL平台并创建一个上下文    context = CreateContext();    // 二、 创建设备并创建命令队列    commandQueue = CreateCommandQueue(context, &device);    //创建和构建程序对象    program = CreateProgram(context, device, filename);//"HelloWorld.cl");    // 四、 创建OpenCL内核并分配内存空间    kernel = clCreateKernel(program, "hello_kernel", NULL);    //创建要处理的数据    int result[widthB*heightA]{0};    int a[midle*heightA];    int b[widthB*midle];    for (int i = 0; i < heightA; i++)    {        for(int j = 0;j < midle;j++)        {            a[i*midle+j]=2;//10.0f * ((int) rand() / (int) RAND_MAX);        }    }    for (int k = 0; k < midle; k++)    {        for(int m = 0;m < widthB;m++)        {            b[k*widthB+m]=3;//10.0f * ((int) rand() / (int) RAND_MAX);        }    }    t1 = clock();  //mach_absolute_time();    printf("t1 = %.8f\n",(double)t1);    //cpu串行处理代码    for(int l=0;l<heightA;l++){        for(int n = 0;n<widthB;n++){            for(int q=0;q<midle;q++){                result[l*widthB+n] +=a [l*midle+q]*b[q*widthB+n];            }            //std::cout<<"r = "<<result[l*widthB+n]<<std::endl;        }    }    t2 = clock(); //mach_absolute_time();    printf("t2 = %.8f\n",(double)t2);    //创建内存对象    if (!CreateMemObjects(context, memObjects, a, b))    {        Cleanup(context, commandQueue, program, kernel, memObjects);        return 1;    }    // 五、 设置内核数据并执行内核    errNum = clSetKernelArg(kernel, 0, sizeof(cl_mem), &memObjects[0]);    errNum |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &memObjects[1]);    errNum |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &memObjects[2]);    errNum |= clSetKernelArg(kernel, 3, sizeof(int), &heightA);    errNum |= clSetKernelArg(kernel, 4, sizeof(int), &widthB);    errNum |= clSetKernelArg(kernel, 5, sizeof(int), &midle);    size_t globalWorkSize[2];    globalWorkSize[0]= heightA;    globalWorkSize[1]=widthB;   // size_t localWorkSize[2] = { 1,1 };    errNum = clEnqueueNDRangeKernel(commandQueue, kernel, 2, NULL,                                    globalWorkSize, NULL,                                    0, NULL, NULL);    // 六、 读取执行结果并释放OpenCL资源    errNum = clEnqueueReadBuffer(commandQueue, memObjects[2], CL_TRUE,                                 0, widthB*heightA * sizeof(int), result,                                 0, NULL, NULL);//    for(int p=0;p<20;p++){//        cout<<"new ="<<result[p];//    }    t3 = clock();  //mach_absolute_time();    printf("cpu t = %.8f\n",(float)(t2-t1)/CLOCKS_PER_SEC);    printf("gpu t = %.8f \n",(double)(t3-t2)/CLOCKS_PER_SEC);    std::cout << std::endl;    std::cout << "Executed program succesfully." << std::endl;    getchar();    Cleanup(context, commandQueue, program, kernel, memObjects);    return 0;}

下面是具体运行效果：
耗时：

维度 cpu gpu 10*10 0.00001900 0.00064400 100*100 0.00279000 0.00073800 500*500 0.38360301 0.00235700

可以看出随着维度的增加cpu的虐势就体现出来了