Space Filling Curves is All You Need: Communication-Avoiding Matrix Multiplication Made Simple

Evangelos Georganas; Alexander Heinecke; Pradeep Dubey

cs.DCJan 22, 2026

Space Filling Curves is All You Need: Communication-Avoiding Matrix Multiplication Made Simple

Authors: Evangelos Georganas, Alexander Heinecke, Pradeep Dubey

Organizations: Intel Corporation

Abstract

General Matrix Multiplication (GEMM) is the cornerstone of HPC workloads and Deep Learning. State-of-the-art (SOTA) vendor libraries tune tensor layouts, parallelization schemes and cache blocking to minimize data movement across the memory hierarchy and maximize throughput. However, optimal settings for these parameters depend on the target platform and matrix shapes, making exhaustive tuning infeasible. In this work, we address this cumbersome scheduling search and tuning using space-filling curves (SFC). We partition the matrix multiplication using advancements in SFC, and obtain platform-oblivious and shape-oblivious matrix multiplication schemes with a high degree of data locality. We extend the SFC-based work partitioning to implement Communication-Avoiding (CA) algorithms with replication techniques in a seamless fashion. The resulting SFC-CA GEMM achieves provable asymptotic communication optimality for both square and rectangular matrix regimes. Across four x86 and Arm platforms, SFC-CA GEMM outperforms vendor libraries by up to 5.5 $\times$ per shape and 1.8 $\times$ in weighted harmonic mean (WHM) throughput. Last, we show the impact of our work on two real-world applications by leveraging our SFC-CA GEMM as a compute backend: i) prefill of LLM inference with speedups up to 1.85 $\times$ over SOTA inference runtimes, and ii) distributed-memory matrix multiplication with speedups up to 2.3 $\times$ over the SOTA distributed-memory GEMM framework with vendor-optimized compute backend.

Figures & tables

Fig. 1: BF16 GEMM performance on 128-core GNR: SFC-CA GEMM (this work), oneDNN, CARMA and OpenBLAS. The x-axis shows

M\times N\times K

dimensions of GEMM.

Fig. 2: Gilbert traversal and recursion [ 13 ] . (a) An

8\times 8

traversal, with endpoints 0 and 63. (b) Standard rule on a square: matching colors identify the four consecutive quadrants and their index ranges in (a), and the dashed cut separates the two half-square runs, 0–31 and 32–63. (c) Long rule: visit the two halves consecutively. The major side

w

is vertical and

h

horizontal.

Fig. 3: SFC-based partitioning of a 4096 $\times$ 4096 $C$ matrix (128 $\times$ 128 blocks of size 32 $\times$ 32) among 64 cores, one color per core. Left: 2D ( $K\_layers=1$ ). Middle: 2.5D ( $K\_layers=2$ ). Right: 3D ( $K\_layers=4$ ).

Matrix shape / regime	Dimension/thread conditions and layer choice	SFC-CA GEMM words $W$ Matched lower-bound type
Square: 2D / 2.5D $M=N=K=n$ (Cor. 3 (i))	$1\leq c\leq T^{1/3}$ , $bk=b$ , and $FastMem=\Theta(cn^{2}/T)$ , in addition to the shared cache-fit condition. $c=1$ gives 2D.	$\Theta\!\big(n^{2}/\sqrt{cT}\big)$ Memory-dependent
1D output bands Any $K$ (Cor. 3 (ii))	$1\leq T\leq M/N$ with $c=1$ . Each core owns $M/T$ full rows of $C$ .	$\Theta\!\big((MK+MN)/T\big)$ Compulsory input/output
Short inner dimension $K\leq N\leq M$ (Cor. 3 (iii))	$M/N\leq T\leq MN/K^{2}$ with $c=1$ . Writing $C$ dominates the total traffic.	$\Theta(MN/T)$ Compulsory output
Large inner dimension $K>N$ , $M\geq N$ (Cor. 3 (iv))	$M/N\leq T\leq MK/N^{2}$ . The exact choice $c=\min(K/N,T)$ must be feasible.	$\Theta(MK/T)$ Compulsory input $A$
Three large dimensions Any ordering of $K$ (Cor. 3 (v))	$T\geq\max(MN/K^{2},MK/N^{2})$ . Assume $bk=b$ and $T,K_{b}$ are powers of two. Choose $c=2^{\lfloor\log_{2}c_{*}\rfloor}$ and $f=K_{b}/c$ .	$\Theta\!\big((MNK/T)^{2/3}\big)$ Memory-independent

