native cpu + gpu inference support

[Lazarus-931]

This began when I realized that cpu + gpu inference is possible with llama.cpp but insanely slow using mlx-lm on recent mlx main branch, all due to how slow it is for cpu.

This pr has several cpu paths to support cpu inference on par with llama.cpp and even cpu and gpu simultaneous inference is on par with llama.cpp.

For starters, this introduces a bf16 SIMD support to cpu, in neon_bf16_simd.h. This was because for models in this numerical percision(bf16), cpu didn't have native bf16 arithmetic, and prior to this, cpu was doing a scalar loop one element at a time, converting to float and back. So this does upcast to fp32, by widening 8 values to compute normal float SIMD, round back with same rounding as the scalar code, making it bit identical.

Secondly, while profiling quantized matmuls, I saw that against llama.cpp mlx uses single cores on cpu as against to 8, so parallel.h handles this, allowing for multiple workers to be launched(8 by default), created and reused for every matmul. If matmuls under 2m ops exist, they aren't routed to this and instead run on the calling thread.

Third, while profiling against llama.cpp, mlx unpacked 4-bit weights into floats and did float multiply-adds; llama.cpp rounds the activations to int8 and uses the CPU's SDOT instruction, so this PR also includes that support on mlx's own quantization format. Activations are quantized to int8 per 64-element group right before the matmul (with the scale kept on the side so results come out correct), and the 4-bit weights are recentered to signed values so the whole dot product runs on vdotq_s32 — 16 multiply-accumulates per instruction, about 3.5x over the float path. The rounding this introduces (~0.5%) is smaller than the error 4-bit weights already carry, and it's the same tradeoff llama.cpp makes on every CPU inference.

Lastly, two extensions to the int8 kernel. It now also covers 8-bit quantized weights, which matters because a bf16 checkpoint converted to 8-bit goes from 7.4 to 59 tok/s with near-lossless quality, and prefill used to re-read the entire weight matrix once per prompt token, now processing 4 prompt rows per weight pass.

All numbers are decode tokens/sec on a single M4 Mac mini (16GB, macOS), generating ~100 tokens from the same prompt per run, using mlx-community 4-bit checkpoints via mlx-lm (CPU-after runs with MLX_CPU_INT8_QMM=1, threading at the default MLX_CPU_THREADS=8) against llama.cpp build 9110 with matching Q4_K_M GGUFs (llama-bench -t 8 for CPU, -ngl 99 for GPU), and simultaneous rows run the same model concurrently on both devices in two separate processes.

Table 1 — CPU: before vs after this PR

model	CPU before	CPU after	speedup
gemma-3-270m	0.93	173.5	187×
Qwen2.5-0.5B	15.5	139.7	9×
Qwen3-0.6B	0.40	93.0	233×
Llama-3.2-1B	6.55	89.7	14×
gemma-3-1b	7.68	70.8	9×
LFM2-1.2B	0.20	94.5	472×
Qwen3-1.7B	0.14	56.0	400×
SmolLM3-3B	0.08	34.3	429×
Llama-3.2-3B	2.46	38.3	16×
Qwen2.5-3B	2.55	38.6	15×
Phi-4-mini	2.06	35.9	17×
Qwen3-4B	0.06	29.1	486×
gemma-3-4b	0.06	26.7	444×
Qwen2.5-7B	1.13	20.6	18×
Qwen2.5-14B	0.58	10.6	18×

Table 2 — simultaneous CPU+GPU: before vs after this PR

model	CPU+GPU before	CPU+GPU after
gemma-3-270m	1.0 + 432.7 = 433.7	155.3 + 354.1 = 509.4
Qwen2.5-0.5B	16.1 + 289.1 = 305.2	105.6 + 223.4 = 329.0
Qwen3-0.6B	0.4 + 248.4 = 248.8	76.5 + 194.7 = 271.2
Llama-3.2-1B	6.7 + 136.0 = 142.7	58.9 + 95.9 = 154.8
gemma-3-1b	8.0 + 150.5 = 158.5	59.2 + 113.2 = 172.4
LFM2-1.2B	0.2 + 151.6 = 151.8	63.6 + 102.7 = 166.3
Qwen3-1.7B	0.2 + 100.8 = 101.0	39.4 + 71.0 = 110.4
SmolLM3-3B	0.1 + 58.3 = 58.4	26.2 + 43.7 = 69.9
Llama-3.2-3B	2.6 + 54.7 = 57.3	24.1 + 36.7 = 60.8
Qwen2.5-3B	2.7 + 57.2 = 59.9	24.8 + 38.9 = 63.7
Phi-4-mini	2.1 + 44.7 = 46.8	21.2 + 29.8 = 51.0
Qwen3-4B	0.1 + 44.1 = 44.2	18.4 + 29.5 = 47.9
gemma-3-4b	0.1 + 45.1 = 45.2	18.3 + 31.2 = 49.5
Qwen2.5-7B	1.2 + 24.9 = 26.1	12.1 + 16.4 = 28.5
Qwen2.5-14B	0.6 + 12.4 = 13.0	1.3 + 12.1 = 13.4 ⚠

