VexCL中的密集矩阵向量乘法

我是VexCL库的开发人员.

我必须承认密集线性代数操作不在我的优先级列表中.我认为很难以VexCL支持的各种设备(即OpenCL/CUDA)的性能可移植方式实现这些.这个任务最好留给供应商BLAS实现(但欢迎补丁!).

您可能还想查看opensource ViennaCL库,它提供密集矩阵操作并支持OpenCL,CUDA和OpenMP后端.他们的自动调整框架允许他们获得与供应商调优的库接近的便携性能.

话虽如此,您还有一些选项(除了提供自定义内核)用于VexCL中的密集矩阵 - 矢量产品.首先,您可以使用基于矩阵向量乘积定义的直接实现:

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一起使用,则可以获得指向Ve​​xCL分配的内存区域的原始指针并调用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()调用的结果是一个向量表达式,原则上你可以将其中的几个组合成一个融合的计算内核.例如,您可以在单个内核调用中计算Ax + By,或sin(Ax) + cos(By),(A + B)(x - y)或任何其他向量表达式:

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调用更有效.我有一些例子,其中VexCL因此而优于cuBLAS 1.5倍.