c++ - VexCL 中的密集矩阵 vector 乘法

Question

VexCL 似乎是一个非常有吸引力的 gpu 编程库。

不幸的是，这是一个非常年轻的图书馆，那里的信息很少。我一直在寻找如何执行矩阵 vector 乘法，但我发现的唯一矩阵表示是 vex::SpMat，它包含一个稀疏矩阵。

如果矩阵是稠密的，那么通常情况下，稀疏表示的计算效率较低。

我所有的矩阵都是密集的，我想知道如何在 VexCL 中有效地执行它。

Answer 1

我是 VexCL 的开发者图书馆。

我不得不承认密集的线性代数运算不在我的优先列表中。我认为很难以一种在 VexCL(即 OpenCL/CUDA)支持的各种设备之间实现性能可移植的方式来实现它们。此任务可能最好留给供应商 BLAS 实现(但欢迎打补丁!)。

您可能还想查看开源 ViennaCL库，它确实提供密集矩阵运算并支持 OpenCL、CUDA 和 OpenMP 后端。他们的autotuning framework允许他们获得接近供应商调优库的可移植性能。

话虽如此，对于 VexCL 中的稠密矩阵 - vector 积，您有几个选项(除了提供自定义内核外)。首先，您可以使用基于矩阵 vector 乘积定义的直接实现:

using namespace vex;
Context ctx(Filter::Env && Filter::Count(1));

// The n x m matrix stored row-wise.
vector<double> A(ctx, n * m);
// The LHS and RHS vectors.
vector<double> x(ctx, m);
vector<double> y(ctx, n);

// Make an n x m matrix from vector x by replicating it along the first
// dimension (reshape), multiply it elementwise by A, and reduce the result
// along the second dimension.
// In other words, y_i = sum_j (A_ij * x_j)
y = reduce<SUM>(
        extents[n][m],  // Shape of the expression to reduce,
        A * reshape(
                x,
                extents[n][m], // (We need an n x m matrix...
                extents[1]     // ... but we only have vector of size m).
            ),          // the expression,
        1               // and the dimension to reduce along.
        );

使用 C++14，这可以很容易地隐藏到函数调用中:

template <class M, class V>
auto prod(size_t n, size_t m, M &&A, V &&x) {
    using namespace vex;
    auto NxM = extents[n][m];
    return reduce<SUM>(NxM, A * reshape(x, NxM, extents[1]), 1);
}

其次，您可以只使用供应商特定的库。例如，如果您将 CUDA 后端与 VexCL 一起使用，您可以获得指向 VexCL 分配的内存区域的原始指针并调用 cuBLAS gemv :

double one  = 1;
double zero = 0;
cublasDgemv(
        cublas_handle, CUBPLAS_OP_N, n, m,
        &zero,
        A(0).raw_ptr(), m,
        x(0).raw_ptr(), 1
        &one,
        y(0).raw_ptr(), 1
        );

第一种方法应该不如调用 cuBLAS 有效。它的优点是 reduce() 调用的结果是一个 vector 表达式，原则上您可以将其中的几个组合成一个单一的融合计算内核。例如，您可以计算 Ax + By，或 sin(Ax) + cos(By)，或 (A + B)(x - y)，或单个内核调用中的任何其他 vector 表达式:

z = prod(n, m, A, x) + prod(n, m, B, y);
z = sin(prod(n, m, A, x)) + cos(prod(n, m, B, y));
z = prod(n, m, A + B, x - y);

这可能比几个链接的 cuBLAS 调用更有效。我有 examples因此，VexCL 的性能优于 cuBLAS 1.5 倍。