Tailored DGEMM Micro-Kernel for Phytium FT-2000/4

A hand-written NEON assembly micro-kernel for double-precision general matrix multiplication (DGEMM), achieving 89% of OpenBLAS single-thread throughput and 90% of OpenBLAS 4-thread throughput on the Phytium FT-2000/4.

This is not a general-purpose BLAS replacement. It is a case study in targeted optimization — identifying the exact microarchitectural characteristics of a specific CPU and designing instruction scheduling around them.

Background

The FT-2000/4 uses four FTC663 cores (ARMv8-A, 2.6 GHz, 16 nm). Each core has a 4-wide out-of-order pipeline with two FP64 FMA units, giving a theoretical peak of 10.4 GFLOPS per core (2 FMA/cycle x 2 FLOP/FMA x 2.6 GHz).

OpenBLAS 0.3.34 ships an optimized DGEMM for this platform and achieves 9.3 GFLOPS single-thread (89% of peak). This project started with a simple question: how close can a from-scratch kernel get?

The answer: 8.3 GFLOPS (80% of peak, 89% of OpenBLAS) — arrived at through 22 iterations of measurement, hypothesis, and correction.

Results

All benchmarks on FT-2000/4 at 2.6 GHz, Ubuntu Linux, GCC 11, -O3 -mcpu=cortex-a72.

Single-thread (512x512x512)

Implementation          GFLOPS   % of Peak   % of OpenBLAS
Naive C (gcc -O3)        3.5       34%          38%
This kernel              8.3       80%          89%
OpenBLAS 0.3.34          9.3       89%         100%
Theoretical peak        10.4      100%           -

4-thread (2000x2000x2000)

Implementation          BEST       AVG        % of OpenBLAS
This kernel             31.1       30.4          90%
OpenBLAS 0.3.34         34.7       34.5         100%

Parallel scaling: 3.70x on 4 cores (92.5% efficiency).

Correctness

Verified against exact reference computation for matrix sizes 64, 128, 257, 512, 1024, 2000. Max relative error below 1e-12 for all cases.

Design

Micro-Kernel: 8x4 Tile

Registers (32 NEON):
  v0-v3:   B data (double-buffered: step N in v0-v1, step N+1 in v2-v3)
  v8-v15:  A data (double-buffered: step N in v8-v11, step N+1 in v12-v15)
  v16-v31: 16 accumulators (8 rows x 2 vectors = 8x4 tile)

Per k-step operations:
  3 x ldp                  (load 8 FP64 of A + 4 FP64 of B)
  16 x fmla .d[0]/.d[1]   (indexed fused multiply-add)
  2 x prfm                 (prefetch A and B)
  Total: 21 instructions per k-step

Key Optimization: FFL Instruction Scheduling

The central finding of this project is that instruction interleaving pattern matters more than instruction count on the FTC663 pipeline.

The kernel interleaves loads and FMAs in a repeating pattern of L(oad) F(mla) F(mla) groups:

Phase 1 (compute step N, load step N+1):
  ldp | fmla fmla | ldp | fmla fmla | ldp | fmla fmla | prfm | fmla fmla ...
  L       FF        L       FF        L       FF        P      FF

Phase 2 (compute step N+1, load step N+2):
  ldp | fmla fmla | ldp | fmla fmla | ldp | fmla fmla | prfm | fmla fmla ...

This fills all 3 issue slots per cycle (2 FMA + 1 load/prefetch) on a 3-wide front-end, and doubles the effective load-to-use distance compared to grouping loads separately.

Other Techniques

Technique                           Effect
ldp post-update addressing          1 insn loads 2 FP64 + updates pointer
Indexed fmla (.d[0] / .d[1])        Single-instruction broadcast-MAC
Double-buffered registers           Compute step N while loading step N+1
Direct C accumulate                 No staging buffer needed
2x k-unroll                         Halves loop branch overhead
Packing (A: 8xKC, B: KCx4)          Converts strided to sequential access
schedule(static)                    Stable multi-thread performance

What Did NOT Work

This is equally important. The following hypotheses were tested and rejected by measurement:

Hypothesis                                     Result
.d[0] indexed FMA faster than .d[1]            No measurable difference
Grouping all loads first                       12% slower (FMA dependency stalls)
Copying OpenBLAS register allocation           5% slower
Shared B-panel packing across threads          5% slower (sync overhead)
4x k-unroll for multi-thread                   8% slower (I-cache pressure)
fmul start (skip accumulator zeroing)          Incompatible with direct C (bug)
KC=256 matching OpenBLAS parameter             3% slower (L2 thrashing)

Repository Structure

src/ft2000_dgemm.h       Public API header
src/ft2000_dgemm.c       Complete implementation (kernel + packing + driver)
bench/bench_dgemm.c      Benchmark harness (vs OpenBLAS, with verification)
docs/optimization_log.md Full optimization history (22 iterations)

Building

gcc -O3 -mcpu=cortex-a72 -ffast-math -funroll-loops \
    -fopenmp -o bench_dgemm bench/bench_dgemm.c \
    src/ft2000_dgemm.c -lopenblas -lm

The kernel itself does not depend on OpenBLAS. It is only linked in the benchmark binary for comparison.

Usage

#include "ft2000_dgemm.h"

// C = A * B  (alpha=1, beta=0)
ft2000_dgemm(m, n, k, A, lda, B, ldb, C, ldc);

Tunable parameters (compile-time): MC=256, NC=256, KC=128, MR=8, NR=4. These were determined empirically for the FT-2000/4 cache hierarchy.

Limitations

Optimized specifically for FT-2000/4 (FTC663 core). On other ARM cores, the FFL scheduling pattern may not be optimal.
Single micro-kernel shape (8x4). No runtime dispatch to other shapes.
Alpha=1, beta=0 is the fast path. General alpha/beta adds overhead.
Single-thread gap (11% vs OpenBLAS) is attributed to pipeline port contention specific to the FTC663 FMA issue logic, which would require cycle-level tuning to close.

Optimization Methodology

Each iteration followed a strict protocol:

Measure — 11-run benchmark with best/avg statistics
Hypothesize — identify one suspected bottleneck
Change — modify exactly one thing
Verify — correctness check + performance comparison
Record — document the result (improvement, regression, or no change)

No iteration changed more than one variable. Regressions were reverted immediately. The complete log is in docs/optimization_log.md.

License

MIT