标签: python opencl
Author | shaniadolphin |
---|---|
349948204@qq.com |
矩阵乘法
对于以下是一个常见的线性方程组,
用矩阵表示就是:
推导出矩阵乘法的计算:
进一步的,使用一般的表现形式(A为m*n,B为n*p):
表示了矩阵A和矩阵B相乘后,新矩阵各元素等于该元素所在A的列与B的行对应相乘后进行累加。
OpenCL的程序
OpenCL是一个为异构平台编写程序的框架,支持大量不同的应用,提供了两种层面的并行机制:任务并行与数据并行,用来加快数据的处理。
以下是OpenCL编程中常见的词条,规范了编程所涉及对象的名称和概念:
- Platform (平台):主机加上OpenCL框架管理下的若干设备,通过平台,主机应用程序可以与设备共享资源并在设备上执行kernel。基本上一个厂商对应一个Platform,比如Intel, AMD,Mali,PowerVR,一个主机可以有多个平台,比如有独显的笔记本同时有核显和独显。
- Device(设备):计算单元(Compute Units)的集合,GPU是典型的device。但是Intel和AMD的多核CPU也提供OpenCL接口,所以也可以作为Device,所以一个平台上可以有多个设备,比如PC上有两个英伟达的显卡。
- Context(上下文):OpenCL的Platform上共享和使用资源的环境,包括kernel、device、memory objects、command queue等。一般一个Platform对应一个Context。
- Command Queue(指令队列):在指定设备上管理多个指令(Command)。队列里指令执行可以顺序也可以乱序。一个设备可以对应多个指令队列。
- Program:OpenCL程序,由kernel函数、其他函数、声明和变量等组成。
- Kernel(内核函数):运行在设备端的函数。
- Memory Object(内存对象):在主机和设备之间传递数据的对象,一般映射到OpenCL程序中的global memory。有两种具体的类型:Buffer Object(缓存对象)和Image Object(图像对象)。
- NDRange:主机端运行设备端kernel函数的主要接口。
- WaitForEvents(同步):在将一个 OpenCL 命令提交到命令队列的时候,用来标识是以阻塞还是非阻塞的方式执行,如果后续处理依赖该内核函数返回的数据就要使用阻塞的方式,否则线程可以马上处理其它的事情,以提高并行效率。
面向异构平台的应用一般按以下步骤实现:
- 查找构成异构系统的平台(clGetPlatformIDs);
- 获得平台特征,使kernel函数能够适应不同硬件单元的特定特性以获得最优性能(clGetDeviceIDs);
- 创建上下文,建立平台设备在主机端的调用接口(clCreateContext);
- 创建将在平台上运行的程序(clCreateProgramWithSource、clBuildProgram、clCreateKernel);
- 创建指令队列(clCreateCommandQueue);
- 建立和管理涉及的内存对象(clCreateBuffer、clEnqueueWriteBuffer);
- 传递kernel函数的参数(clSetKernelArg);
- 执行kernel函数(clEnqueueNDRangeKernel);
- 主机同步(clWaitForEvents、clReleaseEvent);
- 获取执行结果(clEnqueueReadBuffer);
- 释放平台资源(clReleaseProgram、clReleaseContext、clReleaseCommandQueue、clReleaseDevice、clReleaseKernel)。
下图表示了系统中CPU和GPU的互联关系,包括用CPU为GPU创建程序,传递数据,而GPU则执行程序,运算数据,前文介绍了如何在CPU端构建应用,那要如何编程实现运行在GPU的内核函数呢?
首先,来理解一下OpenCL里的工作组、工作项和处理单元的概念:
- 工作项:一个循环中最里面的一次运算,称为一个工作项。
- 工作组:是由访问相同处理资源的工作项组成,包括工作项访问高速内存(也叫局部内存)的同一块内存和同步。
- 处理单元:支持工作组的处理资源被称为处理单元,各个工作组在单个处理单元上的执行,各个处理单元一次只能够执行一个工作组。
处理单元的数量由硬件决定,OpenCL的内核函数应根据处理单元数量合理设置工作组和工作项,而工作项往往是一个具体的运算,比如矩阵乘法里的乘加运算。
使用openCL实现矩阵乘法运算的逻辑如下:
在OpenCL程序中,我们为每个工作项分配一个要计算的乘法矩阵的元素。将i,j的外层循环替换为函数调用工作项ID(get_global_id),并保证得到的工作项ID在矩阵C的范围内计算每个a[i][r]和每个b[r][j]的乘积,再对所有的乘积项求和,代码如下,其中多了参数有效性的判断。
int i = get_global_id(0);
int j = get_global_id(1);
int k;
float tmp;
if ((i < iHeightA) && (j < iWidthB)) // && (iWidthA == iHeightB)
{
tmp = 0.0;
for (k = 0; k < iWidthA; k++)
tmp += pInMatA[i * iWidthA + k] * pInMatB[k * iWidthB + j];
pOutMat[i * iHeightA + j] = tmp;
}
下面用C语言来实现这个算法,对比可知,最外层循环的i和j即上文中传入的工作项ID,它们表示这个在C语言里顺序执行的运算被内核函数并行处理了,比如在这个示例里两个相乘矩阵的第一个矩阵的Height和第二个矩阵的Width(两个矩阵相乘要求第一个矩阵的Width等于第二个矩阵的Height)共同构成了并行的。
#define M 800
#define P 500
#define N 800
void RunAsCpu(const float* pInMatA,const float* pInMatB,float* pOutMat){
for (int i = 0; i < M; i++){
for (int j = 0; j < N; j++){
pOutMat[i*N + j] = 0.0;
for (int k = 0; k < P; k++){
pOutMat[i*N + j] += pInMatA[i*P + k] * pInMatB[k*N + j];
}
}
}
}
矩阵乘法的核心是一个乘加计算(pOutMat[i*N + j] += pInMatA[i*P + k] * pInMatB[k*N + j]
)。所以优化矩阵运算需要尽量减少数据移动,以保证这个计算可以接近峰值性能运行,参考文章中的过程和方法。
C语言实现
下载Intel SDK for OpenCL applications,选择平台后注册帐号后即可下载。
从github的intel_OpenCL_MulMatrix下载工程,用VS打开CapsBasic中的sln工程即编译运行,下文内容中的代码引用自该工程。
以下代码用于找到用于运算的指定的GPU平台,本文所用的是Intel平台GPU。
cl_uint num_of_platforms = 0;
error = clGetPlatformIDs(0, 0, &num_of_platforms);
cl_platform_id* platforms = new cl_platform_id[num_of_platforms];
error = clGetPlatformIDs(num_of_platforms, platforms, 0);
cl_uint selected_platform_index = num_of_platforms;
for (cl_uint i = 0; i < num_of_platforms; ++i){
size_t platform_name_length = 0;
size_t maxComputeUnits = 0;
error = clGetDeviceInfo(deviceIds[0],CL_DEVICE_MAX_COMPUTE_UNITS,sizeof(cl_uint),&maxComputeUnits,NULL);
error = clGetPlatformInfo(platforms[i],CL_PLATFORM_NAME,0,0,&platform_name_length);
char* platform_name = new char[platform_name_length];
error = clGetPlatformInfo(platforms[i],CL_PLATFORM_NAME,platform_name_length,platform_name,0);
if (strstr(platform_name, "Intel") &&selected_platform_index == num_of_platforms)
selected_platform_index = i;
delete[] platform_name;
}
if (selected_platform_index == num_of_platforms) return 1;
然后就可以构建用于运算的queue,由平台platform获得context,再由clCreateCommandQueue将其构成queue:
cl_platform_id platform = platforms[selected_platform_index];
error = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
context = clCreateContext(0, 1, &device, NULL, NULL, &error);
queue = clCreateCommandQueue(context, device, CL_QUEUE_PROFILING_ENABLE, &error);
生成用于测试的数据,实际工作中,可以为图像数据:
float* A_h = new float[M*P];
float* B_h = new float[P*N];
float* C_h = new float[M*N];
srand(100);
for (int i = 0; i < M*P; i++)A_h[i] = (float)(rand() % 50);
for (int i = 0; i < P*N; i++)B_h[i] = (float)(rand() % 50);
将生成的数据传入显存中,使用了clCreateBuffer函数,可以理解为将大小为sizeof(float)MP的A_h数组传入context的A_d显存中:
cl_mem A_d = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float)*M*P, A_h, &error);
cl_mem B_d = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float)*P*N, B_h, &error);
cl_mem C_d = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(float)*M*N, NULL, &error);
载入openCL的源程序,通过fopen打开文件,用fread将文件中的内容读到source中:
FILE* fp = fopen("OpenCLMulMatrix.cl", "rb");
fseek(fp, 0, SEEK_END);
size_t src_size = ftell(fp);
fseek(fp, 0, SEEK_SET);
const char* source = new char[src_size];
fread((void*)source, 1, src_size, fp);
fclose(fp);
编译openCL的源程序,使用 clCreateProgramWithSource函数将程序所在的source传给context,再通过clBuildProgram函数进行编译,最终从编译的GPU可执行文件中找到指定的功能函数,如“RunAsGpu_1”,指定给cl_kernel所定义的run_as_gpu_1进行数据运算:
cl_program program = clCreateProgramWithSource(context, 1, &source, &src_size, &error);
delete[] source;
error = clBuildProgram(program, 1, &device, NULL, NULL, NULL);
char* build_log;
size_t log_size;
clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
build_log = new char[log_size + 1];
clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
build_log[log_size] = '\0';
printf("build_log:%s", build_log);
delete[] build_log;
cl_kernel run_as_gpu_1 = clCreateKernel(program, "RunAsGpu_1", &error);
传入kernel程序运行的参数,通过clSetKernelArg,按*.cl文件中RunAsGpu_1函数所对应的参数顺序进行设定:
cl_int M_d = M;
cl_int P_d = P;
cl_int N_d = N;
error = clSetKernelArg(run_as_gpu_1, 0, sizeof(cl_mem), &A_d);
error |= clSetKernelArg(run_as_gpu_1, 1, sizeof(cl_mem), &B_d);
error |= clSetKernelArg(run_as_gpu_1, 2, sizeof(int), &M_d);
error |= clSetKernelArg(run_as_gpu_1, 3, sizeof(int), &N_d);
error |= clSetKernelArg(run_as_gpu_1, 4, sizeof(int), &P_d);
error |= clSetKernelArg(run_as_gpu_1, 5, sizeof(cl_mem), &C_d);
设定工作组和工作项,clEnqueueNDRangeKernel将在设备上对数据进行划分,运行kernel程序,对数据进行运算:
size_t globalws_1[2] = { M,N };
cl_event ev;
error = clEnqueueNDRangeKernel(queue, run_as_gpu_1, 2, NULL, globalws_1, NULL, 0, NULL, &ev);
clFinish(queue);
clEnqueueNDRangeKernel(cl_command_queue queue, cl_kernel kernel,
cl_uint work_dims, const size_t *global_work_offset,
const size_t *global_work_size, const size_t *local_work_size,
cl_uint num_events, const cl_event *wait_list,cl_event *event)
-
work_dims
:数据的维数,其限定范围一般是1-3; -
global_work_offset
:各个维度上全局ID偏移量,即数据的起始点,比如你希望二维数据的起始点是(1,3),就可以通过这个参数设置; -
global_work_size
:工作项的总体数量。 -
local_work_size
:一个工作组中工作项的数量,如果参数local_work_size的取值被设置成NULL,opencl将分析决定如何在设备上的处理单元间分配工作项。
比如,对于一个12x12的矩阵,将其划分为9个工作组,每个工作组是一个4x4的矩阵,即16个工作项,相应的设置:
cl_uint work_dims =2; //2维数据
size_t global_work_offset[2] = {0,0}; //从(0,0)开始
size_t global_work_size[2] = {12,12};
size_t local_work_size[2] = {4,4};
clEnqueueNDRangeKernel(command_q,kernel,work_dim,global_work_offset,global_work_size,local_work_size,0,NULL,NULL);
kernel运行的起始时间由kernel运行时自动生成,通过clGetEventProfilingInfo读取:
cl_ulong startTime, endTime;
clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_START,sizeof(cl_ulong), &startTime, NULL);
clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_END,sizeof(cl_ulong), &endTime, NULL);
cl_ulong kernelExecTimeNs = endTime - startTime;
从GPU内存中取得运算结果,从queue中将显存gpu_C_1中的内容读到内存C_d中:
float* gpu_C_1 = new float[M*N];
clEnqueueReadBuffer(queue, C_d, CL_TRUE, 0, M*N*sizeof(float), gpu_C_1, 0, NULL, NULL);
释放资源:
delete[] A_h;
delete[] B_h;
delete[] C_h;
delete[] gpu_C_1;
delete[] platforms;
clReleaseKernel(run_as_gpu_1);
clReleaseCommandQueue(queue);
clReleaseContext(context);
clReleaseMemObject(A_d);
clReleaseMemObject(B_d);
clReleaseMemObject(C_d);
运行时可以看到资源管理器中GPU的资源占用有明显上升。
python实现
使用opencl4py
python中的操作与C可以完全对应。需要安装用于PC的opencl4py模块,在github上有这个模块的源码:
pip install opencl4py
找到用于运算的指定的GPU平台:
os.environ["PYOPENCL_CTX"] = "0:0"
platforms = cl.Platforms()
print("OpenCL devices:\n%s"%platforms.dump_devices())
然后就可以构建用于运算的queue,由平台platform获得context,再由clCreateCommandQueue将其构成queue:
ctx = platforms.create_some_context()
prg = ctx.create_program(testopencl.readoclfile("test.cl"))
print(prg.kernel_names)
krn = prg.get_kernel("MatrixMul")
print(krn.attributes)
queue = ctx.create_queue(ctx.devices[0])
生成用于测试的数据,实际工作中,可以为图像数据:
iHeightA = np.array([800], dtype=np.int32)
iWidthA = np.array([500], dtype=np.int32)
pInMatA = np.arange(iHeightA[0] * iWidthA[0], dtype=np.float32)
iHeightB = np.array([500], dtype=np.int32)
iWidthB = np.array([800], dtype=np.int32)
pInMatB = np.arange(iHeightB[0] * iWidthB[0], dtype=np.float32)
pOutMat = np.empty(iHeightA[0] * iWidthB[0], dtype=np.float32)
将生成的数据传入显存中,使用了clCreateBuffer函数,可以理解为将大小为sizeof(float)*M*P的A_h数组传入context的A_d显存中:
pInMatA_buf = ctx.create_buffer(cl.CL_MEM_READ_ONLY | cl.CL_MEM_COPY_HOST_PTR, pInMatA)
pInMatB_buf = ctx.create_buffer(cl.CL_MEM_READ_ONLY | cl.CL_MEM_COPY_HOST_PTR, pInMatB)
pOutMat_buf = ctx.create_buffer(cl.CL_MEM_READ_WRITE | cl.CL_MEM_ALLOC_HOST_PTR, size=pOutMat.nbytes)
创建一个readoclfile函数,载入openCL的源程序:
def readoclfile(filename):
file_object = open(filename, 'r')
oclfiledata = ""
try:
file_context = file_object.read()
oclfiledata = file_context
finally:
file_object.close()
return oclfiledata
通过调用创建运行于内核的程序:
prg = ctx.create_program(testopencl.readoclfile("test.cl"))
print(prg.kernel_names)
指定用于内核程序中的具体功能函数:
krn = prg.get_kernel("MatrixMul")
print(krn.attributes)
传入kernel程序运行的参数,通过 krn.set_args按*.cl文件中功能函数所对应的参数顺序进行设定:
krn.set_args(iHeightA[0:1], iWidthA[0:1], pInMatA_buf, iHeightB[0:1], iWidthB[0:1], pInMatB_buf, pOutMat_buf)
也可以像C程序一样逐个设定:
krn.set_arg(0, iHeightA[0:1])
krn.set_arg(1, iWidthA[0:1])
krn.set_arg(2, pInMatA_buf)
krn.set_arg(3, iHeightB[0:1])
krn.set_arg(4, iWidthB[0:1])
krn.set_arg(5, pInMatB_buf)
krn.set_arg(5, pOutMat_buf)
运行kernel程序:
start = time.time()
ev = queue.execute_kernel(krn, global_size, local_size, need_event=True)
t1 = time.time() - start
从GPU内存中取得运算结果:
queue.read_buffer(pOutMat_buf, pOutMat)
data1 = np.reshape(pOutMat,(iHeightA[0] , iWidthB[0]))
print(data1[0][1:5])
打印由numpy运算的矩阵乘法数据,进行查看和比较:
data2 = np.dot(np.reshape(pInMatA, (iHeightA[0],iWidthA[0])),np.reshape(pInMatB, (iHeightB[0],iWidthB[0])))
print(data2[0][1:5])
释放资源:
del queue
del ctx
del krn
del prg
gc.collect()
以下是程序的终端输出:
OpenCL devices:
Platform 0: NVIDIA CUDA
Device 0: GeForce GTX 970 (4096 Mb, 4096 align, OpenCL C 1.2)
['test', 'matmul', 'testadd', 'MatrixMul']
[3.3233502e+10 3.3233648e+10 3.3233791e+10 3.3233891e+10]
[3.3233508e+10 3.3233648e+10 3.3233781e+10 3.3233904e+10]
0.015639781951904297 0.0
[6.6467004e+10 6.6467295e+10 6.6467582e+10 6.6467783e+10]
[6.6467017e+10 6.6467295e+10 6.6467561e+10 6.6467807e+10]
0.0 0.0
输出的信息包括主机所用到的OPENCL平台和设备型号,比如上面为英伟达的GTX970,所编译的用于GPU的功能名称包括“test”,“matmul”,“testadd”和“matrixmul”。
运行时可以看到资源管理器中GPU的资源占用有明显上升。
使用pyopencl
pyopencl可以通过pip来安装:
如果希望在小机端使用,可以通过源码编译安装:
pi@NanoPi-NEO4:~$ git clone https://github.com/inducer/pyopencl.git
pi@NanoPi-NEO4:~/pyopencl$ git submodule update --init
编译安装过程中可能会报错,可以根据提示安装缺失的软件包:
pi@NanoPi-NEO4:~/pyopencl$ pip3 install pybind11
pi@NanoPi-NEO4:~/pyopencl$ pip3 install mako
然后就可以顺利编译通过了:
pi@NanoPi-NEO4:~/pyopencl$ python3 setup.py bdist_wheel
pi@NanoPi-NEO4:~/pyopencl$ cd dist/
pi@NanoPi-NEO4:~/pyopencl/dist$ ls
pyopencl-2018.2.3-cp36-cp36m-linux_aarch64.whl
将编译得到的文件复制到服务器,以便以后可以直接安装,然后就可以通过pip可安装pyopencl了:
pi@NanoPi-NEO4:~/pyopencl/dist$ cp pyopencl-2018.2.3-cp36-cp36m-linux_aarch64.whl ../../smb/
pi@NanoPi-NEO4:~/smb$ sudo pip3 install pyopencl-2018.2.3-cp36-cp36m-linux_aarch64.whl
接下来就是运行和验证,新建一个“test_ocl.py”文件,将以下内容复制到该文件中:
import pyopencl as cl
print("CL_VERSION:",cl.VERSION)
print("CL_HEADER_VERSION:",cl.get_cl_header_version())
print()
platforms = cl.get_platforms()
print("Platform num:",len(platforms))
for plat in platforms:
print("--Platform Name:",plat.get_info(cl.platform_info.NAME))
# print("--Platform Extensions:",plat.get_info(cl.platform_info.EXTENSIONS))
print("--Platform Profile:",plat.get_info(cl.platform_info.PROFILE))
print("--Platform Vendor:",plat.get_info(cl.platform_info.VENDOR))
print("--Platform Version:",plat.get_info(cl.platform_info.VERSION))
devices = plat.get_devices(cl.device_type.ALL)
print("--device num:",len(devices))
for device in devices:
print("----Name:",device.get_info(cl.device_info.NAME))
print("----OpenCL_C_Version:",device.get_info(cl.device_info.OPENCL_C_VERSION))
print("----Vendor:",device.get_info(cl.device_info.VENDOR))
print("----Version:",device.get_info(cl.device_info.VERSION))
print("----Driver Version:",device.get_info(cl.device_info.DRIVER_VERSION))
print("----MAX_WORK_GROUP_SIZE:",device.get_info(cl.device_info.MAX_WORK_GROUP_SIZE))
print("----MAX_COMPUTE_UNITS:",device.get_info(cl.device_info.MAX_COMPUTE_UNITS))
print("----MAX_WORK_ITEM_SIZES:",device.get_info(cl.device_info.MAX_WORK_ITEM_SIZES))
print("----LOCAL_MEM_SIZE:",device.get_info(cl.device_info.LOCAL_MEM_SIZE))
ctx = cl.Context(dev_type=cl.device_type.ALL,
properties=[(cl.context_properties.PLATFORM,platforms[0])])
通过命令运行后可以得到GPU的相关信息:
pi@NanoPi-NEO4:~/smb$ python3 test_ocl.py
CL_VERSION: (2018, 2, 3)
CL_HEADER_VERSION: (1, 2)
Platform num: 1
--Platform Name: ARM Platform
--Platform Profile: FULL_PROFILE
--Platform Vendor: ARM
--Platform Version: OpenCL 1.2 v1.r13p0-00rel0-git(a4271c9).31ba04af2d3c01618138bef3aed66c2c
--device num: 1
----Name: Mali-T860
----OpenCL_C_Version: OpenCL C 1.2 v1.r13p0-00rel0-git(a4271c9).31ba04af2d3c01618138bef3aed66c2c
----Vendor: ARM
----Version: OpenCL 1.2 v1.r13p0-00rel0-git(a4271c9).31ba04af2d3c01618138bef3aed66c2c
----Driver Version: 1.2
----MAX_WORK_GROUP_SIZE: 256
----MAX_COMPUTE_UNITS: 4
----MAX_WORK_ITEM_SIZES: [256, 256, 256]
----LOCAL_MEM_SIZE: 32768
我将编译的安装文件放到了github上:
https://github.com/shaniadolphin/images/blob/master/python/pyopencl-2018.2.3-cp36-cp36m-linux_aarch64.whl
参考文档
# | 链接地址 | 文档名称 |
---|---|---|
1 | https://github.com/Samsung/opencl4py |
opencl4py |
2 | https://www.cnblogs.com/Reyzal/p/7389993.html |
Intel核心显卡OpenCL环境搭建 |
3 | https://www.khronos.org/registry/OpenCL/sdk/1.2/docs/man/xhtml/ |
OpenCL参考手册 |
4 | https://blog.csdn.net/cloud_desktop/article/details/19822025 |
优化OpenCL矩阵运算 |
5 | https://blog.csdn.net/LIYUAN123ZHOUHUI/article/details/52850282 |
opencl中工作组和工作项 |