为何小规模SAXPY测试有效带宽远超系统DRAM带宽？

测试背景与设置 #

开展「冷」微基准测试：调用SAXPY目标函数50次，每次运行重新分配数据，各工作线程将分配给自己的数据的每个内存页首个数据置零，以此排除页错误影响。测试基于Intel(R) Core(TM) i5-8400 CPU @2.80GHz（6核，无SMT），缓存配置为32KiB L1、1.5MiB L2、9MiB L3，支持AVX2指令集。

SAXPY微基准测试代码 #

/// File This file contains a micro-benchmark program for the SAXPY loop.

#include <string> // std::stoul
#include <utility> // std::abort
#include <iostream> // std::cerr
#include <cstdint> // std::uint64_t
#include <algorithm> // std::min
#include <limits> // std::numeric_limits

#include <omp.h>

#include <stdlib.h> // posix_memalign
#include <unistd.h> // sysconf

void saxpy(const float& a, const float* vec_x, float* vec_y, const std::size_t& n)
{
    #pragma omp parallel for schedule(static)
    for (auto ii = std::size_t{}; ii < n; ++ii)
    {
        vec_y[ii] += vec_x[ii] * a; // fp_count: 2, traffic: 2+1
    }
}

int main(int argc, char** argv)
{
    // extract the problem size
    if (argc < 2)
    {
        std::cerr << "Please provide the problem size as command line argument." << std::endl;
        return 1;
    }
    const auto n = static_cast<std::size_t>(std::stoul(argv[1]));
    if (n < 1)
    {
        std::cerr << "Zero valued problem size provided. The program will now be aborted." << std::endl;
        return 1;
    }
    if (n * sizeof(float) / (1024 * 1024 * 1024) > 40) // let's assume there's only 64 GiB of RAM
    {
        std::cerr << "Problem size is too large. The program will now be aborted." << std::endl;
        return 1;
    }
    
    // report
    std::cout << "Starting runs with problem size n=" << n << ".\nThread count: " << omp_get_max_threads() << "."
        << std::endl;
    
    // details
    const auto page_size = sysconf(_SC_PAGESIZE);
    const auto page_size_float = page_size / sizeof(float);
    
    // experiment loop
    const auto experiment_count = 50;
    const auto warm_up_count = 10;
    const auto run_count = experiment_count + warm_up_count;
    auto durations = std::vector(experiment_count, std::numeric_limits<std::uint64_t>::min());
    const auto a = 10.f;
    float* vec_x = nullptr;
    float* vec_y = nullptr;
    for (auto run_index = std::size_t{}; run_index < run_count; ++run_index)
    {
        // allocate
        const auto alloc_status0 = posix_memalign(reinterpret_cast<void**>(&vec_x), page_size, n * sizeof(float));
        const auto alloc_status1 = posix_memalign(reinterpret_cast<void**>(&vec_y), page_size, n * sizeof(float));
        if (alloc_status0 != 0 || alloc_status1 != 0 || vec_x == nullptr || vec_y == nullptr)
        {
            std::cerr << "Fatal error, failed to allocate memory." << std::endl;
            std::abort();
        }
        
        // "first touch"
        #pragma omp parallel for schedule(static)
        for (auto ii = std::size_t{}; ii < n; ii += page_size_float)
        {
            vec_x[ii] = 0.f;
            vec_y[ii] = 0.f;
        }
        
        // run experiment
        const auto t1 = omp_get_wtime();
        saxpy(a, vec_x, vec_y, n);
        const auto t2 = omp_get_wtime();
        const auto duration_in_us = static_cast<std::int64_t>((t2 - t1) * 1E+6);
        if (duration_in_us <= 0)
        {
            std::cerr << "Fatal error, no time elapsed in the test function." << std::endl;
            std::abort();
        }
        if (run_index + 1 > warm_up_count)
        {
            durations[run_index - warm_up_count] = static_cast<std::uint64_t>(duration_in_us);
        }
        
        // deallocate
        std::free(vec_x);
        std::free(vec_y);
        vec_x = nullptr;
        vec_y = nullptr;
    }
    
    // statistics
    auto min = std::numeric_limits<std::uint64_t>::max();
    auto max = std::uint64_t{};
    auto mean = std::uint64_t{};
    for (const auto& duration : durations)
    {
        min = std::min(min, duration);
        max = std::max(max, duration);
        mean += duration;
    }
    mean /= experiment_count;
    
    // report duration
    std::cout << "Mean duration:      " << mean << " us\n"
        << "Min. duration:      " << min << " us\n"
        << "Max. duration:      " << max << " us.\n";
    
    // compute effective B/W
    const auto traffic = 3 * n * sizeof(float);
    constexpr auto inv_gigi = 1.0 / static_cast<double>(1024 * 1024 * 1024);
    const auto traffic_in_gib = static_cast<double>(traffic) * inv_gigi;
    std::cout << "Traffic per run:    " << traffic << " B (" << traffic_in_gib << " GiB)\n" 
        << "Mean effective B/W: " << static_cast<double>(traffic_in_gib) / (static_cast<double>(mean) * 1E-6) << " GiB/s\n"
        << "Min. effective B/W: " << static_cast<double>(traffic_in_gib) / (static_cast<double>(max) * 1E-6) << " GiB/s\n"
        << "Max. effective B/W: " << static_cast<double>(traffic_in_gib) / (static_cast<double>(min) * 1E-6) << " GiB/s\n"
        << std::endl;
    
    return 0;
}

