Accelerating Multiple-Precision LU Decomposition with Ozaki Scheme II

1 Introduction

• Every scalar multiply-add in the O(n^3) update remains multiple precision in standard approaches, leaving low-precision matrix-engine throughput unused.The cited approaches include GMP/MPFR and multi-component arithmetic.Source: §1 · p. 2

1.1 Prior work: Strassen-type algorithms and their limits

• Strassen reduces operation count, but expensive additions, recursive-copy traffic, and errors of hundreds to thousands of ulps limit its appeal at arbitrary precision.The reported stability penalty is relative to a naive implementation.Source: §1.1 · pp. 2–3
• Ozaki methods replace multiple-precision operations with low-precision GEMMs; scheme II uses O(p) products, while modular reduction and CRT cost O(p^2).The scheme II advantage in product count applies to the GEMM term; which cost dominates depends on n and p.Source: §1.1 · p. 3

1.2 Contributions

• The implementation adds blocked, partially pivoted LU for multi-component and arbitrary precision, compares CPU and GPU at matched precision and size, and tests the ill-conditioned Lotkin matrix.For LU, the cubic Schur-complement update is replaced while the other listed operations remain ordinary multiple-precision arithmetic.Source: §1.2 · pp. 3–4
• New FP16, FP8, and binary64 GPU back-ends broaden the choice; INT8 is fastest on Arm/GB10, while binary64 is fastest for QD on x86/H100.FP8 uses balanced base-17 two-digit decomposition to retain INT8's 362.8-bit CRT capacity; the recommended choice maximizes bits per modulus times engine throughput.Source: §1.2 · p. 4

2.1 dtq (CPU, multi-component precision)

• dtq provides DD, TD, and QD as non-overlapping expansions of two, three, or four binary64 components, which can be passed directly as arrays to the C implementation.The implementation converts these expansions directly to Ozaki scheme II's fixed-point representation.Source: §2.1 · p. 4

2.2 MPFR (CPU, arbitrary precision)

• MPFR provides run-time variable precision with correctly rounded operations and serves as the CPU input and output type for Ozaki GEMM and LU.Its adjustable significand supports precision selection according to the condition number.Source: §2.2 · p. 4

2.3 gdtq (GPU, multi-component precision)

• gdtq provides GPU DD, TD, and QD types with the same memory layout as dtq, and the Ozaki path handles their components directly as integers.This path avoids calling gdtq arithmetic, whose multiplication is device-only.Source: §2.3 · p. 4

2.4 MPC_CUDA (GPU, fixed precision)

• MPC_CUDA's cu_freal is a register-resident GPU fixed-precision type that rounds compatibly with MPFR's RNDN; the work implements and tests its MPFR conversions.A bit-exact round trip was verified over 2000 trials for PB values 128, 256, 512, and 1024.Source: §2.4 · p. 5

3.1 Ozaki scheme II for GEMM

• Ozaki scheme II makes the integer products exact and rounds only once after CRT reconstruction, unlike naive multiple-precision GEMM, which rounds at every multiply-add.The single rounding occurs at the end of the entire inner product.Source: §3.1 · p. 5

3.2 Arbitrary precision: splitting into several slices and CRT

• Scheme II enables wider slices than scheme I by limiting chunk width through CRT capacity rather than accumulator exactness, reducing the surviving digit-pair count to a practical range.The single-pass limit is about 155 bits for K = 4096 and g = 8; TD/QD precision and above therefore require slicing. The slice count balances conversion/CRT cost against GEMM cost, with an interior model optimum of O(p^1.5) rather than O(p^2) at either endpoint.Source: §3.2 · pp. 5–7
• Under cancellation at p = 2048 and K = 256, S = 4 attains 2^-519 and S = 8 attains 2^-194; each extra group recovers about 272 bits here, and four restore exactness.For S = 8, one, two, and three additional groups yield 2^-468, 2^-738, and then four extra groups restore exactness.Source: §3.2 · p. 7

3.3 Direct conversion from multi-component precision

• Direct conversion maps non-overlapping DD/TD/QD expansions to fixed-point integers and back without MPFR, avoiding per-element heap allocation during Schur updates.At n = 1024 and b = 512, the output matrix alone contains 262,144 elements, and MPFR allocation costs exceed the GEMM cost.Source: §3.3 · p. 8

3.4 Incorporation into LU decomposition

• The blocked LU leaves panel factorization, pivot propagation, and forward elimination in multiple precision, replacing only the cubic-cost trailing-submatrix update with Ozaki scheme II GEMM.The modulus and CRT context is built once for each panel width and reused across steps.Source: §3.4 · p. 9

3.5 Implementation

