Problem on fma data sync or thread load unbalance

StgSharp.Common/Mathematics/Linear/MatrixParallel.Task.cs

@@ -41,21 +41,37 @@ using System.Threading;
 
 namespace StgSharp.Mathematics
 {
-    [StructLayout(LayoutKind.Explicit)]
+    [StructLayout(LayoutKind.Explicit, Size = 32, Pack = 8)]
+    public unsafe struct MatrixParallelEnumeration
+    {
+
+        [FieldOffset(00)] public MatrixKernel* kernel;
+
+        [FieldOffset(08)] public int begin;
+        [FieldOffset(12)] public int end;
+        [FieldOffset(16)] public int primCount;
+        [FieldOffset(20)] public int primStride;
+        [FieldOffset(24)] public int secCount;
+        [FieldOffset(28)] public int secStride;
+
+    }
+
+    [StructLayout(LayoutKind.Explicit, Size = 128, Pack = 8)]
     public unsafe struct MatrixParallelTaskPackage
     {
 
-        [FieldOffset(0)] public MatrixKernel* left;
+        [FieldOffset(72)] public MatrixKernel* kernel;
-        [FieldOffset(8)] public MatrixKernel* result;
+        [FieldOffset(0)] public MatrixKernel* Left;
-        [FieldOffset(16)] public MatrixKernel* right;
+        [FieldOffset(32)] public MatrixKernel* Result;
-        [FieldOffset(40)] public ScalarDataPacket* scalerPacket;
+        [FieldOffset(64)] public MatrixKernel* Right;
-        [FieldOffset(24)] public int ComputeMode;
+        [FieldOffset(96)] public ScalarDataPacket* scalerPacket;
-        [FieldOffset(60)] public int count;
+        [FieldOffset(100)] public int ComputeMode;
+        [FieldOffset(104)] public int count;
+        [FieldOffset(108)] public IntPtr executionHandle;
 
-        [FieldOffset(48)] public int leftStride;
+        [FieldOffset(0)] public MatrixParallelEnumeration LeftEnum;
-        [FieldOffset(56)] public int resultStride;
+        [FieldOffset(32)] public MatrixParallelEnumeration ResultEnum;
-        [FieldOffset(52)] public int rightStride;
+        [FieldOffset(64)] public MatrixParallelEnumeration RightEnum;
-        [FieldOffset(32)] public IntPtr executionHandle;
 
     }
 
@@ -87,43 +103,4 @@ namespace StgSharp.Mathematics
         }
 
     }
-
-    public unsafe class MatrixParallelWrap
-    {
-
-        public Queue<IntPtr>[] cache;
-
-        public MatrixParallelWrap(int wrapSize)
-        {
-            cache = new Queue<IntPtr>[wrapSize];
-            for (int i = 0; i < cache.Length; i++) {
-                cache[i] = new Queue<IntPtr>(8);
-            }
-        }
-
-        public Span<IntPtr> EnqueuePackage(Span<IntPtr> packageSpan)
-        {
-            int idx = 0;
-            for (; idx < cache.Length; idx++) {
-                cache[idx].Enqueue(packageSpan[idx]);
-            }
-            return packageSpan[idx..];
-        }
-
-        public void EnqueuePackage(MatrixParallelTaskPackage* task)
-        {
-            Span<IntPtr> span = stackalloc IntPtr[cache.Length];
-            MatrixParallelFactory.
-
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public bool TryGetNextWrapTask(int idInWrap, out IntPtr handle)
-        {
-            return cache[idInWrap].TryDequeue(out handle);
-        }
-
-        public record class Handle(int InWrapIndex, MatrixParallelWrap Wrap);
-
-    }
 }

StgSharp.Common/Mathematics/Linear/MatrixParallel.cs

