如何使用例如 C 语言风格的并行矩阵向量乘法实现现代化`std::transform`？

Question

我正在尝试使用更现代的方法来调整我的编程实践，这些方法利用 STL 容器、算法、执行策略等。

我开发了一个小测试来试验基本矩阵向量乘法的 C 风格实现的现代化：

C[N] = A[N,M] * B[N]

，它与 OpenMP 并行。在这里，遵循高级主程序，我设置问题的参数，用一些随机值初始化向量，然后以两种不同的方法执行矩阵向量乘法。

int main(){

    size_t N(9), M(5);

    std::vector<double> A(N*M), B(M), C(N);

    fill_random(A, 1.0, 1.0); // if min==max the vector is filled with that particular value, otherwise, it's randomly sampled with values in the range [min,max]
    fill_random(B, 1.0, 1.0);

// Run for-loop based routine (parallel OpenMP)
    std::fill(C.begin(), C.end(), 0.0);
    bool err_for = run_for_loop_based_solution(A, B, C, N, M);
    if(!err_for)
        std::cout << "For-loop solution executed successfully!" << std::endl;

    std::for_each(C.begin(), C.end(), [](double x){
        std::cout << x << " ";
    });
    std::cout << std::endl;

// Run transform-based solution 
    std::fill(C.begin(), C.end(), 0.0);
    bool err_task = run_transform_based_based_solution(A, B, C, N, M);
    if(!err_task)
        std::cout << "Transform-based solution executed successfully!" << std::endl;

    std::for_each(C.begin(), C.end(), [](double x){
    std::cout << x << " ";
    });
    std::cout << std::endl;

    return 0;
}

现在，我们来看看这两个实现！

参考实现，即

run_for_loop_based_solution(A, B, C, N, M);

是使用基于 for 循环的 C 风格实现开发的。该计算与 OpenMP 线程并行，每个线程计算输出向量

C[N]

的一个或多个元素的值。简单！

template<typename T>
bool run_for_loop_based_solution(std::vector<T> &A, std::vector<T> &B, std::vector<T> &C, size_t N, size_t M)
{
    if (A.size() != N*M) return 1;
    if (B.size() != M) return 1;
    if (C.size() != N) return 1;

    #pragma omp parallel for
    for(size_t i=0; i<N; ++i)
        for(size_t j=0; j<M; ++j)
            C[i] += A[i*M + j] * B[j];
    
    return 0;
}

现在，我想使用

std::transform

开发一个实现，即

run_transform_based_based_solution(A, B, C, N, M);

，它提供了选择执行策略以启用并行性的选项，例如

std::execution::par

。像这样的东西：

template<typename T>
bool run_transform_based_based_solution(std::vector<T> &A, std::vector<T> &B, std::vector<T> &C, size_t N, size_t M)
{
    if (A.size() != N*M) return 1;
    if (B.size() != M) return 1;
    if (C.size() != N) return 1;

    // std::transform(std::execution::par, C.begin(), C.end(), C.begin(), [&](double x){
    //     return ...  
    // });

    return 0;
}

我在这一点上很挣扎，因为我无法理解如何在

std::transform

体内传播与每个元素的计算相关的索引或范围。是否可以按照这种方法来做到这一点（？），或者除了

std::transform

方法之外还有其他选择吗？

注意： 要运行代码，只需将上述代码片段包含在单个

main.cpp

文件中并使用

-std==c++17

进行编译。另外，在开头包含以下代码片段：

#include <iostream>
#include <algorithm>
#include <vector>
#include <random>
#include <execution>

// fill a std::vector with random values in the range [min, max]. if min==max the vector is filled with that particular value.
template<typename T>
void fill_random(std::vector<T> &v, T min, T max){
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<T> dist(min, max);
    std::transform(v.begin(), v.end(), v.begin(), [&](T x){
        return dist(gen);
    });
}

Answer 1

获取索引可以使用这个answer

的方法

    std::transform(std::execution::par, C.begin(), C.end(), C.begin(),
        [&A, &B, &C, M](T& i) -> T {
            auto index = &i - &C[0]; // Calculate the row index
            return std::inner_product(A.begin() + index * M, A.begin() + (index + 1) * M, B.begin(), i);
    });

演示