• The implementation provides INT8, FP16, FP8, and binary64 backends because the best format depends on bits per modulus multiplied by engine throughput.INT8 is the default, while binary64 can benefit machines with strong FP64 units; the measured fastest backend varies by machine.Source: §3.5 · pp. 9–10

3.5.1 CPU implementation: mapping onto instruction sets

• CPU INT8 GEMM maps to AMX-INT8 on x86-64 and SVE2 i8mm on aarch64, while the binary64 path uses OpenBLAS cblas_dgemm.The batched AVX-512 CRT is x86-only; aarch64 falls back to scalar CRT using GMP's mpn_addmul_1.Source: §3.5.1 · p. 11
• At n = 1024, MPFR arbitrary-precision GEMM speedup is 14.5–27.8× with AMX-INT8 versus 6.9–14.2× with SVE2 i8mm, a 1.6–1.7-times gap (Table 3).The algorithm, parameters, and comparison targets are shared; only the realization of (S4)/(S5) changes.Source: §3.5.1 · p. 11

3.5.2 GPU implementation

• On the GPU, device kernels implement (S1)–(S6), with cuBLASLt GEMM back-ends selectable as INT8, FP16, FP8, or binary64.INT8 is the default; CRT reconstruction obtains the centered representative by round-to-nearest, unlike the CPU procedure.Source: §3.5.2 · p. 11
• The GPU DD/TD/QD path passes contiguous double2/double3/double4 components directly to device kernels and uses expansion-native kernels, avoiding cu_freal and MPFR.When multiple digit groups contribute, their sum is normalized in the expansion format by Priest-style distillation using sweeps of two-sum.Source: §3.5.2 · p. 12
• In LU, one thread per element parallelizes (P1)–(P3), while Ozaki scheme II GEMM replaces only (P4); a blocked naive multiple-precision inner product isolates the update method.The blocked naive version shares the panel and pivoting code with the Ozaki version.Source: §3.5.2 · p. 12

• Arm/GB10 runs used 20 OpenMP threads and SVE2 i8mm, while x86/H100 used 32 threads and Intel AMX-INT8.Source: §4.1 · p. 12, Table 2 · p. 12
• The comparisons target naive dtq/gdtq GEMM and LU, and naive MPFR/cu_freal implementations for arbitrary precision.Source: §4.1 · p. 12
• The earlier arbitrary-precision evaluation concluded that Ozaki scheme I could not accelerate LU decomposition sufficiently, so it is not repeated here.Source: §4.1 · p. 12
• At p = 8192, the binary64 Ozaki scheme II estimate matches the measured count of 716 moduli.Source: §4.1 · pp. 12–13
• At p = 3136–12544, Ozaki scheme I would need 37–164 times as many products for the Lotkin-matrix LU result.Source: §4.1 · p. 13
• Published fast Ozaki scheme I implementations, including ozIMMU, target binary64 emulation rather than DD/TD or several-thousand-bit precision.Source: §4.1 · p. 13
• For DD/TD-equivalent precision, the two schemes’ product counts are comparable at 1.3–2.1×, so the paper does not dismiss Ozaki scheme I in that regime.Source: §4.1 · p. 13

4.2.1 CPU: DD/TD/QD and MPFR arbitrary precision

• At n = 1024, scheme II reaches 2.67/5.37/8.74× speedups on Arm/GB10 and 1.81/5.56/10.52× on x86/H100 for DD/TD/QD.Its relative errors are five to seven orders of magnitude smaller than naive dtq; the speed advantage grows with precision.Source: Table 3 · p. 13, §4.2.1 · p. 14
• The MPFR-mediated conversion path takes 2.4–4.7× as long as direct conversion at n = 1024, reducing or reversing some speedups over naive dtq.Relative errors agree between the paths, so the measured difference is conversion cost.Source: §4.2.1 · p. 14

4.2.2 GPU: cu_freal (MPC_CUDA)

• On GPUs, INT8 is 40–66× faster than naive GEMM, while FP4 is 3–8× slower than INT8 within scheme II.FP4's raw dense GEMM is faster, but its modulus capacity is 168.1 bits versus INT8's 362.8 bits and it needs four sub-GEMMs per modulus.Source: §4.2.2 · p. 14, Table 4 · p. 14

4.2.3 Implementation of the FP16 and FP8 back-ends

• Balanced base-17 FP8 encoding enables the full 362.8-bit CRT capacity of INT8; FP8 matches INT8's correct-bit count for DD/TD/QD and runs within about 1.35–1.44× of its speed.On Arm/GB10, FP16 reaches 84.7 TFLOP/s versus 593.0 TFLOP/s on x86/H100.Source: §4.2.3 · pp. 16–17

4.2.4 Why we do not adopt FP32 or TF32

