Skip to content

Commit e08a79a

Browse files
hbadihbadi
committed
Add Refactorize() to reuse the symbolic analysis for numeric refactorization
SparseLU / SparseCholesky / SparseLDL (Double and Complex) gain a public Refactorize(matrix[, tol]) that re-runs only the numeric factorization, reusing the cached symbolic analysis (fill-reducing ordering, elimination tree). Intended for Newton iterations and time-stepping loops where the matrix values change but the sparsity pattern does not: it skips the ordering/symbolic phase and keeps only the (dominant) numeric work. - LU: a different pattern still factorizes correctly (column ordering is a permutation), only the fill may be sub-optimal. - Cholesky / LDL: the elimination tree is reused, so the pattern must match. - Dimension guard throws ArgumentException on a size mismatch. Adds TestRefactorize for all six factorizations (refactorize B = 2*A on the same pattern reproduces a correct solve; dimension guard throws). All 40 factorization tests pass.
1 parent c736f87 commit e08a79a

12 files changed

Lines changed: 344 additions & 0 deletions

File tree

‎CSparse.Tests/Complex/Factorization/SparseCholeskyTest.cs‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,5 +54,34 @@ public void TestEmptyFactorize()
5454
Assert.That(chol, Is.Not.Null);
5555
Assert.That(chol.NonZerosCount == 0, Is.True);
5656
}
57+
58+
[Test]
59+
public void TestRefactorize()
60+
{
61+
var A = ResourceLoader.Get<Complex>("hermitian-40-spd.mat");
62+
63+
var chol = SparseCholesky.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
64+
65+
// Same sparsity pattern, different values (B = 2*A stays Hermitian SPD) :
66+
// the numeric refactorization must reuse the cached symbolic analysis.
67+
var B = A.Clone();
68+
var bv = B.Values;
69+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
70+
71+
chol.Refactorize(B);
72+
73+
var x = Helper.CreateTestVector(B.ColumnCount);
74+
var b = Helper.Multiply(B, x);
75+
var r = Vector.Clone(b);
76+
77+
chol.Solve(b, x);
78+
B.Multiply(-1.0, x, 1.0, r);
79+
80+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
81+
82+
// Dimension guard.
83+
var small = new SparseMatrix(3, 3, 0);
84+
Assert.Throws<ArgumentException>(() => chol.Refactorize(small));
85+
}
5786
}
5887
}

‎CSparse.Tests/Complex/Factorization/SparseLDLTest.cs‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
namespace CSparse.Tests.Complex.Factorization
22
{
33

4+
using System;
45
using CSparse.Complex;
56
using CSparse.Complex.Factorization;
67
using NUnit.Framework;
@@ -54,5 +55,34 @@ public void TestSolveNonSpd()
5455

5556
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
5657
}
58+
59+
[Test]
60+
public void TestRefactorize()
61+
{
62+
var A = ResourceLoader.Get<Complex>("hermitian-40-spd.mat");
63+
64+
var ldl = SparseLDL.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
65+
66+
// Same sparsity pattern, different values (B = 2*A) : the numeric
67+
// refactorization must reuse the cached symbolic analysis.
68+
var B = A.Clone();
69+
var bv = B.Values;
70+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
71+
72+
ldl.Refactorize(B);
73+
74+
var x = Helper.CreateTestVector(B.ColumnCount);
75+
var b = Helper.Multiply(B, x);
76+
var r = Vector.Clone(b);
77+
78+
ldl.Solve(b, x);
79+
B.Multiply(-1.0, x, 1.0, r);
80+
81+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
82+
83+
// Dimension guard.
84+
var small = new SparseMatrix(3, 3, 0);
85+
Assert.Throws<ArgumentException>(() => ldl.Refactorize(small));
86+
}
5787
}
5888
}

‎CSparse.Tests/Complex/Factorization/SparseLUTest.cs‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,5 +76,34 @@ public void TestEmptyFactorize()
7676
Assert.That(lu, Is.Not.Null);
7777
Assert.That(lu.NonZerosCount == 0, Is.True);
7878
}
79+
80+
[Test]
81+
public void TestRefactorize()
82+
{
83+
var A = ResourceLoader.Get<Complex>("general-40x40.mat");
84+
85+
var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);
86+
87+
// Same sparsity pattern, different values (B = 2*A) : numeric refactorization
88+
// must reuse the cached symbolic ordering and still solve correctly.
89+
var B = A.Clone();
90+
var bv = B.Values;
91+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
92+
93+
lu.Refactorize(B, 1.0);
94+
95+
var x = Helper.CreateTestVector(B.ColumnCount);
96+
var b = Helper.Multiply(B, x);
97+
var r = Vector.Clone(b);
98+
99+
lu.Solve(b, x);
100+
B.Multiply(-1.0, x, 1.0, r);
101+
102+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
103+
104+
// Dimension guard.
105+
var small = new SparseMatrix(3, 3, 0);
106+
Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
107+
}
79108
}
80109
}

