如何矢量化 3x3 2D 卷积？

Question

我正在尝试为 1280x720 图像编写优化的 3x3 2D 图像卷积。

为简单起见，通过将输入填充到 1284*724 来接近边缘条件。

这是我的内核代码：

__kernel 
    __attribute__((vec_type_hint(float4)))
    void conv2d3x3(                                                     
       __global const float* restrict input,                         
       __global float* restrict output,
       __constant const float4* restrict hL, 
       /* 3x3 kernel, padded with 3 zeros on the right, used to calculate
        "left" output samples*/
       __constant const float4* restrict hR
       /*same 3x3 kernel, padded with 3 samples on the left*/)           
    {                                        
       int j = get_global_id(0)*2;  //[0,639]  
       int i = get_global_id(1)*2;  //[0,359]                

       /* load a 4x4 block, note stride is 1284 because input is padded by 4*/   
       float4 data0=vload4(0,input+1284*(i+0)+j);
       float4 data1=vload4(0,input+1284*(i+1)+j);
       float4 data2=vload4(0,input+1284*(i+2)+j);
       float4 data3=vload4(0,input+1284*(i+3)+j);

       /* sum(data[0:2,0:2].* h)*/
       float prodTL=dot(data0,hL[0])+dot(data1,hL[1])+dot(data2,hL[2]);
       /* sum(data[0:2,1:3].* h)*/
       float prodTR=dot(data0,hR[0])+dot(data1,hR[1])+dot(data2,hR[2]);
       /* sum(data[1:3,0:2].* h)*/
       float prodBL=dot(data1,hL[0])+dot(data2,hL[1])+dot(data3,hL[2]);
       /* sum(data[1:3,1:3].* h)*/
       float prodBR=dot(data1,hR[0])+dot(data2,hR[1])+dot(data3,hR[2]);

       output[1280*(i+0)+j]=prodTL;
       output[1280*(i+0)+j+1]=prodTR;
       output[1280*(i+1)+j]=prodBL;
       output[1280*(i+1)+j+1]=prodBR;
    } 

这种设计的合理性是，加载一个 4x4 的数据块，进行四个 3x4 卷积并生成 4 个输出样本。

这段代码有几个明显的问题：

1) 矢量负载未与矢量边界对齐。

2) 输出的存储不是向量化的

3) 性能很差：在 Intel XEON 1245v3 上使用 P4600（使用 Beignet OpenCL 实现）为 3ms，在 Freescale IMX6Q 上使用 GC2000（使用 Freescale OpenCL libOpenCL）为 27ms。

问题：

1) 我做错了什么，为什么这么慢？

2) 就原始 FLOPS 的百分比而言，我应该期待什么样的表现？ （p4600 能够在 350MHz 和 1.2GHz 之间运行 20EU * 2PFU/EU * SIMD8 = 320FLOPS/cycle，而 GC2000 至少能够达到 14GFLOPS）

3）一般来说，如何在不产生过多内存流量和缓存冲突的情况下对固定大小的不可分离二维卷积进行向量化？

Answer 1

首先，我未优化的结果：

Amd FX 8150 @3.3 GHz（32 fp 元素 => 1 add + 1 mul = 64 FLOPS/cycle）：

3.71ms 包括单独的 opencl 缓冲区和 C# 数组之间的复制时间。

2.05ms 不包括数组副本。

使用 1-D ndrange 内核执行而不是 2D。 [0,640x360]

__kernel 
    __attribute__((vec_type_hint(float4)))
    void bench(                                                     
       __global const float* restrict input,                         
       __global float* restrict output,
       __constant const float4* restrict hL, 
       
       __constant const float4* restrict hR
      )           
    {
                int gli=get_global_id(0); 
                int j = (gli%640) * 2 ;
                int i = (gli/640) * 2; 

                /* load a 4x4 block*/
                float4 data0 = vload4(0, input + 1280 * (i + 0) + j);
                float4 data1 = vload4(0, input + 1280 * (i + 1) + j);
                float4 data2 = vload4(0, input + 1280 * (i + 2) + j);
                float4 data3 = vload4(0, input + 1280 * (i + 3) + j);

               
                float prodTL = dot(data0, hL[0]) + dot(data1, hL[1]) + dot(data2, hL[2]);
                
                float prodTR = dot(data0, hR[0]) + dot(data1, hR[1]) + dot(data2, hR[2]);
               
                float prodBL = dot(data1, hL[0]) + dot(data2, hL[1]) + dot(data3, hL[2]);
                
                float prodBR = dot(data1, hR[0]) + dot(data2, hR[1]) + dot(data3, hR[2]);

                output[1280 * (i + 0) + j] = prodTL;
                output[1280 * (i + 0) + j + 1] = prodTR;
                output[1280 * (i + 1) + j] = prodBL;
                output[1280 * (i + 1) + j + 1] = prodBR;
            }

主机端（C# 数组）：

        float[] inp = new float[1280*720*2];
        float[] outp = new float[1280*720*2];
        float[] hL = new float[1024];
        float[] hR = new float[1024];

预取到私有寄存器（我只能希望驱动程序使用 cpu 寄存器）：

2ms

优化部分：

        float4 hl2=hL[2];
        float4 hl1=hL[1];
        float4 hl0=hL[0];

        float4 hr2=hR[2];
        float4 hr1=hR[1];
        float4 hr0=hR[0];

        float prodTL = dot(data0, hl0) + dot(data1, hl1) + dot(data2, hl2);
        
        float prodTR = dot(data0, hr0) + dot(data1, hr1) + dot(data2,  hr2);
       
        float prodBL = dot(data1, hl0) + dot(data2, hl1) + dot(data3, hl2);
        
        float prodBR = dot(data1, hr0) + dot(data2, hr1) + dot(data3,  hr2);

现在增加了点积的并行度：

三个点的总和等于一个大点。

  float16 prodhl    =(float16)(hl0,  hl1,  hl2,  (float4)(0.0f,0.0f,0.0f,0.0f));                    
                float16 prodhl    =(float16)(hr0,  hr1,  hr2,  (float4)(0.0f,0.0f,0.0f,0.0f));   
                float16 prodTdata =(float16)(data0,data1,data2,(float4)(0.0f,0.0f,0.0f,0.0f));    

                float16 prodBdata=(float16)(data1,data2,data3,(float4)(0.0f,0.0f,0.0f,0.0f));    

                float prodTL = dot(prodTdata, prodhl);
                float prodTR = dot(prodTdata, prodhr);

                float prodBL = dot(prodBdata, prodhl);
                float prodBR = dot(prodBdata, prodhr);

在没有任何数组副本的情况下执行：

0.5412 毫秒

也许它只是 CPU 的 AVX 功能。如果没有，那么应该有一些指令级并行发生。

这部分（float16的最新float4部分）浪费了1/4的计算能力，所以必须有办法达到0.4 ms。

注意：线程组大小为 256。我没有尝试增加到 1024，因为它不适合所有设备，例如 amd gpu。

您可以尝试任务进程级别的并行性以提高吞吐量并击败单个 opencl 上下文（如果您已经这样做了）。