为何printf调用能提升OpenMP稀疏矩阵多向量乘法的执行效率？

我在基于OpenMP实现稀疏矩阵-多向量乘法（SpMM）时，发现一个诡异现象：在未并行化的预处理函数的循环中加入若干printf调试调用后，后续并行乘积计算的GFLOPS值从约2大幅提升至10。我尝试将计时方式从clock_gettime()改为omp_get_wtime()，结果几乎一致。

核心主函数代码 #

int main() {

    // some stuff

    // 串行乘积计算
    clock_gettime(CLOCK_MONOTONIC, &t1);
    serial_product(a, x, k, y_s);
    clock_gettime(CLOCK_MONOTONIC, &t2);

    gflops_s = flop / ((t2.tv_sec - t1.tv_sec) * 1.e9 + (t2.tv_nsec - t1.tv_nsec));

    pre_processing(...) // 这里是用于负载均衡的预处理函数

    clock_gettime(CLOCK_MONOTONIC, &t1);
    openmp_spmm(a, rows_idx, num_threads, x, k, y_p);
    clock_gettime(CLOCK_MONOTONIC, &t2);

    gflops_p = flop / ((t2.tv_sec - t1.tv_sec) * 1.e9 + (t2.tv_nsec - t1.tv_nsec));

    // some other stuff
}

如上代码所示，若不在预处理函数的循环中加入调试用的printf，gflops_p约为2 GFLOPS；加入后则跃升至10 GFLOPS，而这些printf与OpenMP并无直接关联。

预处理函数代码 #

void pre_processing(int threads, int tot_nz, const int* irp, int tot_rows, int* rows_idx){
    int j, nz, nz_prev = 0, nz_curr, start_row = 0, r_prev, r_curr = 0;

    for (int i = 0; i < threads; i++) {
        rows_idx[i] = start_row; 

        printf("."); // 为什么这个打印能大幅加速OpenMP乘积计算？

        nz_curr = ((i + 1) * tot_nz) / threads;
        nz = nz_curr - nz_prev;
        nz_prev = nz_curr;

        for (j = start_row; j < tot_rows; j++) {
            r_curr += irp[j + 1] - irp[j]; 

            if (r_curr < nz) { 
                r_prev = r_curr; 
            } else {
                start_row = ((r_curr - nz) < (nz - r_prev)) ? j + 1 : j;
                break;
            }
        }

        r_curr = 0;
    }
    rows_idx[threads] = tot_rows;
    printf("\n");
}

我完全不清楚该现象的原因，猜测可能与stdout刷新、时钟周期或CPU利用率有关？

运行环境（EDIT 1） #

程序运行于搭载CentOS Stream 8（Linux内核4.18.0-448.el8.x86_64）的服务器，CPU参数如下：

Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              40
On-line CPU(s) list: 0-39
Thread(s) per core:  2
Core(s) per socket:  10
Socket(s):           2
NUMA node(s):        2
Vendor ID:           GenuineIntel
CPU family:          6
Model:               85
Model name:          Intel(R) Xeon(R) Silver 4210 CPU @ 2.20GHz
Stepping:            7
CPU MHz:             2200.000
CPU max MHz:         3200.0000
CPU min MHz:         1000.0000
BogoMIPS:            4400.00
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            1024K
L3 cache:            14080K
NUMA node0 CPU(s):   0-9,20-29
NUMA node1 CPU(s):   10-19,30-39

SpMM核心实现与猜测（EDIT 2） #

我猜测该现象可能与SpMM例程中使用的AVX指令有关，有观点指出问题可能源于编译器的向量化与内联优化。这是我首次使用这类指令，对其执行影响理解不足。SpMM代码如下：

void spmm_csr_64(CSR *mat, const int* rows_load, int threads, const Type* x, int k, Type* y){

    int *irp = mat->IRP;
    int *ja = mat->JA;
    Type *as = mat->AS;

    const __m256i scale = _mm256_set1_epi32(k);

    #pragma omp parallel for num_threads(threads) shared(threads, rows_load, irp, k, as, ja, x, y, scale) default(none)
    for (int tid = 0; tid < threads; tid++) {   // 按线程ID并行化
        // 私有参数
        int j, z, iter, lim, r_y, r_x;
        Type val;
        __m256i cols;
        __m512d vals, x_r;

        __m512d t[k];
        #pragma omp unroll partial
        for (z = 0; z < k; z++) {
            t[z] = _mm512_setzero_pd(); // 初始化t向量
        }

        for (int i = rows_load[tid]; i < rows_load[tid + 1]; i++) { // 线程处理分配到的行
            j = irp[i];
            lim = (irp[i+1] - j) / 8;

            for (iter = 0; iter < lim; iter++) {

                vals = _mm512_loadu_pd(&as[j]);                     // 加载8个64位元素
                cols = _mm256_loadu_si256((__m256i*)&ja[j]);        // 加载8个32位列索引
                cols = _mm256_mullo_epi32(cols, scale);             // 将列索引缩放至x的首列位置
                x_r = _mm512_i32gather_pd(cols, x, sizeof(Type));   // 根据列索引从x中收集元素
                t[0] = _mm512_fmadd_pd(vals, x_r, t[0]);            // 执行融合乘加操作

                for (z = 1; z < k; z++) {
                    cols = _mm256_add_epi32(cols, _MM8_1);
                    x_r = _mm512_i32gather_pd(cols, x, sizeof(Type));      // 根据列索引从x中收集元素
                    t[z] = _mm512_fmadd_pd(vals, x_r, t[z]);    // 执行融合乘加操作
                }

                j += 8;
            }

            r_y = i * k;

            // 处理剩余元素（当元素数不是8的倍数时）
            for (; j < irp[i+1]; j++) {
                val = as[j];
                r_x = ja[j] * k;

                #pragma omp unroll partial
                for (z = 0; z < k; z++) {
                    y[r_y + z] += val * x[r_x + z];
                }
            }

            // 对t中的所有64位元素求和归约
            #pragma omp unroll partial
            for (z = 0; z < k; z++) {
                y[r_y + z] += _mm512_reduce_add_pd(t[z]);
                t[z] = _mm512_setzero_pd();
            }
        }
    }
}

我使用Type作为double的占位符，便于后续切换为float，目前使用双精度。同时欢迎提供代码优化建议。

内容的提问来源于stack exchange，提问作者lilith