GEMM/SYRK: per-call malloc of the packing buffers page-faults on every call at some sizes under glibc, up to 2.2x slower

When the GEBP packing buffers of a product exceed EIGEN_STACK_ALLOCATION_LIMIT (128 KiB), general_matrix_matrix_product and the SYRK path allocate them with scratch_malloc on every call and free them on return (GeneralMatrixMatrix.h:313-314, GeneralMatrixMatrixTriangular.h:98-99). Under glibc, at some sizes that memory goes back to the operating system on every free, so every call page-faults its packing buffers in again: a double GEMM at n = 192 takes 112 minor faults per call and runs at 0.6 to 0.78 of its rate at n = 184 or 200. Pinning glibc's thresholds (MALLOC_MMAP_THRESHOLD_ and MALLOC_TRIM_THRESHOLD_) removes the faults and the slowdown completely. Which sizes are hit depends on the blocking sizes and on the process's earlier allocations, which glibc's dynamic mmap and trim thresholds follow, so the same product can be fast in one program and slow in another.

Measurements

Square C.noalias() += A * B and C.selfadjointView<Lower>().rankUpdate(A) on MatrixXd, best of five rounds; minor page faults per call from getrusage over 20 calls. Default glibc malloc against MALLOC_MMAP_THRESHOLD_=33554432 MALLOC_TRIM_THRESHOLD_=1073741824 (fixed thresholds also turn off glibc's dynamic adjustment). Master 210e0a609, one core, taskset-pinned.

Host, build	n	kc, mc, nc	default: GEMM / SYRK GFLOP/s	faults per call	fixed thresholds: GEMM / SYRK
Arm Cortex-X925, GCC 13.3, NEON	144	144, 144, 144	55.7 / 39.3	49	73.8 / 59.4
	160	160, 84, 160	74.0 / 59.7	0	72.2 / 60.1
	192	192, 192, 192	57.5 / 41.7	112	74.7 / 62.3
	200	200, 102, 200	73.8 / 63.4	0	75.1 / 63.1
	240	240, 240, 240	59.3 / 44.8	193	75.1 / 64.8
	448	448, 448, 224	75.1 / 52.5	752 (SYRK)	75.7 / 69.0
	512	512, 512, 256	74.1 / 53.3	992 (SYRK)	74.4 / 68.4
	640	328, 640, 320	75.3 / 71.8	0	75.2 / 71.3
Arm Cortex-X925, GCC 13.3, SVE at 128 bits	192	192, 192, 192	60.3 / 43.5	112	78.6 / 65.8
	512	512, 512, 256	78.4 / 53.6	992 (SYRK)	77.6 / 69.2
AMD Ryzen 7 7800X3D, GCC 15.2, `-mavx2 -mfma`	144	144, 144, 144	37.6 / 24.8	49	64.4 / 56.0
	176	176, 96, 176	43.0 / 31.0	61	64.7 / 57.7
	192	192, 192, 192	38.8 / 26.4	112	65.0 / 57.5
	200	200, 108, 200	65.0 / 58.4	0	65.1 / 58.3
	240	240, 240, 240	41.1 / 28.7	193	64.3 / 58.7

Every row with faults recovers fully under fixed thresholds; every row without faults is unchanged. No float cell faulted on the Cortex-X925, and an AVX-512 build on the x86 host showed no step.

The published cells

The comparison harness binaries behind the DGX Spark pages, re-run on the Cortex-X925 on the cells that dipped there; median of three repetitions, GFLOP/s:

Cell	Binary and cells in the process	default	fixed thresholds
GEMM, double, n = 192	SVE build, Eigen cells only (as in the SVE passes)	59.1	79.1
GEMM, double, n = 192	NEON build, Eigen and Arm Performance Libraries cells, n = 160, 192, 200 only	56.0	73.5
SYRK, float, n = 384	SVE / NEON	108.2 / 104.5	148.7 / 143.2
SYRK, float, n = 768	SVE / NEON	126.3 / 119.1	157.0 / 151.3
SYRK, double, n = 512	SVE / NEON	54.7 / 51.7	71.9 / 67.4

The neighboring Eigen cells (GEMM at n = 160 and 200; SYRK at n = 256, 500 and 1000 in float and 384, 500 and 768 in double) and the library's cells move by 1.2 % or less.

Where it shows on the website

Zen 5 AVX2: GEMM in double at n = 192 at 41 GFLOP/s against 75 at n = 160 and 200; SYRK at 27 against 62 and 68 at n = 128 and 200.
DGX Spark (libeigen/libeigen.gitlab.io!17): the n = 192 step of GEMM and SYRK on the SVE page and the SYRK dips at n = 384 and 768 (float) and 512 (double) on both pages are this effect (table above). In the full campaign runs the SVE processes hit the n = 192 step and the NEON processes did not; the NEON binary hits it as well when it runs those cells alone, so the difference was allocation history, not the instruction set.
NVIDIA Grace runs with 64 KiB pages, a sixteenth of the faults, and shows the n = 192 step only faintly (GEMM in double, NEON build: 38 against 40).

Possible directions

Keep the packing buffers alive between calls: a per-thread, grow-only workspace that gemm_blocking_space borrows when it has no buffer of its own. The cost is memory held per thread (mc * kc + kc * nc scalars, which grows with the problem), and it has to stay correct for nested and parallel products.
Or allocate through a size-class cache in scratch_malloc for buffers above the stack limit.
Until then, users who call mid-size products in a loop can avoid it with mallopt(M_MMAP_THRESHOLD, ...) and mallopt(M_TRIM_THRESHOLD, ...), or the environment variables above.

Reproducer

#include <Eigen/Dense>
#include <algorithm>
#include <chrono>
#include <cstdio>
#include <sys/resource.h>
using namespace Eigen;

long minorFaults() {
  rusage u;
  getrusage(RUSAGE_SELF, &u);
  return u.ru_minflt;
}

int main() {
  for (int n : {184, 192, 200}) {
    MatrixXd A = MatrixXd::Random(n, n), B = MatrixXd::Random(n, n), C = MatrixXd::Zero(n, n);
    C.noalias() += A * B;  // warm-up
    const int calls = 200;
    const long f0 = minorFaults();
    const auto t0 = std::chrono::steady_clock::now();
    for (int i = 0; i < calls; ++i) C.noalias() += A * B;
    const double t = std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count() / calls;
    std::printf("n = %d  %.1f GFLOP/s  %.1f minor faults per call\n", n, 2.0 * n * n * n / t * 1e-9,
                double(minorFaults() - f0) / calls);
  }
}

Run it as is and again with MALLOC_MMAP_THRESHOLD_=33554432 MALLOC_TRIM_THRESHOLD_=1073741824. On an NVIDIA Grace core (64 KiB pages, GCC 13.3) it reports 0, 7 and 6 faults per call at n = 184, 192 and 200 (36.5 and 37.5 GFLOP/s at 192 and 200 against 38.6 and 39.6 with fixed thresholds): with this program's allocation history n = 200 is hit as well.