‎CSparse.Tests/Double/Factorization/SparseCholeskyTest.cs‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,5 +53,34 @@ public void TestEmptyFactorize()
5353
Assert.That(chol, Is.Not.Null);
5454
Assert.That(chol.NonZerosCount == 0, Is.True);
5555
}
56+
57+
[Test]
58+
public void TestRefactorize()
59+
{
60+
var A = ResourceLoader.Get<double>("symmetric-40-spd.mat");
61+
62+
var chol = SparseCholesky.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
63+
64+
// Same sparsity pattern, different values (B = 2*A stays SPD) : the numeric
65+
// refactorization must reuse the cached symbolic analysis and still solve.
66+
var B = A.Clone();
67+
var bv = B.Values;
68+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
69+
70+
chol.Refactorize(B);
71+
72+
var x = Helper.CreateTestVector(B.ColumnCount);
73+
var b = Helper.Multiply(B, x);
74+
var r = Vector.Clone(b);
75+
76+
chol.Solve(b, x);
77+
B.Multiply(-1.0, x, 1.0, r);
78+
79+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
80+
81+
// Dimension guard.
82+
var small = new SparseMatrix(3, 3, 0);
83+
Assert.Throws<ArgumentException>(() => chol.Refactorize(small));
84+
}
5685
}
5786
}

‎CSparse.Tests/Double/Factorization/SparseLDLTest.cs‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,5 +53,34 @@ public void TestSolveNonSpd()
5353

5454
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
5555
}
56+
57+
[Test]
58+
public void TestRefactorize()
59+
{
60+
var A = ResourceLoader.Get<double>("symmetric-40-spd.mat");
61+
62+
var ldl = SparseLDL.Create(A, ColumnOrdering.MinimumDegreeAtPlusA);
63+
64+
// Same sparsity pattern, different values (B = 2*A) : the numeric
65+
// refactorization must reuse the cached symbolic analysis and still solve.
66+
var B = A.Clone();
67+
var bv = B.Values;
68+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
69+
70+
ldl.Refactorize(B);
71+
72+
var x = Helper.CreateTestVector(B.ColumnCount);
73+
var b = Helper.Multiply(B, x);
74+
var r = Vector.Clone(b);
75+
76+
ldl.Solve(b, x);
77+
B.Multiply(-1.0, x, 1.0, r);
78+
79+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
80+
81+
// Dimension guard.
82+
var small = new SparseMatrix(3, 3, 0);
83+
Assert.Throws<ArgumentException>(() => ldl.Refactorize(small));
84+
}
5685
}
5786
}

‎CSparse.Tests/Double/Factorization/SparseLUTest.cs‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,5 +75,39 @@ public void TestEmptyFactorize()
7575
Assert.That(lu, Is.Not.Null);
7676
Assert.That(lu.NonZerosCount == 0, Is.True);
7777
}
78+
79+
[Test]
80+
public void TestRefactorize()
81+
{
82+
// Load matrix from a file.
83+
var A = ResourceLoader.Get<double>("general-40x40.mat");
84+
85+
// Symbolic + numeric factorization of A.
86+
var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);
87+
88+
// Same sparsity pattern, different values (B = 2*A) : numeric refactorization
89+
// must reuse the cached symbolic ordering and still solve correctly.
90+
var B = A.Clone();
91+
var bv = B.Values;
92+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
93+
94+
lu.Refactorize(B, 1.0);
95+
96+
var x = Helper.CreateTestVector(B.ColumnCount);
97+
var b = Helper.Multiply(B, x);
98+
var r = Vector.Clone(b);
99+
100+
// Solve Bx = b with the refactorized decomposition.
101+
lu.Solve(b, x);
102+
103+
// Residual r = b - Bx.
104+
B.Multiply(-1.0, x, 1.0, r);
105+
106+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
107+
108+
// Dimension guard.
109+
var small = new SparseMatrix(3, 3, 0);
110+
Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
111+
}
78112
}
79113
}

‎CSparse/Complex/Factorization/SparseCholesky.cs‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,32 @@ public int NonZerosCount
106106
get { return L.NonZerosCount; }
107107
}
108108

109+
/// <summary>
110+
/// Numerically re-factorizes a Hermitian positive definite matrix that has the
111+
/// <b>same sparsity pattern</b> as the one passed to <c>Create</c>, reusing the
112+
/// cached symbolic analysis (ordering, elimination tree and column counts).
113+
/// </summary>
114+
/// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
115+
/// <remarks>
116+
/// Intended for Newton iterations and time-stepping loops where the matrix
117+
/// values change but its pattern does not: the fill-reducing ordering and the
118+
/// symbolic analysis are skipped, keeping only the numeric work. The matrix
119+
/// must keep the pattern used in <c>Create</c> (the symbolic structure is
120+
/// reused) and stay positive definite.
121+
/// </remarks>
122+
public void Refactorize(CompressedColumnStorage<Complex> A)
123+
{
124+
Check.NotNull(A, "A");
125+
Check.SquareMatrix(A, "A");
126+
127+
if (A.ColumnCount != n)
128+
{
129+
throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
130+
}
131+
132+
Factorize(A, null);
133+
}
134+
109135
/// <summary>
110136
/// Solves a system of linear equations, <c>Ax = b</c>.
111137
/// </summary>