• FP32 offers only 0.5 bit per modulus more than INT8 but has over 20 times lower x86/H100 throughput, while TF32 cannot meet scheme II's exact-GEMM requirement.FP32 has the same modulus limit as FP16 and is 15 times slower; the back-end choice is framed as maximizing bits per modulus times engine throughput.Source: §4.2.4 · p. 17

4.2.5 Precision ceiling of the INT8 back-end

• For arbitrary-precision LU, INT8 is faster on Arm/GB10, whereas binary64 is faster on x86/H100; the GPU difference is larger because conversion is only 3–5% of GPU time.The CPU INT8/binary64 ratio is 0.72–0.88 on Arm/GB10 and 1.24–2.19 on x86/H100; the GPU ratio is 0.07–0.09 and 1.73–1.88, respectively.Source: §4.2.5 · pp. 18–19, Table 8 · p. 18
• With unsplit INT8, doubling precision multiplies GEMM time by 1.80–2.04 but conversion time by 2.87–3.61; optimal slice splitting brings both toward p^1.5 growth.At n = 512 and p = 8192, splitting reduces runtime from 13.32 s to 3.89 s; it also makes p = 16384 feasible where unsplit INT8 exceeds capacity.Source: §4.2.5 · pp. 19–20
• The two capacity routes trade more GEMMs with lighter conversion against fewer GEMMs with heavier conversion, and slice splitting automates the balance through its chosen slice count.The count of low-precision GEMMs grows as O(p) for modulus-table extension and O(p^2) for slice splitting.Source: §4.2.5 · p. 20

4.3.1 CPU: DD/TD/QD

• On x86 CPUs, Ozaki scheme II wins over the best existing implementation at every precision, while on Arm/GB10 it loses at DD (0.78×) and only narrowly wins at TD/QD.The speedup ratio increases with precision on both machines.Source: §4.3.1 · p. 21
• At TD/QD, unblocked native dtq is faster than blocked dtq: 0.277 s versus 0.361 s at TD and 0.518 s versus 0.594 s at QD.At n = 512 on Arm/GB10 in QD, unblocked dtq takes 0.070 s versus 0.144 s for Ozaki scheme II-dtq.Source: §4.3.1 · p. 24
• The direct conversion path shortens the Schur update from 0.056 s to 0.024 s at DD and from 0.078 s to 0.042 s at QD, while avoiding per-element MPFR heap allocation.The MPFR-mediated version is 2.7–15 times slower than dtq in total time.Source: §4.3.1 · p. 24

4.3.2 GPU: gdtq (multi-component)

• On Arm/GB10, Ozaki scheme II/INT8 is 7.4/17.9/29.4× faster in the Schur update and 1.20/1.36/1.45× faster in total; on x86/H100, gdtq is faster at DD/TD and comparable at QD.The x86/H100 ratios are 0.34× for DD, 0.71× for TD, and 1.06× for QD.Source: §4.3.2 · p. 24
• At QD on x86/H100, Ozaki scheme II/FP64 takes 1.83 ms for the Schur update, beating gdtq native at 2.84 ms and INT8 at 2.78 ms; total times are 0.0553 s and 0.0562 s, respectively.At TD, FP64’s 1.46 ms is 1.45× faster than INT8’s 2.12 ms and essentially matches gdtq at 1.42 ms; native gdtq is fastest only at DD.Source: §4.3.2 · p. 26
• FP8 is slowest even on x86/H100, where its raw GEMM reaches 1525.7 TFLOP/s, because its two-digit encoding requires four sub-GEMMs per modulus.At DD on x86/H100, INT8/FP16/FP8 updates take 1.41/1.56/2.14 ms; at QD they take 2.68/2.64/4.16 ms. At DD/TD/QD, the back-ends yield the same correct bits as INT8.Source: §4.3.2 · p. 26

4.4 LU decomposition of the Lotkin matrix and the relative error of the numerical solution

