diff --git a/CompPoly.lean b/CompPoly.lean
index 58d6f008..4c5463fe 100644
--- a/CompPoly.lean
+++ b/CompPoly.lean
@@ -109,6 +109,7 @@ import CompPoly.LinearAlgebra.Dense
 import CompPoly.LinearAlgebra.Dense.Basic
 import CompPoly.LinearAlgebra.Dense.Kernel
 import CompPoly.LinearAlgebra.Dense.KernelCorrectness
+import CompPoly.LinearAlgebra.Dense.KernelInPlaceCorrectness
 import CompPoly.LinearAlgebra.Dense.RowOps
 import CompPoly.LinearAlgebra.Dense.RowOpsCorrectness
 import CompPoly.LinearAlgebra.Dense.RrefSemantics
diff --git a/CompPoly/Bivariate/GuruswamiSudan/Context.lean b/CompPoly/Bivariate/GuruswamiSudan/Context.lean
index e8ef4b01..66fe50c2 100644
--- a/CompPoly/Bivariate/GuruswamiSudan/Context.lean
+++ b/CompPoly/Bivariate/GuruswamiSudan/Context.lean
@@ -80,21 +80,31 @@ structure LinearKernelContext (F : Type*) [Field F] [BEq F] [LawfulBEq F] where
               DenseMatrix.IsHomogeneousSolution M v →
                 ¬ (DenseMatrix.NonzeroVector v)
 
-/-- The dense Gaussian-elimination homogeneous-kernel backend. -/
+/-- The dense Gaussian-elimination homogeneous-kernel backend.
+
+The executable witness is computed by the in-place reduction
+`DenseMatrix.homogeneousWitnessInPlace`, which mutates the matrix array instead of
+copying it on every row operation. It returns the same witness as the copying
+`DenseMatrix.homogeneousWitness` (`DenseMatrix.homogeneousWitnessInPlace_eq`), so the
+backend contract is discharged by the copying-kernel correctness lemmas. -/
 def denseLinearKernelContext (F : Type*) [Field F] [BEq F] [LawfulBEq F] :
     LinearKernelContext F where
-  homogeneousWitness := DenseMatrix.homogeneousWitness
+  homogeneousWitness := DenseMatrix.homogeneousWitnessInPlace
   witness_width := by
     intro M v h
+    rw [DenseMatrix.homogeneousWitnessInPlace_eq] at h
     exact DenseMatrix.homogeneousWitness_width h
   witness_sound := by
     intro M v hM h
+    rw [DenseMatrix.homogeneousWitnessInPlace_eq] at h
     exact (DenseMatrix.homogeneousWitness_sound hM h).2.1
   witness_nonzero := by
     intro M v h
+    rw [DenseMatrix.homogeneousWitnessInPlace_eq] at h
     exact DenseMatrix.homogeneousWitness_nonzero h
   witness_complete := by
     intro M hM h v hvw hv
+    rw [DenseMatrix.homogeneousWitnessInPlace_eq] at h
     exact DenseMatrix.homogeneousWitness_none_complete hM h v hvw hv
 
 /-- Guruswami-Sudan-facing interpolation backend.
diff --git a/CompPoly/Bivariate/GuruswamiSudan/Interpolation/Dense/Correctness.lean b/CompPoly/Bivariate/GuruswamiSudan/Interpolation/Dense/Correctness.lean
index b89e16cc..aba86a90 100644
--- a/CompPoly/Bivariate/GuruswamiSudan/Interpolation/Dense/Correctness.lean
+++ b/CompPoly/Bivariate/GuruswamiSudan/Interpolation/Dense/Correctness.lean
@@ -901,7 +901,7 @@ theorem denseInterpolateWithBasis_exists_of_dimension_slack {F : Type*}
   have hnz := DenseMatrix.homogeneousWitness_nonzero hw
   rcases normalizeVector?_some_of_nonzero hnz with ⟨norm, hn⟩
   refine ⟨interpolationPolynomialOnBasis basis norm, ?_⟩
-  simp [denseLinearKernelContext, hw]
+  simp [denseLinearKernelContext, DenseMatrix.homogeneousWitnessInPlace_eq, hw]
   change normalizeInterpolationPolynomialOnBasis? basis coeffs =
     some (interpolationPolynomialOnBasis basis norm)
   unfold normalizeInterpolationPolynomialOnBasis?
