做多矩阵的矩阵乘法在一个操作

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在他们的文档中演示，但它的速度是Cublas的4倍左右。有谁知道CUBLAS是如何工作的？如果代码在某处可用？

Answer 1

你看过最新的CUBLAS (version 4.1)？它包括专门用于大批量小型矩阵 - 矩阵乘法的新批次GEMM模式。我建议像Jonathan Dursi在他的回答中所建议的那样使用成对的乘法树，使用CUBLAS批处理API加速它，而不是像他所建议的那样编写自己的定制内核。

CUBLAS 4.1包含在CUDA Toolkit v4.1中。

Answer 2

问题是，cublas等被设计为使用所有的SM乘以大的矩阵。这不是你想要的;你想做很多小矩阵乘法。

可能有一些方法可以将这种情况转化为CUBLAS可以为你做的好事，但我没有看到它。我的建议如下：

编写一个内核，使用一个线程块来乘以两个小矩阵，并输出结果。

然后用块吨启动内核日志 N和解决乘法成对：

• 步骤1：乘中号×M个，男×M个 ... M _N-2 x M _N-1，输出M'，M' .. M ' _N/2
• 步骤2：乘M' X M '，M' X M ' ... M' _{N/2 - 2}×M个 _N/2-1，输出M '' ，M '' .. M '' _{N/4 ...}

等

会有50％的内存开销的因素，但我认为你会更好地利用你的核心这种方式。

更新

好吧，如果你真的不希望这样做的阶段，你可以做这种方式，但它会需要更多的编码，并与你的表现可能会加重病情可以得到像cuBLAS和异步传输。我假设你使用的是费米，并且你已经关闭了L1缓存，所以每个SM有48K共享内存。

以2×2块形式存储100个矩阵，每个块在内存中保持连续。所以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乙 _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]从全球
• myblock = B _2,1 [0] = myblock [3] +（myblock [ 1] * myblock [2]）（矩阵mult和加法，全部在共享内存中）。

现在，您可以对链中其余的矩阵序列重复此操作，仅在完成时输出。

完成后，最终会在全局内存中生成（#SMs）matricies，但仍然需要相乘，但在全局内存中不会有任何额外的临时存储空间，不得不将数据复制到全球存储器中，而不是原始的矩阵和要处理的列表。

再说一次，除了不会因为分阶段将数据发送到GPU而烦恼之外，没有真正的理由这样做，性能几乎肯定会变差;全球内存写入较少，但您可能会用手动GEMM支付。好消息是50不是8的倍数，所以你可能不会在共享内存银行冲突方面有太多的东西。

同样，对于奖励积分，您可以预先计算所有成对矩阵产品的所有块，然后再计算列表长度的一半。

Answer 3

LIBXSMM - 库瞄准英特尔架构为小型，密集或稀疏矩阵乘法，和小回旋正好旨在利用对小矩阵乘法最佳性能。

与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）。根据非许可许可（BSD-3条款），LIBXSMM为available。

注意：这不是关于CUBLAS（最初要求）。