diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 4c35e6911..0daf95d5b 100644 --- a/Sources/.editorconfig +++ b/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 diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index b33ebbf8e..3d70dfe8f 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -5453,7 +5453,13 @@ public static Matrix Concat(Direction dir, params Matrix[] matrices) /// /// 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; + } /// /// Creates an instance of that is a matrix. diff --git a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs b/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs index d204e401c..b1d2467bc 100644 --- a/Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs +++ b/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 inverse; @@ -329,13 +336,16 @@ public Entity AsScalar() /// and then applies inner simplification /// 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); + } /// /// The Subtract operator. Performs an active operation @@ -343,13 +353,14 @@ public Entity AsScalar() /// and then applies inner simplification /// 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); + } /// /// The Multiply operator. Performs an active operation @@ -485,7 +496,11 @@ public IEnumerator 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); diff --git a/Sources/AngouriMath/Functions/GenTensorGuard.cs b/Sources/AngouriMath/Functions/GenTensorGuard.cs new file mode 100644 index 000000000..194e71742 --- /dev/null +++ b/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; + +namespace AngouriMath.Functions +{ + /// + /// Makes the matrix operations that go through GenericTensor safe to call from more than one + /// thread, which documents as supported. + /// + /// + /// 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: + /// + /// The scratch-matrix pool behind DeterminantLaplace and Adjoint 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 + /// . 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. + /// + /// The compiled-operation cache behind the piecewise operators is an unsynchronised + /// Dictionary, 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: 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 , leaves the dictionary read-only from then on. + /// + internal static class GenTensorGuard + { + /// + /// Held across any call that reaches GenericTensor's shared scratch-matrix pool -- which + /// is DeterminantLaplace, Adjoint, and InvertMatrix, since inverting + /// goes through the adjugate. DeterminantGaussianSafeDivision and + /// MatrixMultiply do not touch the pool and are deliberately not covered. + /// + [ConstantField] internal static readonly object ScratchPool = new object(); + + [ConstantField] private static readonly Lazy piecewiseCache = + new Lazy(WarmPiecewiseCache, LazyThreadSafetyMode.ExecutionAndPublication); + + /// + /// 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 + /// . + /// + 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; + } + } +} diff --git a/Sources/Tests/UnitTests/Core/MatrixConcurrencyTest.cs b/Sources/Tests/UnitTests/Core/MatrixConcurrencyTest.cs new file mode 100644 index 000000000..3512fe9fa --- /dev/null +++ b/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 +{ + /// + /// A matrix operation asked for on one thread and then on many has to give the same answer. + /// 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 + /// Functions.GenTensorGuard. 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. + /// + public sealed class MatrixConcurrencyTest + { + private const int Matrices = 40; + + /// + /// The entries are deliberately not polynomial. Determinant tries + /// PolynomialDeterminant.Of 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. + /// + 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)); + } + + /// + /// Runs 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. + /// + private static int Disagreements(Func operation) + { + var expected = new Entity?[Matrices]; + for (var i = 0; i < Matrices; i++) + expected[i] = operation(Build(i)); + + var differed = new ConcurrentBag(); + 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)); + + /// + /// 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. + /// + [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]); + } + } + } +}