Add Refactorize() for numeric refactorization reusing the symbolic analysis

CSparse.Tests/Complex/Factorization/SparseCholeskyTest.cs

@@ -54,5 +54,34 @@ public void TestEmptyFactorize()
             Assert.That(chol, Is.Not.Null);
             Assert.That(chol.NonZerosCount == 0, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            var A = ResourceLoader.Get<Complex>("hermitian-40-spd.mat");
+
+            var chol = SparseCholesky.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
+
+            // Same sparsity pattern, different values (B = 2*A stays Hermitian SPD) :
+            // the numeric refactorization must reuse the cached symbolic analysis.
+            var B = A.Clone();
+            var bv = B.Values;
+            for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
+
+            chol.Refactorize(B);
+
+            var x = Helper.CreateTestVector(B.ColumnCount);
+            var b = Helper.Multiply(B, x);
+            var r = Vector.Clone(b);
+
+            chol.Solve(b, x);
+            B.Multiply(-1.0, x, 1.0, r);
+
+            Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => chol.Refactorize(small));
+        }
     }
 }

CSparse.Tests/Complex/Factorization/SparseLDLTest.cs

@@ -1,6 +1,7 @@
 namespace CSparse.Tests.Complex.Factorization
 {
 
+    using System;
     using CSparse.Complex;
     using CSparse.Complex.Factorization;
     using NUnit.Framework;
@@ -54,5 +55,34 @@ public void TestSolveNonSpd()
 
             Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            var A = ResourceLoader.Get<Complex>("hermitian-40-spd.mat");
+
+            var ldl = SparseLDL.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
+
+            // Same sparsity pattern, different values (B = 2*A) : the numeric
+            // refactorization must reuse the cached symbolic analysis.
+            var B = A.Clone();
+            var bv = B.Values;
+            for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
+
+            ldl.Refactorize(B);
+
+            var x = Helper.CreateTestVector(B.ColumnCount);
+            var b = Helper.Multiply(B, x);
+            var r = Vector.Clone(b);
+
+            ldl.Solve(b, x);
+            B.Multiply(-1.0, x, 1.0, r);
+
+            Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => ldl.Refactorize(small));
+        }
     }
 }

CSparse.Tests/Complex/Factorization/SparseLUTest.cs

@@ -76,5 +76,34 @@ public void TestEmptyFactorize()
             Assert.That(lu, Is.Not.Null);
             Assert.That(lu.NonZerosCount == 0, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            var A = ResourceLoader.Get<Complex>("general-40x40.mat");
+
+            var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);
+
+            // Same sparsity pattern, different values (B = 2*A) : numeric refactorization
+            // must reuse the cached symbolic ordering and still solve correctly.
+            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);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
+        }
     }
 }

CSparse.Tests/Double/Factorization/SparseCholeskyTest.cs

@@ -53,5 +53,34 @@ public void TestEmptyFactorize()
             Assert.That(chol, Is.Not.Null);
             Assert.That(chol.NonZerosCount == 0, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            var A = ResourceLoader.Get<double>("symmetric-40-spd.mat");
+
+            var chol = SparseCholesky.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
+
+            // Same sparsity pattern, different values (B = 2*A stays SPD) : the numeric
+            // refactorization must reuse the cached symbolic analysis and still solve.
+            var B = A.Clone();
+            var bv = B.Values;
+            for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
+
+            chol.Refactorize(B);
+
+            var x = Helper.CreateTestVector(B.ColumnCount);
+            var b = Helper.Multiply(B, x);
+            var r = Vector.Clone(b);
+
+            chol.Solve(b, x);
+            B.Multiply(-1.0, x, 1.0, r);
+
+            Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => chol.Refactorize(small));
+        }
     }
 }

CSparse.Tests/Double/Factorization/SparseLDLTest.cs

@@ -53,5 +53,34 @@ public void TestSolveNonSpd()
 
             Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            var A = ResourceLoader.Get<double>("symmetric-40-spd.mat");