多组测试结果 #

Starting run for n=1000000.
Starting runs with problem size n=1000000.
Thread count: 6.
Mean duration:      148 us
Min. duration:      117 us
Max. duration:      417 us.
Traffic per run:    12000000 B (0.0111759 GiB)
Mean effective B/W: 75.5126 GiB/s
Min. effective B/W: 26.8006 GiB/s
Max. effective B/W: 95.5203 GiB/s

----------------------

Starting run for n=10000000.
Starting runs with problem size n=10000000.
Thread count: 6.
Mean duration:      3311 us
Min. duration:      3262 us
Max. duration:      3382 us.
Traffic per run:    120000000 B (0.111759 GiB)
Mean effective B/W: 33.7538 GiB/s
Min. effective B/W: 33.0452 GiB/s
Max. effective B/W: 34.2608 GiB/s

----------------------

Starting run for n=100000000.
Starting runs with problem size n=100000000.
Thread count: 6.
Mean duration:      32481 us
Min. duration:      32137 us
Max. duration:      36431 us.
Traffic per run:    1200000000 B (1.11759 GiB)
Mean effective B/W: 34.4074 GiB/s
Min. effective B/W: 30.6768 GiB/s
Max. effective B/W: 34.7757 GiB/s

----------------------

Intel MLC带宽测试结果 #

Intel(R) Memory Latency Checker - v3.9a
*** Unable to modify prefetchers (try executing 'modprobe msr')
*** So, enabling random access for latency measurements
Measuring idle latencies (in ns)...
        Numa node
Numa node        0  
       0      57.1  

Measuring Peak Injection Memory Bandwidths for the system
Bandwidths are in MB/sec (1 MB/sec = 1,000,000 Bytes/sec)
Using all the threads from each core if Hyper-threading is enabled
Using traffic with the following read-write ratios
ALL Reads        :  41302.4 
3:1 Reads-Writes :  37185.7 
2:1 Reads-Writes :  36779.7 
1:1 Reads-Writes :  36790.3 
Stream-triad like:  37241.2

核心疑问 #

当问题规模n=1000000时，测试得到的有效带宽（约75 GiB/s）为何会高于Intel MLC测得的系统最大读写带宽（约37 GiB/s）？

异常原因分析 #

1. 数据完全命中L3缓存
   当n=1000000时，vec_x和vec_y的总数据量为1000000 * 4B * 2 = 8MB，远小于i5-8400的9MiB三级缓存（LLC）。此时SAXPY运算的所有数据都被缓存到LLC中，测试计算的是缓存带宽而非内存带宽。LLC的带宽远高于内存，i5-8400的LLC带宽通常在70-100GiB/s区间，与测试结果完全匹配。

2. MLC测试的是内存峰值带宽
   MLC给出的"Stream-triad like"结果为37241.2 MB/s（约35.5 GiB/s），这是系统内存的峰值带宽。当数据量超过缓存容量时，运算必须从内存读写，此时带宽就会降到与MLC测试结果接近，这也对应了n=10000000和n=100000000时的测试数据。

3. "冷"测试的缓存命中误区
   虽然每次测试都重新分配数据并执行first touch，但first touch仅确保物理页被分配并绑定到NUMA节点，并未将数据从缓存中驱逐。由于数据量远小于LLC容量，这些数据会被保留在缓存中，后续的SAXPY运行直接命中缓存，并没有真正访问内存，因此所谓的"冷"测试在小数据量下并未达到预期的内存访问效果。

验证方法 #

可以将问题规模修改为刚好超过LLC容量的数值（比如n=2400000，总数据量为2400000*4B*2=19.2MB>9MiB），此时测试得到的有效带宽会降至与MLC测得的内存带宽一致，从而验证缓存命中的结论。

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