From ef2d4fd971e26adb7644cdede9925daf4d093332 Mon Sep 17 00:00:00 2001 From: "Dr. H. BADI" Date: Mon, 8 Jun 2026 13:59:16 +0200 Subject: [PATCH] Allocate factor buffers only on the first Factorize; reuse them on refactorization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SparseLU / SparseCholesky / SparseLDL (Double and Complex) allocate L/U (and pinv, D) only on the first Factorize call, and Clear() + reuse the existing buffers on subsequent calls (same sparsity pattern), avoiding re-allocation on every Refactorize. For SparseLU the end-of-factorization trim (Resize(0)) is now gated on the existing CompressedColumnStorage.AutoTrimStorage flag: disabling it keeps the L/U buffers allocated across re-factorizations (iterative solvers, time-stepping) — fully allocation-free reuse. AutoTrimStorage's doc is updated accordingly. Cholesky/LDL allocate the exact factor size, so no trim is involved there. Adds TestRefactorizeNoTrim (refactorize with AutoTrimStorage = false stays correct). All 293 tests pass. --- .../Double/Factorization/SparseLUTest.cs | 35 +++++++++++++++++++ .../Complex/Factorization/SparseCholesky.cs | 10 +++++- CSparse/Complex/Factorization/SparseLDL.cs | 12 +++++-- CSparse/Complex/Factorization/SparseLU.cs | 25 +++++++++---- .../Double/Factorization/SparseCholesky.cs | 10 +++++- CSparse/Double/Factorization/SparseLDL.cs | 12 +++++-- CSparse/Double/Factorization/SparseLU.cs | 25 +++++++++---- CSparse/Storage/CompressedColumnStorage.cs | 5 ++- 8 files changed, 115 insertions(+), 19 deletions(-) diff --git a/CSparse.Tests/Double/Factorization/SparseLUTest.cs b/CSparse.Tests/Double/Factorization/SparseLUTest.cs index 889d3ab..477a1e3 100644 --- a/CSparse.Tests/Double/Factorization/SparseLUTest.cs +++ b/CSparse.Tests/Double/Factorization/SparseLUTest.cs @@ -2,6 +2,7 @@ namespace CSparse.Tests.Double.Factorization { using CSparse.Double; using CSparse.Double.Factorization; + using CSparse.Storage; using NUnit.Framework; using System; @@ -109,5 +110,39 @@ public void TestRefactorize() var small = new SparseMatrix(3, 3, 0); Assert.Throws(() => lu.Refactorize(small, 1.0)); } + + [Test] + public void TestRefactorizeNoTrim() + { + // With AutoTrimStorage disabled, the L/U buffers are kept across + // re-factorizations (no Resize(0) trim) — the result must stay correct. + var trim = CompressedColumnStorage.AutoTrimStorage; + try + { + CompressedColumnStorage.AutoTrimStorage = false; + + var A = ResourceLoader.Get("general-40x40.mat"); + var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0); + + var B = A.Clone(); + var bv = B.Values; + for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0; + + lu.Refactorize(B, 1.0); + + var x = Helper.CreateTestVector(B.ColumnCount); + var b = Helper.Multiply(B, x); + var r = Vector.Clone(b); + + lu.Solve(b, x); + B.Multiply(-1.0, x, 1.0, r); + + Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True); + } + finally + { + CompressedColumnStorage.AutoTrimStorage = trim; + } + } } } diff --git a/CSparse/Complex/Factorization/SparseCholesky.cs b/CSparse/Complex/Factorization/SparseCholesky.cs index 5629060..23b9efd 100644 --- a/CSparse/Complex/Factorization/SparseCholesky.cs +++ b/CSparse/Complex/Factorization/SparseCholesky.cs @@ -286,7 +286,15 @@ private void Factorize(CompressedColumnStorage A, IProgress pro var ci = C.RowIndices; var cx = C.Values; - this.L = CompressedColumnStorage.Create(n, n, colp[n]); + if (this.L is null) + { + this.L = CompressedColumnStorage.Create(n, n, colp[n]); + } + else + { + // Re-factorization (same pattern) : reuse the existing buffer. + this.L.Clear(); + } var lp = L.ColumnPointers; var li = L.RowIndices; diff --git a/CSparse/Complex/Factorization/SparseLDL.cs b/CSparse/Complex/Factorization/SparseLDL.cs index 8f6a1e4..3d3778f 100644 --- a/CSparse/Complex/Factorization/SparseLDL.cs +++ b/CSparse/Complex/Factorization/SparseLDL.cs @@ -276,8 +276,16 @@ void Factorize(CompressedColumnStorage A, IProgress progress) int[] P = S.q; int[] Pinv = S.pinv; - this.D = new double[n]; - this.L = CompressedColumnStorage.Create(n, n, S.cp[n]); + if (this.L is null) + { + this.D = new double[n]; + this.L = CompressedColumnStorage.Create(n, n, S.cp[n]); + } + else + { + // Re-factorization (same pattern) : reuse the existing buffers. + this.L.Clear(); + } Array.Copy(S.cp, L.ColumnPointers, n + 1); diff --git a/CSparse/Complex/Factorization/SparseLU.cs b/CSparse/Complex/Factorization/SparseLU.cs index 513e11a..4a7fe00 100644 --- a/CSparse/Complex/Factorization/SparseLU.cs +++ b/CSparse/Complex/Factorization/SparseLU.cs @@ -210,9 +210,18 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress int lnz = S.lnz; int unz = S.unz; - this.L = CompressedColumnStorage.Create(n, n, lnz); - this.U = CompressedColumnStorage.Create(n, n, unz); - this.pinv = new int[n]; + if (this.L is null) + { + this.L = CompressedColumnStorage.Create(n, n, lnz); + this.U = CompressedColumnStorage.Create(n, n, unz); + this.pinv = new int[n]; + } + else + { + // Re-factorization (same pattern) : reuse the existing buffers. + this.L.Clear(); + this.U.Clear(); + } // Workspace var x = this.temp; @@ -325,9 +334,13 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress li[p] = pinv[li[p]]; } - // Remove extra space from L and U - L.Resize(0); - U.Resize(0); + // Remove extra space from L and U (skipped when trimming is disabled, + // e.g. to reuse the buffers across re-factorizations in an iterative solver). + if (CompressedColumnStorage.AutoTrimStorage) + { + L.Resize(0); + U.Resize(0); + } } /// diff --git a/CSparse/Double/Factorization/SparseCholesky.cs b/CSparse/Double/Factorization/SparseCholesky.cs index e9f8563..a135928 100644 --- a/CSparse/Double/Factorization/SparseCholesky.cs +++ b/CSparse/Double/Factorization/SparseCholesky.cs @@ -281,7 +281,15 @@ private void Factorize(CompressedColumnStorage A, IProgress prog var ci = C.RowIndices; var cx = C.Values; - this.L = CompressedColumnStorage.Create(n, n, colp[n]); + if (this.L is null) + { + this.L = CompressedColumnStorage.Create(n, n, colp[n]); + } + else + { + // Re-factorization (same pattern) : reuse the existing buffer. + this.L.Clear(); + } var lp = L.ColumnPointers; var li = L.RowIndices; diff --git a/CSparse/Double/Factorization/SparseLDL.cs b/CSparse/Double/Factorization/SparseLDL.cs index d9df329..88a041e 100644 --- a/CSparse/Double/Factorization/SparseLDL.cs +++ b/CSparse/Double/Factorization/SparseLDL.cs @@ -275,8 +275,16 @@ void Factorize(CompressedColumnStorage A, IProgress progress) int[] P = S.q; int[] Pinv = S.pinv; - this.D = new double[n]; - this.L = CompressedColumnStorage.Create(n, n, S.cp[n]); + if (this.L is null) + { + this.D = new double[n]; + this.L = CompressedColumnStorage.Create(n, n, S.cp[n]); + } + else + { + // Re-factorization (same pattern) : reuse the existing buffers. + this.L.Clear(); + } Array.Copy(S.cp, L.ColumnPointers, n + 1); diff --git a/CSparse/Double/Factorization/SparseLU.cs b/CSparse/Double/Factorization/SparseLU.cs index 2e08801..0854bf7 100644 --- a/CSparse/Double/Factorization/SparseLU.cs +++ b/CSparse/Double/Factorization/SparseLU.cs @@ -209,9 +209,18 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress< int lnz = S.lnz; int unz = S.unz; - this.L = CompressedColumnStorage.Create(n, n, lnz); - this.U = CompressedColumnStorage.Create(n, n, unz); - this.pinv = new int[n]; + if (this.L is null) + { + this.L = CompressedColumnStorage.Create(n, n, lnz); + this.U = CompressedColumnStorage.Create(n, n, unz); + this.pinv = new int[n]; + } + else + { + // Re-factorization (same pattern) : reuse the existing buffers. + this.L.Clear(); + this.U.Clear(); + } // Workspace var x = this.temp; @@ -324,9 +333,13 @@ private void Factorize(CompressedColumnStorage A, double tol, IProgress< li[p] = pinv[li[p]]; } - // Remove extra space from L and U - L.Resize(0); - U.Resize(0); + // Remove extra space from L and U (skipped when trimming is disabled, + // e.g. to reuse the buffers across re-factorizations in an iterative solver). + if (CompressedColumnStorage.AutoTrimStorage) + { + L.Resize(0); + U.Resize(0); + } } /// diff --git a/CSparse/Storage/CompressedColumnStorage.cs b/CSparse/Storage/CompressedColumnStorage.cs index feb8169..d25a970 100644 --- a/CSparse/Storage/CompressedColumnStorage.cs +++ b/CSparse/Storage/CompressedColumnStorage.cs @@ -18,7 +18,10 @@ public abstract class CompressedColumnStorage : Matrix /// automatically resized to non-zeros count. Defaults to true. /// /// - /// Affects only sparse matrix addition and multiplication. + /// Affects sparse matrix addition and multiplication, and the trimming of + /// the L/U factors at the end of SparseLU.Factorize. Set to + /// false to keep the factor buffers allocated across repeated + /// Refactorize calls (iterative solvers, time-stepping). /// public static bool AutoTrimStorage { get; set; } = true;