From 6ebf7fd02487a66ae8fc7f25cd986af1f0629367 Mon Sep 17 00:00:00 2001 From: Christian Woltering Date: Sun, 23 Aug 2026 13:10:31 +0200 Subject: [PATCH 1/2] Clear SparseLU workspace variables before refactorize. Fixes issue #58 --- .../Double/Factorization/SparseLUTest.cs | 21 +++++++++++++++++++ CSparse/Complex/Factorization/SparseLU.cs | 3 +++ CSparse/Double/Factorization/SparseLU.cs | 3 +++ 3 files changed, 27 insertions(+) diff --git a/CSparse.Tests/Double/Factorization/SparseLUTest.cs b/CSparse.Tests/Double/Factorization/SparseLUTest.cs index 477a1e3..c940af3 100644 --- a/CSparse.Tests/Double/Factorization/SparseLUTest.cs +++ b/CSparse.Tests/Double/Factorization/SparseLUTest.cs @@ -111,6 +111,27 @@ public void TestRefactorize() Assert.Throws(() => lu.Refactorize(small, 1.0)); } + [Test] + public void TestRefactorize2() + { + const int N = 3; + + var A = CompressedColumnStorage.OfRowMajor(N, N, [0.001, -0.001, 1, -0.001, 10.001, 0, 1, 0, 0]); + + double[] b = { 0.0, 0.0049995, 5.0 }; + + double[] x1 = new double[N]; + double[] x2 = new double[N]; + + var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1e-8); + + lu.Solve(b, x1); + lu.Refactorize(A, 1e-8); + lu.Solve(b, x2); + + Assert.AreEqual(x1[0], x2[0]); + } + [Test] public void TestRefactorizeNoTrim() { diff --git a/CSparse/Complex/Factorization/SparseLU.cs b/CSparse/Complex/Factorization/SparseLU.cs index 4a7fe00..b203316 100644 --- a/CSparse/Complex/Factorization/SparseLU.cs +++ b/CSparse/Complex/Factorization/SparseLU.cs @@ -137,6 +137,9 @@ public void Refactorize(CompressedColumnStorage A, double tol = 1.0) // Ensure tol is in range. tol = Math.Min(Math.Max(tol, 0.0), 1.0); + // Reset workspace + Array.Clear(temp, 0, n); + // Reuse the cached symbolic ordering (S); recompute L, U and the pivoting. Factorize(A, tol, null); } diff --git a/CSparse/Double/Factorization/SparseLU.cs b/CSparse/Double/Factorization/SparseLU.cs index 0854bf7..131e232 100644 --- a/CSparse/Double/Factorization/SparseLU.cs +++ b/CSparse/Double/Factorization/SparseLU.cs @@ -136,6 +136,9 @@ public void Refactorize(CompressedColumnStorage A, double tol = 1.0) // Ensure tol is in range. tol = Math.Min(Math.Max(tol, 0.0), 1.0); + // Reset workspace + Array.Clear(temp, 0, n); + // Reuse the cached symbolic ordering (S); recompute L, U and the pivoting. Factorize(A, tol, null); } From 219d4c7a628f0e6a1cd7aa7c891b3198cee97180 Mon Sep 17 00:00:00 2001 From: Christian Woltering Date: Sun, 23 Aug 2026 13:22:22 +0200 Subject: [PATCH 2/2] SparseLU (re)factorization should not allocate temporary workspace. --- CSparse/Complex/Factorization/SparseLU.cs | 12 ++++++++---- CSparse/Double/Factorization/SparseLU.cs | 14 +++++++++----- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/CSparse/Complex/Factorization/SparseLU.cs b/CSparse/Complex/Factorization/SparseLU.cs index b203316..132cdd4 100644 --- a/CSparse/Complex/Factorization/SparseLU.cs +++ b/CSparse/Complex/Factorization/SparseLU.cs @@ -21,7 +21,8 @@ public class SparseLU : ISparseFactorization CompressedColumnStorage L, U; int[] pinv; // partial pivoting - Complex[] temp; // workspace + readonly Complex[] temp; // workspace (used for factorization and solve) + readonly int[] temp2; // workspace (used for factorization) #region Static methods @@ -99,7 +100,9 @@ public static SparseLU Create(CompressedColumnStorage A, int[] p, doubl private SparseLU(int n) { this.n = n; - this.temp = new Complex[n]; + + temp = new Complex[n]; + temp2 = new int[2 * n]; } /// @@ -139,6 +142,7 @@ public void Refactorize(CompressedColumnStorage A, double tol = 1.0) // Reset workspace Array.Clear(temp, 0, n); + Array.Clear(temp2, 0, 2 * n); // Reuse the cached symbolic ordering (S); recompute L, U and the pivoting. Factorize(A, tol, null); @@ -227,8 +231,8 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress } // Workspace - var x = this.temp; - var xi = new int[2 * n]; + var x = temp; + var xi = temp2; for (i = 0; i < n; i++) { diff --git a/CSparse/Double/Factorization/SparseLU.cs b/CSparse/Double/Factorization/SparseLU.cs index 131e232..d568d65 100644 --- a/CSparse/Double/Factorization/SparseLU.cs +++ b/CSparse/Double/Factorization/SparseLU.cs @@ -20,8 +20,9 @@ public class SparseLU : ISparseFactorization CompressedColumnStorage L, U; int[] pinv; // partial pivoting - double[] temp; // workspace - + readonly double[] temp; // workspace (used for factorization and solve) + readonly int[] temp2; // workspace (used for factorization) + #region Static methods /// @@ -98,7 +99,9 @@ public static SparseLU Create(CompressedColumnStorage A, int[] p, double private SparseLU(int n) { this.n = n; - this.temp = new double[n]; + + temp = new double[n]; + temp2 = new int[2 * n]; } /// @@ -138,6 +141,7 @@ public void Refactorize(CompressedColumnStorage A, double tol = 1.0) // Reset workspace Array.Clear(temp, 0, n); + Array.Clear(temp2, 0, 2 * n); // Reuse the cached symbolic ordering (S); recompute L, U and the pivoting. Factorize(A, tol, null); @@ -226,8 +230,8 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress< } // Workspace - var x = this.temp; - var xi = new int[2 * n]; + var x = temp; + var xi = temp2; for (i = 0; i < n; i++) {