手动矢量化性能差异较大

问题描述 投票：0回答：1

我正在尝试手动向量化两个向量的点积的计算。请注意，我这样做是为了练习，并且我知道使用 BLAS 库会更合适。问题是，在对我的手动矢量化代码进行基准测试时，我所获得的性能差异非常大（我运行基准测试 10,20,...,90 次，并获取每组运行的统计信息）。如果它很重要，我会在 Windows 10 中的 WSL2 上运行所有基准测试，并在配备 AMD Ryzen 9 5900HS 的笔记本电脑上运行。

我的手动矢量化代码如下，

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

#define MATRIX_SIZE 1000

double matrix[MATRIX_SIZE][MATRIX_SIZE];
double vector[MATRIX_SIZE];
double result[MATRIX_SIZE];

#pragma GCC push_options
#pragma GCC optimize("O0")
void initialize() {
    // Initialize matrix and vector with random values
    srand(time(NULL));
    for (int i = 0; i < MATRIX_SIZE; i++) {
        for (int j = 0; j < MATRIX_SIZE; j++) {
            matrix[i][j] = (double)rand() / RAND_MAX;
        }
        vector[i] = (double)rand() / RAND_MAX;
    }
}
#pragma GCC pop_options

// https://stackoverflow.com/questions/49941645/get-sum-of-values-stored-in-m256d-with-sse-avx
double hsum_double_avx(__m256d v) {
    __m128d vlow  = _mm256_castpd256_pd128(v);
    __m128d vhigh = _mm256_extractf128_pd(v, 1); // high 128
            vlow  = _mm_add_pd(vlow, vhigh);     // reduce down to 128

    __m128d high64 = _mm_unpackhi_pd(vlow, vlow);
    return  _mm_cvtsd_f64(_mm_add_sd(vlow, high64));  // reduce to scalar
}

double dot_product(double* vec1, double* vec2) {
    __m256d accum1 = _mm256_setzero_pd();
    __m256d accum2 = _mm256_setzero_pd();
    __m256d accum3 = _mm256_setzero_pd();
    __m256d accum4 = _mm256_setzero_pd();
    __m256d accum5 = _mm256_setzero_pd();
    __m256d accum6 = _mm256_setzero_pd();
    __m256d accum7 = _mm256_setzero_pd();
    __m256d accum8 = _mm256_setzero_pd();

    double* stop_vec = vec1 + MATRIX_SIZE;

    while (vec1 + 32 <= stop_vec)
    {
        __m256d v11 = _mm256_loadu_pd(vec1);
        __m256d v12 = _mm256_loadu_pd(vec1+4);
        __m256d v13 = _mm256_loadu_pd(vec1+8);
        __m256d v14 = _mm256_loadu_pd(vec1+12);
        __m256d v15 = _mm256_loadu_pd(vec1+16);
        __m256d v16 = _mm256_loadu_pd(vec1+20);
        __m256d v17 = _mm256_loadu_pd(vec1+24);
        __m256d v18 = _mm256_loadu_pd(vec1+28);

        __m256d v21 = _mm256_loadu_pd(vec2);
        __m256d v22 = _mm256_loadu_pd(vec2+4);
        __m256d v23 = _mm256_loadu_pd(vec2+8);
        __m256d v24 = _mm256_loadu_pd(vec2+12);
        __m256d v25 = _mm256_loadu_pd(vec2+16);
        __m256d v26 = _mm256_loadu_pd(vec2+20);
        __m256d v27 = _mm256_loadu_pd(vec2+24);
        __m256d v28 = _mm256_loadu_pd(vec2+28);

        accum1 = _mm256_fmadd_pd(v11, v21, accum1);
        accum2 = _mm256_fmadd_pd(v12, v22, accum2);
        accum3 = _mm256_fmadd_pd(v13, v23, accum3);
        accum4 = _mm256_fmadd_pd(v14, v24, accum4);
        accum5 = _mm256_fmadd_pd(v11, v21, accum5);
        accum6 = _mm256_fmadd_pd(v12, v22, accum6);
        accum7 = _mm256_fmadd_pd(v13, v23, accum7);
        accum8 = _mm256_fmadd_pd(v14, v24, accum8);

        vec1 += 32;
        vec2 += 32;
    }

    long long diff = stop_vec-vec1;
    diff -= 1;

    __m256i mask = _mm256_setr_epi64x(diff, diff-1, diff-2, diff-3);
    __m256d v11 = _mm256_maskload_pd(vec1, mask);
    __m256d v21 = _mm256_maskload_pd(vec2, mask);
    accum1 = _mm256_fmadd_pd(v11, v21, accum1);

    mask = _mm256_setr_epi64x(diff-4, diff-5, diff-6, diff-7);
    __m256d v12 = _mm256_maskload_pd(vec1+4, mask);
    __m256d v22 = _mm256_maskload_pd(vec2+4, mask);
    accum2 = _mm256_fmadd_pd(v12, v22, accum2);

    mask = _mm256_setr_epi64x(diff-8, diff-9, diff-10, diff-11);
    __m256d v13 = _mm256_maskload_pd(vec1+8, mask);
    __m256d v23 = _mm256_maskload_pd(vec2+8, mask);
    accum3 = _mm256_fmadd_pd(v13, v23, accum3);

    mask = _mm256_setr_epi64x(diff-12, diff-13, diff-14, diff-15);
    __m256d v14 = _mm256_maskload_pd(vec1+12, mask);
    __m256d v24 = _mm256_maskload_pd(vec2+12, mask);
    accum4 = _mm256_fmadd_pd(v14, v24, accum4);

    mask = _mm256_setr_epi64x(diff-16, diff-17, diff-18, diff-19);
    __m256d v15 = _mm256_maskload_pd(vec1+16, mask);
    __m256d v25 = _mm256_maskload_pd(vec2+16, mask);
    accum5 = _mm256_fmadd_pd(v15, v25, accum5);

    mask = _mm256_setr_epi64x(diff-20, diff-21, diff-22, diff-23);
    __m256d v16 = _mm256_maskload_pd(vec1+20, mask);
    __m256d v26 = _mm256_maskload_pd(vec2+20, mask);
    accum6 = _mm256_fmadd_pd(v16, v26, accum6);

    mask = _mm256_setr_epi64x(diff-24, diff-25, diff-26, diff-27);
    __m256d v17 = _mm256_maskload_pd(vec1+24, mask);
    __m256d v27 = _mm256_maskload_pd(vec2+24, mask);
    accum7 = _mm256_fmadd_pd(v17, v27, accum7);

    mask = _mm256_setr_epi64x(diff-28, diff-29, diff-30, diff-31);
    __m256d v18 = _mm256_maskload_pd(vec1+28, mask);
    __m256d v28 = _mm256_maskload_pd(vec2+28, mask);
    accum8 = _mm256_fmadd_pd(v18, v28, accum8);

    accum1 = _mm256_add_pd(accum1, accum2);
    accum3 = _mm256_add_pd(accum3, accum4);
    accum5 = _mm256_add_pd(accum5, accum6);
    accum7 = _mm256_add_pd(accum7, accum8);
    accum1 = _mm256_add_pd(accum1, accum3);
    accum5 = _mm256_add_pd(accum5, accum7);
    accum1 = _mm256_add_pd(accum1, accum5);
    // accum1 = _mm256_add_pd(accum1, accum3);

    return hsum_double_avx(accum1);
}

