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

Question

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

 
Res = M1.M2.M3. ... .Mn 

Мои матрицы действительно маленькие 100x100 поплавки, но последовательность действительно длинная, порядка миллиарда.

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

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

 
Res1 = M1.M2.M3. ... .Mn/1000-1 
Res1 = M1+n/1000.M2+n/1000.M3+n/1000. ... .M2(n/1000)-1 
... 
Res1 = M1+999*n/1000.M2+999*n/1000.M3+999*n/1000. ... .M1000*(n/1000)-1 
Res = Res1*Res2* ... *Res999 

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

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

Answer 1

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

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

Answer 2

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

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

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

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

• Шаг 1: умножить M х М , М х М ... М _{Н - 2} х М _Н-1, выводя М ', М' .. М ' _N/2
• Шаг 2: умножить М' х М ', М' х М ' ... М' _{N/2 - 2} х М _Н/2-1, вывод M '' , М '' .. М '' _N/4 ...

т.д.

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

Update

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

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

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

Предположим, вы являетесь файловым блоком 1 из 4, и первой из матриц, которые вы должны умножить, является AxB. Вы несете ответственность за (1,1) результата - (AB) _1,1 = A _1,1 B _1,1 + A _1,2 * B _2,1. Вы предварительно загрузите в 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]) (матрица mult и дополнение, все в общей памяти).

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

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

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

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

Answer 3

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

В отличие от NVIDIA CUBLAS (или 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 составляет available по лицензии без разрешения (предложение BSD-3).

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