Skip to content

Commit c1c7acb

Browse files
wo80wo80
authored
Merge pull request #56 from hbadi/feature/factorize-buffer-reuse
Reuse factor buffers across re-factorizations (follow-up to #55)
2 parents 8898899 + ef2d4fd commit c1c7acb

8 files changed

Lines changed: 115 additions & 19 deletions

File tree

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

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ namespace CSparse.Tests.Double.Factorization
22
{
33
using CSparse.Double;
44
using CSparse.Double.Factorization;
5+
using CSparse.Storage;
56
using NUnit.Framework;
67
using System;
78

@@ -109,5 +110,39 @@ public void TestRefactorize()
109110
var small = new SparseMatrix(3, 3, 0);
110111
Assert.Throws<ArgumentException>(() => lu.Refactorize(small, 1.0));
111112
}
113+
114+
[Test]
115+
public void TestRefactorizeNoTrim()
116+
{
117+
// With AutoTrimStorage disabled, the L/U buffers are kept across
118+
// re-factorizations (no Resize(0) trim) — the result must stay correct.
119+
var trim = CompressedColumnStorage<double>.AutoTrimStorage;
120+
try
121+
{
122+
CompressedColumnStorage<double>.AutoTrimStorage = false;
123+
124+
var A = ResourceLoader.Get<double>("general-40x40.mat");
125+
var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);
126+
127+
var B = A.Clone();
128+
var bv = B.Values;
129+
for (int i = 0; i < bv.Length; i++) bv[i] *= 2.0;
130+
131+
lu.Refactorize(B, 1.0);
132+
133+
var x = Helper.CreateTestVector(B.ColumnCount);
134+
var b = Helper.Multiply(B, x);
135+
var r = Vector.Clone(b);
136+
137+
lu.Solve(b, x);
138+
B.Multiply(-1.0, x, 1.0, r);
139+
140+
Assert.That(Vector.Norm(r.Length, r) < EPS, Is.True);
141+
}
142+
finally
143+
{
144+
CompressedColumnStorage<double>.AutoTrimStorage = trim;
145+
}
146+
}
112147
}
113148
}

