From ab9179737f59582eaf2ce95b55352f960ab93aad Mon Sep 17 00:00:00 2001 From: godofecht Date: Mon, 28 Sep 2026 22:44:03 +0100 Subject: [PATCH 1/3] Take the hot loops off the shape that blocks vectorization The gate passed at 188 of 188, and affinity propagation passed it at 1.003x on one run and failed it at 0.947x on the next, on code that had not changed. Our own time moved by two tenths of a percent across those runs and scikit-learn's moved by fifteen, so the row was a coin toss rather than a win. Looking for room in it found something that is not about that row at all. Flow emits a range loop as for (int32_t k = 0; (0 <= n) ? k < n : k > n; k += (0 <= n) ? 1 : -1) so the stride is a runtime select and LLVM declines to vectorize the body. The disassembly of affinity_propagation_fit held no vector instruction at all. Writing the two inner loops as while loops with a constant step took the fit from 5.23 ms to 2.97 ms and the function from zero vector instructions to 38. That is Flow issue #1004, and until it is fixed the workaround is a while loop in anything hot. The same treatment on the kernel methods, plus two things worth doing anyway: The kernel matrix is symmetric under every kernel the module has, so it is computed over one triangle and mirrored, where it used to evaluate all n^2 entries and write each one through matrix_set. The kernel between two rows now takes pointers and loops over the features with a constant step, which reaches every caller including prediction. SVR's inner update, f = f + delta * row, is a saxpy, so it calls one. Local timings at -O3, against the last CI reading of the same rows: svr 10.82 ms -> 1.22 nu_svr 6.61 -> 1.52 affinity_propagation 6.47 -> 2.97 svc 0.90 -> 0.15 linear_svr 0.44 -> 0.20 --- lib/scikit/cluster.flow | 67 +++++++++++++++++----------- lib/scikit/svm.flow | 96 ++++++++++++++++++++++------------------- 2 files changed, 94 insertions(+), 69 deletions(-) diff --git a/lib/scikit/cluster.flow b/lib/scikit/cluster.flow index 3aaae37..3731f80 100644 --- a/lib/scikit/cluster.flow +++ b/lib/scikit/cluster.flow @@ -5,6 +5,13 @@ import "lib/scikit/matrix.flow" import "lib/scikit/prng.flow" import "lib/scikit/blas.flow" +extern { + # clang lowers both to a single instruction, so a loop written with them + # holds no branch and the vectorizer can take it. + function fmaxf(a: f32, b: f32) -> f32 + function fminf(a: f32, b: f32) -> f32 +} + export struct KMeans { centroids: ptr >, n_clusters: i32, @@ -2007,17 +2014,34 @@ export function affinity_propagation_fit(X: Matrix, damping: f32, max_iter: i32, } } - for k in 0 to n { - let mut max_other: f32 = best - if k == best_k { max_other = second } - let v: f32 = take * (si3[k] - max_other) + keep * ri2[k] + # Every element of the row takes the row maximum, and the two + # elements that do not are corrected afterwards. Deciding it per + # element put a branch on k against best_k and another on k + # against i inside a loop that is otherwise a fused multiply add + # and a max, which is the difference between a vector loop and a + # scalar one. + let old_best: f32 = ri2[best_k] + # Written as a while loop with a constant step. A `for k in 0 to n` + # emits an increment chosen by a runtime test of the bounds, which + # is not a stride the LLVM vectorizer will accept, and this loop + # was scalar with it. + let mut k: i32 = 0 + while k < n { + let v: f32 = take * (si3[k] - best) + keep * ri2[k] ri2[k] = v - if k == i { - rkk[i] = v - } elif v > 0.0 { - col_pos[k] = col_pos[k] + v - } + col_pos[k] = col_pos[k] + fmaxf(v, 0.0) + k = k + 1 } + + # The largest entry of the row competes against the runner up. + let fixed: f32 = take * (si3[best_k] - second) + keep * old_best + col_pos[best_k] = col_pos[best_k] + fmaxf(fixed, 0.0) - fmaxf(ri2[best_k], 0.0) + ri2[best_k] = fixed + + # The diagonal is the self responsibility, and it is not part of + # the column sum. + rkk[i] = ri2[i] + col_pos[i] = col_pos[i] - fmaxf(ri2[i], 0.0) } # The availability update. The term that depends only on the column is @@ -2029,21 +2053,16 @@ export function affinity_propagation_fit(X: Matrix, damping: f32, max_iter: i32, for i in 0 to n { let ai2: ptr = A + i * n let ri4: ptr = R + i * n - for k in 0 to i { - let mut val: f32 = base[k] - let r_ik: f32 = ri4[k] - if r_ik > 0.0 { val = val - r_ik } - if val > 0.0 { val = 0.0 } - ai2[k] = take * val + keep * ai2[k] - } - ai2[i] = take * base[i] + keep * ai2[i] - for k in i + 1 to n { - let mut val2: f32 = base[k] - let r_ik2: f32 = ri4[k] - if r_ik2 > 0.0 { val2 = val2 - r_ik2 } - if val2 > 0.0 { val2 = 0.0 } - ai2[k] = take * val2 + keep * ai2[k] - } + let old_diag: f32 = ai2[i] + let mut k2: i32 = 0 + while k2 < n { + let val: f32 = fminf(base[k2] - fmaxf(ri4[k2], 0.0), 0.0) + ai2[k2] = take * val + keep * ai2[k2] + k2 = k2 + 1 + } + # The diagonal keeps the whole column sum rather than the clipped + # one, and is written after the loop for the reason above. + ai2[i] = take * base[i] + keep * old_diag } # Convergence is the set of exemplars holding still, which is what diff --git a/lib/scikit/svm.flow b/lib/scikit/svm.flow index aa9dec8..f98b54a 100644 --- a/lib/scikit/svm.flow +++ b/lib/scikit/svm.flow @@ -2202,61 +2202,68 @@ function _svm_tanh_f32(x: f32) -> f32 { return val } -function _svm_kernel_rows(X1: Matrix, i1: i32, X2: Matrix, i2: i32, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { - # The two rows are addressed once. matrix_at copies the Matrix struct by - # value for every element, and this is called once per pair of rows in - # every kernel prediction, so it was two struct copies per feature per - # pair. - let r1: ptr = X1.data + i1 * X1.cols - let r2: ptr = X2.data + i2 * X2.cols - if kernel == KERNEL_LINEAR { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + r1[k] * r2[k] - } - return s - } - if kernel == KERNEL_POLY { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + r1[k] * r2[k] - } - let base: f32 = gamma * s + coef0 - let val: f32 = _svm_pow_f32(base, degree) - return val - } +# The kernel between two rows, taken as pointers. +# +# Every loop here is a while loop with a constant step. Flow's `for i in a to b` +# emits an increment chosen by a runtime test of the bounds, which LLVM will not +# accept as a stride, so a `for` loop over an array is scalar code. Flow issue +# #1004. +function _svm_kernel_vec(r1: ptr, r2: ptr, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { if kernel == KERNEL_RBF { let mut s: f32 = 0.0 - for k in 0 to n_features { + let mut k: i32 = 0 + while k < n_features { let d: f32 = r1[k] - r2[k] s = s + d * d + k = k + 1 } - let val: f32 = exp((0.0 - gamma * s) as f64) as f32 - return val - } - if kernel == KERNEL_SIGMOID { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + matrix_at(X1, i1, k) * matrix_at(X2, i2, k) - } - let val: f32 = _svm_tanh_f32(gamma * s + coef0) - return val + return exp((0.0 - gamma * s) as f64) as f32 } - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + matrix_at(X1, i1, k) * matrix_at(X2, i2, k) + + let mut dot: f32 = 0.0 + let mut k2: i32 = 0 + while k2 < n_features { + dot = dot + r1[k2] * r2[k2] + k2 = k2 + 1 } - return s + if kernel == KERNEL_LINEAR { return dot } + if kernel == KERNEL_POLY { return _svm_pow_f32(gamma * dot + coef0, degree) } + if kernel == KERNEL_SIGMOID { return _svm_tanh_f32(gamma * dot + coef0) } + return dot } +function _svm_kernel_rows(X1: Matrix, i1: i32, X2: Matrix, i2: i32, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { + # The two rows are addressed once. matrix_at copies the Matrix struct by + # value for every element, and this is called once per pair of rows in + # every kernel prediction, so it was two struct copies per feature per + # pair. + let r1: ptr = X1.data + i1 * X1.cols + let r2: ptr = X2.data + i2 * X2.cols + return _svm_kernel_vec(r1, r2, n_features, kernel, gamma, degree, coef0) +} + +# Every kernel here is symmetric, so one triangle is computed and mirrored. +# The entries are written through row pointers, where matrix_set copied the +# Matrix struct by value for each of the n^2 of them. function _svm_kernel_matrix(X: Matrix, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> Matrix { let n: i32 = X.rows + let p: i32 = X.cols let K: Matrix = matrix_new(n, n) - for i in 0 to n { - for j in 0 to n { - let val: f32 = _svm_kernel_rows(X, i, X, j, X.cols, kernel, gamma, degree, coef0) - matrix_set(K, i, j, val) + let kdata: ptr = K.data + let xdata: ptr = X.data + let mut i: i32 = 0 + while i < n { + let ri: ptr = xdata + i * p + let ki: ptr = kdata + i * n + let mut j: i32 = i + while j < n { + let rj: ptr = xdata + j * p + let val: f32 = _svm_kernel_vec(ri, rj, p, kernel, gamma, degree, coef0) + ki[j] = val + kdata[j * n + i] = val + j = j + 1 } + i = i + 1 } return K } @@ -2682,9 +2689,8 @@ export function svr_fit(X: Matrix, y: ptr, C: f32, epsilon: f32, kernel: i3 let adelta: f32 = fabs((delta) as f64) as f32 if adelta > 0.0 { alpha_diff[i] = a_new - for k in 0 to n { - f[k] = f[k] + row[k] * delta - } + # f = f + delta * row is a saxpy, and BLAS has one. + blas_axpy(n, delta, row, f) if adelta > max_change { max_change = adelta } } } From 8ba0c12e6b3e6e20a1570e9278dec33225fdb740 Mon Sep 17 00:00:00 2001 From: godofecht Date: Mon, 28 Sep 2026 23:02:44 +0100 Subject: [PATCH 2/3] Take four more rows off the loops that were holding them back The gate passed again at 188 of 188, and the row at the bottom of it moved from affinity propagation to label spreading at 1.049x. The whole distribution moved down that run, on rows nobody had touched, which is the runner rather than the code. A row at 1.05x is the next coin toss, so these are the four narrowest. label_spreading propagated its distributions with a scalar triple loop that walked Y down a column of an array of row allocations, a cache miss per element, and allocated a fresh set of rows on every iteration. The propagation is a matrix product, so it is one sgemm against contiguous storage. 1.04 ms to 0.12. knn_regressor read every element of both designs through matrix_at, and ran over all k neighbours for every training point to find the slot to replace when it only has to look when a point actually displaces one. 1.24 ms to 0.53. isomap held its distances, geodesics and centering in arrays of row allocations. Floyd-Warshall now walks contiguous rows with the term that depends only on i and k lifted out, and the comparison written as a minimum so the inner loop holds no branch. The power iteration is one sgemm per pass over the whole component block. 6.22 ms to 4.04. local_outlier_factor and the label spreading affinity get their distance loops as while loops, for the reason the last commit gives: Flow emits a range loop with a stride chosen at runtime, and LLVM will not vectorize that. Flow issue #1004. --- lib/scikit/manifold.flow | 199 +++++++++++++++----------------- lib/scikit/neighbors.flow | 59 +++++++--- lib/scikit/semi_supervised.flow | 89 ++++++++------ 3 files changed, 187 insertions(+), 160 deletions(-) diff --git a/lib/scikit/manifold.flow b/lib/scikit/manifold.flow index d323c36..1a5bd5c 100644 --- a/lib/scikit/manifold.flow +++ b/lib/scikit/manifold.flow @@ -8,6 +8,9 @@ import "lib/scikit/blas.flow" extern { function rand() -> i32 function srand(seed: u32) -> void + # clang lowers this to a single instruction, so a loop written with it + # holds no branch and the vectorizer can take it. + function fminf(a: f32, b: f32) -> f32 } export struct TSNE { @@ -306,125 +309,132 @@ export struct Isomap { fitted: bool } +# Modified Gram-Schmidt on the columns of a row-major n by k block. +# A column that collapses is left at zero, and the remaining columns stay +# orthonormal. +function _orthonormalize_block(V: ptr, n: i32, k: i32) -> void { + for c in 0 to k { + for prev in 0 to c { + let mut dot: f32 = 0.0 + for i in 0 to n { dot = dot + V[i * k + c] * V[i * k + prev] } + for i in 0 to n { V[i * k + c] = V[i * k + c] - dot * V[i * k + prev] } + } + let mut norm: f32 = 0.0 + for i in 0 to n { + let v: f32 = V[i * k + c] + norm = norm + v * v + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + for i in 0 to n { V[i * k + c] = V[i * k + c] * inv } + } + } +} + export function isomap_fit(X: Matrix, n_components: i32, n_neighbors: i32) -> Isomap { let n: i32 = X.rows let p: i32 = X.cols + let kc: i32 = n_components - let D: ptr > = malloc((n as i64) * 8) as ptr > + # Distances, then geodesics, then the double centering, all in contiguous + # blocks. Every inner loop below is a while loop with a constant step, + # because Flow emits a range loop with a stride chosen at runtime and LLVM + # will not vectorize that. Flow issue #1004. + let xdata: ptr = X.data + let geo: ptr = array_new_f32(n * n) for i in 0 to n { - D[i] = array_new_f32(n) - for j in 0 to n { + let ri: ptr = xdata + i * p + let gi: ptr = geo + i * n + for j in i + 1 to n { + let rj: ptr = xdata + j * p let mut s: f32 = 0.0 - for k in 0 to p { - let d: f32 = matrix_at(X, i, k) - matrix_at(X, j, k) + let mut m: i32 = 0 + while m < p { + let d: f32 = ri[m] - rj[m] s = s + d * d + m = m + 1 } - D[i][j] = sqrt((s as f64)) as f32 + let dist: f32 = sqrt((s) as f64) as f32 + gi[j] = dist + geo[j * n + i] = dist } } - let geo: ptr > = malloc((n as i64) * 8) as ptr > - for i in 0 to n { - geo[i] = array_copy_f32(D[i], n) - } - + # Floyd-Warshall. The term that depends only on i and k is lifted out of + # the inner loop, and the comparison is a minimum, so the loop holds no + # branch and vectorizes. for k in 0 to n { + let rowk: ptr = geo + k * n for i in 0 to n { - for j in 0 to n { - if geo[i][k] + geo[k][j] < geo[i][j] { - geo[i][j] = geo[i][k] + geo[k][j] - } + let rowi: ptr = geo + i * n + let dik: f32 = rowi[k] + let mut j: i32 = 0 + while j < n { + rowi[j] = fminf(rowi[j], dik + rowk[j]) + j = j + 1 } } } let mut mean_d: f32 = 0.0 for i in 0 to n { - for j in 0 to n { - geo[i][j] = geo[i][j] * geo[i][j] - mean_d = mean_d + geo[i][j] + let gi2: ptr = geo + i * n + let mut j2: i32 = 0 + while j2 < n { + let v: f32 = gi2[j2] * gi2[j2] + gi2[j2] = v + mean_d = mean_d + v + j2 = j2 + 1 } } mean_d = mean_d / ((n * n) as f32) - let B: ptr > = malloc((n as i64) * 8) as ptr > + let B: ptr = array_new_f32(n * n) for i in 0 to n { - B[i] = array_new_f32(n) - for j in 0 to n { - B[i][j] = -0.5 * (geo[i][j] - mean_d) + let gi3: ptr = geo + i * n + let bi: ptr = B + i * n + let mut j3: i32 = 0 + while j3 < n { + bi[j3] = -0.5 * (gi3[j3] - mean_d) + j3 = j3 + 1 } } + array_free_f32(geo) - let embedding: ptr = array_new_f32(n * n_components) - # Track found eigenvectors for deflation - let found_vecs: ptr > = malloc((n_components as i64) * 8) as ptr > - for c in 0 to n_components { - let v: ptr = array_new_f32(n) - srand(42 + c) - for i in 0 to n { v[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } - - # Orthogonalize against previously found eigenvectors - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + v[i] * found_vecs[prev][i] } - for i in 0 to n { v[i] = v[i] - dot * found_vecs[prev][i] } - } - - let mut iter: i32 = 0 - while iter < 200 { - let Bv: ptr = array_new_f32(n) - for i in 0 to n { - let mut s: f32 = 0.0 - for j in 0 to n { s = s + B[i][j] * v[j] } - Bv[i] = s - } - - # Deflate against previous eigenvectors - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + Bv[i] * found_vecs[prev][i] } - for i in 0 to n { Bv[i] = Bv[i] - dot * found_vecs[prev][i] } - } - - let mut norm: f32 = 0.0 - for i in 0 to n { norm = norm + Bv[i] * Bv[i] } - norm = sqrt((norm as f64)) as f32 - if norm < 0.0000000001 { norm = 0.0000000001 } - - # Convergence check - let mut max_change: f32 = 0.0 - for i in 0 to n { - let new_v: f32 = Bv[i] / norm - let diff: f32 = new_v - v[i] - if diff < 0.0 { diff = -diff } - if diff > max_change { max_change = diff } - v[i] = new_v - } - array_free_f32(Bv) + # Block power iteration. One sgemm per iteration advances every component, + # where the old loop ran a scalar matvec per component per iteration and + # allocated its result buffer inside the iteration. + let V: ptr = array_new_f32(n * kc) + let Z: ptr = array_new_f32(n * kc) + srand(42) + for i in 0 to n * kc { V[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } + _orthonormalize_block(V, n, kc) - if max_change < 0.000001 { break } - iter = iter + 1 - } + let mut iter: i32 = 0 + while iter < 200 { + blas_matmat(B, V, Z, n, kc, n, 1.0, 0.0) + _orthonormalize_block(Z, n, kc) - for i in 0 to n { - embedding[i * n_components + c] = v[i] + let mut max_change: f32 = 0.0 + for i in 0 to n * kc { + let mut diff: f32 = Z[i] - V[i] + if diff < 0.0 { diff = 0.0 - diff } + if diff > max_change { max_change = diff } + V[i] = Z[i] } - found_vecs[c] = v + if max_change < 0.000001 { break } + iter = iter + 1 } - for c in 0 to n_components { array_free_f32(found_vecs[c]) } - free(found_vecs as ptr) - - for i in 0 to n { array_free_f32(D[i]); array_free_f32(geo[i]); array_free_f32(B[i]) } - free(D as ptr) - free(geo as ptr) - free(B as ptr) + array_free_f32(B) + array_free_f32(Z) return Isomap { - embedding: embedding, + embedding: V, n_samples: n, n_features: p, - n_components: n_components, + n_components: kc, n_neighbors: n_neighbors, fitted: true } @@ -456,29 +466,6 @@ export struct LLE { fitted: bool } -# Modified Gram-Schmidt on the columns of a row-major n by k block. -# A column that collapses is left at zero, and the remaining columns stay -# orthonormal. -function _orthonormalize_block(V: ptr, n: i32, k: i32) -> void { - for c in 0 to k { - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + V[i * k + c] * V[i * k + prev] } - for i in 0 to n { V[i * k + c] = V[i * k + c] - dot * V[i * k + prev] } - } - let mut norm: f32 = 0.0 - for i in 0 to n { - let v: f32 = V[i * k + c] - norm = norm + v * v - } - norm = sqrt((norm) as f64) as f32 - if norm > 0.0000000001 { - let inv: f32 = 1.0 / norm - for i in 0 to n { V[i * k + c] = V[i * k + c] * inv } - } - } -} - export function lle_fit(X: Matrix, n_components: i32, n_neighbors: i32) -> LLE { let n: i32 = X.rows let p: i32 = X.cols diff --git a/lib/scikit/neighbors.flow b/lib/scikit/neighbors.flow index ceeb059..b0444da 100644 --- a/lib/scikit/neighbors.flow +++ b/lib/scikit/neighbors.flow @@ -212,39 +212,61 @@ export function knn_regressor_predict(model: KNNRegressor, X: Matrix) -> ptr = array_new_f32(X.rows) let neighbors: ptr = malloc((model.k as i64) * 12) as ptr + # Row pointers into both designs, and the worst of the k held in hand. + # The distance loop read every element through matrix_at, which copies the + # Matrix struct per element, and the loop that finds the slot to replace + # ran over all k for every training point when it only has to run when a + # point actually displaces one. The feature loop is a while loop with a + # constant step, because Flow emits a range loop with a stride chosen at + # runtime and LLVM will not vectorize that. Flow issue #1004. + let cols: i32 = X.cols + let xdata: ptr = X.data + let tdata: ptr = model.X_train.data + let n_train: i32 = model.X_train.rows + let k: i32 = model.k + for i in 0 to X.rows { - for j in 0 to model.k { + let xrow: ptr = xdata + i * cols + for j in 0 to k { neighbors[j].distance = 9999999999.0 neighbors[j].label = 0.0 } + let mut worst_idx: i32 = 0 + let mut worst_dist: f32 = 9999999999.0 - for j in 0 to model.X_train.rows { - # Use f64 for distance to match sklearn precision + for j in 0 to n_train { + let trow: ptr = tdata + j * cols + # f64 for the distance, to match scikit-learn's precision. let mut dist: f64 = 0.0 - for f in 0 to X.cols { - let d: f64 = (matrix_at(X, i, f) - matrix_at(model.X_train, j, f)) as f64 + let mut f: i32 = 0 + while f < cols { + let d: f64 = (xrow[f] - trow[f]) as f64 dist = dist + d * d + f = f + 1 } let dist_f32: f32 = dist as f32 - let mut worst_idx: i32 = 0 - for n in 1 to model.k { - if neighbors[n].distance > neighbors[worst_idx].distance { - worst_idx = n - } - } - - if dist_f32 < neighbors[worst_idx].distance { + if dist_f32 < worst_dist { neighbors[worst_idx].distance = dist_f32 neighbors[worst_idx].label = model.y_train[j] + let mut w: i32 = 0 + let mut wd: f32 = neighbors[0].distance + for m in 1 to k { + if neighbors[m].distance > wd { + wd = neighbors[m].distance + w = m + } + } + worst_idx = w + worst_dist = wd } } let mut sum: f32 = 0.0 - for n in 0 to model.k { + for n in 0 to k { sum = sum + neighbors[n].label } - preds[i] = sum / (model.k as f32) + preds[i] = sum / (k as f32) } free(neighbors as ptr) @@ -597,10 +619,15 @@ export function local_outlier_factor_fit(X: Matrix, n_neighbors: i32) -> LocalOu di[i] = 0.0 for j in i + 1 to n { let rj: ptr = xdata + j * p + # A while loop with a constant step, because Flow emits a range + # loop with a stride chosen at runtime and LLVM will not vectorize + # that. Flow issue #1004. let mut acc: f32 = 0.0 - for m in 0 to p { + let mut m: i32 = 0 + while m < p { let d: f32 = ri[m] - rj[m] acc = acc + d * d + m = m + 1 } let dist: f32 = sqrt((acc as f64)) as f32 di[j] = dist diff --git a/lib/scikit/semi_supervised.flow b/lib/scikit/semi_supervised.flow index 376d6f7..05f618c 100644 --- a/lib/scikit/semi_supervised.flow +++ b/lib/scikit/semi_supervised.flow @@ -2,6 +2,7 @@ # Implements LabelPropagation, LabelSpreading, SelfTrainingClassifier. import "lib/scikit/matrix.flow" +import "lib/scikit/blas.flow" extern { function malloc(size: i64) -> ptr @@ -50,10 +51,15 @@ function _build_affinity(X: Matrix, gamma: f32) -> ptr > { row_i[i] = self_affinity for j in i + 1 to n { let base_j: i32 = j * p + # A while loop with a constant step, because Flow emits a range + # loop with a stride chosen at runtime and LLVM will not vectorize + # that. Flow issue #1004. let mut s: f32 = 0.0 - for k in 0 to p { + let mut k: i32 = 0 + while k < p { let d: f32 = X.data[base_i + k] - X.data[base_j + k] s = s + d * d + k = k + 1 } let w: f32 = exp((-gamma * s) as f64) as f32 row_i[j] = w @@ -259,58 +265,69 @@ export function label_spreading_fit(X: Matrix, y: ptr, n_classes: i32, gamm } } - let S: ptr > = malloc((n as i64) * 8) as ptr > + # The propagation is new_Y = S * Y, a matrix product, so it is one sgemm + # against contiguous storage. It used to be a scalar triple loop that + # walked Y down a column of an array of row allocations, which is a cache + # miss per element, and it allocated a fresh n by n_classes set of rows on + # every iteration. + let Sc: ptr = array_new_f32(n * n) for i in 0 to n { - S[i] = array_new_f32(n) - for j in 0 to n { - S[i][j] = alpha * W[i][j] / deg[i] + let srow: ptr = Sc + i * n + let wrow: ptr = W[i] + let scale: f32 = alpha / deg[i] + let mut j: i32 = 0 + while j < n { + srow[j] = scale * wrow[j] + j = j + 1 } } - let Y: ptr > = malloc((n as i64) * 8) as ptr > - let clamp: ptr > = malloc((n as i64) * 8) as ptr > + let nc: i32 = n_classes + let Yc: ptr = array_new_f32(n * nc) + let Ync: ptr = array_new_f32(n * nc) + let clampc: ptr = array_new_f32(n * nc) for i in 0 to n { - Y[i] = array_new_f32(n_classes) - clamp[i] = array_new_f32(n_classes) - if y[i] >= 0 && y[i] < n_classes { - clamp[i][y[i]] = 1.0 - Y[i][y[i]] = 1.0 + if y[i] >= 0 && y[i] < nc { + clampc[i * nc + y[i]] = 1.0 + Yc[i * nc + y[i]] = 1.0 } } + let hold: f32 = 1.0 - alpha + let total: i32 = n * nc let mut n_iter: i32 = 0 for iter in 0 to max_iter { n_iter = iter + 1 let mut max_change: f32 = 0.0 - let new_Y: ptr > = malloc((n as i64) * 8) as ptr > - for i in 0 to n { - new_Y[i] = array_new_f32(n_classes) - } - - for i in 0 to n { - for c in 0 to n_classes { - let mut s: f32 = 0.0 - for j in 0 to n { - s = s + S[i][j] * Y[j][c] - } - new_Y[i][c] = s + (1.0 - alpha) * clamp[i][c] - } - } + blas_matmat(Sc, Yc, Ync, n, nc, n, 1.0, 0.0) - for i in 0 to n { - for c in 0 to n_classes { - let change: f32 = fabs((new_Y[i][c] - Y[i][c]) as f64) as f32 - if change > max_change { max_change = change } - Y[i][c] = new_Y[i][c] - } - array_free_f32(new_Y[i]) + let mut idx: i32 = 0 + while idx < total { + let v: f32 = Ync[idx] + hold * clampc[idx] + let change: f32 = fabs((v - Yc[idx]) as f64) as f32 + if change > max_change { max_change = change } + Yc[idx] = v + idx = idx + 1 } - free(new_Y as ptr) if max_change < tol { break } } + # The model hands back its distributions a row at a time, so the rows are + # built once here. + let Y: ptr > = malloc((n as i64) * 8) as ptr > + for i in 0 to n { + let row: ptr = array_new_f32(nc) + let src: ptr = Yc + i * nc + for c in 0 to nc { row[c] = src[c] } + Y[i] = row + } + array_free_f32(Sc) + array_free_f32(Yc) + array_free_f32(Ync) + array_free_f32(clampc) + let labels: ptr = malloc((n as i64) * 4) as ptr for i in 0 to n { let mut best_c: i32 = 0 @@ -326,12 +343,8 @@ export function label_spreading_fit(X: Matrix, y: ptr, n_classes: i32, gamm for i in 0 to n { array_free_f32(W[i]) - array_free_f32(S[i]) - array_free_f32(clamp[i]) } free(W as ptr) - free(S as ptr) - free(clamp as ptr) array_free_f32(deg) return LabelSpreading { From 6f4d538b9d26b2bd477d3f402415ea2e7dddb197 Mon Sep 17 00:00:00 2001 From: godofecht Date: Mon, 28 Sep 2026 23:17:09 +0100 Subject: [PATCH 3/3] Take the three generalized linear models off it too gamma_regressor was among the narrowest rows left at 1.43x. Its gradient loop and its weight update are while loops with a constant step now, for the reason the two commits before this one give. 0.68 ms to 0.49 on this machine, and poisson and tweedie share the loop. --- lib/scikit/linear.flow | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/lib/scikit/linear.flow b/lib/scikit/linear.flow index 976f312..c8b5193 100644 --- a/lib/scikit/linear.flow +++ b/lib/scikit/linear.flow @@ -7252,8 +7252,12 @@ function _glm_gradient_fit(X: Matrix, y: ptr, weights: ptr, bias_init: for iter in 0 to max_iter { blas_matvec(xdata, n, p, weights, eta, 1.0, 0.0) + # While loops with a constant step. Flow emits a range loop with a + # stride chosen at runtime, which LLVM will not vectorize. Flow issue + # #1004. let mut grad_b: f32 = 0.0 - for i in 0 to n { + let mut i: i32 = 0 + while i < n { let e: f32 = eta[i] + bias let mut mu: f32 = 0.0 if kind == 0 { mu = _poisson_inverse_link(e) } @@ -7262,11 +7266,14 @@ function _glm_gradient_fit(X: Matrix, y: ptr, weights: ptr, bias_init: let r: f32 = y[i] - mu resid[i] = r grad_b = grad_b + r + i = i + 1 } blas_matvec_trans(xdata, n, p, resid, grad_w, 1.0, 0.0) - for j in 0 to p { + let mut j: i32 = 0 + while j < p { let g: f32 = grad_w[j] * inv_n - alpha * weights[j] weights[j] = weights[j] + lr * g + j = j + 1 } bias = bias + lr * grad_b * inv_n }