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

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

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