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; + } + } +}