Ozaki-Scheme Low-Precision GEMM¶

Intercepts DGEMM, SGEMM, ZGEMM, and CGEMM at link time, replacing them with low-precision Ozaki schemes (mantissa slicing or CRT). Real GEMM calls go through the selected scheme; complex GEMM (ZGEMM/CGEMM) is mapped to real GEMM internally.

Build¶

cd samples/ozaki
make -j $(nproc)

LIBXS must be built first from the repository root. When a sibling libxstream directory is detected, the optional OpenCL/GPU path is compiled in (see OZAKI_OCL below).

Link-Time Interception¶

Two link-time variants are built per precision:

• dgemm-blas.x / sgemm-blas.x — dynamically linked against BLAS
• dgemm-wrap.x / sgemm-wrap.x — linked with --wrap= (GNU ld)
• zgemm-wrap.x / cgemm-wrap.x — complex GEMM wrappers

Static (libwrap.a) and shared (libwrap.so) wrapper libraries are built by default. The shared library can intercept any application:

LD_PRELOAD=/path/to/libwrap.so ./application

The static library requires explicit --wrap flags:

-Wl,--wrap=dgemm_ -Wl,--wrap=sgemm_ -Wl,--wrap=zgemm_ -Wl,--wrap=cgemm_

Test scripts: test-wrap.sh exercises all variants, test-check.sh validates both schemes, test-mhd.sh tests MHD file-based input.

Test Driver¶

dgemm-wrap.x [A.mhd|M [B.mhd|N [K [TA [TB [ALPHA [BETA [LDA [LDB [LDC]]]]]]]]]]

Defaults: M=N=K=257, TA=TB=0, ALPHA=BETA=1.0. The first argument can be 0 followed by MHD filenames to load A and B matrices from files. The zgemm-wrap.x and cgemm-wrap.x drivers call ZGEMM and CGEMM.

Environment Variables¶

Scheme Selection¶

OZAKI=3 (adaptive) starts with Scheme 1 on the first GPU call to learn the effective cutoff from preprocessing occupancy. Subsequent calls compare the Scheme-1 pair count (at the cached cutoff) against the Scheme-2 prime count and pick the cheaper path. On CPU, adaptive falls back to Scheme 2.

Unset means CRT on the CPU and adaptive on the GPU: the GPU default is LIBXSTREAM's, so that it stays the same whichever driver asks for it. Setting OZAKI applies to both.

Accuracy¶

GPU Path (requires LIBXSTREAM)¶

GPU-specific kernel tuning variables (OZAKI_RTM, OZAKI_RTN, OZAKI_WG, OZAKI_SG, OZAKI_KU, OZAKI_RC, OZAKI_PB, OZAKI_HIER, OZAKI_PREFETCH, OZAKI_SCALAR_ACC, OZAKI_CACHE, OZAKI_ARENA) are documented in the LIBXSTREAM Ozaki README.

Monitoring and Diagnostics¶

Benchmark¶

GEMM_HOSTMEM selects how the operands are allocated, which is what an application would do differently rather than a property of the scheme: the default is plain malloc, declared to LIBXSTREAM so that transfers can still take the faster route, whereas a non-zero value allocates page-locked memory through LIBXSTREAM. On a PCIe part the difference is worth several times the transfer time, so a comparison against another library has to state which of the two it used. The default keeps a CPU-only run off the GPU entirely: no device is opened unless the GPU path is enabled (OZAKI_OCL=1), while GEMM_HOSTMEM needs a device to allocate from and therefore brings one up.

Adversarial Input (EVIL)¶

The EVIL variable initializes A and B with controlled exponent structure for stress-testing accuracy and adaptive slice reduction. The magnitude sets the exponent span in bits; the sign selects the distribution:

EVIL=N (N>0) Per-column. Column j of A is scaled by 2^(N*j/(ncols-1)), column j of B by the inverse. Product A*B is well-conditioned. Uniform exponents within each column.

EVIL=-N (N>0) Per-element. Each element gets a pseudorandom exponent in [0,N] via coprime shuffle, with opposite sign for B. Every row of A spans the full exponent range — worst case for row-wise alignment and adaptive cutoff.

EVIL=0 Default shuffle mode (no exponent structure).

The per-column mode (EVIL>0) matches the NVIDIA emulation grading test (diagonal scaling with D and D^-1). The per-element mode (EVIL<0) is adversarial for the adaptive slice-pair reduction: it forces all slices to be populated in every row.

Decay Diagnostic and K-Permutation (OZAKI_DECAY)¶

Controls K-dimension row reordering for smoothness and the forward-difference decay diagnostic. Values map to libxs_sort_t:

Values 2-4 compute a K-permutation from B (via libxs_sort_smooth) before Ozaki-1 slicing. Both A and B are sliced using the permuted K-order. Since C = A*B = (A*P^T)*(P*B) for any permutation matrix P, the result is mathematically identical — only the int8 digit structure along K changes.

When nonzero, reports the forward-difference decay of int8 slice buffers along K, M, and N axes (first K-group only, single- threaded). Uses libxs_fprint per-axis mode internally. Output goes to stderr:

OZ1[MxNxK] Delta-K: d1=... d2=... ... OZ1[MxNxK] Delta-M: d1=... d2=... ... OZ1[MxNxK] Delta-N: d1=... d2=... ...

Decaying values indicate exploitable smoothness; growing values (~2x per order) indicate unstructured data where data-independent schemes (Ozaki-1, Ozaki-2) are well-matched.

Profiling¶

Set LIBXSTREAM_PROFILE=1 for per-kernel device timings, reported at program exit:

PROF ACC/OpenCL: ID=171071 gemm_crt_fused=2.987 ms 5751.0 GFLOPS/s

There is no separate host-side profile: on the CPU the intercepted call and the kernel it spends its time in agree to within a few percent, so the OZAKI GEMM line already reports the kernel. It is a mean over the whole timed loop, so use NREPEAT to average several calls — host timings vary between runs, the more so on a many-core part.

Example¶

./dgemm-wrap.x 256                          # CRT scheme (default)
OZAKI=1 ./dgemm-wrap.x 256                  # mantissa slicing
OZAKI_TRIM=4 OZAKI=1 ./dgemm-wrap.x 256     # drop 4 least significant diagonals
OZAKI=2 OZAKI_GROUPS=4 ./dgemm-wrap.x 4096  # CRT with K-grouping
OZAKI=3 ./dgemm-wrap.x 4096                 # adaptive scheme selection
EVIL=512 ./dgemm-wrap.x 1024                # accuracy grading (wide exponent span)
EVIL=1 ./dgemm-wrap.x 1024                  # narrow span (cutoff drops pairs)
EVIL=-52 ./dgemm-wrap.x 1024                # per-element span (worst for cutoff)
OZAKI_DECAY=1 OZAKI=1 ./dgemm-wrap.x 256    # decay diagnostic (no reorder)
OZAKI_DECAY=2 OZAKI=1 ./dgemm-wrap.x 256    # sort by norm + decay diagnostic
OZAKI_DECAY=4 OZAKI=1 ./dgemm-wrap.x 256    # greedy NN + decay diagnostic