基于CUDA的向量匹配任务并行化可行性与实现方法问询

现有两个整数向量A与B，元素数量从数千到数百万不等：

• A的元素可重复且无序
• B的元素有序且无重复

需要为A中的每个元素v，在B中查找区间[v, v+10]内首个大于等于v的元素，若未找到则结果为0。

当前CPU实现逻辑是：对A排序并记录原索引，单线程遍历排序后的A，同时用计数器遍历B；最后将结果按A的原索引顺序整理。现希望通过CUDA加速该任务，询问：

1. 该任务是否可在CUDA上并行化？
2. 具体的并行化实现方法？

附带的CPU示例代码（Julia版本）：

A = Int32[2, 3, 3, 15, 4, 4, 8, 32, 9, 9, 9, 13, 14]
B = Int32[1, 4, 5, 6, 7, 9, 10, 15, 23, 54]

limit = 10

sortA = sort(A)
inds = sortperm(sortperm(A))

sortR = zeros(Int32, length(A))

b = 1
for i = 1 : length(sortA)
    while B[b] < sortA[i]
        b += 1
        if b == length(B)
            break
        end
    end
    if B[b] <= (sortA[i] + limit)
        sortR[i] = B[b]
    end
end
R = sortR[inds]

# 结果验证
@assert R == [4, 4, 4, 15, 4, 4, 9, 0, 9, 9, 9, 15, 15]

1. 任务是否可并行化？ #

完全可以。原CPU实现的单线程遍历依赖A排序后的顺序复用B的遍历位置，但每个A元素的查询本质是独立的二分查找问题——因为B是有序的，每个v的查询可以独立执行，天然适配CUDA的SIMT并行模型。

2. 具体并行化实现方法 #

核心思路 #

利用B的有序性，为每个A中的元素v独立执行两次二分查找：

• 第一次：在B中找到首个大于等于v的元素位置（记为pos_low）
• 第二次：在B中找到首个大于v+10的元素位置（记为pos_high）
  如果pos_low < pos_high且pos_low未超出B的范围，则结果为B[pos_low]；否则结果为0。

这种方式无需对A排序，直接为每个A元素分配一个CUDA线程并行处理，最后直接得到对应原索引的结果数组，省去CPU版本中排序和索引还原的步骤。

CUDA实现步骤 #

步骤1：数据准备

• 将向量A和B从主机内存（CPU）拷贝到设备内存（GPU）
• 分配设备内存存储结果数组R

步骤2：核函数实现

每个线程负责处理A中的一个元素v = A[threadIdx.x + blockIdx.x * blockDim.x]，执行以下逻辑：

1. 二分查找B中首个≥v的位置pos_low
2. 二分查找B中首个>v+10的位置pos_high
3. 若pos_low < pos_high且pos_low < length(B)，则R[对应索引] = B[pos_low]；否则设为0

步骤3：结果拷贝与验证

将设备内存中的结果数组R拷贝回主机内存，即可得到与CPU版本一致的结果。

示例CUDA C++代码片段 #

#include <cuda_runtime.h>
#include <vector>
#include <algorithm>
#include <cassert>

__device__ int binary_search_low(const int* B, int len, int v) {
    int left = 0, right = len;
    while (left < right) {
        int mid = left + (right - left) / 2;
        if (B[mid] >= v) {
            right = mid;
        } else {
            left = mid + 1;
        }
    }
    return left;
}

__device__ int binary_search_high(const int* B, int len, int v_plus) {
    int left = 0, right = len;
    while (left < right) {
        int mid = left + (right - left) / 2;
        if (B[mid] > v_plus) {
            right = mid;
        } else {
            left = mid + 1;
        }
    }
    return left;
}

__global__ void compute_result(const int* A, const int* B, int* R, int lenA, int lenB, int limit) {
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx >= lenA) return;
    
    int v = A[idx];
    int v_plus = v + limit;
    
    int pos_low = binary_search_low(B, lenB, v);
    int pos_high = binary_search_high(B, lenB, v_plus);
    
    if (pos_low < pos_high && pos_low < lenB) {
        R[idx] = B[pos_low];
    } else {
        R[idx] = 0;
    }
}

int main() {
    std::vector<int> A = {2, 3, 3, 15, 4, 4, 8, 32, 9, 9, 9, 13, 14};
    std::vector<int> B = {1, 4, 5, 6, 7, 9, 10, 15, 23, 54};
    int limit = 10;
    int lenA = A.size();
    int lenB = B.size();
    
    // 设备内存分配
    int *d_A, *d_B, *d_R;
    cudaMalloc(&d_A, lenA * sizeof(int));
    cudaMalloc(&d_B, lenB * sizeof(int));
    cudaMalloc(&d_R, lenA * sizeof(int));
    
    // 数据拷贝到设备
    cudaMemcpy(d_A, A.data(), lenA * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B.data(), lenB * sizeof(int), cudaMemcpyHostToDevice);
    
    // 启动核函数
    int blockSize = 256;
    int gridSize = (lenA + blockSize - 1) / blockSize;
    compute_result<<<gridSize, blockSize>>>(d_A, d_B, d_R, lenA, lenB, limit);
    
    // 结果拷贝回主机
    std::vector<int> R(lenA);
    cudaMemcpy(R.data(), d_R, lenA * sizeof(int), cudaMemcpyDeviceToHost);
    
    // 验证结果
    std::vector<int> expected = {4, 4, 4, 15, 4, 4, 9, 0, 9, 9, 9, 15, 15};
    assert(R == expected);
    
    // 释放内存
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_R);
    
    return 0;
}

性能优化建议 #

• 共享内存复用：如果B的大小适合放入共享内存，可以将B拷贝到每个block的共享内存中，减少全局内存访问延迟，提升查找速度。
• 二分查找优化：使用CUDA Thrust库中的lower_bound或upper_bound函数替代手动实现，利用厂商优化的底层逻辑。
• 批量处理：当A元素数量极大时，可分批次处理，避免占用过多设备内存。

内容的提问来源于stack exchange，提问作者Андрей