Huge progress on Matrix Parallel, a parallel pool has its basic structure.

StgSharp.Common/Mathematics/Linear/Matrix.Basic.cs

@@ -45,9 +45,18 @@ namespace StgSharp.Mathematics
             if (left.ColumnLength != right.ColumnLength || left.RowLength != right.RowLength) {
                 throw new IndexOutOfRangeException("Matrix dimensions must agree.");
             }
-            Matrix<float> result = GCMatrix<float>.FromDefault(left.ColumnLength, left.RowLength);
+            MatrixParallel.LeadParallel();
-            MatrixParallel.PublishTask();
+            try
-            return result;
+            {
+                Matrix<float> result = GCMatrix<float>.FromDefault(left.ColumnLength, left.RowLength);
+                MatrixParallel.PublishTask();
+                MatrixParallel.LaunchParallel(SleepMode.DeepSleep);
+                return result;
+            }
+            finally
+            {
+                MatrixParallel.EndParallel();
+            }
         }
 
     }

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

@@ -31,16 +31,18 @@
 using StgSharp.HighPerformance;
 
 using System;
+using System.Collections.Generic;
 using System.Drawing;
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Runtime.ExceptionServices;
 using System.Runtime.InteropServices;
 using System.Threading;
 
 namespace StgSharp.Mathematics
 {
     [StructLayout(LayoutKind.Explicit)]
-    public unsafe struct MatrixParallelTaskWrap
+    public unsafe struct MatrixParallelTaskPackage
     {
 
         [FieldOffset(0)] public MatrixKernel* left;
@@ -85,4 +87,43 @@ 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.TaskWrapFactory.cs

@@ -42,7 +42,7 @@ namespace StgSharp.Mathematics
 
         private static readonly SlabAllocator<ScalarDataPacket> _scalarAllocator = SlabAllocator<ScalarDataPacket>.Create(64, SlabBufferLayout.Chunked, true);
 
-        private static readonly SlabAllocator<MatrixParallelTaskWrap> _wrapAllocator = SlabAllocator<MatrixParallelTaskWrap>.Create(64, SlabBufferLayout.Chunked, true);
+        private static readonly SlabAllocator<MatrixParallelTaskPackage> _wrapAllocator = SlabAllocator<MatrixParallelTaskPackage>.Create(64, SlabBufferLayout.Chunked, true);
 
         public static ScalarDataPacket* CreateScalarPacket()
         {
@@ -50,9 +50,9 @@ namespace StgSharp.Mathematics
             return p;
         }
 
-        public static MatrixParallelTaskWrap* FromTaskWrap(MatrixParallelTaskWrap* source, int offset)
+        public static MatrixParallelTaskPackage* FromTaskWrap(MatrixParallelTaskPackage* source, int offset)
         {
-            MatrixParallelTaskWrap* p = (MatrixParallelTaskWrap*)_wrapAllocator.Allocate();
+            MatrixParallelTaskPackage* p = (MatrixParallelTaskPackage*)_wrapAllocator.Allocate();
             *p = *source;
             p->left += offset * p->leftStride;
             p->right += offset * p->rightStride;
@@ -65,7 +65,7 @@ namespace StgSharp.Mathematics
             _scalarAllocator.Free((nuint)source);
         }
 
-        public static void Release(MatrixParallelTaskWrap* source)
+        public static void Release(MatrixParallelTaskPackage* source)
         {
             if (source->scalerPacket != null) {
                 _scalarAllocator.Free((nuint)source->scalerPacket);
@@ -73,5 +73,29 @@ namespace StgSharp.Mathematics
             _wrapAllocator.Free((nuint)source);
         }
 
+        public static void SliceTask(MatrixParallelTaskPackage* source, Span<IntPtr> group)
+        {
+            int count = group.Length;
+            group[0] = (nint)source;
+            int l = source->leftStride, r = source->rightStride, res = source->resultStride, c = source->count / count;
+            source->leftStride *= count;
+            source->rightStride *= count;
+            source->resultStride *= count;
+            int remain = source->count - (c * count);
+            for (int i = 1; i < count; i++)
+            {
+                nint handle = (nint)_wrapAllocator.Allocate();
+                MatrixParallelTaskPackage* p = (MatrixParallelTaskPackage*)handle;
+                *p = *source;
+                p->left += l;
+                p->right += r;
+                p->result += res;
+                group[i] = handle;
+            }
+            for (int i = 0; i < remain; i++) {
+                ((MatrixParallelTaskPackage*)group[^i])->count++;
+            }
+        }
+
     }
 }

StgSharp.Common/Mathematics/Linear/MatrixParallel.cs

@@ -28,6 +28,7 @@
 //     
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
+using StgSharp.HighPerformance;
 using StgSharp.HighPerformance.Memory;
 
 using System;
@@ -35,6 +36,8 @@ using System.Collections.Generic;
 using System.IO.Pipes;
 using System.Linq;
 using System.Reflection.Metadata;
+using System.Runtime.CompilerServices;
+using System.Runtime.ExceptionServices;
 using System.Text;
 using System.Threading;
 using System.Threading.Channels;
@@ -45,40 +48,220 @@ namespace StgSharp.Mathematics
     internal static class MatrixParallel
     {
 
-        private static readonly Channel<IntPtr> _taskChannel = Channel.CreateUnbounded<IntPtr>();
+        // TODO 解决多线程协调问题，barrier性能不好，MRES时机不好控制，可能可以使用双MRES
-        private static int _paraState;
+        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),
+         (0, 3), (3, 1), (2, 2), (1, 3), (4, 1),
+         (3, 2), (2, 3), (1, 4), (0, 5),];
+        private static volatile int RemainThreadCount = 0;
 
