基于NVIDIA GPU与CUDA8优化std::vector<double>的L1距离并行计算问询

首先，咱们先从CPU层面的优化说起，这些改动不需要依赖GPU，就能快速提升现有代码的性能：

CPU端优化 #

1. 消除分支，改用无分支的绝对值计算 #

你原来代码里的if (val < 0) val = 0 - val;可以直接替换成val = fabs(val);——编译器会把fabs优化成无分支的SIMD指令，比手动判断分支快得多，尤其是循环次数多的时候。优化后的基础CPU代码：

double compute_l1_distance(const std::vector<double> &v1, const std::vector<double> &v2) {
    if (v1.size() != v2.size()) {
        return -1;
    }
    double result = 0;
    const size_t n = v1.size();
    for (size_t i = 0 ; i < n; i++) {
        result += fabs(v1[i] - v2[i]);
    }
    return result;
}

2. 开启编译器优化选项 #

编译时加上-O3 -mavx2（GCC/Clang）或者/O2 /arch:AVX2（MSVC），让编译器自动做循环展开、SIMD向量化等优化，这一步能带来显著的性能提升，几乎不需要改代码。

3. CPU多核并行（OpenMP） #

如果你的CPU是多核的，用OpenMP给循环加并行化指令，让多个核心同时计算：

#include <omp.h>

double compute_l1_distance_omp(const std::vector<double> &v1, const std::vector<double> &v2) {
    if (v1.size() != v2.size()) {
        return -1;
    }
    double result = 0;
    const size_t n = v1.size();
    #pragma omp parallel for reduction(+:result)
    for (size_t i = 0 ; i < n; i++) {
        result += fabs(v1[i] - v2[i]);
    }
    return result;
}

编译时要加上-fopenmp（GCC/Clang）或者开启OpenMP支持（MSVC）。

GPU加速（基于CUDA 8） #

如果你的向量规模足够大（比如百万级以上，能抵消CPU-GPU数据传输的开销），用CUDA加速会带来质的提升。下面给两种实现方式：

1. 用CUDA自带的Thrust库（最简单高效） #

Thrust是CUDA的标准并行算法库，CUDA 8已经内置了它，不用额外安装。它的transform_reduce可以一步完成“计算元素差的绝对值”+“求和”的操作，代码非常简洁：

#include <thrust/device_vector.h>
#include <thrust/transform_reduce.h>
#include <thrust/functional.h>

// 自定义 functor：计算两个元素差的绝对值（兼容CUDA 8对lambda的有限支持）
struct abs_diff {
    __host__ __device__
    double operator()(double a, double b) const {
        return fabs(a - b);
    }
};

double compute_l1_distance_cuda_thrust(const std::vector<double> &v1, const std::vector<double> &v2) {
    if (v1.size() != v2.size()) {
        return -1;
    }
    const size_t n = v1.size();
    
    // 把CPU数据拷贝到GPU设备向量
    thrust::device_vector<double> d_v1(v1.begin(), v1.end());
    thrust::device_vector<double> d_v2(v2.begin(), v2.end());
    
    // transform_reduce：对每个元素应用abs_diff，然后求和
    double result = thrust::transform_reduce(
        thrust::make_zip_iterator(thrust::make_tuple(d_v1.begin(), d_v2.begin())),
        thrust::make_zip_iterator(thrust::make_tuple(d_v1.end(), d_v2.end())),
        abs_diff(),
        0.0,
        thrust::plus<double>()
    );
    
    return result;
}

2. 手动编写CUDA核函数（更灵活） #

如果需要更精细的控制，可以手动写核函数计算每个元素的绝对值差，再做归约求和：

#include <cuda_runtime.h>

// 核函数：计算每个元素的绝对值差
__global__ void compute_abs_diff(const double* v1, const double* v2, double* diff, size_t n) {
    const size_t i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        diff[i] = fabs(v1[i] - v2[i]);
    }
}

// 简单的归约求和核函数（适合演示，实际优先用Thrust/cuBLAS）
__global__ void reduce_sum(const double* input, double* output, size_t n) {
    __shared__ double sdata[256];
    const size_t tid = threadIdx.x;
    const size_t i = blockIdx.x * blockDim.x + tid;
    
    sdata[tid] = (i < n) ? input[i] : 0.0;
    __syncthreads();
    
    for (size_t s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    
    if (tid == 0) {
        output[blockIdx.x] = sdata[0];
    }
}

double compute_l1_distance_cuda(const std::vector<double> &v1, const std::vector<double> &v2) {
    if (v1.size() != v2.size()) {
        return -1;
    }
    const size_t n = v1.size();
    double result = 0.0;
    
    // 分配设备内存
    double *d_v1, *d_v2, *d_diff, *d_temp;
    cudaMalloc(&d_v1, n * sizeof(double));
    cudaMalloc(&d_v2, n * sizeof(double));
    cudaMalloc(&d_diff, n * sizeof(double));
    
    // 拷贝数据到GPU
    cudaMemcpy(d_v1, v1.data(), n * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_v2, v2.data(), n * sizeof(double), cudaMemcpyHostToDevice);
    
    // 启动计算绝对值差的核函数
    const dim3 block_size(256);
    const dim3 grid_size((n + block_size.x - 1) / block_size.x);
    compute_abs_diff<<<grid_size, block_size>>>(d_v1, d_v2, d_diff, n);
    
    // 归约求和
    const size_t temp_size = grid_size.x;
    cudaMalloc(&d_temp, temp_size * sizeof(double));
    reduce_sum<<<temp_size, block_size>>>(d_diff, d_temp, n);
    
    // 最后对中间结果求和（简化处理，大向量可递归归约）
    std::vector<double> temp(temp_size);
    cudaMemcpy(temp.data(), d_temp, temp_size * sizeof(double), cudaMemcpyDeviceToHost);
    for (double val : temp) {
        result += val;
    }
    
    // 释放内存
    cudaFree(d_v1);
    cudaFree(d_v2);
    cudaFree(d_diff);
    cudaFree(d_temp);
    
    return result;
}

推荐的优化库 #

• Thrust：CUDA自带，上手快，支持各种并行算法（归约、排序、变换等），完全适配CUDA 8，是GPU加速的首选。
• cuBLAS：NVIDIA的线性代数库，可以用cublasDaxpy计算两个向量的差，再用cublasDasum计算差向量的L1范数，适合熟悉线性代数库的开发者。
• OpenMP：CPU多核并行的标准库，代码改动极小，适合不想碰GPU的场景。
• Intel MKL：如果是Intel CPU，MKL的向量运算库可以提供极致的CPU性能，支持SIMD和多核并行。

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