c - Schönauer Triad 基准 - L1D 缓存不可见

Question

我们是两名参与著名的 Schönauer Triad Benchmark 的 HPC 学生，此处报告了其 C 代码及其简短解释:

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

#include <sys/time.h>

#define DEFAULT_NMAX 10000000
#define DEFAULT_NR DEFAULT_NMAX
#define DEFAULT_INC 10
#define DEFAULT_XIDX 0

#define MAX_PATH_LENGTH 1024

// #define WINOS
#define STACKALLOC

#ifdef WINOS 
    #include <windows.h>
#endif

static void dummy(double A[], double B[], double C[], double D[])
{
    return;
}

static double simulation(int N, int R)
{
    int i, j;

    #ifdef STACKALLOC
        double A[N];
        double B[N];
        double C[N];
        double D[N];
    #else
        double * A = malloc(N*sizeof(double));
        double * B = malloc(N*sizeof(double));
        double * C = malloc(N*sizeof(double));
        double * D = malloc(N*sizeof(double));
    #endif

    double elaps;

    for (i = 0; i < N; ++i)
    {
        A[i] = 0.00;
        B[i] = 1.00;
        C[i] = 2.00;
        D[i] = 3.00;
    }

    #ifdef WINOS
        FILETIME tp;
        GetSystemTimePreciseAsFileTime(&tp);
        elaps = - (double)(((ULONGLONG)tp.dwHighDateTime << 32) | (ULONGLONG)tp.dwLowDateTime)/10000000.0;
    #else
        struct timeval tp;
        gettimeofday(&tp, NULL);
        elaps = -(double)(tp.tv_sec + tp.tv_usec/1000000.0);
    #endif

    for(j=0; j<R; ++j)
    {
        for(i=0; i<N; ++i)
            A[i] = B[i] + C[i]*D[i];

        if(A[2] < 0) dummy(A, B, C, D);
    }

    #ifndef STACKALLOC
        free(A);
        free(B); 
        free(C);
        free(D);
    #endif

    #ifdef WINOS
        GetSystemTimePreciseAsFileTime(&tp);
        return elaps + (double)(((ULONGLONG)tp.dwHighDateTime << 32) | (ULONGLONG)tp.dwLowDateTime)/10000000.0;
    #else
        gettimeofday(&tp, NULL);
        return elaps + ((double)(tp.tv_sec + tp.tv_usec/1000000.0));
    #endif
}

int main(int argc, char *argv[])
{
    const int NR = argc > 1 ? atoi(argv[1]) : DEFAULT_NR;
    const int NMAX = argc > 2 ? atoi(argv[2]) : DEFAULT_NMAX;
    const int inc = argc > 3 ? atoi(argv[3]) : DEFAULT_INC;
    const int xidx = argc > 4 ? atoi(argv[4]) : DEFAULT_XIDX;

    int i, j, k;
    FILE * fp;

    printf("\n*** Schonauer Triad benchmark ***\n");

    char csvname[MAX_PATH_LENGTH];
    sprintf(csvname, "data%d.csv", xidx);

    if(!(fp = fopen(csvname, "a+")))
    {
        printf("\nError whilst writing to file\n");
        return 1;
    }

    int R, N;
    double MFLOPS;
    double elaps;

    for(N=1; N<=NMAX; N += inc)
    {
        R = NR/N;
        elaps = simulation(N, R);
        MFLOPS = ((R*N)<<1)/(elaps*1000000);
        fprintf(fp, "%d,%lf\n", N, MFLOPS);
        printf("N = %d, R = %d\n", N, R);
        printf("Elapsed time: %lf\n", elaps);
        printf("MFLOPS: %lf\n", MFLOPS);
    }

    fclose(fp);
    (void) getchar();
    return 0;
}

代码简单地循环 N 并且对于每个 N，它执行 NR 浮点运算，其中 NR 是一个常量，代表在每次最外层迭代中要执行的常量操作数，以便即使对于太短的 N 值也能进行准确的时间测量。要分析的内核显然是模拟子程序。

我们得到了一些奇怪的结果:

我们开始在 E4 E9220 服务器 2U 上对内核进行基准测试，该服务器由 8 个节点组成，每个节点都配备了双插槽 Intel Xeon E5-2697 V2(Ivy Bridge)@ 2.7 GHz，12 核。代码已使用 gcc (GCC) 4.8.2 编译，并已在 Linux CentOS release 6 上运行。下面列出了单个图像中的结果图:

N versus MFlops plots: -Ofast (above) and -Ofast along -march=native (below)

很容易看出 L2 和 L3 下坡非常明显，通过做一些简单的计算并考虑到多道程序问题以及 L2-L3 是 UNIFIED 和 L3 也在所有 12 个中共享的事实，它们在数值上是可以的核心。在第一个图中，L1 是不可见的，而在第二个图中，它是可见的，并且它以 N 值开始，因此根据每个内核的 L1D 大小，得到的 L1D 饱和值正好是 32 KB。第一个问题是:为什么我们没有看到没有 -march=native 架构特化标志的​​ L1 下坡？

经过一些棘手的(显然是错误的) self 解释后，我们决定在配备单路 Intel Core i7-3632QM (Ivy Bridge) 的 Lenovo Z500 上进行基准测试@ 2.2 GHz。这次我们使用了 gcc (Ubuntu 6.3.0-12ubuntu2) 6.3.0 20170406(来自 gcc --version)，结果如下:

N versus MFlops plots: -Ofast (above) and -Ofast along -march=native (below)

第二个问题有点自发:为什么这次我们看到没有 -march=native- 的 L1D 下坡？