-        private static IReadOnlyList<Thread> _threads;
+        private static readonly ManualResetEventSlim _resetEvent = new(false);
+        private static readonly object _lock = new();
+        private static ThreadPriority _AfterWork = ThreadPriority.Lowest;
+
+        public static bool SupportParallel { get; } = Environment.ProcessorCount > 6;
+
+        public static void EndParallel()
+        {
+            Monitor.Exit(_lock);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static (int a, int b) FastDecompose(int c)
+        {
+            if (c < 7) {
+                return (0, 0);
+            }
+            c -= 7;
+            int c_prime = c / 14;
+            int d = c % 14;
+
+            (int a_prime, int b_prime) = LookupTable[d];
+
+            return (a_prime + c_prime, b_prime + c_prime);
+        }
 
         public static void Init()
         {
             int count = Environment.ProcessorCount - 4;
-            count = count < 0 ? 1 : count;
+            count = count == 5 ? 4 : count;
-            List<Thread> _t = new(count);
+            int _3wrap = 0, _4wrap = 0;
-            _threads = _t;
+            _threads = new Thread[count];
-            for (int i = 0; i < count; i++)
+            if (count < 6)
             {
-                _t.Add(new(MatrixParallelCycle)
+                _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
+                        IsBackground = true,
-                    });
+                        Priority = ThreadPriority.Highest,
-                _t[i].Start();
+                        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();
+            RemainThreadCount = _threads.Length;
+            _AfterWork = (ThreadPriority)afterWork;
+            _resetEvent.Reset();
+        }
+
+        public static void LeadParallel()
+        {
+            Monitor.Enter(_lock);
+        }
+
         public static void PublishTask() { }
 
-        public static void Terminate()
+        private static unsafe void MatrixParallelCycle(object obj)
         {
-            _ = Interlocked.Exchange(ref _paraState, 0);
+            MatrixParallelWrap.Handle handle = (obj as MatrixParallelWrap.Handle)!;
-        }
+            MatrixParallelWrap wrap = handle.Wrap;
-
+            int inWrapIndex = handle.InWrapIndex;
-        private static unsafe void MatrixParallelCycle()
+            Thread self = Thread.CurrentThread;
-        {
+            while (true)
-            Thread t_cur = Thread.CurrentThread;
+            {
-            int index = _threads.ToList().IndexOf(t_cur);
+                while (wrap.TryGetNextWrapTask(inWrapIndex, out IntPtr ptr))
-            while (_taskChannel.Reader.TryRead(out IntPtr ptr)) { }
+                {
+                    MatrixParallelTaskPackage* p = (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();
+                self.Priority = ThreadPriority.Highest;
+            }
         }
 
     }
+
+    public enum SleepMode
+    {
+
+        QuickWakeup = ThreadPriority.AboveNormal,
+        NormalSleep = ThreadPriority.BelowNormal,
+        DeepSleep = ThreadPriority.Lowest,
+
+    }
 }

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -73,24 +73,34 @@ typedef uint64_t(SN_DECL *FACTORIALROC)(int n);
 
 INTERNAL void normalize(__m128 *source, __m128 *target);
 
-INTERNAL void SN_DECL kernel_add_sse(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+INTERNAL void SN_DECL kernel_add_sse(mat_kernel const *left, mat_kernel const *right,
-INTERNAL void SN_DECL kernel_add_avx(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+                                     mat_kernel *ans);
-INTERNAL void SN_DECL kernel_add_512(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+INTERNAL void SN_DECL kernel_add_avx(mat_kernel const *left, mat_kernel const *right,
-INTERNAL void SN_DECL kernel_sub_sse(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+                                     mat_kernel *ans);
-INTERNAL void SN_DECL kernel_sub_avx(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+INTERNAL void SN_DECL kernel_add_512(mat_kernel const *left, mat_kernel const *right,
-INTERNAL void SN_DECL kernel_sub_512(mat_kernel const *left, mat_kernel const *right, mat_kernel *ans);
+                                     mat_kernel *ans);
+INTERNAL void SN_DECL kernel_sub_sse(mat_kernel const *left, mat_kernel const *right,
+                                     mat_kernel *ans);
+INTERNAL void SN_DECL kernel_sub_avx(mat_kernel const *left, mat_kernel const *right,
+                                     mat_kernel *ans);
+INTERNAL void SN_DECL kernel_sub_512(mat_kernel const *left, mat_kernel const *right,
+                                     mat_kernel *ans);
 INTERNAL void SN_DECL kernel_transpose_sse(mat_kernel const *source, mat_kernel *target);
 INTERNAL void SN_DECL kernel_transpose_avx(mat_kernel const *source, mat_kernel *target);
 INTERNAL void SN_DECL kernel_transpose_512(mat_kernel const *source, mat_kernel *target);
-INTERNAL void SN_DECL kernel_scalar_mul_sse(mat_kernel const *matrix, float const scalar, mat_kernel *ans);
+INTERNAL void SN_DECL kernel_scalar_mul_sse(mat_kernel const *matrix, float const scalar,
-INTERNAL void SN_DECL kernel_scalar_mul_avx(mat_kernel const *matrix, float const scalar, mat_kernel *ans);
+                                            mat_kernel *ans);
-INTERNAL void SN_DECL kernel_scalar_mul_512(mat_kernel const *matrix, float const scalar, mat_kernel *ans);
+INTERNAL void SN_DECL kernel_scalar_mul_avx(mat_kernel const *matrix, float const scalar,
+                                            mat_kernel *ans);
+INTERNAL void SN_DECL kernel_scalar_mul_512(mat_kernel const *matrix, float const scalar,
+                                            mat_kernel *ans);
 
 #pragma endregion
 
 #pragma region fma
 
-INTERNAL void SN_DECL kernel_fma_sse(mat_kernel *left, mat_kernel *right);
+INTERNAL void SN_DECL kernel_fma_sse(mat_kernel const *restrict left,
+                                     mat_kernel const *restrict right, mat_kernel *restrict ans);
 
 #pragma endregion
 

StgSharp.Native/sn-src/k_fma.c

@@ -1,6 +1,7 @@
 #include "sn_intrinsic.h"
 
-void kernel_fma_sse(mat_kernel *const left, mat_kernel *const right, mat_kernel *ans)
+INTERNAL void SN_DECL kernel_fma_sse(mat_kernel const *restrict left,
+                                     mat_kernel const *restrict right, mat_kernel *restrict ans)
 {
         register __m128 source0 = _mm_load_ps(&left->xmm[0]);
         register __m128 source1 = _mm_load_ps(&left->xmm[1]);