c++ - CUDA 'dot product' 内核能否加速批量 RMS 计算？

Question

这里是 CUDA 的新手，但希望它的使用可以减少单线程 CPU 代码计算均方根所需的时间(具体来说，每个 u 连续不相交的 RMS长度为 v 的子数组包含在长度为 u*v 的数组 A 中，在 1 和 -1 之间 float 。

我曾希望使用 this示例来完成绝大多数工作*，但我发现即使是将 A 的子数组与 0 分开的第一步也是为了内存合并目的(这我正在通过单线程 CPU 代码进行)比使用“ballpark”值 u=200、v=5000 的整个基于 CPU 的 RMS 计算花费的时间更长!

我隐约意识到可能有一种方法可以同时应用填充 A 被复制到设备内存，因为我认为也许可以使用第二个内核来执行填充，但我不确定探索这些方法是否值得。我也知道 Thrust 和 cuBLAS，但在我简单的头脑中，链接的样本似乎更有可能提供所需的加速(假设必须小心准备内核输入)。

上面的“大概”值(类似于示例中的值)是否太小以至于无法让 GPU 发挥作用？唉，这两个数字都不太可能提高到 10 的下一个幂。非常感谢那些比我更熟悉 GPU 计算的人提供的任何意见。 FWIW，这是我试图改进的单线程 CPU 代码(rms 函数)和它的一些上下文:

const size_t num     = 5000; // might be able to increase by a factor of <2
const size_t numSegs =  200; // might be able to increase by a factor of <5

float rms(const float a[], const size_t origin, const size_t length)
{
    float sumOfSquares = 0.0f;

    for (size_t i = 0; i < length; ++i)
        sumOfSquares += a[origin + i] * a[origin + i];

    return std::sqrt(sumOfSquares / float(length));
}

int main()
{
    ...

    float* array = (float*)malloc(num * numSegs * sizeof(float));
    float* RMSes = (float*)malloc(numSegs * sizeof(float));

    // array fill omitted; A[i] lies between -1 and 1 inclusive

    for (size_t segNum = 0; segNum < numSegs; ++segNum)
    {
        RMSes[segNum] = rms(array, segNum * num, num);
    }
    ...
}

*假设 RMS(A)=sqrt(B/C)，其中 B 是“A 点 A”，C 是 A 的长度

编辑:基于 CUDA 的方法确实有效，但目前比主机代码慢很多

编辑(再次):

objective-c PU:x86

目标 GPU:RTX3070(计算能力为 8.6 的“Ampere”)

执行时间:

在 main() 中循环(上图):2227-2259us(6 次运行)

使用 dot_product_avx2 在 main() 中循环 @Michael Roy 的回答:445-491us(3 次运行)

填充 array[]:2279-3748us(6 次运行)

this 的第 119-134 行(不是通过 &hTimer 的行)+ 使用 h_C_GPU 计算 RMS:2050-2306us(3 次运行)

this 的第 119-134 行(不是通过 &hTimer 的行)+ 使用 h_C_GPU without padding array[] 计算 RMS:2125-2382us (3运行)

看到填充对基于 GPU 的方法的影响如此之小，我感到非常惊讶(尽管示例建议“ElementN 强烈推荐为 warp 大小的倍数以满足内存合并的对齐约束”)

并行代码只是链接示例的最小修改版本。 VECTOR_N 对应于 u，ELEMENT_N 对应于大于 v 且可被 warp 大小整除的最小数字， cudaMemcpyHostToDevice只出现一次，A[i]可能为负数，StopWatchInterface相关代码省略，h_C_GPU中的每个点积> 用于计算 RMS 值。

Answer 1

这是一个使用 Thrust/CUB 的实现。它应该作为您提出的任何 CUDA 解决方案的性能引用。由于 CUB 算法不知道问题的规律性，因此理论上完全可以编写比这个更快的 CUDA 实现。但在实践中，通过此引用获得显着加速可能非常重要。

这个问题完全适合 GPU 计算，但您的问题规模可能仍然太小，GPU 无法发挥作用。

我决定在这里使用 cub::DeviceSegmentedReduce 而不是 thrust::reduce_by_key (在后端使用 CUB)以获得更好的性能，因为它更容易获得围绕测量分配临时存储的开销。

#include <cmath>

#include <iostream>

#include <cub/cub.cuh>

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/transform_output_iterator.h>
#include <thrust/random.h>
#include <thrust/random/uniform_real_distribution.h>

constexpr size_t num     = 5000; // might be able to increase by a factor of <2
constexpr size_t numSegs =  200; // might be able to increase by a factor of <5

template <typename T>
class ScaleBy
{
    T factor_;

    public:
    ScaleBy(T factor) noexcept : factor_{factor} {}

    __host__ __device__
    T operator()(T val) const noexcept { return factor_ * val; }
};

template <typename T>
struct Square
{
    __host__ __device__
    T operator()(T val) const noexcept { return val * val; }
};

template <typename T>
class RootMean
{
    T norm_;

    public:
    RootMean(T norm) noexcept : norm_{norm} {}

    __host__ __device__
    T operator()(T sum) const noexcept { return sqrt(sum / norm_); }
};

void segmented_rms_device(thrust::device_vector<float> const &d_array,
                          thrust::device_vector<float> &d_RMSes,
                          uint8_t *d_temp_storage,
                          size_t &temp_storage_bytes)
{
    auto seg_size = d_array.size() / d_RMSes.size();
    auto origin_iter = thrust::make_transform_iterator(
        thrust::make_counting_iterator(0ull),
        ScaleBy<size_t>{seg_size});
    auto input_iter = thrust::make_transform_iterator(
        d_array.cbegin(),
        Square<float>{});
    auto output_iter = thrust::make_transform_output_iterator(
        d_RMSes.begin(),
        RootMean<float>{static_cast<float>(seg_size)});

    cub::DeviceSegmentedReduce::Sum(d_temp_storage,
                                    temp_storage_bytes,
                                    input_iter,
                                    output_iter, numSegs,
                                    origin_iter, origin_iter + 1);
}

float rms(thrust::host_vector<float> const &a, const size_t origin, const size_t length)
{
    float sumOfSquares = 0.0f;

    for (size_t i = 0; i < length; ++i)
        sumOfSquares += a[origin + i] * a[origin + i];

    return std::sqrt(sumOfSquares / float(length));
}

void segmented_rms_host(thrust::host_vector<float> const &array,
                        thrust::host_vector<float> &RMSes)
{
    for (size_t segNum = 0; segNum < numSegs; ++segNum)
    {
        RMSes[segNum] = rms(array, segNum * num, num);
    }
}

int main()
{
    thrust::default_random_engine rng(123456789);
    thrust::uniform_real_distribution<float> dist(-1.0f, 1.0f); // excludes 1.0f, but ok for testing

    thrust::host_vector<float> array(num * numSegs);
    thrust::host_vector<float> RMSes_ref(numSegs);

    for (size_t i = 0ull; i < array.size(); ++i)
    {
        array[i] = dist(rng);
    }

    segmented_rms_host(array, RMSes_ref);

    thrust::device_vector<float> d_array(array);
    thrust::device_vector<float> d_RMSes(numSegs);

    // Determine temporary device storage requirements
    size_t temp_storage_bytes = 0;
    segmented_rms_device(d_array, d_RMSes, nullptr, temp_storage_bytes);
    // Allocate temporary storage
    thrust::device_vector<uint8_t> d_temp_storage(temp_storage_bytes);

    segmented_rms_device(d_array, d_RMSes,
                         thrust::raw_pointer_cast(d_temp_storage.data()), temp_storage_bytes);

    thrust::host_vector<float> RMSes(d_RMSes);
    for (size_t i = 0ull; i < numSegs; ++i)
    {
        if (std::abs(RMSes_ref[i] - RMSes[i]) / RMSes_ref[i] > 1.0e-4f)
        {
            std::cout << "Big deviation detected at i = " << i
                      << ": RMS_ref = " << RMSes_ref[i]
                      << " while RMS = " << RMSes[i] << '\n';
        }
    }

    return 0;
}