Выполнение нескольких матрично-матричных умножений за одну операцию

Question

Я реализую алгоритм, который, по сути, представляет собой серию матрично-матричных умножений:порядка миллиардов.

Я пытался использовать CUBLAS для умножения матриц, но это было медленно, однако я заметил кое-что интересное.

умножение 100x100 на матрицу 100x100 было медленным, но умножение 1.000.000x100 на 100x100 было относительно быстрым, это заставило меня задуматься. Если бы вместо сканирования слева направо было 10.000 сканирований параллельно,Это должно быть довольно быстро, и если бы я умножил свои матрицы, когда закончил с этим, я бы получил тот же результат - просто быстрее.

Res<sub>1</sub> = M<sub>1</sub>.M<sub>2</sub>.M<sub>3</sub>. ... .M<sub>n/1000-1</sub>
Res<sub>1</sub> = M<sub>1+n/1000</sub>.M<sub>2+n/1000</sub>.M<sub>3+n/1000</sub>. ... .M<sub>2(n/1000)-1</sub>
...
Res<sub>1</sub>  = M<sub>1+999*n/1000</sub>.M<sub>2+999*n/1000</sub>.M<sub>3+999*n/1000</sub>. ... .M<sub>1000*(n/1000)-1</sub>
Res = Res<sub>1</sub>*Res<sub>2</sub>* ... *Res<sub>999</sub> 

Ничего не стоит, что M_1 ... M_n находятся в наборе из примерно 100 различных матриц, поэтому использование пространства на самом деле не проблема, все, что мне нужно, это сделать несколько умноженийодна операция.

Теперь вот моя проблема.Я сделал матрично-матричную (sgemm) реализацию, вдохновленную тем, что демонстрирует nvidia в своей документации, но он примерно в 4 раза медленнее, чем cublas.Кто-нибудь знает, как работает CUBLAS?А если код где-то доступен?

Answer 1

Вы смотрели последнюю версию CUBLAS (версия 4.1) ? Он включает в себя новый пакетный режим GEMM, специально предназначенный для больших партий небольших матрично-матричных умножений. Я бы предложил сделать попарно дерево умножения, как предложил Джонатан Дурси в своем ответе, используя пакетный API CUBLAS для его ускорения, вместо того, чтобы писать собственное ядро, как он предлагает.

CUBLAS 4.1 входит в комплект CUDA Toolkit v4.1 .

Answer 2

Проблема в том, что cublas и т. Д. Предназначены для использования всех SM для умножения больших матриц. Это не то, что вы хотите; Вы хотите сделать много маленьких матричных умножений.

Может быть, есть какой-то способ превратить это во что-то, что CUBLAS может сделать для вас хорошо, но я этого не вижу. Мое предложение будет следующим:

Напишите ядро, которое использует один блок потока для умножения двух ваших маленьких матриц, и выведите результат.

Затем запустите журнал ядра ₂ N с тоннами блоков и выполните умножение попарно:

• Шаг 1: умножить M ₁ x M ₂ , M ₃ x M ₄ ... M _{N - 2} x M _N-1 , вывод M ' ₁ , M' ₂ .. M ' _{N / 2}
• Шаг 2: умножить M ' ₁ x M' ₂ , M ' ₃ x M' ₄ ... M ' _{N / 2 - 2} x M _{N / 2-1} , с выводом M' ' ₁ , M' ' ₂ .. M '' _{N / 4} ...

и т.д.

Там будет коэффициент загрузки памяти в 50%, но я думаю, что вы будете лучше использовать ваши ядра таким образом.

Обновление

Хорошо, если вы действительно не хотите делать это поэтапно, вы можете сделать это таким образом, но это потребует большего количества кодирования, и производительность, вероятно, будет хуже по сравнению с тем, что вы могли бы получить с чем-то вроде cuBLAS и асинхронным перечислить. Я предполагаю, что вы используете Fermi, и вы отключили кэш L1, чтобы у вас было 48 Кбайт памяти на SM.

