Параллельное умножение многих малых матриц фиксированным вектором

Question

Ситуация такова: у меня есть число (1000) элементов, которые задаются малыми матрицами размеров 4x2, 9x3 ... вы получаете идею. Все матрицы имеют одинаковую размерность.

Я хочу умножить каждую из этих матриц на фиксированный вектор предварительно просчитанных значений. Короче говоря:

Каков наилучший подход для параллельного использования с помощью Thrust? Как я могу разместить свои данные в памяти?

NB: На графических процессорах могут быть специализированные, более подходящие библиотеки. Меня интересует Thrust, потому что он позволяет мне разворачиваться на разные серверы, а не только CUDA.

Answer 1

Одним из возможных подходов:

1. сплющить массивы (матрицы) в один вектор данных. Это предпочтительный шаг для обеспечения общей обработки тяги в любом случае.

2. использовать механизм strided range принять ваш вектор масштабирования и расширения его общей длину вашего уплощенного вектора данных

3. использование thrust::transform с thrust::multiplies перемножить два вектор вместе.

Если вам необходимо получить доступ матрицы позже из вашего сплющенный вектора данных (или вектора результатов), вы можете сделать это с указателем арифметика, или комбинации fancy iterators.

Если вам нужно повторно использовать расширенный вектор масштабирования, вы можете использовать метод, описанный на шаге 2, точно (т. Е. Создать фактический вектор с использованием этого метода, длина = N матриц, повторяется). Если вы делаете это только один раз, вы можете добиться такого же эффекта с помощью счетного итератора, за которым следует итератор преобразования (по модулю длины вашей матрицы в элементах), за которым следует итератор перестановок, для индексации в исходный вектор масштабирования (длина = 1).

Следующий пример реализует выше, без использования strided метод Диапазон итератора:

#include <iostream> 
#include <stdlib.h> 
#include <thrust/device_vector.h> 
#include <thrust/host_vector.h> 
#include <thrust/functional.h> 
#include <thrust/iterator/permutation_iterator.h> 
#include <thrust/iterator/counting_iterator.h> 
#include <thrust/iterator/transform_iterator.h> 
#include <thrust/transform.h> 

#define N_MAT 1000 
#define H_MAT 4 
#define W_MAT 3 
#define RANGE 1024 

struct my_modulo_functor : public thrust::unary_function<int, int> 
{ 
    __host__ __device__ 
    int operator() (int idx) { 
    return idx%(H_MAT*W_MAT);} 
}; 

int main(){ 

    thrust::host_vector<int> data(N_MAT*H_MAT*W_MAT); 
    thrust::host_vector<int> scale(H_MAT*W_MAT); 
    // synthetic; instead flatten/copy matrices into data vector 
    for (int i = 0; i < N_MAT*H_MAT*W_MAT; i++) data[i] = rand()%RANGE; 
    for (int i = 0; i < H_MAT*W_MAT; i++) scale[i] = rand()%RANGE; 

    thrust::device_vector<int> d_data = data; 
    thrust::device_vector<int> d_scale = scale; 
    thrust::device_vector<int> d_result(N_MAT*H_MAT*W_MAT); 

    thrust::transform(d_data.begin(), d_data.end(), thrust::make_permutation_iterator(d_scale.begin(), thrust::make_transform_iterator(thrust::counting_iterator<int>(0), my_modulo_functor())) ,d_result.begin(), thrust::multiplies<int>()); 

    thrust::host_vector<int> result = d_result; 

    for (int i = 0; i < N_MAT*H_MAT*W_MAT; i++) 
    if (result[i] != data[i] * scale[i%(H_MAT*W_MAT)]) {std::cout << "Mismatch at: " << i << " cpu result: " << (data[i] * scale[i%(H_MAT*W_MAT)]) << " gpu result: " << result[i] << std::endl; return 1;} 
    std::cout << "Success!" << std::endl; 
    return 0; 
} 

РЕДАКТИРОВАТЬ: В ответ на вопрос ниже:

Преимущество фантазии итераторы (т.е. transform(numbers, iterator)) является что они часто позволяют исключить дополнительные копии данных/перемещение данных по сравнению с сборкой other number (что требует дополнительных шагов и перемещения данных), а затем передает его на transform(numbers, other numbers). Если вы только собираетесь использовать other numbers один раз, то фантастические итераторы, как правило, будут лучше. Если вы снова собираетесь использовать other numbers, вы можете захотеть его собрать явно. This preso поучителен, в частности «Fusion».

Для одноразового использования other numbers накладные расходы на сборку на лету с использованием причудливых итераторов и функтора, как правило, ниже, чем явно создание нового вектора, а затем передача этого нового вектора в процедуру transform.

Answer 2

При поиске библиотеки программного обеспечения, которая в сжатом виде предназначена для умножения малых матриц, можно взглянуть на https://github.com/hfp/libxsmm. Ниже код запрашивает специализированное ядро ​​матрицы в соответствии с типичными параметрами GEMM (обратите внимание, что применяются некоторые limitations).

double alpha = 1, beta = 1; 
const char transa = 'N', transb = 'N'; 
int flags = LIBXSMM_GEMM_FLAGS(transa, transb); 
int prefetch = LIBXSMM_PREFETCH_AUTO; 
libxsmm_blasint m = 23, n = 23, k = 23; 
libxsmm_dmmfunction xmm = NULL; 

xmm = libxsmm_dmmdispatch(m, n, k, 
    &m/*lda*/, &k/*ldb*/, &m/*ldc*/, 
    &alpha, &beta, &flags, &prefetch); 

Учитывая выше код, можно продолжить и запустить «XMM» для целого ряда (маленьких) матриц без конкретной структуры данных (ниже кода также использует «местоположение предварительной выборки»).

if (0 < n) { /* check that n is at least 1 */ 
    # pragma parallel omp private(i) 
    for (i = 0; i < (n - 1); ++i) { 
    const double *const ai = a + i * asize; 
    const double *const bi = b + i * bsize; 
    double *const ci = c + i * csize; 
    xmm(ai, bi, ci, ai + asize, bi + bsize, ci + csize); 
    } 
    xmm(a + (n - 1) * asize, b + (n - 1) * bsize, c + (n - 1) * csize, 
    /* pseudo prefetch for last element of batch (avoids page fault) */ 
     a + (n - 1) * asize, b + (n - 1) * bsize, c + (n - 1) * csize); 
} 

В дополнение к контролю ручного контура, как показано выше, libxsmm_gemm_batch (или libxsmm_gemm_batch_omp) могут быть использованы (см ReadTheDocs). Последнее полезно, если существуют структуры данных, описывающие серию операндов (матрицы A, B и C).

Существует две причины, почему эта библиотека обеспечивает превосходную производительность: (1) специализированная специализация по коду с использованием технологии генерации кода в памяти и (2) загрузка следующих матричных операндов при вычислении текущего продукта.

(Учитывая один ищет что-то, что хорошо сочетается с C/C++, эта библиотека поддерживает его. Тем не менее, он не стремится к CUDA/Thrust.)