优化批量矩阵乘法opencl代码

Question

下面是一个opencl内核,它为多个独立矩阵执行阻塞矩阵乘法.selectMatrixA和selectMatrixB以行主要顺序存储多个矩阵(相同大小和方形矩阵).

// Matrix multiplication: C = A * B.


#define BLOCK_SIZE 20
#define MATRIX_SIZE 100 * 100

#define BLOCK_DIMX 5 // Number of blocks in the x dimension

__kernel void
batchedMatrixMul(__global float *selectMatrixC, __global float *selectMatrixA, __global   
float *selectMatrixB, int wA, int wB)
{
    // Block index
    int bx = get_group_id(0);
    int by = get_group_id(1);


    __global float *C = selectMatrixC + (bx/BLOCK_DIMX) * MATRIX_SIZE;
    __global float *A = selectMatrixA + (bx/BLOCK_DIMX) * MATRIX_SIZE;
    __global float *B = selectMatrixB + (bx/BLOCK_DIMX) * MATRIX_SIZE;



    int tx = get_local_id(0);
    int ty = get_local_id(1);

    float Csub = 0;

    // Identify the row and column of the C matrix to work on

    int Row = (by * BLOCK_SIZE)  + ty;
    int Col = ((bx %(BLOCK_DIMX)) * BLOCK_SIZE) + tx;

    // Declaration of the local memory array As used to store the sub-matrix of A
    __local float As[BLOCK_SIZE][BLOCK_SIZE];

    // Declaration of the local memory array Bs used to store the sub-matrix of B
    __local float Bs[BLOCK_SIZE][BLOCK_SIZE];

    // Loop over all the sub-matrices of A and B required to compute the block sub-matrix
    for (int m = 0; m < wA / BLOCK_SIZE; ++m) 
    {

        // Load the matrices from global memory to local memory. Each thread loads one   
        //element of each matrix
        As[ty][tx] = A[Row * wA + m * BLOCK_SIZE + tx];
        Bs[ty][tx] = B[(m * BLOCK_SIZE + ty)*wA + Col];

        // Synchronize to make sure the matrices are loaded
        barrier(CLK_LOCAL_MEM_FENCE);

        // Multiply the two matrices together each thread computes one element of the block 
        //sub-matrix
        for (int k = 0; k < BLOCK_SIZE; ++k)
            Csub += As[ty][k] * Bs[k][tx];

        // Synchronize to make sure that the preceding computation is done before loading 
        //two new sub-matrices of A and B in the next iteration
        barrier(CLK_LOCAL_MEM_FENCE);

    }

    // Write the block sub-matrix to device memory each thread writes one element
    C[Row * wA + Col] = Csub;

}

以下是我启动内核的方法:

localWorkSize[0] = BLOCK_SIZE;
localWorkSize[1] = BLOCK_SIZE;

// for a 100 X 100 matrix, MATRIX_DIMX = MATRIX_DIMY = 100
globalWorkSize[0] = MATRIX_DIMX * NUM_MATRICES;
globalWorkSize[1] = MATRIX_DIMY ;

cl_event             event;
errcode = clEnqueueNDRangeKernel(clCommandQueue, 
          clKernel, 2, NULL, globalWorkSize, 
          localWorkSize, 0, NULL, &event);

以下是在NVIDIA Grid K520上运行时的一些性能数字:

1. matrix size:100 X 100 . Number of matrices = 20000. Time taken for multiplication = 
0.262 seconds. As shown in the code, the block size was set to 20. Block size of 10 was 
slower. This calculates to around 152 GFLOPS

2. matrix size: 10000 X 10000. Number of matrices = 1. Time taken for multiplication = 10.6 
seconds. Here also the block size was 20. Using a block size of 50 is not possible due to   
the size of the local memory.

有人可以帮我理解为什么代码运行缓慢,以及为什么2.比2慢得多.我是OpenCL的新手,我想学习如何根据底层架构细节优化代码.

Answer 1

第一次测试之所以快得多，是因为每次测试所做的工作量有所不同。实际上，是 50 倍。

方阵乘法的 Big-O 为 O(n^3)。参见：为什么方阵乘法的时间复杂度定义为O(n^3)？因此，10k 平方矩阵的乘法实际上比单个 100x100 乘法多花费 100 万倍的工作量。执行 20000 次 100x100 乘法并不能弥补大型矩阵乘法一次所需的大量工作。

矩阵乘法只是许多点积。您的算法仅将点积分成几组以便于处理，并且没有使用任何特殊技巧来减少我下面的计算中的数字。

对于您的小矩阵测试：

Total dot products: 10^4
MADs per dot product: 10^2
Total matrix-multiply operations: 20000 = 2*10^4
Total multiply-adds: 2* 10^(4+2+4) = 2*10^10 = 20,000,000,000

200亿。

大矩阵测试：

Total dot products: 10^8
MADs per dot product: 10^4
Total multiply operations: 1 (or 10^0)
Grand total multiply-adds: 10 ^ (8 + 4 + 0) = 10^12 = 1,000,000,000,000  

10000亿。

从技术上讲，您的 10000x10000 测试运行速度更快 - 只需多出 40 倍的运行时间即可处理多出 50 倍的操作。

在这里阅读有关“特殊技巧”的更多信息： http: //en.wikipedia.org/wiki/Strassen_algorithm。尽管这个算法对于 GPU 计算来说并不被认为是实用的。也存在更复杂的算法，但图形硬件上的强力方法似乎是最常用的。

为什么你的内核通常运行缓慢？您可以使用许多不同的优化来加快速度。下面是一些你可以谷歌搜索并亲自尝试的方法。您可能会遇到一些我在这里没有提到的内容。

• 优化工作组和块大小。请参阅 opencl PREFERRED_WORK_GROUP_SIZE
• 使用 float4 数据类型。opencl 包含一个点积函数，用于计算 floatn 数据类型的点积。
• 在运行内核之前转置矩阵 B。您可以使用另一个内核来进行转置。

Answer 2

在我看来，2. 慢得多的原因是矩阵乘法的访问模式不太适合缓存。如果您需要获取第一行的第一个值和第二行的第一个值，它们将存储在彼此相距很远的内存中。如果矩阵大小增加，它们存储的距离会更远。这将导致大量缓存未命中。

我对矩阵乘法没有任何个人经验，但我只是认为可以将数据存储在Z 顺序曲线中以实现更缓存友好的模式。从维基百科的参考资料来看，Valsalam & al 2002似乎已经完成了类似的事情。

在花费大量时间进行 Z 排序之前，我会尝试另一个快速修复方法，即使用私有变量并消除障碍。即使它需要从全局内存中加载更多负载，编译器也可能对该代码进行更好的优化。