Плотное матрично-векторное умножение в 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 backend с 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() call является векторным выражением, и вы можете в принципе объединить несколько из них в одно объединенное вычислительное ядро. Например, вы можете вычислить 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 раза благодаря этому.