Сохраните 100 матриц в виде блоков 2x2, каждый блок должен быть непрерывным в памяти. Итак, matrix[matrixnum,i,j] начинается с matricies[matrixnum*100*100 + i*100*50 + j*50*50]. Обратите внимание, что каждый блок имеет размер 50 * 50 * 4 байта ~ 10 КБ, поэтому 4 удобно помещаются в общей памяти.

Назначьте для каждого 4-х потоковых блоков (Nmatricies / Nblocks) длинную цепочку матриц для умножения, причем один из четырех отвечает за каждый блок умножения.

Допустим, вы являетесь блоком потоков 1 из 4, и первой из умножаемых матриц является AxB. Вы несете ответственность за (1,1) результата - (AB) _1,1 = A _1,1 B _1,1 + A * +1075 * 1,2 * 1 076 ** В _{2,1 * +1078 *. Вы предварительно загружаете A 1,1 в myblock [0] в разделяемой памяти.}

• загрузить в myblock [1] = B _1,1 из глобальной памяти
• myblock [3] = myblock [0] * myblock [1] (матричный мульт, все в общей памяти)
• загрузить в myblock [1] = A _1,2 от глобальных
• загрузить в myblock [2] = B _2,1 от глобальных
• myblock [0] = myblock [3] + (myblock [1] * myblock [2]) (матрица мульт и сложение, все в общей памяти).

Теперь вы можете повторить это для остальной части последовательности матриц в вашей части цепочки, выводя только когда закончите.

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

Опять же, нет реальной причины делать это, за исключением того, что вы не можете потрудиться передавать данные в GPU поэтапно, и производительность почти наверняка будет хуже; меньше записей в глобальную память, но вы, вероятно, заплатите за это с помощью GEMM. Хорошая новость заключается в том, что 50 не кратно 8, так что вам, вероятно, не придется слишком много избегать конфликтов банков с общей памятью.

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

Answer 3

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

В отличие от NVidiaCUBLAS (или Intel MKL), LIBXSMM не использует пакетный интерфейс.Вместо этого можно организовать отдельные вызовы, а также указать «следующие местоположения», т. Е. Где расположены операнды / матрицы следующих умножений (в памяти).Преимущество состоит в том, что явная структура данных или индексный формат, описывающий пакет, не нужны.

#include <libxsmm.h>

int main()
{
  const libxsmm_gemm_prefetch_type prefetch = LIBXSMM_PREFETCH_AUTO;
  const double alpha = 1.0, beta = 1.0; /* accumulate C-matrix */
  const int m = 23, n = 23, k = 23;     /* some problem size */
  libxsmm_dmmfunction xmm = NULL;       /* function pointer */

  xmm = libxsmm_dmmdispatch(23, 23, 23, /* some problem size */
          NULL/*lda*/, NULL/*ldb*/, NULL/*ldc*/,
          &alpha, &beta, NULL/*flags*/,
          NULL/*&prefetch*/);

  if (xmm) { /* JiT'ted code has been generated */
#   pragma omp parallel for
    for (int i = 0; i < nbatch; ++i) {
      const double *const ai = a + i * asize;
      const double *const bi = b + i * bsize;
      /* e.g., matrix C is accumulated (instead of streamed) */
      double *const ci = c /*+ i * csize*/;
      /* optionally provide "next locations" */
      xmm(ai, bi, ci/*,
          ai + 1 * asize,
          bi + 1 * bsize,
          ci + 0 * csize
      */);
    }
  }
}

LIBXSMM создает высокооптимизированный и специализированный код (JiT), который использует последние расширения набора команд (SSE3, AVX, AVX2и AVX-512).LIBXSMM доступен по недопустимой лицензии (пункт BSD-3).

ПРИМЕЧАНИЕ: Это не о CUBLAS (как первоначально запрашивалось).