Reuse factor buffers across re-factorizations (follow-up to #55)

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<ArgumentException>(() => 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<double>.AutoTrimStorage;
+            try
+            {
+                CompressedColumnStorage<double>.AutoTrimStorage = false;
+
+                var A = ResourceLoader.Get<double>("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<double>.AutoTrimStorage = trim;
+            }
+        }
     }
 }

CSparse/Complex/Factorization/SparseCholesky.cs

@@ -286,7 +286,15 @@ private void Factorize(CompressedColumnStorage<Complex> A, IProgress<double> pro
             var ci = C.RowIndices;
             var cx = C.Values;
 
-            this.L = CompressedColumnStorage<Complex>.Create(n, n, colp[n]);
+            if (this.L is null)
+            {
+                this.L = CompressedColumnStorage<Complex>.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;

CSparse/Complex/Factorization/SparseLDL.cs

@@ -276,8 +276,16 @@ void Factorize(CompressedColumnStorage<Complex> A, IProgress<double> progress)
             int[] P = S.q;
             int[] Pinv = S.pinv;
 
-            this.D = new double[n];
-            this.L = CompressedColumnStorage<Complex>.Create(n, n, S.cp[n]);
+            if (this.L is null)
+            {
+                this.D = new double[n];
+                this.L = CompressedColumnStorage<Complex>.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);
 

CSparse/Complex/Factorization/SparseLU.cs

@@ -210,9 +210,18 @@ private void Factorize(CompressedColumnStorage<Complex> A, double tol, IProgress
             int lnz = S.lnz;
             int unz = S.unz;
 
-            this.L = CompressedColumnStorage<Complex>.Create(n, n, lnz);
-            this.U = CompressedColumnStorage<Complex>.Create(n, n, unz);
-            this.pinv = new int[n];
+            if (this.L is null)
+            {
+                this.L = CompressedColumnStorage<Complex>.Create(n, n, lnz);
+                this.U = CompressedColumnStorage<Complex>.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<Complex> 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<Complex>.AutoTrimStorage)
+            {
+                L.Resize(0);
+                U.Resize(0);
+            }
         }
 
         /// <summary>

CSparse/Double/Factorization/SparseCholesky.cs

@@ -281,7 +281,15 @@ private void Factorize(CompressedColumnStorage<double> A, IProgress<double> prog
             var ci = C.RowIndices;
             var cx = C.Values;
 
-            this.L = CompressedColumnStorage<double>.Create(n, n, colp[n]);
+            if (this.L is null)
+            {
+                this.L = CompressedColumnStorage<double>.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;

CSparse/Double/Factorization/SparseLDL.cs

@@ -275,8 +275,16 @@ void Factorize(CompressedColumnStorage<double> A, IProgress<double> progress)
             int[] P = S.q;
             int[] Pinv = S.pinv;
 
-            this.D = new double[n];
-            this.L = CompressedColumnStorage<double>.Create(n, n, S.cp[n]);
+            if (this.L is null)
+            {
+                this.D = new double[n];
+                this.L = CompressedColumnStorage<double>.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);
 

CSparse/Double/Factorization/SparseLU.cs

@@ -209,9 +209,18 @@ private void Factorize(CompressedColumnStorage<double> A, double tol, IProgress<
             int lnz = S.lnz;
             int unz = S.unz;
 
-            this.L = CompressedColumnStorage<double>.Create(n, n, lnz);
-            this.U = CompressedColumnStorage<double>.Create(n, n, unz);
-            this.pinv = new int[n];
+            if (this.L is null)
+            {
+                this.L = CompressedColumnStorage<double>.Create(n, n, lnz);
+                this.U = CompressedColumnStorage<double>.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<double> 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<double>.AutoTrimStorage)
+            {
+                L.Resize(0);
+                U.Resize(0);
+            }
         }
 
         /// <summary>

CSparse/Storage/CompressedColumnStorage.cs

@@ -18,7 +18,10 @@ public abstract class CompressedColumnStorage<T> : Matrix<T>
         /// automatically resized to non-zeros count. Defaults to true.
         /// </summary>
         /// <remarks>
-        /// 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 <c>SparseLU.Factorize</c>. Set to
+        /// <c>false</c> to keep the factor buffers allocated across repeated
+        /// <c>Refactorize</c> calls (iterative solvers, time-stepping).
         /// </remarks>
         public static bool AutoTrimStorage { get; set; } = true;