From 394854a055e5ba890dbf1fd5c29d53f2964efda2 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 9 Sep 2026 10:34:01 +0000 Subject: [PATCH] Matrix operations give wrong answers when called from more than one thread Two process-wide statics are read and written without synchronisation, and both fail silently rather than throwing. The scratch-matrix pool behind DeterminantLaplace and Inverse hands every caller the same GenTensor for a given size, and the caller writes into it -- GetCofactorMatrix fills the matrix it is given. Two threads taking a determinant of the same size are handed the same buffer and overwrite each other's minors between the write and the read. Sixty 5x5 integer matrices, computed once sequentially for truth and then again under Parallel.For: DeterminantLaplace 53 of 60 wrong Adjoint 60 of 60 wrong Nothing throws and no result is malformed. They are simply different numbers. Making the pool [ThreadStatic] fixes it; the buffer is scratch, so there is nothing to share between threads in the first place. The lock inside GetMatrix goes with it, since nothing outside the owning thread can now see the list -- and that lock was in any case incomplete, because the list was indexed after it had been released, so a concurrent Add could reallocate underneath the reader. DeterminantGaussianSafeDivision is unaffected: it works on its own copy and never asks the pool for anything. A test records that, so the boundary stays recorded rather than assumed. The compiled-operation cache behind PiecewiseAdd and its siblings is a plain Dictionary written from whichever thread asks first. Reached cold from several threads it fails with the runtime's own diagnostic -- "Operations that change non-concurrent collections must have exclusive access. A concurrent update was performed on this collection and corrupted its state" -- and then with a lookup missing a key that had been written. A ConcurrentDictionary fixes it, and a duplicate compilation under GetOrAdd is harmless because two compilations of one key produce delegates that do the same thing. The first version of that test computed its expected values sequentially and only then went parallel, which populated every key before the threads started and passed against the broken code. It has no warm-up now. Determinant, single-threaded, best of 15 x 20 iterations, two runs each: before after Laplace 6x6 0.00870 0.00875 0.00838 0.00825 Laplace 9x9 4.430 4.395 4.201 4.197 Adjoint 6x6 0.00808 0.00799 0.00821 0.00783 Laplace comes out about 5% faster. Reading a [ThreadStatic] field costs a thread-local lookup, and doing that at every one of the n! nodes of the Laplace recursion measured +10%, so the pool is taken once at the top and threaded through the recursion; sizes only shrink on the way down. That also replaces the old per-node Count check and shared-list index with one local index, which is where the 5% comes from. 221 of 221 existing tests pass. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Core/Expressions/IteratorCompiler.cs | 20 +- GenericTensor/Functions/Determinant.cs | 14 +- .../Functions/SquareMatrixFactory.cs | 41 +++- UnitTests/Concurrency.cs | 184 ++++++++++++++++++ 4 files changed, 244 insertions(+), 15 deletions(-) create mode 100644 UnitTests/Concurrency.cs diff --git a/GenericTensor/Core/Expressions/IteratorCompiler.cs b/GenericTensor/Core/Expressions/IteratorCompiler.cs index e18f224..2e042a2 100644 --- a/GenericTensor/Core/Expressions/IteratorCompiler.cs +++ b/GenericTensor/Core/Expressions/IteratorCompiler.cs @@ -26,6 +26,7 @@ using System; +using System.Collections.Concurrent; using System.Collections.Generic; using System.ComponentModel; using System.Linq; @@ -345,15 +346,24 @@ public enum OperationType Division } - private static readonly Dictionary<(OperationType opId, int N, bool parallel), Action, GenTensor, GenTensor>> storage - = new Dictionary<(OperationType opId, int N, bool parallel), Action, GenTensor, GenTensor>>(); + // Concurrent, and it has to be: this is written from whichever thread first asks for a + // given operation and rank, and Dictionary<,> is documented as safe for concurrent readers + // only while nobody is writing. An Add that resizes the buckets while another thread is + // walking them does not fail cleanly -- it returns the wrong entry, or spins forever + // inside the lookup. + // + // A duplicate compilation under GetOrAdd is harmless: two compilations of the same + // (operation, rank, parallel) produce delegates that do the same thing, so the loser of + // the race is simply discarded. + private static readonly ConcurrentDictionary<(OperationType opId, int N, bool parallel), Action, GenTensor, GenTensor>> storage + = new ConcurrentDictionary<(OperationType opId, int N, bool parallel), Action, GenTensor, GenTensor>>(); private static Action, GenTensor, GenTensor> GetFunc(int N, Func operation, bool parallel, OperationType ot) { var key = (ot, N, parallel); - if (!storage.ContainsKey(key)) - storage[key] = CompileForNDimensions(N, operation, parallel); - return storage[key]; + if (storage.TryGetValue(key, out var cached)) + return cached; + return storage.GetOrAdd(key, CompileForNDimensions(N, operation, parallel)); } public static GenTensor PiecewiseAdd(GenTensor a, GenTensor b, bool parallel) diff --git a/GenericTensor/Functions/Determinant.cs b/GenericTensor/Functions/Determinant.cs index d67448b..e402850 100644 --- a/GenericTensor/Functions/Determinant.cs +++ b/GenericTensor/Functions/Determinant.cs @@ -25,6 +25,7 @@ #endregion +using System.Collections.Generic; using GenericTensor.Core; namespace GenericTensor.Functions @@ -33,13 +34,22 @@ internal static class Determinant where TWrapper : struct, IOperati { #region Matrix Determinant #region Laplace + // The scratch pool is taken once, here, and threaded through the recursion rather than + // looked up at every level. The pool is [ThreadStatic] -- it has to be, since the caller + // writes into the matrix it is handed -- and a thread-local lookup at every one of the n! + // nodes of this recursion is about 10% of the whole determinant. Sizes only shrink on the + // way down, so growing it to diagLength - 1 once covers every level. internal static T DeterminantLaplace(GenTensor t, int diagLength) + => DeterminantLaplace(t, diagLength, SquareMatrixFactory.GetPool(diagLength - 1)); + + private static T DeterminantLaplace(GenTensor t, int diagLength, + List> pool) { if (diagLength == 1) return t.GetValueNoCheck(0, 0); var det = default(TWrapper).CreateZero(); var sign = default(TWrapper).CreateOne(); - var temp = SquareMatrixFactory.GetMatrix(diagLength - 1); + var temp = pool[diagLength - 2]; for (int i = 0; i < diagLength; i++) { Inversion.GetCofactorMatrix(t, temp, 0, i, diagLength); @@ -48,7 +58,7 @@ internal static T DeterminantLaplace(GenTensor t, int diagLength) sign, default(TWrapper).Multiply( t.GetValueNoCheck(0, i), - DeterminantLaplace(temp, diagLength - 1) + DeterminantLaplace(temp, diagLength - 1, pool) )) ); sign = default(TWrapper).Negate(sign); diff --git a/GenericTensor/Functions/SquareMatrixFactory.cs b/GenericTensor/Functions/SquareMatrixFactory.cs index 2b8f93f..233e193 100644 --- a/GenericTensor/Functions/SquareMatrixFactory.cs +++ b/GenericTensor/Functions/SquareMatrixFactory.cs @@ -25,6 +25,7 @@ #endregion +using System; using System.Collections.Generic; using GenericTensor.Core; using GenericTensor.Functions; @@ -34,17 +35,41 @@ namespace GenericTensor.Core internal static class SquareMatrixFactory where TWrapper : struct, IOperations { // [0] is 1x1 matrix, [1] is 2x2 matrix, etc. - static readonly List> tensorTempFactorySquareMatrices = new List>(); + // + // These are scratch buffers, and the caller writes into the one it is given -- see + // Inversion.GetCofactorMatrix, which fills the matrix returned here. So the pool must be + // per thread. One pool shared by every thread hands two threads computing a determinant + // or an inverse of the same size the same buffer, and each overwrites the other's minor + // between the write and the read. Neither call throws and neither result is malformed; + // they are simply wrong, which is the worst way for this to fail. + // + // [ThreadStatic] cannot carry an initializer -- it would run only on whichever thread + // happened to trigger the static constructor and leave every other thread with null -- + // so the list is created on first use on each thread. + [ThreadStatic] private static List>? tensorTempFactorySquareMatrices; - internal static GenTensor GetMatrix(int diagLength) + /// + /// This thread's pool, grown so that every square matrix up to + /// is in it. Index n - 1 is the n x n one. + /// + /// Callers that ask repeatedly should take the list once and index it, rather than calling + /// in a loop: reading a [ThreadStatic] field costs a thread-local + /// lookup, which is nothing on its own and about 10% of a Laplace determinant when it sits + /// inside that recursion. + /// + internal static List> GetPool(int diagLength) { - if (diagLength >= tensorTempFactorySquareMatrices.Count + 1) - lock (tensorTempFactorySquareMatrices) - if (diagLength >= tensorTempFactorySquareMatrices.Count + 1) - for (int i = tensorTempFactorySquareMatrices.Count + 1; i <= diagLength; i++) - tensorTempFactorySquareMatrices.Add(new GenTensor(i, i)); - return tensorTempFactorySquareMatrices[diagLength - 1]; + // No lock: nothing outside this thread can see the list, which is also why the + // grow-then-index here is safe where the previous shared version was not (it + // indexed after releasing the lock, so a concurrent Add could reallocate underneath). + var matrices = tensorTempFactorySquareMatrices ??= new List>(); + for (var i = matrices.Count + 1; i <= diagLength; i++) + matrices.Add(new GenTensor(i, i)); + return matrices; } + + internal static GenTensor GetMatrix(int diagLength) + => GetPool(diagLength)[diagLength - 1]; } // It is here just for a test to avoid InternalsToVisible diff --git a/UnitTests/Concurrency.cs b/UnitTests/Concurrency.cs new file mode 100644 index 0000000..fd5ceae --- /dev/null +++ b/UnitTests/Concurrency.cs @@ -0,0 +1,184 @@ +#region copyright +/* + * MIT License + * + * Copyright (c) 2020-2021 WhiteBlackGoose + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +#endregion + + +using System; +using System.Collections.Concurrent; +using System.Threading.Tasks; +using GenericTensor.Core; +using GenericTensor.Functions; +using Microsoft.VisualStudio.TestTools.UnitTesting; + +namespace UnitTests +{ + /// + /// The same operation, asked for on one thread and then on many, has to give the same answer. + /// + /// These bind against two shared statics: the scratch-matrix pool behind + /// DeterminantLaplace and Inverse, and the compiled-operation cache behind the + /// piecewise operators. Both used to be process-wide, so the failure is a wrong number rather + /// than an exception, and the assertions below report how many entries disagreed rather than + /// stopping at the first one -- a count is what says whether this is a rare interleaving or + /// the normal case. + /// + [TestClass] + public class Concurrency + { + private const int Matrices = 60; + private const int Size = 5; + + private static GenTensor Build(int seed) + { + var random = new Random(seed); + var m = GenTensor.CreateSquareMatrix(Size); + for (var x = 0; x < Size; x++) + for (var y = 0; y < Size; y++) + m.SetValueNoCheck(random.Next(-4, 5), x, y); + return m; + } + + /// + /// Runs over freshly built matrices sequentially, then again + /// in parallel, and returns the indices whose answers differ. + /// + /// The matrices are rebuilt inside each pass on purpose. Reusing the instances would let + /// any caching inside the library answer the parallel pass from work the sequential pass + /// already did, which is a green test that measured nothing. + /// + private static ConcurrentBag Disagreements( + Func, TResult> operation, + Func same) + { + var expected = new TResult[Matrices]; + for (var i = 0; i < Matrices; i++) + expected[i] = operation(Build(i)); + + var differed = new ConcurrentBag(); + Parallel.For(0, Matrices, i => + { + if (!same(expected[i], operation(Build(i)))) + differed.Add(i); + }); + return differed; + } + + [TestMethod] + public void DeterminantLaplaceIsThreadSafe() + { + var differed = Disagreements( + m => m.DeterminantLaplace(), + (a, b) => a == b); + Assert.AreEqual(0, differed.Count, + $"{differed.Count} of {Matrices} determinants changed value when computed in parallel"); + } + + /// + /// A control rather than a regression test: this passed before the scratch pool was made + /// per thread, and it is here to record that the Gaussian determinant is *not* affected. + /// It works on its own copy of the matrix, so it never asks the pool for anything, which + /// is what bounds the damage of the pool bug to the Laplace determinant and the inverse. + /// + [TestMethod] + public void DeterminantGaussianIsThreadSafe() + { + var differed = Disagreements( + m => m.DeterminantGaussianSafeDivision(), + (a, b) => a == b); + Assert.AreEqual(0, differed.Count, + $"{differed.Count} of {Matrices} determinants changed value when computed in parallel"); + } + + [TestMethod] + public void AdjugateIsThreadSafe() + { + var differed = Disagreements( + m => m.Adjoint(), + (a, b) => a == b); + Assert.AreEqual(0, differed.Count, + $"{differed.Count} of {Matrices} adjugates changed value when computed in parallel"); + } + + /// + /// The compiled-operation cache has to be filled correctly when several threads reach it + /// at once and it is still empty. + /// + /// This one deliberately has no warm-up pass. The first version of this test computed the + /// expected values sequentially and only then went parallel, which populated every key + /// before the threads started -- so the parallel pass did nothing but read a finished + /// dictionary and the test passed against the unsynchronised version too. The cache is + /// keyed on (operation, rank, parallel) and holds a handful of entries, so a warm-up is + /// enough to hide the defect entirely. + /// + /// Cold, it fails against the unsynchronised version with the runtime's own diagnostic: + /// "Operations that change non-concurrent collections must have exclusive access. A + /// concurrent update was performed on this collection and corrupted its state", followed + /// by a lookup missing a key that had been written. Note the failure mode is not + /// guaranteed on any single run -- a torn Dictionary can also return a wrong entry or + /// spin -- so a green result on the old code would not have exonerated it. + /// + [TestMethod] + public void PiecewiseOperationsOnAColdCacheAreThreadSafe() + { + const int Ranks = 5; + var operations = new Func, GenTensor, GenTensor>[] + { + (a, b) => GenTensor.PiecewiseAdd(a, b), + (a, b) => GenTensor.PiecewiseSubtract(a, b), + (a, b) => GenTensor.PiecewiseMultiply(a, b), + }; + + var results = new GenTensor[Ranks * operations.Length]; + Parallel.For(0, results.Length, i => + { + var rank = i / operations.Length + 1; + var left = BuildOfRank(rank, seed: i); + var right = BuildOfRank(rank, seed: i + 1000); + results[i] = operations[i % operations.Length](left, right); + }); + + for (var i = 0; i < results.Length; i++) + { + var rank = i / operations.Length + 1; + var expected = operations[i % operations.Length]( + BuildOfRank(rank, seed: i), BuildOfRank(rank, seed: i + 1000)); + Assert.IsTrue(expected == results[i], + $"the rank-{rank} operation at index {i} gave a different answer on a cold cache"); + } + } + + private static GenTensor BuildOfRank(int rank, int seed) + { + var random = new Random(seed); + var shape = new int[rank]; + for (var i = 0; i < rank; i++) + shape[i] = 2; + var t = new GenTensor(shape); + foreach (var (index, _) in t.Iterate()) + t.SetValueNoCheck(random.Next(-9, 10), index); + return t; + } + } +}