From 9110c48ab7c6d80a229b4085788674b9fc6a2892 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 9 Sep 2026 10:47:08 +0000 Subject: [PATCH] Matrix operations answer wrongly when called from more than one thread MathS.Multithreading documents concurrent use as supported, and matrix operations do not honour it. The cause is not in this repository: GenericTensor 1.0.4 keeps two process-wide mutable statics, reported upstream as https://github.com/asc-community/GenericTensor/issues/40 with a fix in https://github.com/asc-community/GenericTensor/pull/41. This is what stands in until a release carrying that fix exists. Measured on c2c8eb83, forty 4x4 matrices with non-polynomial entries, each computed once sequentially and then rebuilt and recomputed under Parallel.For: Determinant 38 of 40 disagree Inverse 18 of 40 disagree Adjugate 11 of 40 disagree m + m, m - m, PointwiseMultiplication throws, on a corrupted Dictionary Nothing about the first three looks wrong. They are well-formed entities of the right shape carrying another computation's values. The two statics need opposite treatments, because only one 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 three call sites take a lock. That serialises determinants and inverses of matrices the polynomial elimination declines -- PolynomialDeterminant.Of runs first and never reaches the pool, which is why the entries in the test are sines and cosines, and why a polynomial matrix was not affected at all. The compiled-operation cache behind the elementwise operators is an unsynchronised Dictionary, so it is only unsafe while being filled; once an entry is there the reads are pure. That one needs no lock, because the set of entries this library can ever ask for is fixed: Entity.Matrix rejects any tensor that is not rank 2, and every call here uses the default single-threaded mode, so the only reachable keys are addition, subtraction and multiplication at rank 2. Filling all three once under Lazy leaves the dictionary read-only, at no steady-state cost. Answers are unchanged, so there is no BREAKING-CHANGES entry: a single-threaded caller computed the same values before and computes them now. The tests rebuild their matrices inside each pass, which is load-bearing rather than tidiness. Determinant, Inverse and Adjugate are cached lazy properties, so reusing the instances lets the parallel pass read what the sequential pass already computed -- a green test that ran the operation once. The elementwise test has no sequential pass at all for the mirror-image reason: a warm-up populates every cache key before the threads start and passes against the unguarded code. Full suite green: 8499 and 1485. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- Sources/.editorconfig | 6 + Sources/AngouriMath/Convenience/MathS.cs | 8 +- .../Core/Entity/Omni/Entity.Matrix.cs | 49 ++++--- .../AngouriMath/Functions/GenTensorGuard.cs | 75 +++++++++++ .../UnitTests/Core/MatrixConcurrencyTest.cs | 121 ++++++++++++++++++ 5 files changed, 241 insertions(+), 18 deletions(-) create mode 100644 Sources/AngouriMath/Functions/GenTensorGuard.cs create mode 100644 Sources/Tests/UnitTests/Core/MatrixConcurrencyTest.cs 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]); + } + } + } +}