如何使用thrust::inclusive_scan实现CUDA二阶递归计算并行化

你给出的串行递推逻辑为：

for (int i = 2; i<size; i++)
{
    result[i] = oldArray[i] + k * result[i-2];
}

核心实现思路 #

该二阶递推可以拆分为两个完全独立的一阶递推序列，直接复用一阶偏移（i-1）的并行方案即可：

• 偶数索引序列：仅包含i=0,2,4...的计算项，递推仅依赖前一个偶数索引的结果，转换后为标准一阶递推形式：result[2j] = oldArray[2j] + k * result[2(j-1)]，初始值为result[0]
• 奇数索引序列：仅包含i=1,3,5...的计算项，同理递推仅依赖前一个奇数索引的结果，转换后为：result[2j+1] = oldArray[2j+1] + k * result[2(j-1)+1]，初始值为result[1]
  两个序列之间不存在任何依赖关系，可以并行分别处理。

基于thrust::inclusive_scan的实现方案 #

有两种常见实现方式可选：

方案1：拆分奇偶序列分别扫描 #

这种方式兼容性最好，支持任意长度的输入数组：

1. 分别从oldArray中提取偶数、奇数索引的元素，生成两个独立的输入序列
2. 复用你之前实现的一阶递推扫描算子，分别对两个序列执行inclusive_scan，分别传入result[0]和result[1]作为初始值
3. 将两个扫描后的结果序列按奇偶位置回写回result数组即可

方案2：打包双元素直接扫描 #

这种方式性能更高，无需额外拆分数组，适合大部分场景：
自定义二元结合算子，将相邻两个元素打包为一个状态执行扫描：

// 示例为浮点计算，可按需更换为其他数值类型
struct SecondOrderOp {
    const float k;
    SecondOrderOp(float k_) : k(k_) {}
    __host__ __device__
    float2 operator()(const float2& prev_state, const float2& curr_input) const {
        float res_even = curr_input.x + k * prev_state.x;
        float res_odd = curr_input.y + k * prev_state.y;
        return float2{res_even, res_odd};
    }
};

调用逻辑：

1. 将oldArray从索引2开始的元素按相邻两个一组打包为float2数组，每组第一个元素为偶数索引值，第二个为奇数索引值
2. 初始扫描状态设置为float2{result[0], result[1]}
3. 调用thrust::inclusive_scan传入自定义算子执行扫描
4. 将扫描得到的float2结果数组解包，按位置回写回result数组即可。如果输入数组长度为奇数，单独处理最后一个偶数索引的元素即可。

内容的提问来源于stack exchange，提问作者defladamouse