尝试特别使用内在函数时出现分段错误_mm256_storeu_pd()答案

【问题标题】：Segmentation fault when trying to use intrinsics specifically _mm256_storeu_pd()尝试特别使用内在函数时出现分段错误_mm256_storeu_pd()
【发布时间】：2016-09-15 16:55:58
【问题描述】：

似乎自己通过在 mm256 调用中键入 cij2 指针来修复它

所以_mm256_storeu_pd((double *)cij2,vecC);

我不知道为什么这会改变任何东西......

我正在编写一些代码并尝试利用英特尔手动矢量化。但是每当我运行代码时，我都会在尝试使用我的双 *cij2 时遇到分段错误。

if( q == 0)
{
    __m256d vecA;
    __m256d vecB;
    __m256d vecC;
    for (int i = 0; i < M; ++i)
      for (int j = 0; j < N; ++j)
      {
        double cij = C[i+j*lda];
        double *cij2 = (double *)malloc(4*sizeof(double));
        for (int k = 0; k < K; k+=4)
        {
          vecA = _mm256_load_pd(&A[i+k*lda]);
          vecB = _mm256_load_pd(&B[k+j*lda]);
          vecC = _mm256_mul_pd(vecA,vecB);
          _mm256_storeu_pd(cij2, vecC);
          for (int x = 0; x < 4; x++)
          {
            cij += cij2[x];
          }

        }
        C[i+j*lda] = cij;
      }

我已将问题定位到 cij2 指针。如果我注释掉包含该指针的 2 行代码运行良好，它不会像它应该的那样工作，但它实际上会运行。

我的问题是为什么我会在这里遇到分段错误？我知道我已经正确分配了内存，并且内存是 256 个双精度向量，大小为 64 位。

在阅读了 cmets 之后，我来补充一些说明。我做的第一件事是使用 malloc 将 _mm_malloc 更改为普通分配。不应该影响任何一种方式，但理论上会给我更多的喘息空间。

第二个问题不是来自分配的空返回，我添加了几个循环来增加数组并确保我可以修改内存而不会崩溃，所以我相对确定这不是问题。问题似乎源于将数据从 vecC 加载到数组。

最后我不能使用 BLAS 调用。这是一个并行类。我知道调用比我更聪明的方法会简单得多，但不幸的是，如果我尝试这样做，我会得到 0。

【问题讨论】：

我不确定_mm_malloc()，但标准malloc() 和calloc() 在分配失败的情况下返回NULL。假设_mm_malloc() 做同样的事情，您将无法检查这种可能性。如果发生这种情况，当您随后尝试使用指针时出现分段错误也就不足为奇了。
另一方面，我观察到 _mm256_storeu_pd() 似乎需要一个 256 位对齐的指针，但您的分配只确保 64 位对齐。仅从文档来看，您似乎想使用_mm_malloc(4*sizeof(double), 256) 进行分配。
@JohnBollinger：storeu_pd 是未对齐的存储。 store_pd 是对齐存储。 _mm_malloc 以字节而不是位为单位进行对齐 arg，因此 64 是 64 字节对齐的。 _mm_malloc 可能会失败，因为 cij2 从未被释放。使用自动存储会更明智。
@PeterCordes, ... 这就是为什么我没有从我以前不熟悉的文档（显然不够充分）的阅读中做出答案。
此代码即使在您开始工作后也不会正常运行。在矩阵乘法的内部循环内的向量元素上编写标量循环绝对是可怕的。内部循环内的水平添加已经够糟糕了，但这样做很讨厌。特别是使用动态分配的缓冲区。使用 BLAS 库中经过良好调整的 DGEMM 函数，您将获得更好的结果。现代 CPU 的最佳矩阵乘法需要具有缓存感知能力，并且有很多技巧（例如，根据需要将一个数组部分转置）。我不建议你自己写。

标签： c pointers intrinsics

【解决方案1】：

您动态分配double *cij2 = (double *)malloc(4*sizeof(double));，但您从未释放它。这很愚蠢。使用double cij2[4]，特别是如果你不想费心去对齐它。您一次不需要多个暂存缓冲区，而且它的大小固定很小，所以只需使用自动存储即可。

在 C++11 中，您将使用 alignas(32) double cij2[4]，因此您可以使用 _mm256_store_pd 而不是 storeu。（或者只是为了确保 storeu 不会被未对齐的地址拖慢）。

如果你真的想调试你的原始文件，当它出现段错误时，使用调试器来捕获它，并查看指针值。确保它是明智的。

你测试内存是否有效的方法（比如循环它，或者注释掉一些东西）听起来可能会导致你的很多循环被优化掉，所以问题不会发生。

当您的程序崩溃时，您还可以查看 asm 指令。向量内在函数相当直接地映射到 x86 asm（编译器认为更有效的方式除外）。

如果您将水平总和从 k 上的循环中拉出，您的实现会少很多。不是存储每个乘法结果并水平相加，而是使用向量加法到向量累加器中。 hsum 它在 k 的循环之外。

    __m256d cij_vec = _mm256_setzero_pd();
    for (int k = 0; k < K; k+=4) {
      vecA = _mm256_load_pd(&A[i+k*lda]);
      vecB = _mm256_load_pd(&B[k+j*lda]);
      vecC = _mm256_mul_pd(vecA,vecB);
      cij_vec = _mm256_add_pd(cij_vec, vecC);  // TODO: use multiple accumulators to keep multiple VADDPD or VFMAPD instructions in flight.
    }
    C[i+j*lda] = hsum256_pd(cij_vec);  // put the horizontal sum in an inline function

对于良好的 hsum256_pd 实现（除了存储到内存和使用标量循环），请参阅Fastest way to do horizontal float vector sum on x86（我在其中包含了一个 AVX 版本。应该很容易将改组模式适应 256b 双精度。）这个将有助于您的代码很多，因为您仍然有 O(N^2) 水平总和（但没有 O(N^3) 与此更改）。

理想情况下，您可以并行累积 4 个 i 值的结果，而不需要水平求和。

VADDPD 具有 3 到 4 个时钟的延迟，每 1 到 0.5 个时钟的吞吐量为 1，因此您需要 3 到 8 个向量累加器来使执行单元饱和。或者使用 FMA，最多 10 个矢量累加器（例如在 Haswell 上，FMA...PD 具有 5c 延迟和每 0.5c 吞吐量一个）。请参阅Agner Fog's instruction tables and optimization guides 了解更多信息。还有x86 标签wiki。

此外，理想情况下，嵌套循环的方式可以让您连续访问三个数组中的两个，因为缓存访问模式对于 matmul（大量数据重用）至关重要。即使您不喜欢并一次转置适合缓存的小块。即使转置一个输入矩阵也可能是一种胜利，因为这会花费 O(N^2) 并加快 O(N^3) 过程。我看到您的内部循环当前在访问 A[] 时的步幅为 lda。

【讨论】：