Answer 1

有内部“TRIAD”循环的汇编片段 (A[i] = B[i] + C[i]*D[i]: per i迭代 2 次 double_precision 触发器，3 次读取 double，1 次写入 double)。

perf annotate 中的精确百分比不是很有用，因为您将具有不同性能的所有区域分析到单次运行中。而且长的性能报告根本没用，通常只需要在 # 之后的第一行 5-10 行。您可以尝试将测试限制在 4*N*sizeof(double) < sizeof(L1d_cache) 的有趣区域并重新收集 perf 注释并获得 perf stat ./program 和 的结果perf stat -d ./program(并了解英特尔特定的 perf 包装器 ocperf.py - https://github.com/andikleen/pmu-tools 和其他工具)。

来自 gcc-6.3.0 -Ofast - 128 位(2 个 double )XMM registers和 SSE2 movupd/movups使用(SSE2 是 x86_64 cpu 的默认 FPU)，i 的 2 次迭代用于每个汇编程序循环(movupd 从内存中加载 2 个 double)

         :                              A[i] = B[i] + C[i]*D[i];
    0.03 :        d70:       movupd (%r11,%rax,1),%xmm1    # load C[i:i+1] into xmm1
   14.87 :        d76:       add    $0x1,%ecx              # advance 'i/2' loop counter by 1
    0.10 :        d79:       movupd (%r10,%rax,1),%xmm0    # load D[i:i+1] into xmm0
   14.59 :        d7f:       mulpd  %xmm1,%xmm0            # multiply them into xmm0
    2.78 :        d83:       addpd  (%r14,%rax,1),%xmm0    # load B[i:i+1] and add to xmm0
   17.69 :        d89:       movups %xmm0,(%rsi,%rax,1)    # store into A[i:i+1]
    2.71 :        d8d:       add    $0x10,%rax             # advance array pointer by 2 doubles (0x10=16=2*8)
    1.68 :        d91:       cmp    %edi,%ecx              # check for end of loop (edi is N/2)
    0.00 :        d93:       jb     d70 <main+0x4c0>       # if not, jump to 0xd70

来自 gcc-6.3.0 -Ofast -march=native:vmovupd 不仅仅是 vector (SSE2 somethingpd 也是 vector )，它们是 AVX instructions可以使用 2 倍宽的寄存器 YMM(256 位，每个寄存器 4 个 double 值)。有更长的循环，但每个循环迭代处理 4 个 i 迭代

    0.02 :        db6:       vmovupd (%r10,%rdx,1),%xmm0   # load C[i:i+1] into xmm0 (low part of ymm0)
    8.42 :        dbc:       vinsertf128 $0x1,0x10(%r10,%rdx,1),%ymm0,%ymm1  # load C[i+2:i+3] into high part of ymm1 and copy xmm0 into lower part; ymm1 is C[i:i+3]
    7.37 :        dc4:       add    $0x1,%esi              # loop counter ++
    0.06 :        dc7:       vmovupd (%r9,%rdx,1),%xmm0    # load D[i:i+1] -> xmm0
   15.05 :        dcd:       vinsertf128 $0x1,0x10(%r9,%rdx,1),%ymm0,%ymm0  # load D[i+2:i+3] and get D[i:i+3] in ymm0
    0.85 :        dd5:       vmulpd %ymm0,%ymm1,%ymm0      # mul C[i:i+3] and D[i:i+3] into ymm0
    1.65 :        dd9:       vaddpd (%r11,%rdx,1),%ymm0,%ymm0  # soad 4 doubles of B[i:i+3] and add to ymm0
   21.18 :        ddf:       vmovups %xmm0,(%r8,%rdx,1)    # store low 2 doubles to A[i:i+1]
    1.24 :        de5:       vextractf128 $0x1,%ymm0,0x10(%r8,%rdx,1)  # store high 2 doubles to A[i+2:i+3]
    2.04 :        ded:       add    $0x20,%rdx             # advance array pointer by 4 doubles
    0.02 :        df1:       cmp    -0x460(%rbp),%esi      # loop cmp
    0.00 :        df7:       jb     db6 <main+0x506>       # loop jump to 0xdb6

启用 AVX 的代码(使用 -march=native)更好，因为它使用更好的展开，但它使用 2 个 double 的窄负载。随着更多的真实测试，数组将更好地对齐，编译器可能会选择最宽的 256-bit vmovupd进入 ymm，无需插入/提取指令。

您现在拥有的代码可能非常慢以致于它无法完全加载(饱和)到 L1 的接口(interface)数据缓存 在大多数情况下使用短数组。另一种可能性是数组之间对齐不良。

您在 /image/2ovxm.png 的下图中有高带宽的短尖峰- 6 个“GFLOPS”，这很奇怪。进行计算将其转换为 GByte/s 并找到 Ivy Bridge 的 L1d 带宽和负载发布率的限制...类似 https://software.intel.com/en-us/forums/software-tuning-performance-optimization-platform-monitoring/topic/532346 “Haswell 内核每个周期只能发出两个负载，因此它们必须是 256 位 AVX 负载才能有机会达到 64 字节/周期的速率。”(TRIAD 专家的话和 STREAM 的作者 John D. McCalpin，博士“带宽博士”，请搜索他的帖子)和 http://www.overclock.net/t/1541624/how-much-bandwidth-is-in-cpu-cache-and-how-is-it-calculated “L1 带宽取决于每滴答指令和指令步幅(AVX = 256 位，SSE = 128 位等)。IIRC，Sandy Bridge 每滴答有 1 条指令”