Factorize() should clear the SparseLU workspace now that it can be refactorized.

[Mervill]

Issue

int n = 3;

// A: an MNA matrix with an ideal-voltage-source row/column
// (column 2's only structural entry is at row 0 - its own diagonal is absent):
//   [ 0.001   -0.001   1 ]
//   [-0.001   10.001   0 ]
//   [ 1        0       0 ]   <- row 2 encodes the constraint x[0] == 5
var coo = new CoordinateStorage<double>(n, n, 6);
coo.At(0, 0, 0.001);
coo.At(1, 0, -0.001);
coo.At(2, 0, 1.0);
coo.At(0, 1, -0.001);
coo.At(1, 1, 10.001);
coo.At(0, 2, 1.0);

var A = CompressedColumnStorage<double>.OfIndexed(coo, inplace: true);

double[] someNonzeroB = { 1.0, 2.0, 3.0 };
double[] B = { 0.0, 0.0049995, 5.0 };
double[] throwaway = new double[n];
double[] x = new double[n];
double[] xFresh = new double[n];

var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1e-8);
lu.Solve(someNonzeroB, throwaway); // leaves `temp` dirty with a real solution
lu.Refactorize(A, 1e-8);           // reuses the dirty `temp`
lu.Solve(B, x);

var luFresh = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1e-8);
luFresh.Solve(B, xFresh);

Console.WriteLine($"x (via Refactorize) = [{string.Join(", ", x)}]");
Console.WriteLine($"xFresh (via Create) = [{string.Join(", ", xFresh)}]");
Console.WriteLine($"expected x[0] == 5.0 (row 2's constraint); via Refactorize got {x[0]}");

// Output:
//   x (via Refactorize) = [5.015042122304683, 0.0010013540768227859, -0.005014040768227683]
//   xFresh (via Create) = [5, 0.0009998500149985, -0.004999000149985001]
//   expected x[0] == 5.0 (row 2's constraint); via Refactorize got 5.015042122304683

Solution

SparseLU.Factorize(...) should clear the workspace:

// Workspace
var x = this.temp;
Array.Clear(x, 0, n);
var xi = new int[2 * n];

For both Double and Complex versions.

More

In a Factorize call, x is only self-cleaning for indices in the current column's elimination "reach" set and SolveSp clears those before use. The diagonal-preference pivot check, if (pinv[col] < 0 && Math.Abs(x[col]) >= a * tol), reads x[col] unconditionally, without checking whether col is actually in that reach set. When a column's own diagonal is structurally absent from A and not reachable through the partialy built L (e.g. an MNA circuit's ideal voltage-source branch column), x[col] is never written this iteration, so the read returns whatever was last left in temp[col]. At low pivot tolerances that stale value easily passes the >= a * tol check.

[Mervill]

Oh also var xi = new int[2 * n]; should be hoisted to a member alongside temp so then Factorize() wont allocate anything,

[wo80]

Thorough analysis and valid fix, thanks!

Please take a look at #59 .