Function gemm_f64 

pub fn gemm_f64(
    m: usize,
    k: usize,
    n: usize,
    a: &[f64],
    b: &[f64],
) -> Result<Vec<f64>, ForgeError>

Best-path double-precision dense GEMM: row-major C[M×N] = A[M×K] · B[K×N], all f64.

The 3-tier f64 chain (“best f64 path on every machine”)

WGSL has no native f64, so double precision on the GPU is reached two different ways depending on the hardware; this is the whole reason gemm_f64 resolves through three tiers rather than the single accelerator arm of gemm_f32:

tier	path	when
1	native CUDA-f64 ([gemm_f64_cuda], NVIDIA only)	caps().cuda and ≥ GEMM_GPU_THRESHOLD FMAs
2	df64 / double-single WGSL (gemm_f64_df64, any other GPU)	caps().wgpu and ≥ GEMM_GPU_THRESHOLD FMAs
3	CPU floor (gemm_cpu_f64)	otherwise, or if every eligible accelerator errors

Tier 1 is exact double (native double + fma.rn.f64). Tier 2 emulates each f64 as a hi/lo pair of f32 with error-free transforms (~44–48 effective mantissa bits, well beyond a single f32’s 24) — so a non-NVIDIA GPU (AMD, Intel, Apple, mobile) now gets real f64 acceleration instead of dropping straight to the CPU. On any accelerator runtime error the call falls through to the next tier (errors are never propagated), so it is never broken. The CUDA arm is compiled in only under the cuda feature; the df64 arm is always present (it needs only a wgpu adapter).

a must have m * k elements, b must have k * n; both row-major. Returns m * n row-major elements.