Matrix operations answer wrongly when called from more than one thread

Sources/.editorconfig

@@ -273,3 +273,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
 
 [Tests/UnitTests/Core/Transformations/RuleEffectsMeasuredTest.cs]
 file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
+
+[AngouriMath/Functions/GenTensorGuard.cs]
+file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
+
+[Tests/UnitTests/Core/MatrixConcurrencyTest.cs]
+file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n

Sources/AngouriMath/Convenience/MathS.cs

@@ -5453,7 +5453,13 @@ public static Matrix Concat(Direction dir, params Matrix[] matrices)
             /// </code>
             /// </example>
             public static Matrix PointwiseMultiplication(Matrix m1, Matrix m2)
-                => (Matrix)new Matrix(GenTensor.PiecewiseMultiply(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified;
+            {
+                // GenericTensor caches the compiled elementwise loop in an unsynchronised
+                // dictionary, so the entry has to be there before several threads read it.
+                // See Functions.GenTensorGuard.
+                Functions.GenTensorGuard.EnsurePiecewiseCacheWarmed();
+                return (Matrix)new Matrix(GenTensor.PiecewiseMultiply(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified;
+            }
 
             /// <summary>
             /// Creates an instance of <see cref="Entity.Matrix"/> that is a matrix.

Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs

@@ -301,7 +301,11 @@ public Entity AsScalar()
                     if (Functions.PolynomialDeterminant.Of(@this.RowCount, (r, c) => @this[r, c])
                         is { } byElimination)
                         return byElimination.InnerSimplified;
-                    return @this.InnerMatrix.DeterminantLaplace().InnerSimplified;
+                    // GenericTensor's Laplace determinant writes into a process-wide scratch
+                    // matrix, so two threads here silently corrupt each other's minors. See
+                    // Functions.GenTensorGuard.
+                    lock (Functions.GenTensorGuard.ScratchPool)
+                        return @this.InnerMatrix.DeterminantLaplace().InnerSimplified;
                 },
                 this
                 );
@@ -317,7 +321,10 @@ public Entity AsScalar()
                     return null;
                 if (@this.Determinant == 0)
                     return null;
-                cp.InvertMatrix();
+                // Inverting goes through the adjugate, which takes its minors in the same
+                // process-wide scratch matrix the determinant does. See Functions.GenTensorGuard.
+                lock (Functions.GenTensorGuard.ScratchPool)
+                    cp.InvertMatrix();
                 return ToMatrix(new Matrix(cp).InnerSimplified);
             }, this);
             private LazyPropertyA<Matrix?> inverse;
@@ -329,27 +336,31 @@ public Entity AsScalar()
             /// and then applies inner simplification
             /// </summary>
             public static Matrix operator +(Matrix m1, Matrix m2)
