gcc 不自动向量化矩阵向量乘法

Question

我刚刚开始使用我的矢量化代码。我的矩阵向量乘法代码没有被 gcc 自动向量化，我想知道为什么。 This pastebin contains the output from -fopt-info-vec-missed.

我无法理解输出告诉我的内容以及它与我在代码中编写的内容的匹配程度。

例如，我看到很多行说 not enough data-refs in basic block，我在网上搜索 google 找不到更多详细信息。我还看到存在与内存对齐相关的问题，例如Unknown misalignment, naturally aligned 和 vector alignment may not be reachable。我所有的内存分配都是针对使用 malloc 的 double 类型的，我相信它可以保证与该类型对齐。

环境：在 WSL2

上使用 gcc 编译

gcc -v: gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04)

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#define N 4000 // Matrix size will be N x N
#define T 1

//gcc -fopenmp -g vectorisation.c -o main -O3 -march=native -fopt-info-vec-missed=missed.txt

void doParallelComputation(double *A, double *V, double *results, unsigned long matrixSize, int numThreads)
{
    omp_set_num_threads(numThreads);
    unsigned long i, j;

    #pragma omp parallel for simd private(j)
    for (i = 0; i < matrixSize; i++)
    {
        // double *AHead = &A[i * matrixSize];
        // double tmp = 0;

        for (j = 0; j < matrixSize; j++)
        {
            results[i] += A[i * matrixSize + j] * V[j];
            // also tried tmp += A[i * matrixSize + j] * V[j];
        }
        // results[i] = tmp;
    }

}

void genRandVector(double *S, unsigned long size)
{
    srand(time(0));
    unsigned long i;

    for (i = 0; i < size; i++)
    {
        double n = rand() % 5;
        S[i] = n;
    }
}

void genRandMatrix(double *A, unsigned long size)
{
    srand(time(0));
    unsigned long i, j;
        for (i = 0; i < size; i++)
        {
            for (j = 0; j < size; j++)
            {
                double n = rand() % 5;
                A[i*size + j] = n;
            }

        }
    }

int main(int argc, char *argv[])
{

    double *V = (double *)malloc(N * sizeof(double));     // v in our A*v = parV computation
    double *parV = (double *)malloc(N * sizeof(double));  // Parallel computed vector
    double *A = (double *)malloc(N * N * sizeof(double)); // NxN Matrix to multiply by V
    genRandVector(V, N);
    doParallelComputation(A, V, parV, N, T);

    free(parV);
    free(A);
    free(V);
    
    return 0;
}

Answer 1

添加 double *restrict results 以承诺 non-overlapping input/output 有所帮助，没有 OpenMP 但有 -ffast-math。 https://godbolt.org/z/qaPh1v

您需要具体告知 OpenMP 缩减，让它放宽 FP-math 关联性。（-ffast-math 对 OpenMP 矢量器没有帮助）。有了它，我们也得到了你想要的：
#pragma omp simd reduction(+:tmp)

只有 restrict 而没有 -ffast-math 或 -fopenmp，你会得到完全的垃圾：它执行 SIMD FP 乘法，然后将其解压缩为 4x vaddsd标量累加器，根本无助于隐藏 FP 延迟。

使用 restrict 和 -fopenmp（没有 fast-math），它只做标量 FMA。

有 restrict 和 -ffast-math（ 没有 -fopenmp 或 #pragma 评论）它 auto-vectorizes 很好： vfmadd231pd ymm 在循环内部，在外部随机播放/添加水平总和。（但不并行化）。 https://godbolt.org/z/f36oG3

使用 restrict 和 -ffast-math（使用 -fopenmp）它仍然没有 auto-vectorize。 OpenMP 向量化器不同，可能不会利用 fast-math，而是需要您告诉它有关缩减的信息？

另请注意，对于您的数据布局，您要并行化（外部）的循环与您要使用 SIMD（内部）矢量化的循环不同。 两个输入内部 dot-product 循环的“向量”在连续的内存中，因此读取它们最有意义，而不是尝试将来自 4 个不同列的 SIMD 数据洗牌到一个向量中以累积 4 result[i+0..3] 结果1 个矢量。

但是，将外循环展开 4 次以使用来自 4 个不同列的数据的每个 V[j+0..3] 将提高计算强度（每个 FMA 接近 1 个负载，而不是 2 个）

(只要V[] 和矩阵的一行适合L1d缓存，这很好。如果不是，它实际上很糟糕，应该得到cache-blocked。或者实际上，如果你展开外循环，矩阵的 4 行。）

另请注意 double tmp = 0; 是个好主意：您当前的版本添加到 result[i]，请在编写前阅读它。这需要 zero-init 才能将其用作纯输出。

Auto-vec auto-par 版本：

我认为这是正确的； asm 看起来像 auto-par 等位化以及 auto-vectorizing 内部循环。

void doParallelComputation(double *restrict A, double *restrict V, double *restrict results, unsigned long matrixSize, int numThreads)
{
    omp_set_num_threads(numThreads);
    unsigned long i, j;

    #pragma omp parallel for private(j)
    for (i = 0; i < matrixSize; i++)
    {
        // double *AHead = &A[i * matrixSize];
        double tmp = 0;

         // TODO: unroll outer loop and cache-block it.
        #pragma omp simd reduction(+:tmp)
        for (j = 0; j < matrixSize; j++)
        {
            //results[i] += A[i * matrixSize + j] * V[j];
            tmp += A[i * matrixSize + j] * V[j];  // 
        }
        results[i] = tmp;  // write-only to results, not adding to old value.
    }

}

使用 OpenMPified 辅助函数内的矢量化内循环编译 (Godbolt) doParallelComputation._omp_fn.0:

# gcc7.5 -xc -O3 -fopenmp -march=skylake
.L6:
        add     rdx, 1               # loop counter; newer GCC just compares the end-pointer
        vmovupd ymm2, YMMWORD PTR [rcx+rax]            # 32-byte load
        vfmadd231pd     ymm0, ymm2, YMMWORD PTR [rsi+rax]  # 32-byte memory-source FMA
        add     rax, 32                                # pointer increment
        cmp     rdi, rdx
        ja      .L6

然后是循环后效率一般的横向求和；不幸的是，OpenMP 向量化器不如“普通”-ftree-vectorize 向量化器智能，但这需要 -ffast-math 在这里做任何事情。

gcc 不自动向量化矩阵向量乘法

gcc not autovectorising matrix-vector multiplication

c

gcc

simd

openmp

auto-vectorization

Auto-vec auto-par 版本：