diff --git a/CompPoly/LinearAlgebra/Dense.lean b/CompPoly/LinearAlgebra/Dense.lean
index a55c7ee5..f9ae8402 100644
--- a/CompPoly/LinearAlgebra/Dense.lean
+++ b/CompPoly/LinearAlgebra/Dense.lean
@@ -10,7 +10,9 @@ import CompPoly.LinearAlgebra.Dense.RowOpsCorrectness
 import CompPoly.LinearAlgebra.Dense.RrefSemantics
 import CompPoly.LinearAlgebra.Dense.RrefShape
 import CompPoly.LinearAlgebra.Dense.Kernel
+import CompPoly.LinearAlgebra.Dense.KernelInPlace
 import CompPoly.LinearAlgebra.Dense.KernelCorrectness
+import CompPoly.LinearAlgebra.Dense.KernelInPlaceCorrectness
 
 /-!
 # Dense Linear Algebra
diff --git a/CompPoly/LinearAlgebra/Dense/KernelInPlace.lean b/CompPoly/LinearAlgebra/Dense/KernelInPlace.lean
new file mode 100644
index 00000000..0ed96b17
--- /dev/null
+++ b/CompPoly/LinearAlgebra/Dense/KernelInPlace.lean
@@ -0,0 +1,146 @@
+/-
+Copyright (c) 2026 CompPoly Contributors. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Valerii Huhnin
+-/
+
+import CompPoly.LinearAlgebra.Dense.Kernel
+
+/-!
+# In-Place Dense Homogeneous Kernels
+
+A memory-efficient executable variant of the dense Gauss-Jordan homogeneous-kernel
+backend.
+
+The functions in `RowOps.lean` and `Kernel.lean` thread a whole `DenseMatrix`
+structure through every row operation. Because `scaleRow`/`addScaledRow`/`swapRows`
+read the matrix through `M.get` while the structure `M` is still live, the backing
+`Array` is multiply-referenced at the first `setD`, so `Array.setIfInBounds` copies
+the entire array on every row-operation call. Gauss-Jordan elimination performs
+`Θ(rows)` row operations per pivot and `Θ(min rows cols)` pivots, so the copies turn
+an `O(n³)` algorithm into `O(n⁴)` work plus `Θ(n⁴)` words of short-lived allocation.
+
+This module performs the same arithmetic, in the same order, threading a bare
+`Array F` (with the dimensions carried as separate `Nat` arguments) through the
+reduction. The array is uniquely referenced as it flows from one operation to the
+next, so `setIfInBounds` mutates in place after at most one fork from a shared input.
+The reduced matrix and pivot data are bit-for-bit identical to `rref`, so the
+extracted witnesses agree with `homogeneousWitness`.
+
+The free-column witness extraction (`freeColumns`, `basisVectorForFreeColumn`) is
+cheap and shared verbatim with `Kernel.lean`.
+-/
+
+namespace CompPoly
+
+namespace DenseMatrix
+
+variable {F : Type*}
+
+/-- Swap rows `rowA` and `rowB` of a width-`cols` row-major array, in place. -/
+def swapRowsData [Zero F] (cols rowA rowB : Nat) (data : Array F) : Array F :=
+  if rowA = rowB then
+    data
+  else
+    Id.run do
+      let mut data := data
+      for col in [0:cols] do
+        let ia := rowA * cols + col
+        let ib := rowB * cols + col
+        let a := data.getD ia 0
+        let b := data.getD ib 0
+        data := setD data ia b
+        data := setD data ib a
+      pure data
+
+/-- Scale row `row` of a width-`cols` row-major array by `factor`, in place. -/
+def scaleRowData [Field F] (cols row : Nat) (factor : F) (data : Array F) :
+    Array F :=
+  Id.run do
+    let mut data := data
+    for col in [0:cols] do
+      let idx := row * cols + col
+      data := setD data idx (factor * data.getD idx 0)
+    pure data
+
+/-- Add `factor` times row `source` to row `target` (with `target ≠ source`), in
+place. -/
+def addScaledRowData [Field F] (cols target source : Nat) (factor : F)
+    (data : Array F) : Array F :=
+  Id.run do
+    let mut data := data
+    for col in [0:cols] do
+      let ti := target * cols + col
+      let si := source * cols + col
+      data := setD data ti (data.getD ti 0 + factor * data.getD si 0)
+    pure data
+
+/-- Find a pivot row at or below `startRow` in column `col`. -/
+def findPivotRowData [Zero F] [BEq F] (rows cols startRow col : Nat)
+    (data : Array F) : Option Nat :=
+  Id.run do
+    let mut out : Option Nat := none
+    for row in [startRow:rows] do
+      if out.isNone && !(data.getD (row * cols + col) 0 == 0) then
+        out := some row
+    pure out
+
+/-- Normalize the pivot row and clear the pivot column in every other row, in
+place. -/
+def normalizeAndEliminateData [Field F] [BEq F] (rows cols pivotRow pivotCol : Nat)
+    (data : Array F) : Array F :=
+  let pivot := data.getD (pivotRow * cols + pivotCol) 0
+  if pivot == 0 then
+    data
+  else
+    Id.run do
+      let mut data := scaleRowData cols pivotRow pivot⁻¹ data
+      for row in [0:rows] do
+        if row != pivotRow then
+          let factor := -(data.getD (row * cols + pivotCol) 0)
+          data := addScaledRowData cols row pivotRow factor data
+      pure data
+
+/-- Fuel-bounded Gauss-Jordan reduction over a bare row-major array. -/
+def rrefLoopData [Field F] [BEq F] :
+    Nat → Nat → Nat → Nat → Nat → Array F → Array Nat → Array F × Array Nat
+  | 0, _rows, _cols, _col, _row, data, pivots => (data, pivots)
+  | fuel + 1, rows, cols, col, row, data, pivots =>
+      if col < cols && row < rows then
+        match findPivotRowData rows cols row col data with
+        | none => rrefLoopData fuel rows cols (col + 1) row data pivots
+        | some pivotRow =>
+            let data := swapRowsData cols pivotRow row data
+            let data := normalizeAndEliminateData rows cols row col data
+            rrefLoopData fuel rows cols (col + 1) (row + 1) data (pivots.push col)
+      else
+        (data, pivots)
+
+/-- Reduced row-echelon form with pivot-column metadata, computed in place.
+
+This destructures `M` so the backing array is owned by the reduction loop rather
+than aliased through the structure, which is what lets `setIfInBounds` mutate in
+place. The result equals `rref M`. -/
+def rrefInPlace [Field F] [BEq F] (M : DenseMatrix F) : RrefResult F :=
+  match M with
+  | { rows := rows, cols := cols, data := data } =>
+      let (data, pivots) := rrefLoopData (cols + 1) rows cols 0 0 data #[]
+      { matrix := { rows := rows, cols := cols, data := data }, pivots := pivots }
+
+/-- Homogeneous kernel basis extracted from the in-place RREF free columns. -/
+def homogeneousKernelBasisInPlace [Field F] [BEq F] (M : DenseMatrix F) :
+    Array (Array F) :=
+  let R := rrefInPlace M
+  (freeColumns R.matrix.cols R.pivots).map
+    (basisVectorForFreeColumn R.matrix R.pivots)
+
+/-- One nonzero homogeneous-kernel witness via in-place reduction, if a free
+column exists. Output-identical to `homogeneousWitness`. -/
+def homogeneousWitnessInPlace [Field F] [BEq F] (M : DenseMatrix F) :
+    Option (Array F) :=
+  let basis := homogeneousKernelBasisInPlace M
+  if basis.size = 0 then none else some (basis.getD 0 #[])
+
+end DenseMatrix
+
+end CompPoly
diff --git a/CompPoly/LinearAlgebra/Dense/KernelInPlaceCorrectness.lean b/CompPoly/LinearAlgebra/Dense/KernelInPlaceCorrectness.lean
new file mode 100644
index 00000000..82d21b0c
--- /dev/null
+++ b/CompPoly/LinearAlgebra/Dense/KernelInPlaceCorrectness.lean
@@ -0,0 +1,271 @@
+/-
+Copyright (c) 2026 CompPoly Contributors. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Valerii Huhnin
+-/
+
+import CompPoly.LinearAlgebra.Dense.KernelInPlace
+import CompPoly.LinearAlgebra.Dense.RowOpsCorrectness
+import CompPoly.LinearAlgebra.Dense.RrefSemantics
+
+/-!
+# In-Place Dense Kernel Correctness
+
+The in-place homogeneous-kernel backend (`KernelInPlace.lean`) performs the same
+arithmetic, in the same order, as the copying backend (`RowOps.lean`/`Kernel.lean`),
+so it returns the same witness. This file proves that equivalence, culminating in
+`homogeneousWitnessInPlace_eq`, which lets the certified Guruswami-Sudan backend run
+the in-place implementation while reusing every existing correctness theorem.
+-/
+
+namespace CompPoly
+
+namespace DenseMatrix
+
+variable {F : Type*}
+
+/-- `setD` at a different index leaves the queried entry unchanged. -/
+private theorem getD_setD_ne [Zero F] (d : Array F) {i j : Nat} (h : i ≠ j) (v : F) :
+    (setD d i v).getD j 0 = d.getD j 0 := by
+  unfold setD
+  rw [Array.getD_eq_getD_getElem?, Array.getD_eq_getD_getElem?,
+    Array.getElem?_setIfInBounds_ne h]
+
+/-- A self-reading scale fold equals the constant-reading scale fold, as long as
+the accumulator still agrees with the original array on the entries it reads. -/
+private theorem foldl_scale_self_eq_const [Field F]
+    (cols r : Nat) (factor : F) (data₀ : Array F) :
+    ∀ (xs : List Nat), xs.Nodup →
+      ∀ (d : Array F),
+        (∀ c, c ∈ xs → d.getD (r * cols + c) 0 = data₀.getD (r * cols + c) 0) →
+          List.foldl (fun d c ↦ setD d (r * cols + c) (factor * d.getD (r * cols + c) 0)) d xs =
+            List.foldl
+              (fun d c ↦ setD d (r * cols + c) (factor * data₀.getD (r * cols + c) 0)) d xs := by
+  intro xs
+  induction xs with
+  | nil => intro _ d _; rfl
+  | cons c xs ih =>
+      intro hnodup d hagree
+      simp only [List.foldl_cons]
+      rw [hagree c (List.mem_cons_self ..)]
+      apply ih hnodup.tail
+      intro c' hc'
+      have hne : c ≠ c' := fun h ↦ (List.nodup_cons.mp hnodup).1 (h ▸ hc')
+      rw [getD_setD_ne _ (by omega)]
+      exact hagree c' (List.mem_cons_of_mem _ hc')
+
+/-- The in-place scale of a row produces the same backing array as the copying scale. -/
+theorem scaleRowData_eq [Field F] (M : DenseMatrix F) (r : Nat) (factor : F) :
+    scaleRowData M.cols r factor M.data = (scaleRow M r factor).data := by
+  have h := foldl_scale_self_eq_const M.cols r factor M.data (List.range' 0 M.cols)
+    (List.nodup_range' (s := 0) (n := M.cols)) M.data (fun c _ ↦ rfl)
+  simpa [scaleRowData, scaleRow, get, index, Std.Legacy.Range.forIn_eq_forIn_range'] using h
+
+/-- Row-major indices in distinct rows are distinct (within the column range). -/
+private theorem row_col_index_ne {a b cols c c' : Nat} (hrow : a ≠ b)
+    (hc : c < cols) (hc' : c' < cols) : a * cols + c ≠ b * cols + c' := by
+  intro h
+  have hcols : 0 < cols := Nat.lt_of_le_of_lt (Nat.zero_le c) hc
+  apply hrow
+  have hdiv := congrArg (· / cols) h
+  simpa [Nat.mul_comm a cols, Nat.mul_comm b cols, Nat.mul_add_div hcols,
+    Nat.div_eq_of_lt hc, Nat.div_eq_of_lt hc'] using hdiv
+
+/-- A self-reading add-scaled fold equals the constant-reading version, as long as the
+accumulator still agrees with the original array on the target and source entries. -/
+private theorem foldl_addScaled_self_eq_const [Field F]
+    (cols target source : Nat) (factor : F) (data₀ : Array F) (hts : target ≠ source) :
+    ∀ (xs : List Nat), xs.Nodup → (∀ c, c ∈ xs → c < cols) →
+      ∀ (d : Array F),
+        (∀ c, c ∈ xs → d.getD (target * cols + c) 0 = data₀.getD (target * cols + c) 0) →
+        (∀ c, c ∈ xs → d.getD (source * cols + c) 0 = data₀.getD (source * cols + c) 0) →
+          List.foldl (fun d c ↦ setD d (target * cols + c)
+              (d.getD (target * cols + c) 0 + factor * d.getD (source * cols + c) 0)) d xs =
+            List.foldl (fun d c ↦ setD d (target * cols + c)
+              (data₀.getD (target * cols + c) 0 + factor * data₀.getD (source * cols + c) 0))
+              d xs := by
+  intro xs
+  induction xs with
+  | nil => intro _ _ d _ _; rfl
+  | cons c xs ih =>
+      intro hnodup hbounds d hT hS
+      have hcc : c < cols := hbounds c (List.mem_cons_self ..)
+      simp only [List.foldl_cons]
+      rw [hT c (List.mem_cons_self ..), hS c (List.mem_cons_self ..)]
+      apply ih hnodup.tail (fun c' hc' ↦ hbounds c' (List.mem_cons_of_mem _ hc'))
+      · intro c' hc'
+        have hne : c ≠ c' := fun h ↦ (List.nodup_cons.mp hnodup).1 (h ▸ hc')
+        rw [getD_setD_ne _ (by omega)]
+        exact hT c' (List.mem_cons_of_mem _ hc')
+      · intro c' hc'
+        have hcc' : c' < cols := hbounds c' (List.mem_cons_of_mem _ hc')
+        rw [getD_setD_ne _ (row_col_index_ne hts hcc hcc')]
+        exact hS c' (List.mem_cons_of_mem _ hc')
+
+/-- The in-place add-scaled row produces the same backing array as the copying version,
+when the target and source rows are distinct. -/
+theorem addScaledRowData_eq [Field F] (M : DenseMatrix F) (target source : Nat)
+    (factor : F) (hts : target ≠ source) :
+    addScaledRowData M.cols target source factor M.data =
+      (addScaledRow M target source factor).data := by
+  have h := foldl_addScaled_self_eq_const M.cols target source factor M.data hts
+    (List.range' 0 M.cols) (List.nodup_range' (s := 0) (n := M.cols))
+    (fun c hc ↦ by simpa using (List.mem_range'_1.mp hc).2)
+    M.data (fun c _ ↦ rfl) (fun c _ ↦ rfl)
+  simpa [addScaledRowData, addScaledRow, get, index,
+    Std.Legacy.Range.forIn_eq_forIn_range'] using h
+
+/-- A self-reading swap fold equals the constant-reading version, as long as the
+accumulator still agrees with the original array on the two swapped rows. -/
+private theorem foldl_swap_self_eq_const [Zero F]
+    (cols rowA rowB : Nat) (data₀ : Array F) (hAB : rowA ≠ rowB) :
+    ∀ (xs : List Nat), xs.Nodup → (∀ c, c ∈ xs → c < cols) →
+      ∀ (d : Array F),
+        (∀ c, c ∈ xs → d.getD (rowA * cols + c) 0 = data₀.getD (rowA * cols + c) 0) →
+        (∀ c, c ∈ xs → d.getD (rowB * cols + c) 0 = data₀.getD (rowB * cols + c) 0) →
+          List.foldl (fun d c ↦ setD (setD d (rowA * cols + c) (d.getD (rowB * cols + c) 0))
+              (rowB * cols + c) (d.getD (rowA * cols + c) 0)) d xs =
+            List.foldl (fun d c ↦ setD (setD d (rowA * cols + c) (data₀.getD (rowB * cols + c) 0))
+              (rowB * cols + c) (data₀.getD (rowA * cols + c) 0)) d xs := by
+  intro xs
+  induction xs with
+  | nil => intro _ _ d _ _; rfl
+  | cons c xs ih =>
+      intro hnodup hbounds d hA hB
+      have hcc : c < cols := hbounds c (List.mem_cons_self ..)
+      simp only [List.foldl_cons]
+      rw [hA c (List.mem_cons_self ..), hB c (List.mem_cons_self ..)]
+      apply ih hnodup.tail (fun c' hc' ↦ hbounds c' (List.mem_cons_of_mem _ hc'))
+      · intro c' hc'
+        have hne : c ≠ c' := fun h ↦ (List.nodup_cons.mp hnodup).1 (h ▸ hc')
+        have hcc' : c' < cols := hbounds c' (List.mem_cons_of_mem _ hc')
+        rw [getD_setD_ne _ (row_col_index_ne hAB.symm hcc hcc'), getD_setD_ne _ (by omega)]
+        exact hA c' (List.mem_cons_of_mem _ hc')
+      · intro c' hc'
+        have hne : c ≠ c' := fun h ↦ (List.nodup_cons.mp hnodup).1 (h ▸ hc')
+        have hcc' : c' < cols := hbounds c' (List.mem_cons_of_mem _ hc')
+        rw [getD_setD_ne _ (by omega), getD_setD_ne _ (row_col_index_ne hAB hcc hcc')]
+        exact hB c' (List.mem_cons_of_mem _ hc')
+
+/-- The in-place row swap produces the same backing array as the copying swap. -/
+theorem swapRowsData_eq [Zero F] (M : DenseMatrix F) (rowA rowB : Nat) :
+    swapRowsData M.cols rowA rowB M.data = (swapRows M rowA rowB).data := by
+  by_cases hAB : rowA = rowB
+  · simp [swapRowsData, swapRows, hAB]
+  · have h := foldl_swap_self_eq_const M.cols rowA rowB M.data hAB
+      (List.range' 0 M.cols) (List.nodup_range' (s := 0) (n := M.cols))
+      (fun c hc ↦ by simpa using (List.mem_range'_1.mp hc).2)
+      M.data (fun c _ ↦ rfl) (fun c _ ↦ rfl)
+    simpa [swapRowsData, swapRows, hAB, get, index,
+      Std.Legacy.Range.forIn_eq_forIn_range'] using h
+
+/-- The in-place pivot search agrees with the copying one. -/
+theorem findPivotRowData_eq [Zero F] [BEq F] (M : DenseMatrix F) (startRow col : Nat) :
+    findPivotRowData M.rows M.cols startRow col M.data = findPivotRow M startRow col := by
+  simp only [findPivotRowData, findPivotRow, get, index]
+
+/-- The in-place elimination loop tracks the copying one entry-for-entry. -/
+private theorem forIn_eliminate_eq [Field F] (cols pivotRow pivotCol : Nat) :
+    ∀ (xs : List Nat) (out : DenseMatrix F), out.cols = cols →
+      (Id.run (forIn xs out.data fun row r ↦
+          if row = pivotRow then pure (ForInStep.yield r)
+          else pure (ForInStep.yield
+            (addScaledRowData cols row pivotRow (-(r.getD (row * cols + pivotCol) 0)) r))))
+        = (Id.run (forIn xs out fun row r ↦
+            if row = pivotRow then pure (ForInStep.yield r)
+            else pure (ForInStep.yield
+              (addScaledRow r row pivotRow (-(r.get row pivotCol)))))).data := by
+  intro xs
+  induction xs with
+  | nil => intro _ _; rfl
+  | cons row xs ih =>
+      intro out hcols
+      by_cases hr : row = pivotRow
+      · simpa [List.forIn_cons, hr] using ih out hcols
+      · have hget : out.data.getD (row * cols + pivotCol) 0 = out.get row pivotCol := by
+          simp [get, index, hcols]
+        have hbridge := addScaledRowData_eq out row pivotRow (-(out.get row pivotCol)) hr
+        rw [hcols] at hbridge
+        have hcols' : (addScaledRow out row pivotRow (-(out.get row pivotCol))).cols = cols := by
+          rw [← hcols]; rfl
+        have hih := ih (addScaledRow out row pivotRow (-(out.get row pivotCol))) hcols'
+        rw [← hbridge] at hih
+        simpa [List.forIn_cons, hr, ← Array.getD_eq_getD_getElem?, hget] using hih
+
+/-- The in-place normalize-and-eliminate step produces the same backing array. -/
+theorem normalizeAndEliminateData_eq [Field F] [BEq F] (M : DenseMatrix F)
+    (pivotRow pivotCol : Nat) :
+    normalizeAndEliminateData M.rows M.cols pivotRow pivotCol M.data
+      = (normalizeAndEliminate M pivotRow pivotCol).data := by
+  have hpiv : M.data.getD (pivotRow * M.cols + pivotCol) 0 = M.get pivotRow pivotCol := by
+    simp [get, index]
+  unfold normalizeAndEliminateData normalizeAndEliminate
+  rw [hpiv]
+  by_cases hp : M.get pivotRow pivotCol == 0
+  · simp [hp]
+  · simp only [hp, if_false, Bool.false_eq_true]
+    have hstart : scaleRowData M.cols pivotRow (M.get pivotRow pivotCol)⁻¹ M.data
+        = (scaleRow M pivotRow (M.get pivotRow pivotCol)⁻¹).data :=
+      scaleRowData_eq M pivotRow _
+    have hcols : (scaleRow M pivotRow (M.get pivotRow pivotCol)⁻¹).cols = M.cols := by
+      rw [scaleRow]
+    have heli := forIn_eliminate_eq M.cols pivotRow pivotCol
+      (List.range' 0 (Std.Legacy.Range.size [:M.rows]))
+      (scaleRow M pivotRow (M.get pivotRow pivotCol)⁻¹) hcols
+    rw [hstart]
+    simpa [Std.Legacy.Range.forIn_eq_forIn_range'] using heli
+
+/-- The in-place reduction loop tracks the copying reduction loop, returning the
+same reduced backing array and the same pivot columns. -/
+theorem rrefLoopData_eq [Field F] [BEq F] :
+    ∀ (fuel col row : Nat) (M : DenseMatrix F) (pivots : Array Nat),
+      rrefLoopData fuel M.rows M.cols col row M.data pivots
+        = ((rrefLoop fuel col row M pivots).matrix.data,
+           (rrefLoop fuel col row M pivots).pivots) := by
+  intro fuel
+  induction fuel with
+  | zero => intro col row M pivots; rfl
+  | succ fuel ih =>
+      intro col row M pivots
+      rw [rrefLoopData, rrefLoop, findPivotRowData_eq]
+      by_cases hguard : col < M.cols && row < M.rows
+      · simp only [hguard, if_true]
+        cases hpiv : findPivotRow M row col with
+        | none => simpa using ih (col + 1) row M pivots
+        | some pivotRow =>
+            simp only
+            rw [swapRowsData_eq]
+            rw [show M.rows = (swapRows M pivotRow row).rows from (swapRows_rows _ _ _).symm,
+                show M.cols = (swapRows M pivotRow row).cols from (swapRows_cols _ _ _).symm,
+                normalizeAndEliminateData_eq]
+            rw [show (swapRows M pivotRow row).rows
+                  = (normalizeAndEliminate (swapRows M pivotRow row) row col).rows from
+                  (normalizeAndEliminate_rows _ _ _).symm,
+                show (swapRows M pivotRow row).cols
+                  = (normalizeAndEliminate (swapRows M pivotRow row) row col).cols from
+                  (normalizeAndEliminate_cols _ _ _).symm]
+            exact ih (col + 1) (row + 1)
+              (normalizeAndEliminate (swapRows M pivotRow row) row col) (pivots.push col)
+      · simp only [hguard, if_false, Bool.false_eq_true]
+
+/-- The in-place RREF equals the copying RREF. -/
+theorem rrefInPlace_eq [Field F] [BEq F] (M : DenseMatrix F) : rrefInPlace M = rref M := by
+  obtain ⟨rows, cols, data⟩ := M
+  have key := rrefLoopData_eq (cols + 1) 0 0 (⟨rows, cols, data⟩ : DenseMatrix F) #[]
+  simp only at key
+  simp only [rrefInPlace, rref]
+  rw [key]
+  set N := rrefLoop (cols + 1) 0 0 (⟨rows, cols, data⟩ : DenseMatrix F) #[] with hN
+  have hr : N.matrix.rows = rows := rrefLoop_rows _ _ _ _ _
+  have hc : N.matrix.cols = cols := rrefLoop_cols _ _ _ _ _
+  rw [← hr, ← hc]
+
+/-- The in-place homogeneous-kernel witness equals the copying witness. -/
+theorem homogeneousWitnessInPlace_eq [Field F] [BEq F] (M : DenseMatrix F) :
+    homogeneousWitnessInPlace M = homogeneousWitness M := by
+  simp only [homogeneousWitnessInPlace, homogeneousWitness,
+    homogeneousKernelBasisInPlace, homogeneousKernelBasis, rrefInPlace_eq, rref_cols]
+
+end DenseMatrix
+
+end CompPoly
diff --git a/bench/CompPolyBench/Bivariate/GuruswamiSudan.lean b/bench/CompPolyBench/Bivariate/GuruswamiSudan.lean
index 702d6cd6..7f6f9302 100644
--- a/bench/CompPolyBench/Bivariate/GuruswamiSudan.lean
+++ b/bench/CompPolyBench/Bivariate/GuruswamiSudan.lean
@@ -58,27 +58,41 @@ private def runGsInterpolationSolveKoala (preset : BenchPreset) (gen : StdGen) :
   let fastPoints := gsSmallBenchmarkPoints fastMessage
   let matrix := interpolationMatrix points gsSmallParams
   let fastMatrix := interpolationMatrix fastPoints gsSmallParams
-  let kernelContext := denseLinearKernelContext KoalaBear.Field
-  let fastKernelContext := denseLinearKernelContext KoalaBear.Fast.Field
   let warmup := gsWarmupIterations preset
   let measured := preset.selectNat 1 1 1
   let fastMeasured := preset.selectNat 2 1 1
-  let checksumIterations := groupChecksumIterations measured [fastMeasured]
+  let inPlaceMeasured := preset.selectNat 8 2 1
+  let fastInPlaceMeasured := preset.selectNat 16 3 1
+  let checksumIterations := groupChecksumIterations measured
+    [fastMeasured, inPlaceMeasured, fastInPlaceMeasured]
   let row <- runTimed
-    "guruswami-sudan-interp-solve" "DenseMatrix" "Homogeneous interpolation solve"
+    "guruswami-sudan-interp-solve-copying" "DenseMatrix"
+    "Homogeneous interpolation solve, copying"
     "KoalaBear.Field" gsSmallInterpInputShape preset warmup measured
-    (fun _ ↦ kernelContext.homogeneousWitness matrix)
+    (fun _ ↦ DenseMatrix.homogeneousWitness matrix)
+    (checksumOptionArray checksumKoalaBear) checksumIterations
+  let inPlaceRow <- runTimed
+    "guruswami-sudan-interp-solve" "DenseMatrix"
+    "Homogeneous interpolation solve, in-place"
+    "KoalaBear.Field" gsSmallInterpInputShape preset warmup inPlaceMeasured
+    (fun _ ↦ DenseMatrix.homogeneousWitnessInPlace matrix)
     (checksumOptionArray checksumKoalaBear) checksumIterations
   let fastRow <- runTimed
-    "guruswami-sudan-interp-solve-fast" "DenseMatrix"
-    "Homogeneous interpolation solve"
+    "guruswami-sudan-interp-solve-copying-fast" "DenseMatrix"
+    "Homogeneous interpolation solve, copying"
     "KoalaBear.Fast.Field" gsSmallInterpInputShape preset warmup fastMeasured
-    (fun _ ↦ fastKernelContext.homogeneousWitness fastMatrix)
+    (fun _ ↦ DenseMatrix.homogeneousWitness fastMatrix)
+    (checksumOptionArray checksumKoalaBearFast) checksumIterations
+  let fastInPlaceRow <- runTimed
+    "guruswami-sudan-interp-solve-inplace-fast" "DenseMatrix"
+    "Homogeneous interpolation solve, in-place"
+    "KoalaBear.Fast.Field" gsSmallInterpInputShape preset warmup fastInPlaceMeasured
+    (fun _ ↦ DenseMatrix.homogeneousWitnessInPlace fastMatrix)
     (checksumOptionArray checksumKoalaBearFast) checksumIterations
   pure ({
     groupKey := "guruswami-sudan-interp-solve-small-koalabear",
     title := "Guruswami-Sudan dense interpolation solving, small (KoalaBear)",
-    records := #[row, fastRow]
+    records := #[row, inPlaceRow, fastRow, fastInPlaceRow]
   }, gen)
 
 private def runGsInterpolationSmallKoala (preset : BenchPreset) (gen : StdGen) :
diff --git a/bench/CompPolyBench/Common.lean b/bench/CompPolyBench/Common.lean
index 0c885e5b..022cf51c 100644
--- a/bench/CompPolyBench/Common.lean
+++ b/bench/CompPolyBench/Common.lean
@@ -747,6 +747,8 @@ def implementationMethodLabels : List (String × String) := [
   ("Dense linear", "Dense linear"),
   ("Interpolation system construction", "System construction"),
   ("Homogeneous interpolation solve", "Dense solve"),
+  ("Homogeneous interpolation solve, copying", "Dense solve (copying)"),
+  ("Homogeneous interpolation solve, in-place", "Dense solve (in-place)"),
   ("Roth-Ruckenstein root finding with nonlinear field-root equations",
     "Roth-Ruckenstein"),
   ("Roth-Ruckenstein root finding with NTTFast field-root equations",