在 OpenCL 中使用本地工作线程进行大型矩阵计算

//C# host code
...
float[] BigMatrix = new float[161*161*161]; //1-D result array
CLCalc.Program.Variable dev_BigMatrix = new CLCalc.Program.Variable(BigMatrix);
CLCalc.Program.Variable dev_OtherArray = new CLCalc.Program.Variable(otherArray);
//...load some other variables here too.
CLCalc.Program.Variable[] args = new CLCalc.Program.Variable[7] {//stuff...}
//Here, I execute the kernel, with a 2-dimensional worker pool:
BigMatrixCalc.Execute(args, new int[2]{N*N*N,192});
dev_BigMatrix.ReadFromDeviceTo(BigMatrix);

__kernel void MyKernel(
__global float * BigMatrix
__global float * otherArray
//various other variables...
)
{
    int N = 161; //Size of matrix edges
    int pixel_id = get_global_id(0); //The location of the pixel in the 1D array
    int array_id = get_global_id(1); //The location within the otherArray

    //Finding the x,y,z values of the pixel_id.
    float3 p;
    p.x = pixel_id % N;    
    p.y = ((pixel_id % (N*N))-p.x)/N;
    p.z = (pixel_id - p.x - p.y*N)/(N*N);
    float result;
    //...
    //Some long calculation for 'result' involving otherArray and p...
    //...
    BigMatrix[pixel_id] += result;
}

我认为等待内存访问会丢失很多性能。我已经回答了一个类似的SO问题。我希望我的帖子对您有所帮助。请提出任何问题。

优化：

1. 我的内核版本的最大提升来自于将 otherArray 读取到本地内存中。
2. 每个工作项在 BigMatrix 中计算 4 个值。这意味着它们可以在同一缓存行上同时写入。并行度损失最小，因为仍有> 1M 个工作项要执行。

。

#define N 161
#define Nsqr N*N
#define Ncub N*N*N
#define otherSize 192
__kernel void MyKernel(__global float * BigMatrix, __global float * otherArray)
{
    //using 1 quarter of the total size of the matrix
    //this work item will be responsible for computing 4 consecutive values in BigMatrix
    //also reduces global size to (N^3)/4  ~= 1043000 for N=161
    int global_id = get_global_id(0) * 4; //The location of the first pixel in the 1D array
    int pixel_id;
    //array_id won't be used anymore. work items will process BigMatrix[pixel_id] entirely
    int local_id = get_local_id(0); //work item id within the group
    int local_size = get_local_size(0); //size of group

    float result[4]; //result cached for 4 global values
    int i, j;
    float3 p;
    //cache the values in otherArray to local memory
    //now each work item in the group will be able to read the values efficently
    //each element in otherArray will be read a total of N^3 times, so this is important
    //opencl specifies at least 16kb of local memory, so up to 4k floats will work fine
    __local float otherValues[otherSize];
    for(i=local_id; i<otherSize; i+= local_size){
        otherValues[i] = otherArray[i];
    }
    mem_fence(CLK_LOCAL_MEM_FENCE);
    //now this work item can compute the complete result for pixel_id 
    for(j=0;j<4;j++){
        result[j] = 0;
        pixel_id = global_id + j;
        //Finding the x,y,z values of the pixel_id.
        //TODO: optimize the calculation of p.y and p.z
        //they will be the same most of the time for a given work item
        p.x = pixel_id % N;    
        p.y = ((pixel_id % Nsqr)-p.x)/N;
        p.z = (pixel_id - p.x - p.y*N)/Nsqr;
        for(i=0;i<otherSize;i++){
            //...
            //Some long calculation for 'result' involving otherValues[i] and p...
            //...
            //result[j] += ...
        }
    }
    //4 consecutive writes to BigMatrix will fall in the same cacheline (faster)
    BigMatrix[global_id] += result[0];
    BigMatrix[global_id + 1] += result[1];
    BigMatrix[global_id + 2] += result[2];
    BigMatrix[global_id + 3] += result[3];
}

笔记：

1. 全局工作规模需要是 4 的倍数。理想情况下，是 4*工作组大小的倍数。这是因为没有错误检查来查看每个pixel_id是否在以下范围内：0..N^3-1。未处理的元素可以在等待内核执行时由 CPU 处理。
2. 工作组规模应该相当大。这将强制更多地使用缓存的值，并且在LDS中缓存数据的好处将会增加。

3. p.x/y/z 的计算需要进一步优化，以避免太多昂贵的除法和模运算。 请参阅下面的代码。

   __kernel void MyKernel(__global float * BigMatrix, __global float * otherArray)   {
   int global_id = get_global_id(0) * 4; //The location of the first pixel in the 1D array
   int pixel_id = global_id;
   int local_id = get_local_id(0); //work item id within the group
   int local_size = get_local_size(0); //size of group
   
   float result[4]; //result cached for 4 global values
   int i, j;
   float3 p;
   //Finding the initial x,y,z values of the pixel_id.
   p.x = pixel_id % N;    
   p.y = ((pixel_id % Nsqr)-p.x)/N;
   p.z = (pixel_id - p.x - p.y*N)/Nsqr;
   //cache the values here. same as above...
   //now this work item can compute the complete result for pixel_id 
   for(j=0;j<4;j++){
       result[j] = 0;
   //increment the x,y,and z values instead of computing them all from scratch
       p.x += 1;
       if(p.x >= N){
           p.x = 0;
           p.y += 1;
           if(p.y >= N){
               p.y = 0;
               p.z += 1;
           }
       }
       for(i=0;i<otherSize;i++){
           //same i loop as above...
       }
   }
   

我给你一些建议：

1. 我认为您的代码具有竞争条件。 最后一行代码具有由多个不同工作项修改的 BigMatrix 的相同元素。
2. 如果您的矩阵确实是 161x161x161，则此处有很多工作项将这些维度用作唯一的维度。 您已经有> 400 万个工作项，这对于您的计算机来说应该是足够的并行性。 你不需要 192 倍。 此外，如果不将单个像素的计算拆分为多个工作项，则无需同步最终添加。
3. 如果你的全局工作大小不是 2 的大幂的好倍数，你可以尝试将其填充出来，使其成为一个。 即使您将 NULL 作为本地工作大小传递，某些 OpenCL 实现也会为不能很好地划分的全局大小选择低效的本地大小。
4. 如果您的算法不需要本地内存或障碍，则几乎可以跳过本地工作组。

希望这有帮助！