-                =>
-                m1.InnerMatrix.Shape != m2.InnerMatrix.Shape
-                ?
-                throw new InvalidMatrixOperationException(
-                    $"Cannot add matrices or vectors of shapes {m1.InnerMatrix.Shape} and {m2.InnerMatrix.Shape}")
-                :
-                ToMatrix(new Matrix(GenTensor.PiecewiseAdd(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified);
+            {
+                if (m1.InnerMatrix.Shape != m2.InnerMatrix.Shape)
+                    throw new InvalidMatrixOperationException(
+                        $"Cannot add matrices or vectors of shapes {m1.InnerMatrix.Shape} and {m2.InnerMatrix.Shape}");
+                // GenericTensor caches the compiled elementwise loop in an unsynchronised
+                // dictionary, so the entry has to be there before several threads read it.
+                // See Functions.GenTensorGuard.
+                Functions.GenTensorGuard.EnsurePiecewiseCacheWarmed();
+                return ToMatrix(new Matrix(GenTensor.PiecewiseAdd(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified);
+            }
 
             /// <summary>
             /// The Subtract operator. Performs an active operation
             /// (elementwise subtraction of two matrices or vectors)
             /// and then applies inner simplification
             /// </summary>
             public static Matrix operator -(Matrix m1, Matrix m2)
-                =>
-                m1.InnerMatrix.Shape != m2.InnerMatrix.Shape
-                ?
-                throw new InvalidMatrixOperationException(
-                    $"Cannot subtract matrices or vectors of shapes {m1.InnerMatrix.Shape} and {m2.InnerMatrix.Shape}")
-                :
-                ToMatrix(new Matrix(GenTensor.PiecewiseSubtract(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified);
+            {
+                if (m1.InnerMatrix.Shape != m2.InnerMatrix.Shape)
+                    throw new InvalidMatrixOperationException(
+                        $"Cannot subtract matrices or vectors of shapes {m1.InnerMatrix.Shape} and {m2.InnerMatrix.Shape}");
+                // See Functions.GenTensorGuard, as for the addition above.
+                Functions.GenTensorGuard.EnsurePiecewiseCacheWarmed();
+                return ToMatrix(new Matrix(GenTensor.PiecewiseSubtract(m1.InnerMatrix, m2.InnerMatrix)).InnerSimplified);
+            }
 
             /// <summary>
             /// The Multiply operator. Performs an active operation
@@ -485,7 +496,11 @@ public IEnumerator<Entity> GetEnumerator()
                     {
                         if (!@this.IsSquare)
                             return null;
-                        var innerSimplified = new Matrix(@this.InnerMatrix.Adjoint()).InnerSimplified;
+                        // The adjugate takes every minor in one process-wide scratch matrix.
+                        // See Functions.GenTensorGuard.
+                        Entity innerSimplified;
+                        lock (Functions.GenTensorGuard.ScratchPool)
+                            innerSimplified = new Matrix(@this.InnerMatrix.Adjoint()).InnerSimplified;
                         return ToMatrix(innerSimplified);
                     },
                     this);

Sources/AngouriMath/Functions/GenTensorGuard.cs

@@ -0,0 +1,75 @@
+//
+// Copyright (c) 2019-2026 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using System;
+using System.Threading;
+using GenTensor = GenericTensor.Core.GenTensor<AngouriMath.Entity, AngouriMath.Entity.Matrix.EntityTensorWrapperOperations>;
+
+namespace AngouriMath.Functions
+{
+    /// <summary>
+    /// Makes the matrix operations that go through GenericTensor safe to call from more than one
+    /// thread, which <see cref="MathS.Multithreading"/> documents as supported.
+    /// </summary>
+    /// <remarks>
+    /// GenericTensor 1.0.4 keeps two process-wide mutable statics, and neither fails loudly.
+    /// Reported upstream as https://github.com/asc-community/GenericTensor/issues/40 with a fix in
+    /// https://github.com/asc-community/GenericTensor/pull/41; this type is what stands in until a
+    /// release carrying that fix exists, and it should be deleted when one does.
+    ///
+    /// The two need opposite treatments, because only one of them can be closed without a lock:
+    ///
+    /// <b>The scratch-matrix pool</b> behind <c>DeterminantLaplace</c> and <c>Adjoint</c> hands
+    /// every caller the same tensor for a given size, and the caller writes into it. There is
+    /// nothing to warm and no way to avoid it from out here, so those calls take
+    /// <see cref="ScratchPool"/>. Measured on GenericTensor's own suite, sixty 5x5 matrices
+    /// computed sequentially and then again in parallel: 53 of 60 Laplace determinants and 60 of
+    /// 60 adjugates came back with different values. Serialising them is the cost of not doing
+    /// that.
+    ///
+    /// <b>The compiled-operation cache</b> behind the piecewise operators is an unsynchronised
+    /// <c>Dictionary</c>, so it is only unsafe while it is being filled -- once an entry is there
+    /// the reads are pure. That one needs no lock at all, because the set of entries this library
+    /// can ever ask for is fixed and tiny: <see cref="Entity.Matrix"/> rejects any tensor that is
+    /// not rank 2, and every call here uses the default single-threaded mode, so the only keys
+    /// reachable are addition, subtraction and multiplication at rank 2. Filling all three once,
+    /// under <see cref="Lazy{T}"/>, leaves the dictionary read-only from then on.
+    /// </remarks>
+    internal static class GenTensorGuard
+    {
+        /// <summary>
+        /// Held across any call that reaches GenericTensor's shared scratch-matrix pool -- which
+        /// is <c>DeterminantLaplace</c>, <c>Adjoint</c>, and <c>InvertMatrix</c>, since inverting
+        /// goes through the adjugate. <c>DeterminantGaussianSafeDivision</c> and
+        /// <c>MatrixMultiply</c> do not touch the pool and are deliberately not covered.
+        /// </summary>
+        [ConstantField] internal static readonly object ScratchPool = new object();
+
+        [ConstantField] private static readonly Lazy<bool> piecewiseCache =
+            new Lazy<bool>(WarmPiecewiseCache, LazyThreadSafetyMode.ExecutionAndPublication);
+
+        /// <summary>
+        /// Call before any piecewise operator. The first caller fills GenericTensor's compiled-
+        /// operation cache single-threaded; everyone after that pays a read of an already-computed
+        /// <see cref="Lazy{T}"/>.
+        /// </summary>
+        internal static void EnsurePiecewiseCacheWarmed() => _ = piecewiseCache.Value;
+
+        private static bool WarmPiecewiseCache()
+        {
+            // Rank 2, because that is the only rank Entity.Matrix admits, and one element, because
+            // the point is to insert the cache entry rather than to compute anything. The results
+            // are discarded.
+            var seed = new GenTensor(1, 1);
+            seed.SetValueNoCheck(Entity.Number.Integer.Zero, 0, 0);
+            GenTensor.PiecewiseAdd(seed, seed);
+            GenTensor.PiecewiseSubtract(seed, seed);
+            GenTensor.PiecewiseMultiply(seed, seed);
+            return true;
+        }
+    }
+}

Sources/Tests/UnitTests/Core/MatrixConcurrencyTest.cs

@@ -0,0 +1,121 @@
+//
+// Copyright (c) 2019-2026 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using AngouriMath;
+using System;
+using System.Collections.Concurrent;
+using System.Threading.Tasks;
+using Xunit;
+using static AngouriMath.Entity;
+
+namespace AngouriMath.Tests.Core
+{
+    /// <summary>
+    /// A matrix operation asked for on one thread and then on many has to give the same answer.
+    /// <see cref="MathS.Multithreading"/> documents concurrent use as supported, so this is a
+    /// claim about the public surface rather than about the test runner.
+    ///
+    /// These bind against GenericTensor's process-wide statics, reported as
+    /// https://github.com/asc-community/GenericTensor/issues/40 and worked around in
+    /// <c>Functions.GenTensorGuard</c>. Measured on the commit before that guard: 38 of 40
+    /// determinants, 18 of 40 inverses and 11 of 40 adjugates came back with different values,
+    /// and the elementwise operators threw on a corrupted dictionary. The first three fail by
+    /// returning a wrong entity rather than by throwing, so those assertions count disagreements
+    /// instead of stopping at the first one.
+    /// </summary>
+    public sealed class MatrixConcurrencyTest
+    {
+        private const int Matrices = 40;
+
+        /// <summary>
+        /// The entries are deliberately not polynomial. <c>Determinant</c> tries
+        /// <c>PolynomialDeterminant.Of</c> first and only falls back to GenericTensor's Laplace
+        /// where that declines, so a matrix of polynomials never reaches the shared scratch buffer
+        /// -- which is exactly why the first version of this test found nothing.
+        /// </summary>
+        private static Matrix Build(int seed)
+        {
+            var random = new Random(seed);
+            var a = MathS.Var("a");
+            return MathS.Matrices.Matrix(4, 4,
+                MathS.Sin(a) * random.Next(1, 4), random.Next(-3, 4), MathS.Cos(a), random.Next(-3, 4),
+                random.Next(-3, 4), MathS.Cos(a) * random.Next(1, 4), random.Next(-3, 4), MathS.Sin(a),
+                MathS.Sin(a), random.Next(-3, 4), MathS.Sin(a) * random.Next(1, 4), random.Next(-3, 4),
+                random.Next(-3, 4), MathS.Sin(a), random.Next(-3, 4), MathS.Cos(a) * random.Next(1, 4));
+        }
+
+        /// <summary>
+        /// Runs <paramref name="operation"/> over freshly built matrices sequentially, then again
+        /// in parallel, and counts the ones whose answers differ.
+        ///
+        /// The matrices are rebuilt inside each pass, and that is load-bearing. Determinant,
+        /// Inverse and Adjugate are all cached lazy properties, so reusing the instances would let
+        /// the parallel pass read the values the sequential pass had already computed. That test
+        /// passes against broken code because it never runs the operation twice.
+        /// </summary>
+        private static int Disagreements(Func<Matrix, Entity?> operation)
+        {
+            var expected = new Entity?[Matrices];
+            for (var i = 0; i < Matrices; i++)
+                expected[i] = operation(Build(i));
+
+            var differed = new ConcurrentBag<int>();
+            Parallel.For(0, Matrices, i =>
+            {
+                if (operation(Build(i)) != expected[i])
+                    differed.Add(i);
+            });
+            return differed.Count;
+        }
+
+        [Fact]
+        public void DeterminantIsTheSameOnManyThreads()
+            => Assert.Equal(0, Disagreements(m => m.Determinant));
+
+        [Fact]
+        public void AdjugateIsTheSameOnManyThreads()
+            => Assert.Equal(0, Disagreements(m => m.Adjugate));
+
+        [Fact]
+        public void InverseIsTheSameOnManyThreads()
+            => Assert.Equal(0, Disagreements(m => m.Inverse));
+
+        /// <summary>
+        /// The elementwise operators reach a different shared static: an unsynchronised dictionary
+        /// of compiled loops. It is only unsafe while it is being filled, so this one has no
+        /// sequential warm-up -- a warm-up would populate every key before the threads started and
+        /// pass against the unguarded code.
+        /// </summary>
+        [Fact]
+        public void ElementwiseOperatorsAreTheSameOnManyThreads()
+        {
+            var results = new Entity?[Matrices];
+            Parallel.For(0, Matrices, i =>
+            {
+                var m = Build(i);
+                results[i] = (i % 3) switch
+                {
+                    0 => m + m,
+                    1 => m - m,
+                    _ => MathS.Matrices.PointwiseMultiplication(m, m)
+                };
+            });
+
+            for (var i = 0; i < Matrices; i++)
+            {
+                var m = Build(i);
+                Entity expected = (i % 3) switch
+                {
+                    0 => m + m,
+                    1 => m - m,
+                    _ => MathS.Matrices.PointwiseMultiplication(m, m)
+                };
+                Assert.Equal(expected, results[i]);
+            }
+        }
+    }
+}