PyOpenCL - 没有看到预期的加速

PyOpenCL - not seeing expected speedup

在试验 PyOpenCL 时,我注意到我的代码 运行 比预期的要慢。事实证明,它在 CPU 上比在 GPU 上 运行 更快(在这两种情况下,在 PyOpenCL 上 运行,仅达到 1 GFLOP)。

为了调试它,我随后尝试了朴素的矩阵乘法作为比较,并且只看到 GPU 与 CPU 相比有 2 倍的加速(~20 GFLOPs 对~10 GFLOPs)。我的系统是i7 8750H + GTX 1070 Max-Q。

有没有人可以分享我可能做错了什么?我知道下面的代码不是最优的,但我预计随着我的 GPU 的浮点能力和内存带宽的增加,会有更大的差异。

import pyopencl as cl
import pyopencl.array as pycl_array
import numpy as np
import numpy.linalg as la
import time

size = 4000

m1 = np.random.normal(size = [size,size]).astype(np.float32)
m2 = np.random.normal(size = [size,size]).astype(np.float32)

ctx = cl.create_some_context(interactive=True)
queue = cl.CommandQueue(ctx)

a = pycl_array.to_device(queue, m1)
b = pycl_array.to_device(queue, m2)
res = pycl_array.empty_like(a)

prg = cl.Program(ctx, """
    __kernel void multiplymatrices(const unsigned int size, __global const float * a, 
    __global const float * b, __global float * res) {

    int i = get_global_id(0); 
    int j = get_global_id(1);

    res[size * i + j] = 0;

    for (int k = 0; k < size; k++)
    {
        res[size * i + j] += a[k + size * j] * b[i + size * k];
    }

    }
    """).build()

t = time.time()
task = prg.multiplymatrices(queue, m1.shape, None, np.int32(size), a.data, b.data, res.data)

task.wait()
tot_time = time.time()-t
print("gflops", 2*size**3/(tot_time*1000**3))

按照使用本地寄存器累加结果的建议,我修改了我的代码如下,在大约 360 GB/s 内存带宽(这是我的 GPU 能够承受的最大带宽)下获得大约 90 gflops的)。改进 gflops 需要更复杂的矩阵乘法算法,该算法多次重用存储在缓存中的相同数据,但这超出了这个问题的范围。

__kernel void multiplymatrices(const unsigned int size, __global const float * a, 
__global const float * b, __global float * res) {

int i = get_global_id(0); 
int j = get_global_id(1);
float temp = 0;

for (int k = 0; k < size; k++)
{
    temp += a[k + size * j] * b[i + size * k];
}
res[size * i + j] = temp;
}

编辑:对于那些寻找快速矩阵乘法示例的人,它展示了将本地内存与工作组一起使用以及 2D 寄存器平铺,我根据教程 here 创建了以下内容。它在我的 GPU 上获得 1.4 TFLOPs。

prg4 = cl.Program(ctx, """
__kernel void multiplymatrices(const unsigned int size, __global const float * A, 
__global const float * B, __global float * res) {

int ig = get_group_id(0);
int jg = get_group_id(1);
int il = get_local_id(0);
int jl = get_local_id(1);

const int memtile = 64;
const int regtile = 4;
volatile int il2;
volatile int jl2;
int iglob = memtile*ig + regtile*il;
int jglob = memtile*jg + regtile*jl;

__local float Asub[64][64];
__local float Bsub[64][64];
float acc[4][4];
float Areg;
float Breg[4];

for (int k = 0; k < regtile; k++) {
    for (int m = 0; m < regtile; m++) {
        acc[k][m] = 0;
    }
}

for (int l = 0; l < size/memtile; l++) {
    for (int k = 0; k < regtile; k++) {
        for (int m = 0; m < regtile; m++) {
            il2 = il*regtile + k;
            jl2 = jl*regtile + m;
            Asub[il2][jl2] = A[size*(iglob + k) + memtile*l + jl2];
            Bsub[il2][jl2] = B[size*(memtile*l + il2) + jglob + m];
        }
    }
    barrier(CLK_LOCAL_MEM_FENCE);

    for (int k = 0; k < regtile; k++) {
        for (int r = 0; r < regtile; r++) {
            Breg[r] = Bsub[il*regtile+k][jl*regtile+r];
        }
        for (int m = 0; m < regtile; m++) {
            Areg = Asub[il*regtile+m][jl*regtile+k];
            for (int r = 0; r < regtile; r++) {
                acc[k][m] += Areg*Breg[r];
            }
        }
    }
}
for (int k = 0; k < regtile; k++) {
    for (int m = 0; m < regtile; m++) {
        res[size*(iglob+k)+jglob+m] = acc[k][m];
    }
}
}
""").build()

t = time.time()
memtile = 64
regtile = 4
wgsize = int(memtile/regtile)
global_size = int(size/regtile)
task = prg4.multiplymatrices(queue, (global_size,global_size), (wgsize,wgsize), np.int32(size), a.data, b.data, res.data)

queue.finish()
tot_time = time.time()-t
print("gflops", 2*size**3/(tot_time*1000**3))
print("GB/s total", 2*4*size**3/(tot_time*1000**3))
print("GB/s global", 2*4*size**3/(memtile*tot_time*1000**3))