Matrix operations give wrong answers when called from more than one thread

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<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>>> storage 
-                 = new Dictionary<(OperationType opId, int N, bool parallel), Action<GenTensor<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>>>();
+        // 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<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>>> storage
+                 = new ConcurrentDictionary<(OperationType opId, int N, bool parallel), Action<GenTensor<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>>>();
 
         private static Action<GenTensor<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>> GetFunc(int N, Func<Expression, Expression, Expression> 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<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, bool parallel)

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<T, TWrapper> 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, TWrapper> t, int diagLength)
+            => DeterminantLaplace(t, diagLength, SquareMatrixFactory<T, TWrapper>.GetPool(diagLength - 1));
+
+        private static T DeterminantLaplace(GenTensor<T, TWrapper> t, int diagLength,
+            List<GenTensor<T, TWrapper>> pool)
         {
             if (diagLength == 1)
                 return t.GetValueNoCheck(0, 0);
             var det = default(TWrapper).CreateZero();
             var sign = default(TWrapper).CreateOne();
-            var temp = SquareMatrixFactory<T, TWrapper>.GetMatrix(diagLength - 1);
+            var temp = pool[diagLength - 2];
             for (int i = 0; i < diagLength; i++)
             {
                 Inversion<T, TWrapper>.GetCofactorMatrix(t, temp, 0, i, diagLength);
@@ -48,7 +58,7 @@ internal static T DeterminantLaplace(GenTensor<T, TWrapper> 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);

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<T, TWrapper> where TWrapper : struct, IOperations<T>
     {
         // [0] is 1x1 matrix, [1] is 2x2 matrix, etc.
-        static readonly List<GenTensor<T, TWrapper>> tensorTempFactorySquareMatrices = new List<GenTensor<T, TWrapper>>();
+        //
+        // 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<GenTensor<T, TWrapper>>? tensorTempFactorySquareMatrices;
 
-        internal static GenTensor<T, TWrapper> GetMatrix(int diagLength)
+        /// <summary>
+        /// This thread's pool, grown so that every square matrix up to <paramref name="diagLength"/>
+        /// 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
+        /// <see cref="GetMatrix"/> 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.
+        /// </summary>
+        internal static List<GenTensor<T, TWrapper>> 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<T, TWrapper>(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<GenTensor<T, TWrapper>>();
+            for (var i = matrices.Count + 1; i <= diagLength; i++)
+                matrices.Add(new GenTensor<T, TWrapper>(i, i));
+            return matrices;
         }
+
+        internal static GenTensor<T, TWrapper> GetMatrix(int diagLength)
+            => GetPool(diagLength)[diagLength - 1];
     }
 
     // It is here just for a test to avoid InternalsToVisible

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
+{
+    /// <summary>
+    /// 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
+    /// <c>DeterminantLaplace</c> and <c>Inverse</c>, 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.
+    /// </summary>
+    [TestClass]
+    public class Concurrency
+    {
+        private const int Matrices = 60;
+        private const int Size = 5;
+
+        private static GenTensor<int, IntWrapper> Build(int seed)
+        {
+            var random = new Random(seed);
+            var m = GenTensor<int, IntWrapper>.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;
+        }
+
+        /// <summary>
+        /// Runs <paramref name="operation"/> 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.
+        /// </summary>
+        private static ConcurrentBag<int> Disagreements<TResult>(
+            Func<GenTensor<int, IntWrapper>, TResult> operation,
+            Func<TResult, TResult, bool> same)
+        {
+            var expected = new TResult[Matrices];
+            for (var i = 0; i < Matrices; i++)
+                expected[i] = operation(Build(i));
+
+            var differed = new ConcurrentBag<int>();
+            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");
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        [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");
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        [TestMethod]
+        public void PiecewiseOperationsOnAColdCacheAreThreadSafe()
+        {
+            const int Ranks = 5;
+            var operations = new Func<GenTensor<int, IntWrapper>, GenTensor<int, IntWrapper>, GenTensor<int, IntWrapper>>[]
+            {
+                (a, b) => GenTensor<int, IntWrapper>.PiecewiseAdd(a, b),
+                (a, b) => GenTensor<int, IntWrapper>.PiecewiseSubtract(a, b),
+                (a, b) => GenTensor<int, IntWrapper>.PiecewiseMultiply(a, b),
+            };
+
+            var results = new GenTensor<int, IntWrapper>[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<int, IntWrapper> 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<int, IntWrapper>(shape);
+            foreach (var (index, _) in t.Iterate())
+                t.SetValueNoCheck(random.Next(-9, 10), index);
+            return t;
+        }
+    }
+}