TABLE I: Proved asymptotic communication-optimal regimes and coverage of the dimension/thread ranges.

Fig. 4: SFC-based partitioning of $C$ matrices with different $M$ : $N$ aspect ratios using 64 cores. The SFC-based partitioning yields 2D core decompositions with aspect ratios matching the corresponding $M$ : $N$ ratio.

Fig. 5: EMR BF16 GEMM performance over all 125 shapes: SFC-CA GEMM, oneDNN, CARMA, OpenBLAS and the roofline. The $M\times N\times K$ configurations are sorted by operational intensity.

Fig. 6: Total L2 misses on EMR for two problem sizes.

Fig. 7: Performance for 100 random GEMMs on EMR.

Panel of Figure 9 — Fig. 8: BF16 GEMM over 125 $M\times N\times K$ shapes per platform, with a shared configuration order and x-axis labels shown only at the bottom. Top: GNR, comparing SFC-CA GEMM, oneDNN, CARMA, OpenBLAS and the roofline. Middle: ZEN5 and Bottom: GVT4, comparing SFC-CA GEMM, oneDNN (ACL through oneDNN on GVT4) and the roofline.

Fig. 9: LLM inference prefill of BF16 Llama-3 8B for various input lengths and batch sizes: i) PyTorch-TPP+SFC-CA GEMM (blue), ii) PyTorch-TPP+PARLOOPER GEMMs (green), iii) stock PyTorch+oneDNN (orange). Top: EMR, Bottom: GNR.

Fig. 10: Strong scaling of BF16 32k

\times

32k

\times

32k GEMM in COSMA with SFC-CA GEMM (blue) and oneDNN (red), on 1–16 EMR nodes (2 MPI ranks per node). Dark segments show compute, and light segments show communication plus other overhead. Boxed labels give aggregate compute-only TFLOP/s, unboxed labels effective TFLOP/s, and green labels wall-time speedups.

Explore similar work

Self-Supervised Learning for Sparse Matrix Reordering
May 17, 2026Ziwei Li, Tao Yuan, Fangfang Liu +3Sparsity Graph Edge Sparsification
Institute of Software, Chinese Academy of Sciences, Beijing, China
CUDA-L2: Surpassing cuBLAS Performance for Matrix Multiplication through Reinforcement Learning
Dec 2, 2025Songqiao Su, Xiaoya Li, Albert Wang +3Cuda Matrix Multiplication
DeepReinforce Team
FP8 is All You Need (Part 1): Debunking Hardware FP64 as the HPC Holy Grail (Sep 3rd version)
Date pendingSatoshi MatsuokaHigh-Performance Computing Fp8
Director, RIKEN Center for Computational Science (R-CCS) Kobe, Hyogo, Japan

May 17, 2026, cs.LGJ/K move · Enter open · S save

Self-Supervised Learning for Sparse Matrix Reordering

Ziwei Li, Tao Yuan, Fangfang Liu, Shuzi Niu+2

Institute of Software, Chinese Academy of Sciences, Beijing, China

Sparsity Graph Edge Sparsification Loss Function Self-Supervised Learning

Open paper AskPDF arXiv

Rearranging the rows or columns of a sparse matrix using an appropriate ordering can significantly reduce fill-ins, i.e., new nonzeros introduced during matrix factorization, decreasing memory usage and runtime. However, finding an ordering that minimizes fill-ins is NP-complete. Existing approaches, including graph-theoretic and deep learning methods, rely on surrogate objectives without theoretical guarantees. The Fill-Path Theorem reveals a direct and intrinsic relationship between fill-in generation and the sparse structure of the matrix as path triplet inequalities. Here we first employ a multigrid graph network to capture structural information for each vertex. We then derive a triplet sampling strategy based on inequalities. Finally, we introduce an end-max chain loss function to reduce the number of triplets whose predicted scores satisfy these inequalities. Experimental evaluations on the publicly available SuiteSparse matrix collection demonstrate the superiority of the proposed method in terms of both fill-in reduction and speedup in LU factorization time.

Abstract

Explore similar work

Self-Supervised Learning for Sparse Matrix Reordering

CUDA-L2: Surpassing cuBLAS Performance for Matrix Multiplication through Reinforcement Learning

FP8 is All You Need (Part 1): Debunking Hardware FP64 as the HPC Holy Grail (Sep 3rd version)