提高 CUDA 中 Compact/Scatter 的效率答案

【问题标题】：Improving the Efficiency of Compact/Scatter in CUDA提高 CUDA 中 Compact/Scatter 的效率
【发布时间】：2015-05-09 13:20:29
【问题描述】：

总结：

关于如何进一步改进 CUDA 中的基本分散操作有什么想法吗？特别是如果有人知道它只会用于将更大的阵列压缩成更小的阵列？或者为什么以下向量化内存操作和共享内存的方法不起作用？我觉得我可能缺少一些基本的东西，任何帮助都将不胜感激。

编辑 03/09/15：所以我发现了这个Parallel For All Blog post“使用 Warp-Aggregated Atomics 优化过滤”。为此，我曾假设原子本质上会变慢，但我错了——尤其是因为我认为我不关心在模拟过程中维护数组中的元素顺序。我得再考虑一下，然后实施它看看会发生什么！

编辑 01/04/16：我意识到我从来没有写过我的结果。不幸的是，在 Parallel for All 博客文章中，他们将紧凑的全局原子方法与 Thrust 前缀和紧凑方法进行了比较，这实际上非常慢。 CUB 的 Device::IF 比 Thrust 快得多——我使用 CUB 的 Device::Scan + 自定义代码编写的前缀和版本也是如此。 warp-aggregate 全局原子方法的速度仍然快了大约 5-10%，但远不及我根据博客中的结果所希望的 3-4 倍。我仍然使用前缀和方法，因为虽然不需要维护元素顺序，但我更喜欢前缀和结果的一致性，并且原子的优势不是很大。我仍然尝试各种方法来改进紧凑性，但到目前为止，对于显着增加代码复杂性而言，最多只能进行微小的改进（2%）。

详情：

我正在 CUDA 中编写一个模拟，其中我压缩了我不再对每 40-60 个时间步进行模拟感兴趣的元素。从分析看来，分散操作在压缩时占用的时间最多 - 比过滤器内核或前缀总和更多。现在我使用一个非常基本的分散函数：

    __global__ void scatter_arrays(float * new_freq, const float * const freq, const int * const flag, const int * const scan_Index, const int freq_Index){
            int myID =  blockIdx.x*blockDim.x + threadIdx.x;
            for(int id = myID; id < freq_Index; id+= blockDim.x*gridDim.x){
                 if(flag[id]){
                    new_freq[scan_Index[id]] = freq[id];
                 }
             } 
    }

freq_Index 是旧数组中的元素数。标志数组是过滤器的结果。 Scan_ID 是标志数组上前缀和的结果。

我为改进它所做的尝试是首先将标记的频率读入共享内存，然后从共享内存写入全局内存——这个想法是对全局内存的写入将在扭曲之间更加合并（例如，相反线程 0 写入位置 0 和线程 128 写入位置 1，线程 0 将写入 0，线程 1 将写入 1）。我还尝试将读取和写入矢量化——而不是读取和写入浮点数/整数，我尽可能从全局数组中读取/写入 float4/int4，因此一次四个数字。我认为这可能会通过更少的内存操作传输更多的内存来加速分散。具有矢量化内存加载/存储和共享内存的“厨房水槽”代码如下：

    const int compact_threads = 256;
    __global__ void scatter_arrays2(float * new_freq, const float * const freq, const int * const flag, const int * const scan_Index, const int freq_Index){
        int gID =  blockIdx.x*blockDim.x + threadIdx.x; //global ID
        int tID = threadIdx.x; //thread ID within block
        __shared__ float row[4*compact_threads];
        __shared__ int start_index[1];
        __shared__ int end_index[1];
        float4 myResult;
        int st_index;
        int4 myFlag;
        int4 index;
        for(int id = gID; id < freq_Index/4; id+= blockDim.x*gridDim.x){
            if(tID == 0){
                index = reinterpret_cast<const int4*>(scan_Index)[id];
                myFlag = reinterpret_cast<const int4*>(flag)[id];
                start_index[0] = index.x;
                st_index = index.x;
                myResult = reinterpret_cast<const float4*>(freq)[id];
                if(myFlag.x){ row[0] = myResult.x; }
                if(myFlag.y){ row[index.y-st_index] = myResult.y; }
                if(myFlag.z){ row[index.z-st_index] = myResult.z; }
                if(myFlag.w){ row[index.w-st_index] = myResult.w; }
            }
            __syncthreads();
            if(tID > 0){
                myFlag = reinterpret_cast<const int4*>(flag)[id];
                st_index = start_index[0];
                index = reinterpret_cast<const int4*>(scan_Index)[id];
                myResult = reinterpret_cast<const float4*>(freq)[id];
                if(myFlag.x){ row[index.x-st_index] = myResult.x; }
                if(myFlag.y){ row[index.y-st_index] = myResult.y; }
                if(myFlag.z){ row[index.z-st_index] = myResult.z; }
                if(myFlag.w){ row[index.w-st_index] = myResult.w; }
                if(tID == blockDim.x -1 || gID == mutations_Index/4 - 1){ end_index[0] = index.w + myFlag.w; }
            }
            __syncthreads();
            int count = end_index[0] - st_index;

            int rem = st_index & 0x3; //equivalent to modulo 4
            int offset = 0;
            if(rem){ offset = 4 - rem; }

            if(tID < offset && tID < count){
                new_mutations_freq[population*new_array_Length+st_index+tID] = row[tID];
            }

            int tempID = 4*tID+offset;
            if((tempID+3) < count){
                reinterpret_cast<float4*>(new_freq)[tID] = make_float4(row[tempID],row[tempID+1],row[tempID+2],row[tempID+3]);
            }

            tempID = tID + offset + (count-offset)/4*4;
            if(tempID < count){ new_freq[st_index+tempID] = row[tempID]; }
        }
        int id = gID + freq_Index/4 * 4; 
        if(id < freq_Index){
            if(flag[id]){
                new_freq[scan_Index[id]] = freq[id];
            }
        }
    }