‎CSparse/Complex/Factorization/SparseLDL.cs‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,31 @@ public int NonZerosCount
112112
get { return L.NonZerosCount; }
113113
}
114114

115+
/// <summary>
116+
/// Numerically re-factorizes a matrix that has the <b>same sparsity pattern</b>
117+
/// as the one passed to <c>Create</c>, reusing the cached symbolic analysis
118+
/// (ordering and elimination tree).
119+
/// </summary>
120+
/// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
121+
/// <remarks>
122+
/// Intended for Newton iterations and time-stepping loops where the matrix
123+
/// values change but its pattern does not: the fill-reducing ordering and the
124+
/// symbolic analysis are skipped, keeping only the numeric work. The matrix
125+
/// must keep the pattern used in <c>Create</c> (the symbolic structure is reused).
126+
/// </remarks>
127+
public void Refactorize(CompressedColumnStorage<Complex> A)
128+
{
129+
Check.NotNull(A, "A");
130+
Check.SquareMatrix(A, "A");
131+
132+
if (A.ColumnCount != n)
133+
{
134+
throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
135+
}
136+
137+
Factorize(A, null);
138+
}
139+
115140
/// <summary>
116141
/// Solves a linear system Ax=b, where A is symmetric positive definite.
117142
/// </summary>

‎CSparse/Complex/Factorization/SparseLU.cs‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,37 @@ public int NonZerosCount
110110
get { return (L.NonZerosCount + U.NonZerosCount - n); }
111111
}
112112

113+
/// <summary>
114+
/// Numerically re-factorizes a matrix that has the <b>same sparsity pattern</b>
115+
/// as the one passed to <c>Create</c>, reusing the cached symbolic ordering.
116+
/// </summary>
117+
/// <param name="A">Column-compressed matrix, same size as the one used in <c>Create</c>.</param>
118+
/// <param name="tol">Partial pivoting tolerance (from 0.0 to 1.0).</param>
119+
/// <remarks>
120+
/// Intended for Newton iterations and time-stepping loops where the matrix
121+
/// values change but its pattern does not: the fill-reducing column ordering
122+
/// and symbolic analysis are skipped, keeping only the (dominant) numeric work.
123+
/// The column ordering is a permutation, so a different pattern still factorizes
124+
/// correctly, but the fill may then be sub-optimal.
125+
/// </remarks>
126+
public void Refactorize(CompressedColumnStorage<Complex> A, double tol = 1.0)
127+
{
128+
Check.NotNull(A, "A");
129+
Check.SquareMatrix(A, "A");
130+
Check.NotNaN(tol, "tol");
131+
132+
if (A.ColumnCount != n)
133+
{
134+
throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
135+
}
136+
137+
// Ensure tol is in range.
138+
tol = Math.Min(Math.Max(tol, 0.0), 1.0);
139+
140+
// Reuse the cached symbolic ordering (S); recompute L, U and the pivoting.
141+
Factorize(A, tol, null);
142+
}
143+
113144
/// <summary>
114145
/// Solves a system of linear equations, <c>Ax = b</c>.
115146
/// </summary>

‎CSparse/Double/Factorization/SparseCholesky.cs‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,32 @@ public int NonZerosCount
105105
get { return L.NonZerosCount; }
106106
}
107107

108+
/// <summary>
109+
/// Numerically re-factorizes a symmetric positive definite matrix that has the
110+
/// <b>same sparsity pattern</b> as the one passed to <c>Create</c>, reusing the
111+
/// cached symbolic analysis (ordering, elimination tree and column counts).
112+
/// </summary>
113+
/// <param name="A">Column-compressed matrix, same size and pattern as in <c>Create</c>.</param>
114+
/// <remarks>
115+
/// Intended for Newton iterations and time-stepping loops where the matrix
116+
/// values change but its pattern does not: the fill-reducing ordering and the
117+
/// symbolic analysis are skipped, keeping only the numeric work. The matrix
118+
/// must keep the pattern used in <c>Create</c> (the symbolic structure is
119+
/// reused) and stay positive definite.
120+
/// </remarks>
121+
public void Refactorize(CompressedColumnStorage<double> A)
122+
{
123+
Check.NotNull(A, "A");
124+
Check.SquareMatrix(A, "A");
125+
126+
if (A.ColumnCount != n)
127+
{
128+
throw new ArgumentException("Matrix dimensions don't match the factorization.", "A");
129+
}
130+
131+
Factorize(A, null);
132+
}
133+
108134
/// <summary>
109135
/// Solves a system of linear equations, <c>Ax = b</c>.
110136
/// </summary>

0 commit comments

Comments
 (0)