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

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

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

for(i = 1...n)
    X[i] = M[i] . N;

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

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

Источник

2 ответа

Решение

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

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

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

  3. используйте thrust::transform с thrust::multiplies, чтобы умножить два вектора вместе.

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

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

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

#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 опять же, тогда вы можете захотеть собрать его явно. Это пресо поучительно, в частности "Fusion".

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

Источник

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

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.)

Источник
Другие вопросы по тегам