Function dot 

fn dot(a: &[f64], b: &[f64]) -> f64

Inner product with Neumaier-compensated summation, so the sum is accurate to about one rounding whatever the length (a plain left-to-right sum over the millions of background voxels carries ~1e-13 relative error).

cgsolve.m on PDF.m’s normal equations is very sensitive to rounding where the k-space kernel is used (B0 along an axis, or space = kspace): by the 30th iteration it amplifies a 1e-13 difference in α at the first iteration to ~10 % in α, concentrated in a few voxels (QSM-CI’s simulated head: one voxel near (89, 95, 46)). MATLAB’s MKL ddot and numpy’s BLAS sum in blocks and stay within ~1e-16 of the exact sum; a plain sum did not, and was what made the port differ from PDF.m by 0.61 % (sim) and 0.74 % (KKI, kspace) at the 30-iteration cap (astewartau/QSM.rs#405). With this sum: 0.007 % and 0.026 %, inside PDF.m’s own rounding spread (a numpy PDF.m is 0.003 % / 0.03 % from MATLAB, and moves 0.03 % on KKI for a 1-ulp change of its input).