使用 openmp 并行化矩阵乘法并使用 avx2 进行矢量化

Posted 2023-02-16

【中文标题】使用 openmp 并行化矩阵乘法并使用 avx2 进行矢量化

【英文标题】：Matrix multiplication paralellize with openmp and vectorized with avx2

【发布时间】：2020-04-26 19:51:10

【问题描述】：

(1)对于某些大小（矩阵大小）代码工作正常，但对于某些大小，它计算错误的矩阵乘法，虽然我使用 Avx2 指令集仔细但我不知道问题出在哪里。

(2) 当我仅使用 Avx2 指令集对代码进行向量化时，与使用 avx2 对代码进行向量化并使用 Openmp 并行化时相比，执行时间更少。虽然同时使用向量化（Avx2）和并行化（Openmp）时，执行时间应该会更短。

void mat_mul_pl(int size, double **mat1, double **mat2, double **result)

__m256d vec_multi_res = _mm256_setzero_pd(); //Initialize vector to zero
__m256d vec_mat1 = _mm256_setzero_pd(); //Initialize vector to zero
__m256d vec_mat2 = _mm256_setzero_pd();


int i, j, k;

// #pragma omp parallel for schedule(static)
for (i = 0; i < size; i++)

    for (j = 0; j < size; ++j)
    
        //Stores one element in mat1 and use it in all computations needed before proceeding
        //Stores as vector to increase computations per cycle
        vec_mat1 = _mm256_set1_pd(mat1[i][j]);
#pragma omp parallel for
        for (k = 0; k < size; k += 8)
        
            vec_mat2 = _mm256_loadu_pd((void*)&mat2[j][k]); //Stores row of second matrix (eight in each iteration)
            vec_multi_res = _mm256_loadu_pd((void*)&result[i][k]); //Loads the result matrix row as a vector
            vec_multi_res = _mm256_add_pd(vec_multi_res ,_mm256_mul_pd(vec_mat1, vec_mat2));//Multiplies the vectors and adds to th the result vector

            _mm256_storeu_pd((void*)&result[i][k], vec_multi_res); //Stores the result vector into the result array
        
    
    

【问题讨论】：

虽然在同时使用矢量化 (Avx2) 和并行化 (Openmp) 时执行时间应该更短。 并非总是如此。对于足够小的向量，启动[和停止]线程的开销高于非线程情况。此外，[假设x86 cpu]，某些处理器型号上的某些 SMT cpu 内核可能 [必须] 共享 H/W，并且可能必须对所使用的 SIMD H/W 进行序列化访问。更多的线程意味着更多的内存总线争用[过多的线程使内存总线饱和]。 openmp 是否理解 AVX2 等。人。缓存的引用位置呢？

您的乘法代码假定矩阵大小是 8 的倍数。通常 OpenMP 指令应位于最外层循环（本例中为 i 循环）。

@CraigEstey：“序列化”并不是对 SMT / 超线程在物理内核的加载/存储和 FMA 执行单元上竞争周期的准确描述。是的竞争，但仅是争夺前端周期和整数 ALU 端口。但是序列化意味着一些排序。但是，是的，这可能不受延迟或分支未命中的限制，因此 SMT 并不是很有帮助；单个硬件线程可能会使该内核的内存管道保持满。每个 FMA 有 2 次加载和 1 次存储，即使它全部命中 L1d 缓存，这也将成为瓶颈。

@PeterCordes 嗨，彼得。是的，“序列化”有点松散。也许“失速”会更好？如果 [两个] 同级超线程共享相同的 FPU/SIMD H/W（特定于模型？），则一个必须等​​待另一个完成。或者，如果 H/W 有 N 个流水线阶段，并且第 1 个兄弟姐妹已经完成第 0 阶段[在第 1 阶段工作]，那么第 2 个兄弟姐妹可以立即加载到第 0 阶段吗？是否（如何）SIMD H/W 负载平衡以防止兄弟姐妹垄断它？

@CraigEstey：我会说“争夺”或“竞争”对这些执行资源的访问权。 “停顿”通常意味着更大的等待时间，例如高速缓存未命中，这会停顿整个流水线，而不是执行一条指令或在一个 ALU 端口上丢失一个周期。 OoO exec 意味着将一条指令延迟一个周期不会延迟独立工作。顺便说一句，同级超线程竞争一切，包括整数 ALU；这就是超线程/ SMT 的重点。也许您正在考虑 Bulldozer 式 CMT，其中两个弱整数内核共享一个 FPU/SIMD 单元。

【参考方案1】：

我在编译时没有遇到任何错误，但是当我执行程序时，我得到一个分段错误（核心转储）。通过调试，我注意到在 vec_mat2 = _mm256_load_pd((void*)&mat2[j][k]) 这条指令处，核心以一定大小（大于 100）转储。请参阅下面的更新代码按照您的指示。

void mat_mul_pl(int size, double **mat1, double **mat2, double **result)


int i, j, k;

#pragma omp Parallel shared(vec_mat1,vec_mat2,vec_multi_res) private(i,j,k)


#pragma omp for schedule(static)
    for (i = 0; i < size; i++)
           
        for (j = 0; j < size; ++j)
        
            __m256d vec_mat1 = _mm256_setzero_pd(); 
            vec_mat1 = _mm256_set1_pd(mat1[i][j]);
            for (k = 0; k < size; k += 4)
             
                __m256d vec_multi_res = _mm256_setzero_pd(); 
                __m256d vec_mat2 = _mm256_setzero_pd();
                vec_mat2 = _mm256_load_pd((void*)&mat2[j][k]); 
                vec_multi_res = _mm256_load_pd((void*)&result[i][k]); 
                vec_multi_res = _mm256_add_pd(vec_multi_res ,_mm256_mul_pd(vec_mat1, vec_mat2));
                _mm256_store_pd((void*)&result[i][k], vec_multi_res);