Skip to content

Fix SparseLU refactorization #59

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
wo80 merged 2 commits into from
Aug 24, 2026
Merged

Fix SparseLU refactorization #59

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6ebf7fd
Clear SparseLU workspace variables before refactorize.
wo80 Aug 23, 2026
219d4c7
SparseLU (re)factorization should not allocate temporary workspace.
wo80 Aug 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 21 additions & 0 deletions CSparse.Tests/Double/Factorization/SparseLUTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,27 @@ public void TestRefactorize()
Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
}

[Test]
public void TestRefactorize2()
{
const int N = 3;

var A = CompressedColumnStorage<double>.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()
{
Expand Down
15 changes: 11 additions & 4 deletions CSparse/Complex/Factorization/SparseLU.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ public class SparseLU : ISparseFactorization<Complex>
CompressedColumnStorage<Complex> 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

Expand Down Expand Up @@ -99,7 +100,9 @@ public static SparseLU Create(CompressedColumnStorage<Complex> 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];
}

/// <summary>
Expand Down Expand Up @@ -137,6 +140,10 @@ public void Refactorize(CompressedColumnStorage<Complex> 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);
Array.Clear(temp2, 0, 2 * n);

// Reuse the cached symbolic ordering (S); recompute L, U and the pivoting.
Factorize(A, tol, null);
}
Expand Down Expand Up @@ -224,8 +231,8 @@ private void Factorize(CompressedColumnStorage<Complex> 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++)
{
Expand Down
17 changes: 12 additions & 5 deletions CSparse/Double/Factorization/SparseLU.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,9 @@ public class SparseLU : ISparseFactorization<double>
CompressedColumnStorage<double> 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

/// <summary>
Expand Down Expand Up @@ -98,7 +99,9 @@ public static SparseLU Create(CompressedColumnStorage<double> 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];
}

/// <summary>
Expand Down Expand Up @@ -136,6 +139,10 @@ public void Refactorize(CompressedColumnStorage<double> 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);
Array.Clear(temp2, 0, 2 * n);

// Reuse the cached symbolic ordering (S); recompute L, U and the pivoting.
Factorize(A, tol, null);
}
Expand Down Expand Up @@ -223,8 +230,8 @@ private void Factorize(CompressedColumnStorage<double> 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++)
{
Expand Down