perf: reuse delay matrices across FreqResponsePointwise omega points

frequency.go

@@ -222,13 +222,21 @@ func (e frequencyEvaluator) responsePointwise(omega []float64) (*FreqResponseMat
 	pm := e.p * e.m
 	data := make([]complex128, len(omega)*pm)
 	ws := newSSEvalWorkspace(e.n, e.p, e.m)
+	// The effective delay matrices are pure functions of sys, which is
+	// not mutated across the sweep, so hoisting them out of the loop
+	// keeps every per-point value bit-identical while dropping the
+	// per-omega allocations applyIODelayAtS would repeat.
+	delaySS := effectiveIODelayMatrix(e.sys, e.p, e.m, true)
+	var delayTF *mat.Dense
 	var tf *TransferFunc
 	for k, w := range omega {
 		s := e.sAt(w)
 		dst := data[k*pm : (k+1)*pm]
 		if err := evalFrSSInto(ws, e.sys, s, e.n, e.p, e.m); err == nil {
 			copy(dst, ws.g[:pm])
-			applyIODelayAtS(e.sys, s, dst, e.p, e.m, true)
+			if delaySS != nil {
+				applyIODelayMatrixAtS(e.sys, s, dst, e.p, e.m, delaySS)
+			}
 			continue
 		}
 		if tf == nil {
@@ -237,9 +245,12 @@ func (e frequencyEvaluator) responsePointwise(omega []float64) (*FreqResponseMat
 				return nil, err
 			}
 			tf = res.TF
+			delayTF = effectiveIODelayMatrix(e.sys, e.p, e.m, false)
 		}
 		tf.evalInto(s, dst)
-		applyIODelayAtS(e.sys, s, dst, e.p, e.m, false)
+		if delayTF != nil {
+			applyIODelayMatrixAtS(e.sys, s, dst, e.p, e.m, delayTF)
+		}
 	}
 	return e.matrix(data, omega), nil
 }

frequency_pointwise_test.go

@@ -175,6 +175,43 @@ func TestFreqResponsePointwise_PerPointFallback(t *testing.T) {
 	assertBitIdenticalToSinglePoint(t, sys, omega)
 }
 
+// The pointwise sweep reuses one workspace (solve buffers and effective
+// delay matrices) across omega points, so allocations per sweep must stay
+// constant as the sweep grows: only the result buffer scales with
+// len(omega), and that is a single allocation. A per-point allocation
+// regression (like rebuilding the delay matrix inside the loop) makes the
+// long-sweep count scale with n and fails this test.
+func TestFreqResponsePointwise_AllocsAmortizeAcrossSweepLengths(t *testing.T) {
+	sys, err := New(
+		mat.NewDense(2, 2, []float64{-0.125, 0, 0.125, -1.0 / 3.0}),
+		mat.NewDense(2, 1, []float64{0.1875, 0}),
+		mat.NewDense(1, 2, []float64{0, 1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if err := sys.SetInputDelay([]float64{1.5}); err != nil {
+		t.Fatal(err)
+	}
+
+	allocsAt := func(n int) float64 {
+		omega := pointwiseTestGrid(0.001, 10.0, n)
+		return testing.AllocsPerRun(20, func() {
+			if _, err := sys.FreqResponsePointwise(omega); err != nil {
+				t.Fatal(err)
+			}
+		})
+	}
+
+	short := allocsAt(10)
+	long := allocsAt(1000)
+	if long > short {
+		t.Fatalf("allocations scale with sweep length: %v allocs at 10 points, %v at 1000 (want constant)", short, long)
+	}
+}
+
 func TestFreqResponsePointwise_EmptyOmega(t *testing.T) {
 	sys, err := New(
 		mat.NewDense(1, 1, []float64{-1}),