PyOpenCL示例－02.向量与矩阵运算

​“Python是目前较流行的一种科学计算语言。语法简洁,上手快，易于维护的优点。但其运算速度是真心的慢。那我们能否利用Python的简洁+OpenCL的运算能力呢？答案是可以的，那就是PyOpenCL。”

---

>>向量加法

这里我们利用OpenCL的并行特性做一个向量加法的小程序。OpenCL程序由主程序和内核程序组成。下面我们写个计算[1,2,3,4]+[5,6,7,8]的程序。

main.py源代码

import pyopencl as cl
import numpy as np
 
if __name__ == '__main__':
     
    #step 1:选择并创建上下文
    context = cl.create_some_context() 
     
    #step 3:创建命令队列
    queue = cl.CommandQueue(context)
     
    #定义两个数组，每个数组包含4个元素
    buf_A = np.array([1.0,2.0,3.0,4.0]).astype(np.float32)
    buf_B = np.array([5.0,6.0,7.0,8.0]).astype(np.float32)
     
    #获得OpenCL内存标志集
    mf = cl.mem_flags
    a_buf = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=buf_A)
    b_buf = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=buf_B)
    c_buf = cl.Buffer(context, mf.WRITE_ONLY , buf_A.nbytes )
     
    #加载并创建CL程序
    f_cl = open("add.cl","r")
    f_buf = f_cl.read()
    prg = cl.Program(context,f_buf).build()
     
    #运行
    prg.sum(queue,buf_A.shape,None,a_buf,b_buf,c_buf)
     
    #克隆一块相同的缓冲区
    c_ret = np.empty_like(buf_A)
    #复制设备内存中的值
    cl.enqueue_copy(queue,c_ret,c_buf)
    #输出结果
    print(buf_A,"+",buf_B,"=",c_ret)

add.cl内核代码

__kernel void sum(__global const float *a,__global const float *b,__global float *c)
{  
    int cid = get_global_id(0);  
    c[cid] = a[cid] + b[cid];
}

#运行结果

[6,8,10,12]

---

>>矩阵乘法

矩阵相乘的例子，他是把矩阵分成一个二维的4X4数组。计算时，16个数位同时计算，每个数位都进行四次乘法和四次加法的运算。

其实我们不难发现矩阵乘法其实就是行列式的点乘计算。在OpenCL里内置有dot运算符，通过这个运算符我们可以减少每一数位的运算次数，把四次乘法和四次加法合并成一次点乘计算。当然，我们必须把之前的int类型先转变成float数据类型,另外把数组B先转置。下面我们来看代码。

main.py源代码

import pyopencl as cl
import numpy as np
 
if __name__ == '__main__':
     
    #step 1:选择并创建上下文
    context = cl.create_some_context() 
     
    #step 3:创建命令队列
    queue = cl.CommandQueue(context)
     
    #定义两个数组，每个数组包含4个元素
    buf_A = np.array([1,0,0,0,
                      0,1,0,0,
                      0,0,1,0,
                      0,2,3,1,]).astype(np.float32)
                       
    buf_B = np.array([1,0,4,5,
                      0,1,0,0,
                      0,0,1,0,
                      0,0,0,1,]).astype(np.float32)
                       
    #数组B转置
    buf_B = buf_B.reshape(4,4).T
     
    #二维数组转一维
    buf_B = buf_B.flatten()
       
                       
     
    #获得OpenCL内存标志集
    mf = cl.mem_flags
    a_buf = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=buf_A)
    b_buf = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=buf_B)
    c_buf = cl.Buffer(context, mf.WRITE_ONLY , buf_A.nbytes )
     
    #加载并创建CL程序
    f_cl = open("mul.cl","r")
    f_buf = f_cl.read()
    prg = cl.Program(context,f_buf).build()
     
    #运行
    prg.matrix_mult(queue,buf_A.shape,None,a_buf,b_buf,c_buf)
     
    #克隆一块相同的缓冲区
    c_ret = np.empty_like(buf_A)
    #复制设备内存中的值
    cl.enqueue_copy(queue,c_ret,c_buf)
    #输出结果
    print(c_ret)

mul.cl内核代码

__kernel void matrix_mult(__global float4* A,
                    __global float4* B,
                    __global float* C)
{
    #获取索引号，这里是一维的
    除以4取整得到组号（0-3组）
    int idx = get_global_id(0)/4;
       
    //一共四组数字
    int num_rows = 4;   
    int i=0;
       
    for(i=0;i<num_rows;i++)
    {
        C[idx*num_rows+i] = dot(A[idx],B[i]);
    }   
}

#运行结果：

[1. 0. 4. 5. 0. 1. 0. 0. 0. 0. 1. 0. 0. 2. 3. 1.]