@@ -45,24 +45,49 @@ using System.Threading.Tasks;
 
 namespace StgSharp.Mathematics
 {
+    /*
+    * Thread wrapping to accelerate matrix operations by dividing threads into groups (wraps).
+    * For M wraps, each containing an average of N threads, we satisfy the condition M * N = T - 4,
+    * where T is the total thread count of the CPU. Typically, for modern consumer CPUs, T is around 16.
+    * For simplicity, let M * N = 16.
+    *
+    * The total time complexity for publishing tasks is O(M + N):
+    * - O(M) accounts for distributing tasks across M wraps.
+    * - O(N) accounts for threads within each wrap fetching tasks from the wrap queue.
+    *
+    * Using the inequality Sqrt(M * N) <= (M + N) / 2 (derived from the arithmetic-geometric mean inequality),
+    * the optimal configuration occurs when M = N. For M * N = 16, this gives M = N = 4.
+    * 
+    * To ensure compatibility across a wide range of CPU architectures, 
+    * we allow each wrap to contain either 3 or 4 threads, providing flexibility while maintaining efficiency.
+    */
+
     internal static class MatrixParallel
     {
 
-        // TODO 解决多线程协调问题，barrier性能不好，MRES时机不好控制，可能可以使用双MRES
         private static Thread[] _threads;
         private static MatrixParallelWrap[] _wraps;
         private static readonly (int a, int b)[] LookupTable =
-        [(1, 1), (0, 2), (3, 0), (2, 1), (1, 2),
+        [(1, 1), (0, 2), (3, 0), (2, 1),
-         (0, 3), (3, 1), (2, 2), (1, 3), (4, 1),
+         (1, 2), (0, 3), (3, 1), (2, 2), (1, 3),
-         (3, 2), (2, 3), (1, 4), (0, 5),];
+         (0, 4), (3, 2), (2, 3), (1, 4), (0, 5),];
         private static volatile int RemainThreadCount = 0;
 
-        private static readonly ManualResetEventSlim _resetEvent = new(false);
+        private static readonly ManualResetEventSlim _restartEvent = new(false);
+        private static readonly ManualResetEventSlim _startEvent = new(false);
+
+        private static MatrixParallelTaskPackage _totalTask;
         private static readonly object _lock = new();
         private static ThreadPriority _AfterWork = ThreadPriority.Lowest;
 
+        public static ref MatrixParallelTaskPackage TotalTask => ref _totalTask;
+
         public static bool SupportParallel { get; } = Environment.ProcessorCount > 6;
 
+        public static bool IsParallelMeaningful { get; } = Environment.ProcessorCount > 6;
+
+        internal static MatrixParallelWrap.Handle LeaderHandle { get; private set; }
+
         public static void EndParallel()
         {
             Monitor.Exit(_lock);
@@ -89,99 +114,19 @@ namespace StgSharp.Mathematics
             count = count == 5 ? 4 : count;
             int _3wrap = 0, _4wrap = 0;
             _threads = new Thread[count];
-            if (count < 6)
+            if (count < 6) { } else if (count == 6) { } else { }
-            {
+            LeaderHandle = new MatrixParallelWrap.Handle(0, _wraps[0]);
-                _wraps = new MatrixParallelWrap[1];
-                _wraps[0] = new MatrixParallelWrap(count);
-                for (int i = 1; i < count; i++)                     //first is reverse for leader thread
-                {
-                    _threads[i] = new Thread(MatrixParallelCycle)
-                    {
-                        IsBackground = true,
-                        Priority = ThreadPriority.Highest,
-                        Name = $"StgSharp Matrix Parallel Worker {i}"
-                    };
-                    _threads[i].Start(new MatrixParallelWrap.Handle(i, _wraps[0]));
-                }
-            } else if (count == 6)
-            {
-                _3wrap = 2;
-                _threads = new Thread[6];
-                _wraps = new MatrixParallelWrap[2];
-                _wraps[0] = new(3);
-                _wraps[1] = new(3);
-                for (int i = 1; i < 3; i++)                          //first is reverse for leader thread
-                {
-                    _threads[i] = new Thread(MatrixParallelCycle)
-                    {
-                        IsBackground = true,
-                        Priority = ThreadPriority.Highest,
-                        Name = $"StgSharp Matrix Parallel Worker {i}"
-                    };
-                    _threads[i].Start(new MatrixParallelWrap.Handle(i, _wraps[0]));
-                }
-                for (int i = 3; i < 6; i++)
-                {
-                    _threads[i] = new Thread(MatrixParallelCycle)
-                    {
-                        IsBackground = true,
-                        Priority = ThreadPriority.Highest,
-                        Name = $"StgSharp Matrix Parallel Worker {i}"
-                    };
-                    _threads[i].Start(new MatrixParallelWrap.Handle(i - 3, _wraps[1]));
-                }
-            } else
-            {
-                (_3wrap, _4wrap) = FastDecompose(count);
-                _wraps = new MatrixParallelWrap[_3wrap + _4wrap];
-                int i = 0, wi = 0;
-                for (; i < 3 * _3wrap; i += 3)
-                {
-                    MatrixParallelWrap wrap = new(3);
-                    _wraps[wi] = wrap;
-                    wi++;
-                    for (int j = i; j < i + 3; j++)
-                    {
-                        if (j == 0)
-                        {
-                            continue;                               //first is reverse for leader thread
-                        }
-                        _threads[j] = new Thread(MatrixParallelCycle)
-                        {
-                            IsBackground = true,
-                            Priority = ThreadPriority.Highest,
-                            Name = $"StgSharp Matrix Parallel Worker {j}"
-                        };
-                        _threads[j].Start(new MatrixParallelWrap.Handle(j, wrap));
-                    }
-                }
-                for (; i < count; i += 4)
-                {
-                    MatrixParallelWrap wrap = new(4);
-                    _wraps[wi] = wrap;
-                    wi++;
-                    for (int j = i; j < i + 4; j++)
-                    {
-                        _threads[j] = new Thread(MatrixParallelCycle)
-                        {
-                            IsBackground = true,
-                            Priority = ThreadPriority.Highest,
-                            Name = $"StgSharp Matrix Parallel Worker {j}"
-                        };
-                        _threads[j].Start(new MatrixParallelWrap.Handle(j, wrap));
-                    }
-                }
-            }
         }
 
         public static void LaunchParallel(SleepMode afterWork)
         {
-            _resetEvent.Set();
+            _startEvent.Set();
+            _restartEvent.Reset();
             RemainThreadCount = _threads.Length;
             _AfterWork = (ThreadPriority)afterWork;
-            _resetEvent.Reset();
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static void LeadParallel()
         {
             Monitor.Enter(_lock);
@@ -189,7 +134,67 @@ namespace StgSharp.Mathematics
 
         public static void PublishTask() { }
 
-        private static unsafe void MatrixParallelCycle(object obj)
+        public unsafe class MatrixParallelWrap
+        {
+
+            private Queue<IntPtr>[] _cache;
+            private Thread[] _threads;
+            private MatrixParallelTaskPackage _wrapTask;
+
+            public MatrixParallelWrap(int id, int wrapSize, bool isFirst)
+            {
+                ID = id;
+                _cache = new Queue<IntPtr>[wrapSize];
+                for (int i = 0; i < _cache.Length; i++) {
+                    _cache[i] = new Queue<IntPtr>(8);
+                }
+                _threads = new Thread[wrapSize];
+                int begin = isFirst ? 0 : 1,end = wrapSize - 1;
+                for (int i = 0; i < end; i++)
+                {
+                    _threads[i] = new Thread(CommonMatrixParallelCycle)
+                    {
+                        IsBackground = true,
+                        Priority = ThreadPriority.BelowNormal,
+                        Name = $"MPT#{ID}_{i}"
+                    };
+                    _threads[i].Start(new Handle(i, this));
+                }
+                _threads[end] = new Thread(SchedulerMatrixParallelCycle)
+                {
+                    IsBackground = true,
+                    Priority = ThreadPriority.BelowNormal,
+                    Name = $"MPT#{ID}_{end}*"
+                };
+                _threads[end].Start(new Handle(end, this));
+            }
+
+            public ref MatrixParallelTaskPackage WrapTask => ref _wrapTask;
+
+            public int ID { get; init; }
+
+            public Span<IntPtr> EnqueuePackage(Span<IntPtr> packageSpan)
+            {
+                int idx = 0;
+                for (; idx < _cache.Length; idx++) {
+                    _cache[idx].Enqueue(packageSpan[idx]);
+                }
+                return packageSpan[idx..];
+            }
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public bool TryGetNextWrapTask(int idInWrap, out IntPtr handle)
+            {
+                return _cache[idInWrap].TryDequeue(out handle);
+            }
+
+            public record class Handle(int InWrapIndex, MatrixParallelWrap Wrap);
+
+        }
+
+        #region ThreadOperation
+
+        private static unsafe void CommonMatrixParallelCycle(object obj)
         {
             MatrixParallelWrap.Handle handle = (obj as MatrixParallelWrap.Handle)!;
             MatrixParallelWrap wrap = handle.Wrap;
@@ -197,63 +202,88 @@ namespace StgSharp.Mathematics
             Thread self = Thread.CurrentThread;
             while (true)
             {
-                while (wrap.TryGetNextWrapTask(inWrapIndex, out IntPtr ptr))
+                self.Priority = ThreadPriority.Highest;
-                {
+                while (wrap.TryGetNextWrapTask(inWrapIndex, out IntPtr ptr)) {
-                    MatrixParallelTaskPackage* p = (MatrixParallelTaskPackage*)ptr;
+                    ExecuteTaskPackage((MatrixParallelTaskPackage*)ptr);
-                    MatrixKernel* left = p->left, right = p->left, result = p->result;
-                    int count = p->count;
-                    ScalarDataPacket* scalar = p->scalerPacket;
-                    switch (p->ComputeMode)
-                    {
-                        case 0:
-                            delegate*<MatrixKernel*, MatrixKernel*, void> op_0 = (delegate*<MatrixKernel*, MatrixKernel*, void>)p->executionHandle;
-                            for (int i = 0; i < count; i++)
-                            {
-                                op_0(left, result);
-                                left += p->leftStride;
-                                result += p->resultStride;
-                            }
-                            break;
-                        case 1:
-                            delegate*<MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void> op_1 = (delegate*<MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void>)p->executionHandle;
-                            for (int i = 0; i < count; i++)
-                            {
-                                op_1(left, result, scalar);
-                                left += p->leftStride;
-                                result += p->resultStride;
-                            }
-                            break;
-                        case 2:
-                            delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, void> op_2 = (delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, void>)p->executionHandle;
-                            for (int i = 0; i < count; i++)
-                            {
-                                op_2(left, right, result);
-                                right += p->rightStride;
-                                left += p->leftStride;
-                                result += p->resultStride;
-                            }
-                            break;
-                        case 3:
-                            delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void> op_3 = (delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void>)p->executionHandle;
-                            for (int i = 0; i < count; i++)
-                            {
-                                op_3(left, right, result, scalar);
-                                right += p->rightStride;
-                                left += p->leftStride;
-                                result += p->resultStride;
-                            }
-                            break;
-                        default:
-                            break;
-                    }
                 }
                 _ = Interlocked.Decrement(ref RemainThreadCount);
                 self.Priority = _AfterWork;
-                _resetEvent.Wait();
+                _restartEvent.Wait();
-                self.Priority = ThreadPriority.Highest;
             }
         }
 
+        private static unsafe void ExecuteTaskPackage(MatrixParallelTaskPackage* p)
+        {
+            MatrixKernel* left = p->left, right = p->left, result = p->result;
+            int count = p->count;
+            ScalarDataPacket* scalar = p->scalerPacket;
+            switch (p->ComputeMode)
+            {
+                case 0:
+                    delegate*<MatrixKernel*, MatrixKernel*, void> op_0 = (delegate*<MatrixKernel*, MatrixKernel*, void>)p->executionHandle;
+                    for (int i = 0; i < count; i++)
+                    {
+                        op_0(left, result);
+                        left += p->leftStride;
+                        result += p->resultStride;
+                    }
+                    break;
+                case 1:
+                    delegate*<MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void> op_1 = (delegate*<MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void>)p->executionHandle;
+                    for (int i = 0; i < count; i++)
+                    {
+                        op_1(left, result, scalar);
+                        left += p->leftStride;
+                        result += p->resultStride;
+                    }
+                    break;
+                case 2:
+                    delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, void> op_2 = (delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, void>)p->executionHandle;
+                    for (int i = 0; i < count; i++)
+                    {
+                        op_2(left, right, result);
+                        right += p->rightStride;
+                        left += p->leftStride;
+                        result += p->resultStride;
+                    }
+                    break;
+                case 3:
+                    delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void> op_3 = (delegate*<MatrixKernel*, MatrixKernel*, MatrixKernel*, ScalarDataPacket*, void>)p->executionHandle;
+                    for (int i = 0; i < count; i++)
+                    {
+                        op_3(left, right, result, scalar);
+                        right += p->rightStride;
+                        left += p->leftStride;
+                        result += p->resultStride;
+                    }
+                    break;
+                default:
+                    break;
+            }
+        }
+
+        private static unsafe void SchedulerMatrixParallelCycle(object obj)
+        {
+            MatrixParallelWrap.Handle handle = (obj as MatrixParallelWrap.Handle)!;
+            MatrixParallelWrap wrap = handle.Wrap;
+            int inWrapIndex = handle.InWrapIndex;
+            Thread self = Thread.CurrentThread;
+            ref MatrixParallelTaskPackage wt = ref wrap.WrapTask, tt = ref TotalTask;
+            int l1 = wt
+
+            while (true)
+            {
+                self.Priority = ThreadPriority.Highest;
+                while (wrap.TryGetNextWrapTask(inWrapIndex, out IntPtr ptr)) {
+                    ExecuteTaskPackage((MatrixParallelTaskPackage*)ptr);
+                }
+                _ = Interlocked.Decrement(ref RemainThreadCount);
+                self.Priority = _AfterWork;
+                _restartEvent.Wait();
+            }
+        }
+
+        #endregion
     }
 
     public enum SleepMode

StgSharp.Native/sn-src/k_fma.c

@@ -1,5 +1,6 @@
 #include "sn_intrinsic.h"
 
+// ans = left * right(T) + ans
 INTERNAL void SN_DECL kernel_fma_sse(mat_kernel const *restrict left,
                                      mat_kernel const *restrict right, mat_kernel *restrict ans)
 {
@@ -29,8 +30,13 @@ INTERNAL void SN_DECL kernel_fma_sse(mat_kernel const *restrict left,
                         _mm_add_ps(ans3, _mm_mul_ps(source3, broadcast));
                 }
         }
-        _mm_store_ps(&ans->xmm[0], ans0);
+        source0 = _mm_load_ps(&ans->xmm[0]);
-        _mm_store_ps(&ans->xmm[1], ans1);
+        source1 = _mm_load_ps(&ans->xmm[1]);
-        _mm_store_ps(&ans->xmm[2], ans2);
+        source2 = _mm_load_ps(&ans->xmm[2]);
-        _mm_store_ps(&ans->xmm[3], ans3);
+        source3 = _mm_load_ps(&ans->xmm[3]);
+
+        _mm_store_ps(&ans->xmm[0], _mm_add_ps(source0, ans0));
+        _mm_store_ps(&ans->xmm[1], _mm_add_ps(source1, ans1));
+        _mm_store_ps(&ans->xmm[2], _mm_add_ps(source2, ans2));
+        _mm_store_ps(&ans->xmm[3], _mm_add_ps(source3, ans3));
 }