• At n = 2048, Ozaki scheme II reaches a maximum speedup of 2.84× over the faster naive variant; speedup increases monotonically with n in all four environments.Across sizes, speedups are 1.34→1.80→2.18× on the Arm/GB10 CPU, 1.91→2.02→2.64× on the x86/H100 CPU, 1.26→2.65→2.84× on the Arm/GB10 GPU, and 1.01→1.63→2.06× on the x86/H100 GPU.Source: Table 14 · p. 27, §4.4 · p. 27
• At n = 2048, the relative error reaches 10^-645, and CPU and GPU correct-bit counts agree within 1–3 bits.The methods agree with one another, so the speedup costs no accuracy.Source: §4.4 · p. 27
• At n = 2048, the Arm/GB10 GPU takes 17604.1 s versus 1348.4 s on its CPU, while x86/H100 GPU and CPU times are 919.3 s and 734.5 s.The reported Arm/GB10 GPU time is 13.1 times its CPU time; the passage attributes this to weak binary64 capability affecting panel factorization.Source: §4.4 · p. 29
• Blocking alone changes GPU naive total times by only −4.2% to +5.8%, while Ozaki scheme II’s speedups over the faster naive variant range from 1.26–2.84× on Arm/GB10 and 1.01–2.06× on x86/H100.The passage attributes the small blocking effect to locality gains being roughly cancelled by panel-width overhead.Source: §4.4 · p. 29
• Panel factorization remains the common bottleneck, taking 77% of CPU Lotkin time and 96% of GPU time at n = 1024/p = 8192; at n = 2048 on the Arm/GB10 GPU it accounts for 98%.At that Arm/GB10 GPU size, the Schur update shrinks 49.4×, but the overall speedup remains 2.84×.Source: §4.4 · p. 29

4.5 Summary: profiling the LU decomposition

• Ozaki scheme II accelerates the Schur update by 1.50–3.87× on CPU and 7.3–29.5× on GPU, while total-time improvement is only 0.78–1.11× on CPU and 1.20–1.45× on GPU.The panel factorization grows to 80–98% of Ozaki scheme II total time, leaving the update no longer as the bottleneck.Source: §4.5 · pp. 30–31, Figure 2 · p. 30, Figure 3 · p. 30
• The unblocked b = 1 bars have an almost invisible panel, indicating that blocking cost is not negligible in multiple precision.The figures separate panel factorization and Schur update; their vertical scales differ between panels.Source: Figure 2 · p. 30, Figure 3 · p. 30, §4.5 · p. 31
• At QD, the Schur update is 1.6% of total GPU time (0.0063/0.4020), and for MPFR it is 0.7% (3.7/500.8).The profile identifies panel factorization, rather than further Schur-update optimization, as the next target.Source: §4.5 · p. 31
• For GPU MPFR at n = 1024 and p = 8192, reducing b from 512 to 256 lowers total time from 500.8 s to 300.4 s as panel time falls from 497.1 s to 288.4 s.The preferred panel width balances panel cost, which scales as O(nb^2), against Schur-update preprocessing and CRT overhead.Source: §4.5 · p. 31

5.1 Conclusions

• At n=1024, CPU Arm/GB10 matrix multiplication sped up 3.69/7.01/12.42× for DD/TD/QD, while Schur updates sped up 1.59/2.80/4.12×.Schur-update figures use the same multi-component types directly, without cu_freal or MPFR.Source: §5.1 · p. 31
• Larger size and higher precision consistently favor Ozaki scheme II, but total-LU gains vary by machine: 1.04–1.06× on Arm/GB10 CPU versus 1.27–2.36× on x86/H100 CPU at n=1024.GPU total-LU results range from 1.20–1.45× on Arm/GB10 to 0.34–1.07× on x86/H100; after acceleration, panel factorisation accounts for 80–98% of total time.Source: §5.1 · pp. 31–32
• The fastest back-end depends on the machine: arbitrary-precision Arm/GB10 favors INT8 at every tested precision, whereas x86/H100 favors binary64.The scheme's design principle is to maximize bits per modulus multiplied by engine throughput, rather than always choosing a low-precision engine.Source: §5.1 · p. 32
• FP8 matches INT8's 362.8-bit CRT capacity through balanced base-17 two-digit encoding, but is slowest within Ozaki scheme II because its two-digit decomposition requires four times as many sub-GEMMs.Even FP8's 1.94× raw GEMM throughput advantage on x86/H100 does not recover that sub-GEMM count.Source: §5.1 · p. 32
• At n=2048, arbitrary-precision Lotkin LU reached relative error 10^-645 and was 2.18× faster on Arm/GB10 CPU, 2.64× on x86/H100 CPU, 2.84× on Arm/GB10 GPU, and 2.06× on x86/H100 GPU.The tested precision was p ≈ 1.2 log2 cond; speedups compare with existing implementations.Source: §5.1 · p. 32
• At n=512, CPU slice splitting reduced runtime from 13.32 s to 3.89 s (3.4×) at p=8192 and from 31.28 s to 10.19 s (3.1×) at p=16384.More slices trade additional GEMMs for lower conversion cost; truncating low-order digit groups means severe cancellation requires a guard sized to Γ.Source: §5.1 · p. 32, §5.2 · p. 33

5.2 Future work: choosing the Ozaki scheme II implementation to match the number of digits

• Future work calls for precision- and size-aware implementation selection, faster panel factorisation, larger-n GPU evaluation, and balancing moduli against panel width.Table 15 gives a measurement-based implementation guideline within this work's scope.Source: §5.2 · p. 33, Table 15 · p. 33