当行不是内存连续时，用于矩阵乘法的英特尔 MKL

Question

我们的硬件是英特尔至强融核，因此我们鼓励我们使用英特尔 MKL 替换手写的线性代数运算（例如方阵乘法），以充分利用它。

问题是对于这种情况应该使用哪种正确的 MKL，因为我们的矩阵行在内存中不连续的问题可能会禁止使用某些函数，例如cblas_dgemm.

这是稀疏 BLAS 的用例吗？

具有非连续行的矩阵示例：

#include <iostream>

int main()
{

        // construct this matrix:
        //
        // ( 1 2 3 )
        // ( 4 5 6 )

        const int NCOLS = 3;

        // allocate two perhaps-not-contiguous blocks of memory
        double *row1 = (double*)malloc(NCOLS * sizeof(double));
        double *row2 = (double*)malloc(NCOLS * sizeof(double));

        // initialize them to the desired values, i.e.
        row1[0] = 1;
        row1[1] = 2;
        row1[2] = 3;
        row2[0] = 4;
        row2[1] = 5;
        row2[2] = 6;

        // allocate a block of two memory elements the size of a pointer
        double **matrix = (double**)malloc(2 * sizeof(double*));

        // set them to point to the two (perhaps-not-contiguous) previous blocks
        matrix[0] = &row1[0];
        matrix[1] = &row2[0];

        // print
        for (auto j=0; j<2; j++)
        {
                for (auto i=0; i<3; i++)
                {
                        std::cout << matrix[j][i] << ",";
                }
                std::cout << "\n";
        }
}

Answer 1

密集 BLAS 操作可以对具有给定固定步幅的矩阵进行操作，但在您的情况下，步幅不是恒定的。稀疏矩阵旨在对包含大量零的矩阵进行运算，这显然不是您的情况（至少在提供的示例中不是）。

由于您的矩阵实际上很大 (20k x 20k)，最快的解决方案是 复制矩阵行 成一个大的连续矩阵。事实上，即使支持 BLAS 密集操作，它在数组的数组上也不会有效（无论如何 AFAIK 并非如此）。这是因为数组的数组通常不能有效地存储在缓存中（如果幸运的话，可以几乎连续地分配行）并且需要更复杂的索引。

例如，在像我这样的 x86-64 PC 上，具有 ~40 GB/s RAM 和 i5-9600KF 处理器，能够在此类矩阵上以 ~300 GFlops 的速度运行，O(n**3) 矩阵乘法大约需要 2 * 20_000**3 / 300e9 = 53.3 秒。同时，矩阵的最佳副本大约需要 2 * 8 * 20_000**2 / 40e9 = 0.16 秒。因此，复制的时间与实际的矩阵乘法时间相比可以忽略不计。也就是说，肯定需要三个副本（两个输入矩阵和一个输出矩阵）。此外，快速 BLAS 实现使用具有更好渐近复杂度的 Strassen 算法，在这种情况下应该快大约 2~3 倍。尽管如此，与实际矩阵乘法时间（>17.8 秒）相比，所有副本的时间（~0.5 秒）仍然非常小。

在 KNL Xeon Phi 上，MCDRAM 达到 ~400 GB/s 的吞吐量，主 RAM 达到 ~90 GB/s 的吞吐量，而处理器可以达到 3 TFlops。因此，如果将矩阵存储在 MCDRAM 中，则矩阵乘法的结果应为 1.8~5.3 秒，所有副本的结果应为 0.05 秒。如果矩阵存储在慢速 DRAM 中，则复制时间为 0.21 秒，与计算时间相比要大得多，但仍然没有那么大。如果您没有足够的 space 来将矩阵存储在 MCDRAM 中，那么您可以 将矩阵拆分成大块 （例如 10k x 10k）并分别计算每个块（为每个图块使用副本和 BLAS DGEMM）。

如果你想获得更高的性能，那么你可以使用 Xeon Phi 处理器的几个线程来复制一些块，以便重叠复制时间与计算时间.但是，这肯定会使代码变得更复杂以进行小的改进。

当行不是内存连续时，用于矩阵乘法的英特尔 MKL

Intel MKL for matrix multiplication when rows are not memory-contiguous

c++

blas

intel-mkl