首页 / 技术问答 / C++问答列表 / 问答详情

C++ 手动矢量化性能差异较大

我正在try 手动向量化两个向量的点积的计算.请注意，我这样做是一种练习，我知道使用BLAS库会更合适.问题是，在对我的手动向量化代码进行基准测试时，我获得了非常高的性能差异(我运行基准测试10，20，...，90次，并获得每组运行的统计数据).如果意义重大，我在Windows10的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编译.以下是自动矢量化的统计数据，

Manual Vectorization Statistics

For the automatic vectorization from GCC, I use the same code as above but use the commented out matrix_vector_multiply() function instead. I compile with the command gcc -O3 -march=native -ftree-vectorize -funroll-loops -ffast-math -mavx2 -mfma mat-vec-auto.c -o auto.o. Below are the statistics for the manual vectorization, Automatic Vectorization Statistics

正如您从统计数据中看到的，尽管使用自动向量化的代码的性能有了显著的提高，但方差也有很大的增加.判断生成的汇编代码(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，并再次运行基准，结果没有明显的变化.

我错过了什么？或者这种差异是预期的吗？

推荐答案

您的测试函数matrix_vector_multiply执行得非常快.在我的i9-11950H上，它在1毫秒内完成.为了减少差异，您可能希望在循环中执行函数(比如matrix_vector_multiply0次迭代)并进行测量.

您可能还需要将clock()电话替换为std::chrono::high_resolution_clock电话.在MSVC上，CLOCKS_PER_SEC等于clock()0，这不足以测试一遍函数.使用OP的GCC，这个价值可能会更高.

此外，您可以在测量部分之前进行热身.这将确保您的CPU没有关闭以节省电力.预热至少需要几毫秒.

我用下面的代码重新测试了一下，得到了更一致的结果，在13到17 GFlop之间(不使用调试器运行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();

C++相关问答推荐

C strlen on char array

不同到达时间的轮询实现

C指针算法在函数参数中的应用

当打印字符串时，为什么在c中没有使用常量限定符时我会收到警告？

使用NameSurname扫描到两个单独的字符串

如何将字符串argv[]赋给C中的整型数组？

字符是否必须转换为无符号字符，然后才能与getc家族的返回值进行比较？

如何在C中打印包含扫描字符和整数的语句？

对于C中给定数组中的每个查询，如何正确编码以输出给定索引范围(1到N)中所有数字的总和？

tick.q中的Kdb+键控表语法

变量的作用域是否在C中的循环未定义行为或实现定义行为的参数中初始化？

Wcstok导致分段故障

从C中的函数返回静态字符串是不是一种糟糕的做法？

GETS()在C++中重复它前面的行

即使我在C++中空闲，也肯定会丢失内存

在下面的C程序中，.Ap0是如何解释的？

将char*数组深度复制到 struct 中？

Ubuntu编译：C中的文件格式无法识别错误

为什么实现文件中的自由函数默认没有内部链接？

使用 SDL2 的 C 程序中的内存泄漏

实用课程推荐
相关教程推荐
休闲君评测 |  TOOLFK工具网 |  古诗文网 |  良许Linux教程网 |  Python之禅 |  红色石头 |  入党申请书 | 

© 2020 Hi LearnFk 意见反馈 免责声明 关于我们 桂ICP备11002319号-3 友链&广告位+QQ: 1963612630