+
+            var ldl = SparseLDL.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
+
+            // Same sparsity pattern, different values (B = 2*A) : the numeric
+            // refactorization must reuse the cached symbolic analysis and still solve.
+            var B = A.Clone();
+            var bv = B.Values;
+            for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
+
+            ldl.Refactorize(B);
+
+            var x = Helper.CreateTestVector(B.ColumnCount);
+            var b = Helper.Multiply(B, x);
+            var r = Vector.Clone(b);
+
+            ldl.Solve(b, x);
+            B.Multiply(-1.0, x, 1.0, r);
+
+            Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => ldl.Refactorize(small));
+        }
     }
 }

CSparse.Tests/Double/Factorization/SparseLUTest.cs

@@ -75,5 +75,39 @@ public void TestEmptyFactorize()
             Assert.That(lu, Is.Not.Null);
             Assert.That(lu.NonZerosCount == 0, Is.True);
         }
+
+        [Test]
+        public void TestRefactorize()
+        {
+            // Load matrix from a file.
+            var A = ResourceLoader.Get<double>("general-40x40.mat");
+
+            // Symbolic + numeric factorization of A.
+            var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);
+
+            // Same sparsity pattern, different values (B = 2*A) : numeric refactorization
+            // must reuse the cached symbolic ordering and still solve correctly.
+            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);
+
+            // Solve Bx = b with the refactorized decomposition.
+            lu.Solve(b, x);
+
+            // Residual r = b - Bx.
+            B.Multiply(-1.0, x, 1.0, r);
+
+            Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
+
+            // Dimension guard.
+            var small = new SparseMatrix(3, 3, 0);
+            Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
+        }
     }
 }

CSparse/Complex/Factorization/SparseCholesky.cs

@@ -106,6 +106,32 @@ public int NonZerosCount
             get { return L.NonZerosCount; }
         }
 
+        /// <summary>
+        /// Numerically re-factorizes a Hermitian positive definite matrix that has the
+        /// <b>same sparsity pattern</b> as the one passed to <c>Create</c>, reusing the
+        /// cached symbolic analysis (ordering, elimination tree and column counts).
+        /// </summary>
+        /// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
+        /// <remarks>
+        /// Intended for Newton iterations and time-stepping loops where the matrix
+        /// values change but its pattern does not: the fill-reducing ordering and the
+        /// symbolic analysis are skipped, keeping only the numeric work. The matrix
+        /// must keep the pattern used in <c>Create</c> (the symbolic structure is
+        /// reused) and stay positive definite.
+        /// </remarks>
+        public void Refactorize(CompressedColumnStorage<Complex> A)
+        {
+            Check.NotNull(A, "A");
+            Check.SquareMatrix(A, "A");
+
+            if (A.ColumnCount != n)
+            {
+                throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
+            }
+
+            Factorize(A, null);
+        }
+
         /// <summary>
         /// Solves a system of linear equations, <c>Ax = b</c>.
         /// </summary>

CSparse/Complex/Factorization/SparseLDL.cs

@@ -112,6 +112,31 @@ public int NonZerosCount
             get { return L.NonZerosCount; }
         }
 
+        /// <summary>
+        /// Numerically re-factorizes a matrix that has the <b>same sparsity pattern</b>
+        /// as the one passed to <c>Create</c>, reusing the cached symbolic analysis
+        /// (ordering and elimination tree).
+        /// </summary>
+        /// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
+        /// <remarks>
+        /// Intended for Newton iterations and time-stepping loops where the matrix
+        /// values change but its pattern does not: the fill-reducing ordering and the
+        /// symbolic analysis are skipped, keeping only the numeric work. The matrix
+        /// must keep the pattern used in <c>Create</c> (the symbolic structure is reused).
+        /// </remarks>
+        public void Refactorize(CompressedColumnStorage<Complex> A)
+        {
+            Check.NotNull(A, "A");
+            Check.SquareMatrix(A, "A");
+
+            if (A.ColumnCount != n)
+            {
+                throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
+            }
+
+            Factorize(A, null);
+        }
+
         /// <summary>
         /// Solves a linear system Ax=b, where A is symmetric positive definite.
         /// </summary>

CSparse/Complex/Factorization/SparseLU.cs

