Function std::mps_matmul_f16

int mps_matmul_f16(
    const &float a,
    const &float b,
    &float out,
    unsigned long M,
    unsigned long K,
    unsigned long N,
)

Real half/bfloat compute on the GPU, not just storage: gpu_mps_kernels.metal's matmul_kernel_f16/matmul_kernel_bf16/relu_kernel_f16/relu_kernel_bf16 hold/multiply in the reduced-precision type and accumulate in float (standard mixed-precision-GEMM practice) — verified to compile through this machine's real Metal toolchain (both 'half' and 'bfloat' are natively supported MSL types on this SDK). The functions below take and return ordinary 'float*' host arrays, same shape as mps_matmul_f32/ mps_relu_f32, so a caller doesn't need to touch fp16/bf16 bit patterns directly: conversion (see float16.h's f32_to_fp16_bulk/f32_to_bf16_bulk and their inverses) happens on the CPU side of the upload/readback, invisibly, inside these functions.

This is the naive (TILE_SIZE=16, threadgroup-memory-tiled) kernel tier only — the smma/multi-simdgroup tiers (see mps_matmul_f32's comment) use simdgroup_float8x8, a float-only type on this Metal SDK; porting those to half/bfloat accumulation is future work, not attempted here. Every call here does its own upload/dispatch/readback (no persistent- buffer or batching integration) — the simplest, most directly testable shape, matching mps_matmul_f32's non-batched fallback path.

Expected precision loss is real and by design: fp16's ~3-4 significant decimal digits (10 mantissa bits) and bf16's ~2-3 (7 mantissa bits) are both coarser than float32's ~7 — see float16.h's precision notes. Returns 1 on success, 0 on failure (no Metal device, pipeline compilation failure, or allocation failure) — 'out' is left untouched on failure.