‎CSparse/Complex/Factorization/SparseCholesky.cs‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -286,7 +286,15 @@ private void Factorize(CompressedColumnStorage<Complex> A, IProgress<double> pro
286286
var ci = C.RowIndices;
287287
var cx = C.Values;
288288

289-
this.L = CompressedColumnStorage<Complex>.Create(n, n, colp[n]);
289+
if (this.L is null)
290+
{
291+
this.L = CompressedColumnStorage<Complex>.Create(n, n, colp[n]);
292+
}
293+
else
294+
{
295+
// Re-factorization (same pattern) : reuse the existing buffer.
296+
this.L.Clear();
297+
}
290298

291299
var lp = L.ColumnPointers;
292300
var li = L.RowIndices;

‎CSparse/Complex/Factorization/SparseLDL.cs‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -276,8 +276,16 @@ void Factorize(CompressedColumnStorage<Complex> A, IProgress<double> progress)
276276
int[] P = S.q;
277277
int[] Pinv = S.pinv;
278278

279-
this.D = new double[n];
280-
this.L = CompressedColumnStorage<Complex>.Create(n, n, S.cp[n]);
279+
if (this.L is null)
280+
{
281+
this.D = new double[n];
282+
this.L = CompressedColumnStorage<Complex>.Create(n, n, S.cp[n]);
283+
}
284+
else
285+
{
286+
// Re-factorization (same pattern) : reuse the existing buffers.
287+
this.L.Clear();
288+
}
281289

282290
Array.Copy(S.cp, L.ColumnPointers, n + 1);
283291

‎CSparse/Complex/Factorization/SparseLU.cs‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -210,9 +210,18 @@ private void Factorize(CompressedColumnStorage<Complex> A, double tol, IProgress
210210
int lnz = S.lnz;
211211
int unz = S.unz;
212212

213-
this.L = CompressedColumnStorage<Complex>.Create(n, n, lnz);
214-
this.U = CompressedColumnStorage<Complex>.Create(n, n, unz);
215-
this.pinv = new int[n];
213+
if (this.L is null)
214+
{
215+
this.L = CompressedColumnStorage<Complex>.Create(n, n, lnz);
216+
this.U = CompressedColumnStorage<Complex>.Create(n, n, unz);
217+
this.pinv = new int[n];
218+
}
219+
else
220+
{
221+
// Re-factorization (same pattern) : reuse the existing buffers.
222+
this.L.Clear();
223+
this.U.Clear();
224+
}
216225

217226
// Workspace
218227
var x = this.temp;
@@ -325,9 +334,13 @@ private void Factorize(CompressedColumnStorage<Complex> A, double tol, IProgress
325334
li[p] = pinv[li[p]];
326335
}
327336

328-
// Remove extra space from L and U
329-
L.Resize(0);
330-
U.Resize(0);
337+
// Remove extra space from L and U (skipped when trimming is disabled,
338+
// e.g. to reuse the buffers across re-factorizations in an iterative solver).
339+
if (CompressedColumnStorage<Complex>.AutoTrimStorage)
340+
{
341+
L.Resize(0);
342+
U.Resize(0);
343+
}
331344
}
332345

333346
/// <summary>

‎CSparse/Double/Factorization/SparseCholesky.cs‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -281,7 +281,15 @@ private void Factorize(CompressedColumnStorage<double> A, IProgress<double> prog
281281
var ci = C.RowIndices;
282282
var cx = C.Values;
283283

284-
this.L = CompressedColumnStorage<double>.Create(n, n, colp[n]);
284+
if (this.L is null)
285+
{
286+
this.L = CompressedColumnStorage<double>.Create(n, n, colp[n]);
287+
}
288+
else
289+
{
290+
// Re-factorization (same pattern) : reuse the existing buffer.
291+
this.L.Clear();
292+
}
285293

286294
var lp = L.ColumnPointers;
287295
var li = L.RowIndices;

‎CSparse/Double/Factorization/SparseLDL.cs‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -275,8 +275,16 @@ void Factorize(CompressedColumnStorage<double> A, IProgress<double> progress)
275275
int[] P = S.q;
276276
int[] Pinv = S.pinv;
277277

278-
this.D = new double[n];
279-
this.L = CompressedColumnStorage<double>.Create(n, n, S.cp[n]);
278+
if (this.L is null)
279+
{
280+
this.D = new double[n];
281+
this.L = CompressedColumnStorage<double>.Create(n, n, S.cp[n]);
282+
}
283+
else
284+
{
285+
// Re-factorization (same pattern) : reuse the existing buffers.
286+
this.L.Clear();
287+
}
280288

281289
Array.Copy(S.cp, L.ColumnPointers, n + 1);
282290

‎CSparse/Double/Factorization/SparseLU.cs‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -209,9 +209,18 @@ private void Factorize(CompressedColumnStorage<double> A, double tol, IProgress<
209209
int lnz = S.lnz;
210210
int unz = S.unz;
211211

212-
this.L = CompressedColumnStorage<double>.Create(n, n, lnz);
213-
this.U = CompressedColumnStorage<double>.Create(n, n, unz);
214-
this.pinv = new int[n];
212+
if (this.L is null)
213+
{
214+
this.L = CompressedColumnStorage<double>.Create(n, n, lnz);
215+
this.U = CompressedColumnStorage<double>.Create(n, n, unz);
216+
this.pinv = new int[n];
217+
}
218+
else
219+
{
220+
// Re-factorization (same pattern) : reuse the existing buffers.
221+
this.L.Clear();
222+
this.U.Clear();
223+
}
215224

216225
// Workspace
217226
var x = this.temp;
@@ -324,9 +333,13 @@ private void Factorize(CompressedColumnStorage<double> A, double tol, IProgress<
324333
li[p] = pinv[li[p]];
325334
}
326335

327-
// Remove extra space from L and U
328-
L.Resize(0);
329-
U.Resize(0);
336+
// Remove extra space from L and U (skipped when trimming is disabled,
337+
// e.g. to reuse the buffers across re-factorizations in an iterative solver).
338+
if (CompressedColumnStorage<double>.AutoTrimStorage)
339+
{
340+
L.Resize(0);
341+
U.Resize(0);
342+
}
330343
}
331344

332345
/// <summary>

‎CSparse/Storage/CompressedColumnStorage.cs‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,10 @@ public abstract class CompressedColumnStorage<T> : Matrix<T>
1818
/// automatically resized to non-zeros count. Defaults to true.
1919
/// </summary>
2020
/// <remarks>
21-
/// Affects only sparse matrix addition and multiplication.
21+
/// Affects sparse matrix addition and multiplication, and the trimming of
22+
/// the L/U factors at the end of <c>SparseLU.Factorize</c>. Set to
23+
/// <c>false</c> to keep the factor buffers allocated across repeated
24+
/// <c>Refactorize</c> calls (iterative solvers, time-stepping).
2225
/// </remarks>
2326
public static bool AutoTrimStorage { get; set; } = true;
2427

0 commit comments

Comments
 (0)