@@ -110,6 +110,37 @@ public int NonZerosCount
             get { return (L.NonZerosCount + U.NonZerosCount - n); }
         }
 
+        /// <summary>
+        /// Numerically re-factorizes a matrix that has the <b>same sparsity pattern</b>
+        /// as the one passed to <c>Create</c>, reusing the cached symbolic ordering.
+        /// </summary>
+        /// <param name="A">Column-compressed matrix, same size as the one used in <c>Create</c>.</param>
+        /// <param name="tol">Partial pivoting tolerance (from 0.0 to 1.0).</param>
+        /// <remarks>
+        /// Intended for Newton iterations and time-stepping loops where the matrix
+        /// values change but its pattern does not: the fill-reducing column ordering
+        /// and symbolic analysis are skipped, keeping only the (dominant) numeric work.
+        /// The column ordering is a permutation, so a different pattern still factorizes
+        /// correctly, but the fill may then be sub-optimal.
+        /// </remarks>
+        public void Refactorize(CompressedColumnStorage<Complex> A, double tol = 1.0)
+        {
+            Check.NotNull(A, "A");
+            Check.SquareMatrix(A, "A");
+            Check.NotNaN(tol, "tol");
+
+            if (A.ColumnCount != n)
+            {
+                throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
+            }
+
+            // Ensure tol is in range.
+            tol = Math.Min(Math.Max(tol, 0.0), 1.0);
+
+            // Reuse the cached symbolic ordering (S); recompute L, U and the pivoting.
+            Factorize(A, tol, null);
+        }
+
         /// <summary>
         /// Solves a system of linear equations, <c>Ax = b</c>.
         /// </summary>

CSparse/Double/Factorization/SparseCholesky.cs

@@ -105,6 +105,32 @@ public int NonZerosCount
             get { return L.NonZerosCount; }
         }
 
+        /// <summary>
+        /// Numerically re-factorizes a symmetric positive definite matrix that has the
+        /// <b>same sparsity pattern</b> as the one passed to <c>Create</c>, reusing the
+        /// cached symbolic analysis (ordering, elimination tree and column counts).
+        /// </summary>
+        /// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
+        /// <remarks>
+        /// Intended for Newton iterations and time-stepping loops where the matrix
+        /// values change but its pattern does not: the fill-reducing ordering and the
+        /// symbolic analysis are skipped, keeping only the numeric work. The matrix
+        /// must keep the pattern used in <c>Create</c> (the symbolic structure is
+        /// reused) and stay positive definite.
+        /// </remarks>
+        public void Refactorize(CompressedColumnStorage<double> A)
+        {
+            Check.NotNull(A, "A");
+            Check.SquareMatrix(A, "A");
+
+            if (A.ColumnCount != n)
+            {
+                throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
+            }
+
+            Factorize(A, null);
+        }
+
         /// <summary>
         /// Solves a system of linear equations, <c>Ax = b</c>.
         /// </summary>

CSparse/Double/Factorization/SparseLDL.cs

@@ -111,6 +111,31 @@ public int NonZerosCount
             get { return L.NonZerosCount; }
         }
 
+        /// <summary>
+        /// Numerically re-factorizes a matrix that has the <b>same sparsity pattern</b>
+        /// as the one passed to <c>Create</c>, reusing the cached symbolic analysis
+        /// (ordering and elimination tree).
+        /// </summary>
+        /// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
+        /// <remarks>
+        /// Intended for Newton iterations and time-stepping loops where the matrix
+        /// values change but its pattern does not: the fill-reducing ordering and the
+        /// symbolic analysis are skipped, keeping only the numeric work. The matrix
+        /// must keep the pattern used in <c>Create</c> (the symbolic structure is reused).
+        /// </remarks>
+        public void Refactorize(CompressedColumnStorage<double> A)
+        {
+            Check.NotNull(A, "A");
+            Check.SquareMatrix(A, "A");
+
+            if (A.ColumnCount != n)
+            {
+                throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
+            }
+
+            Factorize(A, null);
+        }
+
         /// <summary>
         /// Solves a linear system Ax=b, where A is symmetric positive definite.
         /// </summary>