显然它变得有点复杂。 :) 虽然当数组中有数十万个元素时，上述内核看起来很稳定，但我注意到当数组数以千万计时出现竞争情况。我仍在尝试追踪错误。

但无论如何，没有一种方法（共享内存或矢量化）一起或单独提高性能。我对向量化内存操作缺乏好处感到特别惊讶。它对我编写的其他函数有所帮助，但现在我想知道它是否有帮助，因为它在其他函数的计算步骤中增加了指令级并行性，而不是减少了内存操作。

【问题讨论】：

我很困惑。在进行流compaction时，你看的不是聚集而不是分散操作吗？
@njuffa 也许我使用了错误的术语，但我遵循的压缩数组的过程是：过滤器>扫描>分散：过滤器确定要保留数组的哪些元素，前缀扫描确定每个保留元素的新索引，然后将旧数组的元素分散到新数组中。我相信这最后一步也可以在聚集操作中完成，但我认为分散更有效。有没有更好的压缩方法？
@njuffa http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html "Horn 在 2005 年实现流压缩的 GPU 没有 scatter 功能，所以 Horn 用一系列收集步骤来模拟 scatter。压缩 n 个元素需要 log n收集步骤，虽然这些步骤可以在一个片段程序中实现，但这种“收集搜索”操作相当昂贵，并且需要更多的内存操作。在最近的 GPU 中添加本机分散使得流压缩大大提高了效率。”跨度>
其性能取决于数据。例如，如果原始数组中标记的数据元素之间的距离大于块尺寸（并且由于网格跨步循环而忽略网格尺寸间距），则共享内存步骤将没有任何好处。我还认为作为回写式缓存的 L2 缓存将显着减轻这种行为。当这些行被驱逐时，它们可能大部分都合并到主内存，即使没有共享内存模块。你没有说你正在运行什么样的设备，也没有提供有人可以测试的完整代码
@cr_dave：据我所知，术语如下：“收集”是指非连续负载。 “分散”是指不连续的商店。 “流压缩”是将非连续输入流转换为连续输出流的操作。因此我的问题。对于访问之间的平均距离较小的聚集操作，读取纹理路径，例如使用 LDG（参见文档），通常会有所帮助。

标签： performance cuda scatter stream-compaction

【解决方案1】：

我发现poster 中提到的算法（paper 中也讨论了类似的算法）效果很好，特别是对于压缩大型数组。它使用更少的内存来执行此操作，并且比我以前的方法（5-10%）略快。我对海报的算法进行了一些调整：1）消除阶段 1 中的最终扭曲 shuffle 减少，可以在计算元素时简单地求和，2）使函数能够处理的不仅仅是大小为 a 的数组1024 的倍数 + 添加网格跨步循环，以及 3) 允许每个线程在第 3 阶段同时加载它们的寄存器，而不是一次加载一个。我还使用 CUB 而不是 Thrust for Inclusive sum 来加快扫描速度。我可以做更多的调整，但现在这很好。

//kernel phase 1
int myID =  blockIdx.x*blockDim.x + threadIdx.x;
//padded_length is nearest multiple of 1024 > true_length
for(int id = myID; id < (padded_length >> 5); id+= blockDim.x*gridDim.x){
    int lnID = threadIdx.x % warp_size;
    int warpID = id >> 5;

    unsigned int mask;
    unsigned int cnt=0;//;//

    for(int j = 0; j < 32; j++){
        int index = (warpID<<10)+(j<<5)+lnID;
        
        bool pred;
        if(index > true_length) pred = false;
        else pred = predicate(input[index]);
        mask = __ballot(pred); 

        if(lnID == 0) {
            flag[(warpID<<5)+j] = mask;
            cnt += __popc(mask);
        }
    }

    if(lnID == 0) counter[warpID] = cnt; //store sum
}

//kernel phase 2 -> CUB Inclusive sum transforms counter array to scan_Index array

//kernel phase 3
int myID =  blockIdx.x*blockDim.x + threadIdx.x;

for(int id = myID; id < (padded_length >> 5); id+= blockDim.x*gridDim.x){
    int lnID = threadIdx.x % warp_size;
    int warpID = id >> 5;

    unsigned int predmask;
    unsigned int cnt;

    predmask = flag[(warpID<<5)+lnID];
    cnt = __popc(predmask);

    //parallel prefix sum
#pragma unroll
    for(int offset = 1; offset < 32; offset<<=1){
        unsigned int n = __shfl_up(cnt, offset);
        if(lnID >= offset) cnt += n;
    }

    unsigned int global_index = 0;
    if(warpID > 0) global_index = scan_Index[warpID - 1];

    for(int i = 0; i < 32; i++){
        unsigned int mask = __shfl(predmask, i); //broadcast from thread i
        unsigned int sub_group_index = 0;
        if(i > 0) sub_group_index = __shfl(cnt, i-1);
        if(mask & (1 << lnID)){
            compacted_array[global_index + sub_group_index + __popc(mask & ((1 << lnID) - 1))] = input[(warpID<<10)+(i<<5)+lnID]; 
        }
    }
}

}

编辑：海报作者的一个子集有一个更新的article，他们在其中检查了比上面写的更快的紧凑变体。但是，他们的新版本不保留订单，所以对我自己没有用，我还没有实现它来测试它。也就是说，如果您的项目不依赖对象顺序，那么他们更新的紧凑版本可能会加快您的算法。

【讨论】：