Table 3 — CPU: mlx (with PR) vs llama.cpp

model	mlx CPU	llama.cpp CPU
gemma-3-270m	173.5	229.8
Qwen2.5-0.5B	139.7	164.6
Qwen3-0.6B	93.0	137.2
Llama-3.2-1B	89.7	89.7
gemma-3-1b	70.8	90.4
LFM2-1.2B	94.5	99.5
Qwen3-1.7B	56.0	59.3
SmolLM3-3B	34.3	36.5
Llama-3.2-3B	38.3	35.3
Qwen2.5-3B	38.6	36.4
Phi-4-mini	35.9	29.4
Qwen3-4B	29.1	27.6
gemma-3-4b	26.7	28.5
Qwen2.5-7B	20.6	17.2
Qwen2.5-14B	10.6	8.8

Table 4 — GPU: mlx vs llama.cpp

model	mlx GPU	llama.cpp GPU
gemma-3-270m	436.6	231.6
Qwen2.5-0.5B	287.8	163.3
Qwen3-0.6B	248.3	165.4
Llama-3.2-1B	137.4	108.2
gemma-3-1b	151.4	98.7
LFM2-1.2B	148.4	120.2
Qwen3-1.7B	97.4	77.3
SmolLM3-3B	57.4	47.6
Llama-3.2-3B	54.3	45.5
Qwen2.5-3B	57.6	46.9
Phi-4-mini	45.6	37.0
Qwen3-4B	43.4	36.3
gemma-3-4b	44.4	37.5
Qwen2.5-7B	25.4	22.6
Qwen2.5-14B	13.0	11.7

Table 5 — simultaneous CPU+GPU: mlx vs llama.cpp

model	mlx CPU+GPU	llama.cpp CPU+GPU
gemma-3-270m	155.3 + 354.1 = 509.4	169.6 + 202.8 = 372.4
Qwen2.5-0.5B	105.6 + 223.4 = 329.0	119.9 + 142.5 = 262.4
Qwen3-0.6B	76.5 + 194.7 = 271.2	105.1 + 139.2 = 244.3
Llama-3.2-1B	58.9 + 95.9 = 154.8	63.2 + 87.2 = 150.4
gemma-3-1b	59.2 + 113.2 = 172.4	62.4 + 82.0 = 144.5
LFM2-1.2B	63.6 + 102.7 = 166.3	70.4 + 96.5 = 166.9
Qwen3-1.7B	39.4 + 71.0 = 110.4	44.2 + 62.0 = 106.3
SmolLM3-3B	26.2 + 43.7 = 69.9	25.5 + 37.0 = 62.5
Llama-3.2-3B	24.1 + 36.7 = 60.8	25.6 + 36.0 = 61.6
Qwen2.5-3B	24.8 + 38.9 = 63.7	26.1 + 37.1 = 63.2
Phi-4-mini	21.2 + 29.8 = 51.0	21.1 + 29.2 = 50.3
Qwen3-4B	18.4 + 29.5 = 47.9	20.3 + 28.8 = 49.1
gemma-3-4b	18.3 + 31.2 = 49.5	20.3 + 29.2 = 49.5
Qwen2.5-7B	12.1 + 16.4 = 28.5	12.4 + 17.4 = 29.8
Qwen2.5-14B	1.3 + 12.1 = 13.4	7.0 + 11.0 = 18.0

Also, argmax over a large vocab compared scores one element at a time in arg_reduce.cpp (~0.4ms per token on a 262k vocab). This adds a SIMD fast path for contiguous float rows, rows containing NaN fall back to the original loop, so results are always identical to before.

All local tests green

The int8 path is opt-in (MLX_CPU_INT8_QMM=1) so default numerics are unchanged and the full test suite passes as-is.