5

我正在实现一个算法，它本质上是一系列矩阵-矩阵乘法，如下所示：

水库 = M 1 .M 2 .M 3。... .M n

我的矩阵是非常小的 100x100 浮点数，但序列非常长，大约数十亿。

我尝试使用 CUBLAS 来进行矩阵乘法，但这很慢，但我确实注意到了一些有趣的事情。

将 100x100 与 100x100 矩阵相乘很慢，但将 1.000.000x100 与 100x100 相乘相对较快，这让我想到。如果我不是从左到右进行扫描，而是并行进行 10.000 次扫描。这应该很快，如果我在完成后乘以我的​​矩阵，我会得到相同的结果——只是更快。

水库1 = M 1 .M 2 .M 3。... .M n/1000-1 分辨率
1 = M 1+n/1000 .M 2+n/1000 .M 3+n/1000。... .M 2(n/1000)-1
...
水库1   = M 1+999*n/1000 .M 2+999*n/1000 .M 3+999*n/1000。... .M 1000*(n/1000)-1
分辨率 = 分辨率1 * 分辨率2 * ... * 分辨率999

M_1 ... M_n 在一组大约 100 个不同的矩阵中毫无价值，所以空间消耗并不是真正的问题，我需要做的就是在一次操作中进行多次乘法运算。

现在这是我的问题。我已经完成了一个矩阵-矩阵（sgemm）实现，灵感来自 nvidia 在他们的文档中展示的一个，但它的速度大约是 cublas 的 4 倍。有人知道 CUBLAS 是如何工作的吗？如果代码在某处可用？

4

3 回答 3

11

你看过最新的CUBLAS（4.1版）吗？它包括一个新的批处理 GEMM 模式，专门用于大批量的小型矩阵矩阵乘法。我建议像 Jonathan Dursi 在他的回答中建议的那样做一个成对的乘法树，使用 CUBLAS 批处理 API 来加速它，而不是像他建议的那样编写自己的自定义内核。

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

CUBLAS BATCHED GEMM API 提高了小矩阵批次的性能

于 2012-02-10T02:36:52.760 回答
2

问题是 cublas 等是为使用所有 SM 来乘以大型矩阵而设计的。那不是你想要的；你想做很多小矩阵乘法。

可能有一些方法可以把它变成 CUBLAS 可以为你做的事情，但我没有看到它。我的建议如下：

编写一个内核，使用一个线程块将两个小矩阵相乘，然后输出结果。

然后启动带有吨块的内核 log 2 N 并成对处理乘法：

  • 第 1 步：乘 M 1 x M 2 , M 3 x M 4 ... M N - 2 x M N-1，输出 M' 1 ,M' 2 .. M' N/2
  • 步骤 2：乘 M' 1 x M' 2 , M' 3 x M' 4 ... M' N/2 - 2 x M N/2-1，输出 M'' 1 ,M'' 2 .. M '' N/4 ...

等等

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

更新

好的，如果您真的不想分阶段执行此操作，您可以这样做，但它需要更多的编码，并且与使用 cuBLAS 和异步传输之类的东西可以获得的性能相比，性能可能会更差。我假设您使用的是 Fermi，并且您已关闭 L1 缓存，因此每个 SM 有 48K 共享内存。

以 2x2 块的形式存储 100 个矩阵，每个块在内存中连续。所以matrix[matrixnum,i,j]从matricies[matrixnum*100*100 + i*100*50 + j*50*50]. 请注意，每个块是 50*50*4 字节 ~ 10K，因此 4 可以轻松放入共享内存中。

为每个 4 个线程块分配一个 (Nmatricies/Nblocks) 长矩阵链以进行相乘，其中四个线程块中的一个负责乘法的每个块。

假设您是 4 个线程块 1，您要相乘的第一个矩阵是 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] (matrix mult, all in shared memory)
  • 从全局加载 myblock[1] = A 1,2
  • 从全局加载 myblock[2] = B 2,1
  • myblock[0] = myblock[3] + (myblock[1] * myblock[2]) （矩阵乘法和加法，都在共享内存中）。

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

完成后，您将在全局内存中得到 (#SMs) 矩阵，这些矩阵仍然需要相乘，但全局内存中不会有任何额外的临时存储，您也不必将数据复制到全局内存中，而不是原始矩阵和要处理的列表。

同样，没有真正的理由这样做，除非您不必费心将数据分阶段传送到 GPU，而且性能几乎肯定会更差；全局内存写入较少，但您可能会使用手动 GEMM 支付。好消息是 50 不是 8 的倍数，因此您可能不会有太多的共享内存库冲突。

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

于 2012-02-10T01:44:07.623 回答
0

LIBXSMM - 一个针对英特尔架构的库，用于小型、密集或稀疏矩阵乘法，而小型卷积正是为了利用小型矩阵乘法的最佳性能。

与 NVidia CUBLAS（或英特尔 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 无关（如最初要求的那样）。

于 2017-02-23T13:22:50.533 回答