// Auto dot product code
//void matrix_vector_multiply() {
    // Perform matrix-vector multiplication
//    for (int i = 0; i < MATRIX_SIZE; i++) {
//        result[i] = 0.0;
//        for (int j = 0; j < MATRIX_SIZE; j++) {
//            result[i] += matrix[i][j] * vector[j];
//        }
//    }
//}

void matrix_vector_multiply() {
    // Perform matrix-vector multiplication
    for (int i = 0; i < MATRIX_SIZE; i++) {
        result[i] += dot_product(matrix[i], vector);
    }
}

int main() {
    initialize();

    clock_t start_time = clock();
    matrix_vector_multiply();
    clock_t end_time = clock();
    double elapsed_time = ((double) (end_time - start_time)) / CLOCKS_PER_SEC;

    // Print time taken
    // printf("Time taken: %f seconds\n", elapsed_time);

    // Calculate total floating point operations
    long long flops = 2 * MATRIX_SIZE * MATRIX_SIZE;

    // Calculate GFLOPS
    double gflops = flops / (elapsed_time * 1e9);
    printf("%f\n", gflops);

    return 0;
}

我使用 

#pragma GCC push_options
来防止 GCC 将点积计算转移到初始化函数中（在我这样做之前，它在 -O3 上这样做）。我尝试使用 4 个累加器和 8 个累加器，8 个累加器对我来说效果更好。我使用命令 
gcc -O3 -march=native -mavx2 -mfma mat-vec-manual.c -o manual.o
进行编译。以下是自动矢量化的统计数据，

对于 GCC 的自动矢量化，我使用与上面相同的代码，但使用注释掉的 

matrix_vector_multiply()
函数。我使用命令 
gcc -O3 -march=native -ftree-vectorize -funroll-loops -ffast-math -mavx2 -mfma mat-vec-auto.c -o auto.o
进行编译。以下是手动矢量化的统计数据，

从统计中可以看到，虽然使用自动向量化后代码的性能有显着的提升，但方差也有很大的提升。检查生成的汇编代码 (gcc v11.4.0) 表明矢量化正在按预期进行，并且也正在使用 FMA 指令。

我认为这可能是由于数组未对齐到 32 字节而发生的，但将数组声明更改为，

double matrix[MATRIX_SIZE][MATRIX_SIZE] __attribute__((aligned(32)));
double vector[MATRIX_SIZE] __attribute__((aligned(32)));
double result[MATRIX_SIZE];

并将每次调用 

_mm256_loadu_pd
更改为 
_mm256_load_pd
并再次运行基准测试，结果没有明显变化。

我在这里缺少什么？或者这种类型的差异是预期的吗？

c performance vectorization avx2 fma
1个回答
0
投票

您的测试函数

matrix_vector_multiply
执行得非常快。在我的 i9-11950H 上，它在 1 毫秒内完成。为了减少方差，您可能需要在循环中执行该函数（例如 1000 次迭代）并进行测量。

您可能还想将 

clock()
替换为 
std::chrono::high_resolution_clock
调用。在我的 Windows 计算机上，
CLOCKS_PER_SEC
为 1000，这不足以测试该函数的一次传递。在 Linux 上，该值可能更高。

此外，您可以在测量部分之前进行热身。这将确保您的 CPU 不会因节能而降频。预热至少需要几毫秒。

我使用以下内容重新测试并获得了更一致的结果，在 13 到 17 GFlops 之间（在没有调试器的情况下运行 Release）。

 constexpr int N = 1000;
 // warmup
 for (int i = 0; i < N; i++)
     matrix_vector_multiply();

 // reset the dest array
 memset(result, 0, sizeof(double) * MATRIX_SIZE);

 auto t0 = std::chrono::high_resolution_clock::now();
 
 for(int i = 0; i < N; i++)
     matrix_vector_multiply();

 auto t1 = std::chrono::high_resolution_clock::now();
 double elapsed_time = 1e-9 * (double)std::chrono::duration_cast<std::chrono::nanoseconds>(t1 - t0).count();
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.