perf: CSR-flatten L/U factors for cache locality and fewer allocs

simplex/ft.go

@@ -23,7 +23,8 @@ func (f *ftLU) clone() *ftLU {
 		rlist: append([]ftR(nil), f.rlist...),
 		nUpd:  f.nUpd,
 		lPr:   append([]int(nil), f.lPr...),
-		lCol:  make([][]int32, f.m), lVal: make([][]float64, f.m),
+		lFR:   append([]int32(nil), f.lFR...), lFV: append([]float64(nil), f.lFV...),
+		lPtr: append([]int32(nil), f.lPtr...),
 		uCol: make([][]int32, f.m), uVal: make([][]float64, f.m),
 		ucRows: make([][]int32, f.m),
 		// transient scratch is shared (single-threaded, sequential solves);
@@ -35,8 +36,6 @@ func (f *ftLU) clone() *ftLU {
 		wColLen: f.wColLen,
 	}
 	for s := range f.m {
-		c.lCol[s] = append([]int32(nil), f.lCol[s]...)
-		c.lVal[s] = append([]float64(nil), f.lVal[s]...)
 		c.uCol[s] = append([]int32(nil), f.uCol[s]...)
 		c.uVal[s] = append([]float64(nil), f.uVal[s]...)
 		c.ucRows[s] = append([]int32(nil), f.ucRows[s]...)
@@ -48,11 +47,12 @@ func (f *ftLU) clone() *ftLU {
 // keyed by stable orig rows, U by stable basis cols; "step" order remaps freely.
 type ftLU struct {
 	m       int
-	lCol    [][]int32   // factorize step s -> orig rows eliminated at s
-	lVal    [][]float64 // matching L multipliers
-	lPr     []int       // factorize step s -> its pivot orig row (immutable)
-	udiag   []float64   // udiag[s] = U[s][s]
-	uCol    [][]int32   // step-row s -> basis cols c with rinvcol[c] > s
+	lFR     []int32   // flat L: step s eliminated rows, range [lPtr[s],lPtr[s+1])
+	lFV     []float64 // matching L multipliers (contiguous, cache-friendly)
+	lPtr    []int32   // step s -> start offset into lFR/lFV (len m+1)
+	lPr     []int     // factorize step s -> its pivot orig row (immutable)
+	udiag   []float64 // udiag[s] = U[s][s]
+	uCol    [][]int32 // step-row s -> basis cols c with rinvcol[c] > s
 	uVal    [][]float64
 	ucRows  [][]int32 // basis col -> orig rows whose U row may hold it (lazy)
 	prow    []int     // prow[step] = orig row
@@ -106,7 +106,6 @@ func newFTLU(m int, a []float64) *ftLU {
 func newFTLUSparse(m int, colRow [][]int32, colVal [][]float64) *ftLU {
 	f := &ftLU{
 		m: m, udiag: make([]float64, m),
-		lCol: make([][]int32, m), lVal: make([][]float64, m),
 		lPr:  make([]int, m),
 		uCol: make([][]int32, m), uVal: make([][]float64, m),
 		ucRows: make([][]int32, m),
@@ -131,13 +130,13 @@ func (f *ftLU) rebuild(colRow [][]int32, colVal [][]float64) bool {
 	m := f.m
 	f.rlist = f.rlist[:0]
 	f.nUpd = 0
+	f.lFR, f.lFV, f.lPtr = f.lFR[:0], f.lFV[:0], f.lPtr[:0]
 	rowCol, rowVal, colRows := f.wRowCol, f.wRowVal, f.wColRows
 	buckets := f.wBuck
 	rowUsed, colUsed, mark, acc := f.wRowUsed, f.wColUsed, f.wMrk, f.wAcc
 	for r := range m {
 		rowCol[r], rowVal[r], colRows[r] = rowCol[r][:0], rowVal[r][:0], colRows[r][:0]
 		rowUsed[r], colUsed[r], mark[r], acc[r] = false, false, false, 0
-		f.lCol[r], f.lVal[r] = f.lCol[r][:0], f.lVal[r][:0]
 		f.uCol[r], f.uVal[r] = f.uCol[r][:0], f.uVal[r][:0]
 		f.ucRows[r] = f.ucRows[r][:0]
 		buckets[r] = buckets[r][:0]
@@ -159,6 +158,7 @@ func (f *ftLU) rebuild(colRow [][]int32, colVal [][]float64) bool {
 	}
 	minL := 0
 	for step := range m {
+		f.lPtr = append(f.lPtr, int32(len(f.lFR)))
 		// pivot row: shortest active row (via buckets)
 		pr := -1
 		for minL <= m {
@@ -238,8 +238,8 @@ func (f *ftLU) rebuild(colRow [][]int32, colVal [][]float64) bool {
 				continue
 			}
 			mult := hit / pv
-			f.lCol[step] = append(f.lCol[step], int32(rr))
-			f.lVal[step] = append(f.lVal[step], mult)
+			f.lFR = append(f.lFR, int32(rr))
+			f.lFV = append(f.lFV, mult)
 			rc, rv := rowCol[rr], rowVal[rr]
 			w := 0
 			for k, c := range rc {
@@ -281,6 +281,7 @@ func (f *ftLU) rebuild(colRow [][]int32, colVal [][]float64) bool {
 			acc[c] = 0
 		}
 	}
+	f.lPtr = append(f.lPtr, int32(len(f.lFR)))
 	copy(f.lPr, f.prow) // L stays keyed by factorize-time pivot rows
 	f.buildURows(rowCol, rowVal)
 	return true
@@ -309,9 +310,8 @@ func (f *ftLU) forwardLR(b []float64) {
 		if v == 0 {
 			continue
 		}
-		lc, lv := f.lCol[s], f.lVal[s]
-		for k, r := range lc {
-			b[r] -= lv[k] * v
+		for k := f.lPtr[s]; k < f.lPtr[s+1]; k++ {
+			b[f.lFR[k]] -= f.lFV[k] * v
 		}
 	}
 	for _, op := range f.rlist {
@@ -362,9 +362,8 @@ func (f *ftLU) btran(c []float64) {
 	for s := m - 1; s >= 0; s-- { // L^-T back by column
 		pr := f.lPr[s]
 		v := z[pr]
-		lc, lv := f.lCol[s], f.lVal[s]
-		for k, r := range lc {
-			v -= lv[k] * z[r]
+		for k := f.lPtr[s]; k < f.lPtr[s+1]; k++ {
+			v -= f.lFV[k] * z[f.lFR[k]]
 		}
 		z[pr] = v
 		c[pr] = v

simplex/ftreal_test.go

@@ -0,0 +1,74 @@
+package simplex
+
+import (
+	"math/rand"
+	"testing"
+)
+
+// BenchmarkFTReplaceReal: per-pivot cost at the engine's refactor cadence
+// (rebuild every maxEtas), unlike BenchmarkFTReplace's unbounded rlist growth.
+func BenchmarkFTReplaceReal(b *testing.B) {
+	const m = 300
+	rng := rand.New(rand.NewSource(9))
+	a := make([]float64, m*m)
+	for i := range a {
+		if rng.Float64() < 0.03 {
+			a[i] = rng.NormFloat64()
+		}
+	}
+	for d := range m {
+		a[d*m+d] += float64(m)
+	}
+	cr := make([][]int32, m)
+	cv := make([][]float64, m)
+	for c := range m {
+		for r := range m {
+			if v := a[c*m+r]; v != 0 {
+				cr[c] = append(cr[c], int32(r))
+				cv[c] = append(cv[c], v)
+			}
+		}
+	}
+	nv := make([]float64, m)
+	for i := range nv {
+		nv[i] = rng.NormFloat64()
+	}
+	rhs := make([]float64, m)
+	for i := range rhs {
+		rhs[i] = rng.NormFloat64()
+	}
+	f := newFTLU(m, a)
+	v := make([]float64, m)
+	n := 0
+	for b.Loop() {
+		if f.nUpd >= maxEtas {
+			f.rebuild(cr, cv)
+		}
+		col := n % m
+		nv[col] += float64(m)
+		f.replaceColumn(col, nv)
+		nv[col] -= float64(m)
+		copy(v, rhs)
+		f.ftran(v)
+		n++
+	}
+}
+
+// BenchmarkFTClone: per-B&B-child factor deep copy.
+func BenchmarkFTClone(b *testing.B) {
+	const m = 300
+	rng := rand.New(rand.NewSource(9))
+	a := make([]float64, m*m)
+	for i := range a {
+		if rng.Float64() < 0.03 {
+			a[i] = rng.NormFloat64()
+		}
+	}
+	for d := range m {
+		a[d*m+d] += float64(m)
+	}
+	f := newFTLU(m, a)
+	for b.Loop() {
+		_ = f.clone()
+	}
+}

simplex/sparselu.go

@@ -8,14 +8,16 @@ import "math"
 // per-step multiplier columns, U upper stored as per-step rows.
 type sparseLU struct {
 	k         int
-	rowPerm   []int32 // step -> kernel-local row
-	colPerm   []int32 // step -> kernel-local col
-	stepOfRow []int32 // kernel-local row -> step
-	lIdx      [][]int32
-	lVal      [][]float64 // L[row, step] multipliers, rows in kernel-local ids
+	rowPerm   []int32   // step -> kernel-local row
+	colPerm   []int32   // step -> kernel-local col
+	stepOfRow []int32   // kernel-local row -> step
+	lFlatIdx  []int32   // flat L: step s multipliers in [lPtr[s],lPtr[s+1])
+	lFlatVal  []float64 // matching L values, kernel-local row ids
+	lPtr      []int32   // step -> start offset into lFlat* (len k+1)
 	uDiag     []float64
-	uIdx      [][]int32 // U[step, step'] entries, step' > step
-	uVal      [][]float64
+	uFlatIdx  []int32 // flat U: step s entries in [uPtr[s],uPtr[s+1])
+	uFlatVal  []float64
+	uPtr      []int32 // step -> start offset into uFlat* (len k+1)
 	work      []float64
 }
 
@@ -88,11 +90,9 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64, w *luWS) *sparseLU
 		rowPerm:   make([]int32, k),
 		colPerm:   make([]int32, k),
 		stepOfRow: make([]int32, k),
-		lIdx:      make([][]int32, k),
-		lVal:      make([][]float64, k),
 		uDiag:     make([]float64, k),
-		uIdx:      make([][]int32, k),
-		uVal:      make([][]float64, k),
+		lPtr:      make([]int32, k+1),
+		uPtr:      make([]int32, k+1),
 		work:      make([]float64, k),
 	}
 
@@ -143,13 +143,12 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64, w *luWS) *sparseLU
 				touched = append(touched, c)
 			}
 		}
-		uIdx := make([]int32, len(touched))
-		uVal := make([]float64, len(touched))
-		for t, c := range touched {
-			uIdx[t] = c // remapped to steps after the loop
-			uVal[t] = acc[c]
+		lu.lPtr[step] = int32(len(lu.lFlatIdx))
+		lu.uPtr[step] = int32(len(lu.uFlatIdx))
+		for _, c := range touched {
+			lu.uFlatIdx = append(lu.uFlatIdx, c) // remapped to steps below
+			lu.uFlatVal = append(lu.uFlatVal, acc[c])
 		}
-		lu.uIdx[step], lu.uVal[step] = uIdx, uVal
 
 		// eliminate the pivot column from every remaining row
 		for rr := range k {
@@ -167,8 +166,8 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64, w *luWS) *sparseLU
 				continue
 			}
 			mult := hit / pv
-			lu.lIdx[step] = append(lu.lIdx[step], int32(rr))
-			lu.lVal[step] = append(lu.lVal[step], mult)
+			lu.lFlatIdx = append(lu.lFlatIdx, int32(rr))
+			lu.lFlatVal = append(lu.lFlatVal, mult)
 
 			// rewrite the row in place: the write index never passes the
 			// read index, and fill-in only appends after the scan
@@ -205,10 +204,10 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64, w *luWS) *sparseLU
 			acc[c] = 0
 		}
 	}
-	for s := range k {
-		for t, c := range lu.uIdx[s] {
-			lu.uIdx[s][t] = colPos[c]
-		}
+	lu.lPtr[k] = int32(len(lu.lFlatIdx))
+	lu.uPtr[k] = int32(len(lu.uFlatIdx))
+	for t := range lu.uFlatIdx {
+		lu.uFlatIdx[t] = colPos[lu.uFlatIdx[t]]
 	}
 	return lu
 }
@@ -222,15 +221,15 @@ func (lu *sparseLU) solve(v []float64) {
 		ys := v[lu.rowPerm[s]]
 		y[s] = ys
 		if ys != 0 {
-			for t, rr := range lu.lIdx[s] {
-				v[rr] -= lu.lVal[s][t] * ys
+			for t := lu.lPtr[s]; t < lu.lPtr[s+1]; t++ {
+				v[lu.lFlatIdx[t]] -= lu.lFlatVal[t] * ys
 			}
 		}
 	}
 	for s := k - 1; s >= 0; s-- {
 		x := y[s]
-		for t, cs := range lu.uIdx[s] {
-			x -= lu.uVal[s][t] * y[cs]
+		for t := lu.uPtr[s]; t < lu.uPtr[s+1]; t++ {
+			x -= lu.uFlatVal[t] * y[lu.uFlatIdx[t]]
 		}
 		y[s] = x / lu.uDiag[s]
 	}
@@ -249,16 +248,16 @@ func (lu *sparseLU) solveT(w []float64) {
 		x := w[lu.colPerm[s]] / lu.uDiag[s]
 		y[s] = x
 		if x != 0 {
-			for t, cs := range lu.uIdx[s] {
-				w[lu.colPerm[cs]] -= lu.uVal[s][t] * x
+			for t := lu.uPtr[s]; t < lu.uPtr[s+1]; t++ {
+				w[lu.colPerm[lu.uFlatIdx[t]]] -= lu.uFlatVal[t] * x
 			}
 		}
 	}
 	// L^T backward, gathering from later steps
 	for s := k - 1; s >= 0; s-- {
 		x := y[s]
-		for t, rr := range lu.lIdx[s] {
-			x -= lu.lVal[s][t] * y[lu.stepOfRow[rr]]
+		for t := lu.lPtr[s]; t < lu.lPtr[s+1]; t++ {
+			x -= lu.lFlatVal[t] * y[lu.stepOfRow[lu.lFlatIdx[t]]]
 		}
 		y[s] = x
 	}

simplex/sparselu_bench_test.go

@@ -2,10 +2,51 @@ package simplex
 
 import (
 	"encoding/gob"
+	"math/rand"
 	"os"
 	"testing"
 )
 
+// synthKernel builds a random diagonally-dominant sparse k*k kernel as sparse
+// columns; representative of a mid-size refactorized basis block.
+func synthKernel(k int, density float64) ([][]int32, [][]float64) {
+	rng := rand.New(rand.NewSource(7))
+	colIdx := make([][]int32, k)
+	colVal := make([][]float64, k)
+	for c := range k {
+		for r := range k {
+			if r == c {
+				colIdx[c] = append(colIdx[c], int32(r))
+				colVal[c] = append(colVal[c], float64(k))
+			} else if rng.Float64() < density {
+				colIdx[c] = append(colIdx[c], int32(r))
+				colVal[c] = append(colVal[c], rng.NormFloat64())
+			}
+		}
+	}
+	return colIdx, colVal
+}
+
+// BenchmarkLUFactorizeSynth: self-contained factor + solve + solveT on a
+// synthetic kernel; -benchmem tracks the per-factorize allocation count.
+func BenchmarkLUFactorizeSynth(b *testing.B) {
+	const k = 300
+	colIdx, colVal := synthKernel(k, 0.03)
+	w := &luWS{}
+	v := make([]float64, k)
+	for b.Loop() {
+		lu := luFactorize(k, colIdx, colVal, w)
+		if lu == nil {
+			b.Fatal("singular")
+		}
+		for j := range v {
+			v[j] = float64(j%17) - 8
+		}
+		lu.solve(v)
+		lu.solveT(v)
+	}
+}
+
 func loadKernel(b *testing.B) (int, [][]int32, [][]float64) {
 	path := os.Getenv("SOLVER_KERNEL_GOB")
 	if path == "" {