VexCL似乎是一个非常适合GPU编程的库。
不幸的是,它是一个非常年轻的库,信息很少。我一直在寻找如何执行矩阵向量乘法,但我发现的唯一矩阵表示是vex::SpMat,它保存了稀疏矩阵。
如果矩阵是密集的,则稀疏表示对计算效率来说通常不太高效。
我的所有矩阵都是密集的,我想知道如何在VexCL中高效地执行此操作。
我是VexCL库的开发者。
我必须承认,密集的线性代数运算并不是我的优先考虑。我认为在各种由VexCL支持的设备(即OpenCL/CUDA)之间实现性能可移植的方式是非常困难的。这项任务可能最好留给供应商的BLAS实现(但欢迎贡献!)。
您可能还想查看开源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);
}
其次,您可以使用供应商特定的库。例如,如果您使用带有VexCL的CUDA后端,则可以获取指向由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 less effective。它的优点是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倍。