working on matrix mul matrix

StgSharp.Common/HighPerformance/MatrixPanel.cs

@@ -34,82 +34,63 @@ using System.Runtime.InteropServices;
 namespace StgSharp.HighPerformance
 {
     public unsafe struct MatrixPanel
-    {
-
-        private byte Head;
-
-    }
-
-    // Single POD-style panel: ZMM-compatible size. Alignment handled by allocator, not by the struct.
-    public unsafe struct MatrixPanel<T> where T : unmanaged, INumber<T>
     {
 
         public const int AlignmentBytes = 64;
 
-        internal T Head; // first element of the panel buffer
+        private byte Head;
-
-        public ref T this[int column, int row]
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => ref Unsafe.Add(ref Head, row * ElementSideLength + column);
-        }
-
-        public static int ElementSideLength
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => Math.Max(8, AlignmentBytes / sizeof(T));
-        }
-
-        public static int MaxElementCount
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => ElementSideLength * ElementSideLength;
-        }
-
-        public static nuint Size
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => (nuint)(MaxElementCount * sizeof(T));
-        }
-
-        public Span<T> Elements
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get
-            {
-                ref T first = ref Head;
-                unsafe
-                {
-                    void* p = Unsafe.AsPointer(ref first);
-                    return new Span<T>(p, MaxElementCount);
-                }
-            }
-        }
 
         // Allocate aligned memory for one panel; caller must free via Destroy.
-        public static unsafe MatrixPanel<T>* Create()
+        public static unsafe MatrixPanel<T>* Create<T>() where T : unmanaged, INumber<T>
         {
-            return (MatrixPanel<T>*)NativeMemory.AlignedAlloc(Size, AlignmentBytes);
+            return (MatrixPanel<T>*)NativeMemory.AlignedAlloc(Size<T>(), AlignmentBytes);
         }
 
-        public static unsafe void Destroy(MatrixPanel<T>* panel)
+        public static unsafe void Destroy<T>(MatrixPanel<T>* panel) where T : unmanaged, INumber<T>
         {
             NativeMemory.AlignedFree(panel);
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe Span<T> Elements<T>(ref T head) where T : unmanaged, INumber<T>
+        {
+            void* p = Unsafe.AsPointer(ref head);
+            return new Span<T>(p, MaxElementCount<T>());
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static int ElementSideLength<T>() where T : unmanaged, INumber<T>
+        {
+            return Math.Max(8, AlignmentBytes / sizeof(T));
+        }
+
         // View an existing pointer; caller guarantees lifetime/alignment.
-        public static unsafe MatrixPanel<T>* FromPointer(void* ptr)
+        public static unsafe MatrixPanel<T>* FromPointer<T>(void* ptr)
+            where T : unmanaged, INumber<T>
         {
             return (MatrixPanel<T>*)Align64(ptr);
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static int MaxElementCount<T>() where T : unmanaged, INumber<T>
+        {
+            int side = ElementSideLength<T>();
+            return side * side;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static nuint Size<T>() where T : unmanaged, INumber<T>
+        {
+            return (nuint)(MaxElementCount<T>() * sizeof(T));
+        }
+
         /// <summary>
         ///   Number of 4x4 kernels that fit in a panel for the given ISA and element T.
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static int KernelCapacity()
+        internal static int KernelCapacity<T>() where T : unmanaged, INumber<T>
         {
-            int perSide = KernelSideCount();
+            int perSide = KernelSideCount<T>();
             return perSide * perSide;
         }
 
@@ -117,13 +98,13 @@ namespace StgSharp.HighPerformance
         ///   How many 4x4 kernels fit along a single panel dimension for the given ISA.
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static int KernelSideCount()
+        internal static int KernelSideCount<T>() where T : unmanaged, INumber<T>
         {
             return NativeIntrinsic.IntrinsicLevel switch
             {
-                IntrinsicLevel.SSE => KernelSideCountSse(),
+                IntrinsicLevel.SSE => KernelSideCountSse<T>(),
-                IntrinsicLevel.AVX2 => KernelSideCountAvx2(),
+                IntrinsicLevel.AVX2 => KernelSideCountAvx2<T>(),
-                IntrinsicLevel.AVX512 => KernelSideCountAvx512(),
+                IntrinsicLevel.AVX512 => KernelSideCountAvx512<T>(),
                 _ => 0
             };
         }
@@ -132,7 +113,7 @@ namespace StgSharp.HighPerformance
         private static unsafe nuint Align64(void* p) => ((nuint)p + (AlignmentBytes - 1)) & ~(nuint)(AlignmentBytes - 1);
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static int KernelSideCountAvx2()
+        private static int KernelSideCountAvx2<T>() where T : unmanaged, INumber<T>
         {
             const int minQ = 8;
             int lanes = 32 / sizeof(T);
@@ -141,7 +122,7 @@ namespace StgSharp.HighPerformance
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static int KernelSideCountAvx512()
+        private static int KernelSideCountAvx512<T>() where T : unmanaged, INumber<T>
         {
             const int minQ = 8;
             int lanes = 64 / sizeof(T);
@@ -150,7 +131,7 @@ namespace StgSharp.HighPerformance
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        private static int KernelSideCountSse()
+        private static int KernelSideCountSse<T>() where T : unmanaged, INumber<T>
         {
             const int minQ = 4;
             int lanes = 16 / sizeof(T);
@@ -159,4 +140,25 @@ namespace StgSharp.HighPerformance
         }
 
     }
+
+    // Single POD-style panel: ZMM-compatible size. Alignment handled by allocator, not by the struct.
+    public unsafe struct MatrixPanel<T> where T : unmanaged, INumber<T>
+    {
+
+        internal T Head; // first element of the panel buffer 
+
+        public ref T this[int column, int row]
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => ref Unsafe.Add(ref Head, row * MatrixPanel.ElementSideLength<T>() + column);
         }
+
+        public Span<T> Elements
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => MatrixPanel.Elements(ref Head);
+        }
+
+    }
+}
+

StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs

@@ -45,6 +45,7 @@ namespace StgSharp.Internal.Intrinsic
         public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> f32_buffer_sub;             //f32_kernel_sub
         public delegate* unmanaged[Cdecl]<void*, void*, void> f32_buffer_transpose;                //f32_kernel_transpose
         public delegate* unmanaged[Cdecl]<MatrixKernel*, MatrixPanel*, int, int, int, int, void> f32_build_panel;              //f32_build_panel
+        public delegate* unmanaged[Cdecl]<MatrixPanel*, void> f32_clear_panel;
         public delegate* unmanaged[Cdecl]<Vector4*, Vector4*, void> f32_normalize_3;                //f32_normalize_3
         public delegate* unmanaged[Cdecl]<void*, void*, void*, void> f32_panel_fma;               //f32_panel_fma
         public delegate* unmanaged[Cdecl]<void*, void*, int, int, int, int, void> f32_store_panel; //f32_store_panel

StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs

@@ -0,0 +1,54 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixCompute.Additional"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using StgSharp.HighPerformance;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace StgSharp.Mathematics.Numeric
+{
+    public partial class MatrixCompute
+    {
+
+        private static unsafe void ComputeFmaPanelF32(
+                                   MatrixKernel* left,
+                                   MatrixKernel* right,
+                                   MatrixKernel* ans,
+                                   ScalarPacket* scalar)
+        {
+            const int eSize = 32;
+            int colum = scalar->Data<int>(0);
+            int row = scalar->Data<int>(1);
+            int k = 
+
+        }
+
+    }
+}

StgSharp.Mathematics/Numeric/MatrixCompute.Basic.cs

@@ -39,13 +39,10 @@ using System.Threading.Tasks;
 
 namespace StgSharp.Mathematics.Numeric
 {
-    public static class MatrixCompute
+    public static partial class MatrixCompute
     {
 
-        public static unsafe Matrix<float> Add(
+        public static unsafe void Add(Matrix<float> left, Matrix<float> right, Matrix<float> ans)
-                                           Matrix<float> left,
-                                           Matrix<float> right,
-                                           Matrix<float> ans)
         {
             if (left.ColumnLength != right.ColumnLength ||
                 left.RowLength != right.RowLength ||
@@ -61,7 +58,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* ansPtr = ans.Buffer;
                 NativeIntrinsic.Intrinsic.f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
 
-                return ans;
+                return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
@@ -72,22 +69,22 @@ namespace StgSharp.Mathematics.Numeric
                 NativeIntrinsic.Intrinsic.f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
-                return ans;
+                return;
             }
             MatrixParallelTaskPackage* package = MatrixParallelFactory.CreateBaseTask<float>();
             package->Left = (nint)left.Buffer;
             package->Right = (nint)right.Buffer;
             package->Result = (nint)ans.Buffer;
             package->ElementSize = sizeof(float);
-            package->PrimCount = left.KernelColumnLength;
+            package->PrimCount = ans.KernelColumnLength;
-            package->SecCount = left.KernelRowLength;
+            package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixOperationStyle.BUFFER_OP, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
             package->ComputeHandle = (IntPtr)NativeIntrinsic.Intrinsic.f32_buffer_add;
             MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
             taskGroup.Dispose();
-            return ans;
+            return;
         }
 
         public static unsafe void Fill(Matrix<float> ans, float value)
@@ -120,10 +117,7 @@ namespace StgSharp.Mathematics.Numeric
             taskGroup.Dispose();
         }
 
-        public static unsafe Matrix<float> Mul(
+        public static unsafe void Mul(Matrix<float> left, Matrix<float> right, Matrix<float> ans)
-                                           Matrix<float> left,
-                                           Matrix<float> right,
-                                           Matrix<float> ans)
         {
             // For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
             // ColumnLength = number of columns (N), RowLength = number of rows (M)
@@ -134,30 +128,68 @@ namespace StgSharp.Mathematics.Numeric
                 throw new IndexOutOfRangeException("Matrix dimensions must agree.");
             }
             long count = (long)left.KernelColumnLength * left.KernelRowLength;
-            if (count <= 256)
+            if (count <= 1024)
             {
-                MatrixPanel<float>* leftPanel = MatrixPanel<float>.Create();
+                MatrixParallelThread[] pool = null;
-                MatrixPanel<float>* rightPanel = MatrixPanel<float>.Create();
+                if (count > 256) {
-                MatrixPanel<float>* ansPanel = MatrixPanel<float>.Create();
+                    pool = MatrixParallel.LeadParallel(1);
-                int size = MatrixPanel<float>.KernelSideCount();
+                }
+                MatrixPanel<float>* leftPanel = MatrixPanel.Create<float>();
+                MatrixPanel<float>* rightPanel = MatrixPanel.Create<float>();
+                MatrixPanel<float>* ansPanel = MatrixPanel.Create<float>();
+                int size = MatrixPanel.KernelSideCount<float>();
                 for (int i = 0; i < left.RowLength; i += size)
                 {
                     for (int j = 0; j < right.ColumnLength; j += size)
                     {
-                        for (int k = 0; k < left.ColumnLength; k += size) {
+                        NativeIntrinsic.Intrinsic.f32_clear_panel((MatrixPanel*)ansPanel);
+                        for (int k = 0; k < left.ColumnLength; k += size)
+                        {
                             NativeIntrinsic.Intrinsic
                                            .f32_build_panel((MatrixKernel*)left.Buffer,
                                                             (MatrixPanel*)leftPanel,
                                                             left.ColumnLength, left.RowLength, i,
                                                             k);
+                            NativeIntrinsic.Intrinsic
+                                           .f32_build_panel((MatrixKernel*)right.Buffer,
+                                                            (MatrixPanel*)rightPanel,
+                                                            right.ColumnLength, right.RowLength, k,
+                                                            j);
+                            NativeIntrinsic.Intrinsic
+                                           .f32_panel_fma((MatrixPanel*)leftPanel,
+                                                            (MatrixPanel*)rightPanel,
+                                                            (MatrixPanel*)ansPanel);
+                        }
+                        NativeIntrinsic.Intrinsic
+                                       .f32_store_panel((MatrixKernel*)ans.Buffer,
+                                                        (MatrixPanel*)ansPanel, ans.ColumnLength,
+                                                        ans.RowLength, i, j);
                     }
                 }
+                MatrixPanel.Destroy(leftPanel);
+                MatrixPanel.Destroy(rightPanel);
+                MatrixPanel.Destroy(ansPanel);
+                if (count > 256) {
+                    MatrixParallel.ReturnParallelResources(ref pool!);
                 }
+                return;
             }
-            throw new NotImplementedException();
+            MatrixParallelTaskPackage* package = MatrixParallelFactory.CreateBaseTask<float>();
+            package->Right = (nint)left.Buffer;
+            package->Result = (nint)ans.Buffer;
+            package->ElementSize = sizeof(float);
+            package->PrimCount = left.KernelColumnLength;
+            package->SecCount = left.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixOperationStyle.BUFFER_OP, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
+            package->ComputeHandle = (IntPtr)NativeIntrinsic.Intrinsic.f32_buffer_scalar_mul;
+            MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
+
+            MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
         }
 
-        public static unsafe Matrix<float> Mul(Matrix<float> left, float right, Matrix<float> ans)
+        public static unsafe void Mul(Matrix<float> left, float right, Matrix<float> ans)
         {
             if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
                 throw new IndexOutOfRangeException("Matrix dimensions must agree.");
@@ -171,7 +203,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* ansPtr = ans.Buffer;
                 NativeIntrinsic.Intrinsic.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
 
-                return ans;
+                return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
@@ -181,7 +213,7 @@ namespace StgSharp.Mathematics.Numeric
                 NativeIntrinsic.Intrinsic.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
-                return ans;
+                return;
             }
             MatrixParallelTaskPackage* package = MatrixParallelFactory.CreateBaseTask<float>();
             package->Right = (nint)left.Buffer;
@@ -197,13 +229,10 @@ namespace StgSharp.Mathematics.Numeric
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
             taskGroup.Dispose();
             MatrixParallelFactory.Release(scalar);
-            return ans;
+            return;
         }
 
-        public static unsafe Matrix<float> Sub(
+        public static unsafe void Sub(Matrix<float> left, Matrix<float> right, Matrix<float> ans)
-                                           Matrix<float> left,
-                                           Matrix<float> right,
-                                           Matrix<float> ans)
         {
             if (left.ColumnLength != right.ColumnLength ||
                 left.RowLength != right.RowLength ||
@@ -219,7 +248,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* ansPtr = ans.Buffer;
                 NativeIntrinsic.Intrinsic.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
 
-                return ans;
+                return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
@@ -230,7 +259,7 @@ namespace StgSharp.Mathematics.Numeric
                 NativeIntrinsic.Intrinsic.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
-                return ans;
+                return;
             }
             MatrixParallelTaskPackage* package = MatrixParallelFactory.CreateBaseTask<float>();
             package->Left = (nint)left.Buffer;
@@ -245,7 +274,7 @@ namespace StgSharp.Mathematics.Numeric
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
             taskGroup.Dispose();
-            return ans;
+            return;
         }
 
     }

StgSharp.Mathematics/Numeric/MatrixCompute.cs

@@ -0,0 +1,37 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixCompute"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace StgSharp.Mathematics.Numeric
+{
+    public static partial class MatrixCompute { }
+}

StgSharp.Mathematics/Numeric/MatrixComputeModel.cs

@@ -1,4 +1,4 @@
-//-----------------------------------------------------------------------
+//-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
 // file="MatrixComputeModel"
 // Project: StgSharp
@@ -76,12 +76,13 @@ namespace StgSharp.Mathematics.Numeric
             long offset = Unsafe.As<int, long>(ref p->PrimTileOffset) * size;
             ((delegate* unmanaged[Cdecl]< void*, void*, ScalarPacket*, long, void>)p->ComputeHandle)((void*)(p->Right + offset),
                                                                                                      (void*)(p->Result + offset),
-                                                                                                     p->Scalar, count);
+                                                                                                     p->Scalar,
+                                                                                                     count);
             MatrixParallelFactory.Release(p);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static unsafe void Compute_Binary(MatrixParallelTaskPackage* p)
+        internal static unsafe void KernelComputeBinary(MatrixParallelTaskPackage* p)
         {
             int eSize = p->ElementSize;
 
@@ -133,7 +134,7 @@ namespace StgSharp.Mathematics.Numeric
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static unsafe void Compute_BinaryScalar(MatrixParallelTaskPackage* p)
+        internal static unsafe void KernelComputeBinaryScalar(MatrixParallelTaskPackage* p)
         {
             int eSize = p->ElementSize;
 
@@ -185,7 +186,7 @@ namespace StgSharp.Mathematics.Numeric
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static unsafe void Compute_Unary(MatrixParallelTaskPackage* p)
+        internal static unsafe void KernelComputeUnary(MatrixParallelTaskPackage* p)
         {
             int eSize = p->ElementSize;
 
@@ -230,7 +231,7 @@ namespace StgSharp.Mathematics.Numeric
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal static unsafe void Compute_UnaryScalar(MatrixParallelTaskPackage* p)
+        internal static unsafe void KernelComputeUnaryScalar(MatrixParallelTaskPackage* p)
         {
             int eSize = p->ElementSize;
 
@@ -274,5 +275,7 @@ namespace StgSharp.Mathematics.Numeric
             MatrixParallelFactory.Release(p);
         }
 
+        internal static unsafe void PanelComputeBinary(MatrixParallelTaskPackage* p) { }
+
     }
 }

StgSharp.Mathematics/Numeric/MatrixParallel.cs

@@ -1,4 +1,4 @@
-//-----------------------------------------------------------------------
+//-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
 // file="MatrixParallel"
 // Project: StgSharp
@@ -27,15 +27,14 @@
 // -----------------------------------------------------------------------
 using StgSharp.Collections;
 using System;
-using System.Collections.Concurrent;
 using System.Collections.Generic;
-using System.Diagnostics.CodeAnalysis;
+using System.Linq;
-using System.Runtime.CompilerServices;
+using System.Text;
-using System.Threading;
+using System.Threading.Tasks;
 
 namespace StgSharp.Mathematics.Numeric
 {
-    public static partial class MatrixParallel
+    public static class MatrixParallel
     {
 
         private static bool _meaningful;

StgSharp.Mathematics/Numeric/MatrixParallelThread.cs

@@ -1,4 +1,4 @@
-//-----------------------------------------------------------------------
+//-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
 // file="MatrixParallelThread"
 // Project: StgSharp
@@ -114,16 +114,16 @@ namespace StgSharp.Mathematics.Numeric
                         switch (p->ComputeMode.ParameterStyle)
                         {
                             case MatrixOperationParam.RIGHT_ANS_PARAM:
-                                MatrixComputeModel.Compute_Unary(p);
+                                MatrixComputeModel.KernelComputeUnary(p);
                                 break;
                             case MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM:
-                                MatrixComputeModel.Compute_UnaryScalar(p);
+                                MatrixComputeModel.KernelComputeUnaryScalar(p);
                                 break;
                             case MatrixOperationParam.LEFT_RIGHT_ANS_PARAM:
-                                MatrixComputeModel.Compute_Binary(p);
+                                MatrixComputeModel.KernelComputeBinary(p);
                                 break;
                             case MatrixOperationParam.LEFT_RIGHT_ANS_SCALAR_PARAM:
-                                MatrixComputeModel.Compute_BinaryScalar(p);
+                                MatrixComputeModel.KernelComputeBinaryScalar(p);
                                 break;
                             default:
                                 break;
@@ -400,7 +400,9 @@ namespace StgSharp.Mathematics.Numeric
         private static extern IntPtr GetCurrentThread();
 
         [DllImport("kernel32.dll")]
-        private static extern nuint SetThreadAffinityMask(IntPtr hThread, nuint dwThreadAffinityMask);
+        private static extern nuint SetThreadAffinityMask(
+                                    IntPtr hThread,
+                                    nuint dwThreadAffinityMask);
 
         [LibraryImport("kernel32.dll")]
         [return: MarshalAs(UnmanagedType.Bool)]

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -140,7 +140,9 @@ typedef void(SN_DECL *PNL_PROC_LRAS)(mat_panel *left, mat_panel *right, mat_pane
 typedef void(SN_DECL *PNL_PROC_RAS)(mat_panel *right, mat_panel *ans, scalar_pack *scalar);
 
 #define PNL_PROC(T_type, T_arch, T_op_name) T_type##_pnl_##T_op_name##_##T_arch
-
+#define DECLARE_PNL_PROC_ANS(T_type, T_arch, T_op_name) \
+        INTERNAL void SN_DECL PNL_PROC(T_type, T_arch,  \
+                                       T_op_name)(MAT_PANEL_##T_arch##(T_type) *restrict ans)
 #define DECLARE_PNL_PROC_LEFT_RIGHT_ANS(T_type, T_arch, T_op_name)                           \
         INTERNAL void SN_DECL PNL_PROC(T_type, T_arch,                                       \
                                        T_op_name)(MAT_PANEL_##T_arch##(T_type) const *left,  \
@@ -223,8 +225,14 @@ DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill);
 
 #pragma region fma
 
+DECLARE_PNL_PROC_ANS(float, sse, clear);
+DECLARE_PNL_PROC_ANS(float, avx, clear);
+DECLARE_PNL_PROC_ANS(float, 512, clear);
+
 DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, sse, fma);
 DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx, fma);
+DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx_fma, fma);
+DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, 512, fma);
 
 #pragma endregion
 
@@ -266,6 +274,7 @@ typedef struct sn_intrinsic {
         BUF_PROC_LRA f32_buffer_sub;
         KER_PROC_RA f32_buffer_transpose;
         PANEL_PROC f32_build_panel;
+        PANEL_PROC f32_clear_panel;
         VECTORNORMALIZEPROC f32_normalize_3;
         PNL_PROC_LRA f32_panel_fma;
         PANEL_PROC f32_store_panel;

StgSharp.Native/sn-head/sn_matkernel.h

@@ -117,7 +117,11 @@ extern __declspec(align(64)) sn_zero_kernel_u64 SN_ZERO_KERNEL;
                                                          int col_index, int row_index)
 
 /* Declarations for AVX float variant */
+DECLARE_BUILD_PANEL(sse, float);
 DECLARE_BUILD_PANEL(avx, float);
+DECLARE_BUILD_PANEL(512, float);
+DECLARE_STORE_PANEL(sse, float);
 DECLARE_STORE_PANEL(avx, float);
+DECLARE_STORE_PANEL(512, float);
 
 #endif /* SN_MATKERNEL */

StgSharp.Native/sn-src/k_basicop.c

@@ -1,5 +1,32 @@
 #include "sn_intrinsic.h"
 
+DECLARE_PNL_PROC_ANS(float, sse, clear)
+{
+        register __m128 zero = _mm_setzero_ps();
+        PANEL_COL_VEC((*ans), 0, 0) = zero;
+        PANEL_COL_VEC((*ans), 1, 0) = zero;
+        PANEL_COL_VEC((*ans), 2, 0) = zero;
+        PANEL_COL_VEC((*ans), 3, 0) = zero;
+}
+
+DECLARE_PNL_PROC_ANS(float, avx, clear)
+{
+        register __m256 zero = _mm256_setzero_ps();
+        PANEL_COL_VEC((*ans), 0, 0) = zero;
+        PANEL_COL_VEC((*ans), 1, 0) = zero;
+        PANEL_COL_VEC((*ans), 2, 0) = zero;
+        PANEL_COL_VEC((*ans), 3, 0) = zero;
+}
+
+DECLARE_PNL_PROC_ANS(float, 512, clear)
+{
+        register __m512 zero = _mm512_setzero_ps();
+        PANEL_COL_VEC((*ans), 0, 0) = zero;
+        PANEL_COL_VEC((*ans), 1, 0) = zero;
+        PANEL_COL_VEC((*ans), 2, 0) = zero;
+        PANEL_COL_VEC((*ans), 3, 0) = zero;
+}
+
 // Matrix transpose implementations
 DECLARE_KER_PROC_RIGHT_ANS(float, sse, transpose)
 {

StgSharp.Native/sn-src/k_panel.c

@@ -41,6 +41,51 @@ DECLARE_BUILD_PANEL(avx, float)
         }
 }
 
+DECLARE_BUILD_PANEL(sse, float)
+{
+        /* Single kernel -> 4 columns, each is one __m128 column vector. */
+        const MAT_KERNEL(float) *src = matrix + (col_index * row_length + row_index);
+
+        for (int pc = 0; pc < 4; ++pc) {
+                PANEL_COL_VEC((*panel), pc, 0) = src->xmm[pc];
+        }
+}
+
+DECLARE_BUILD_PANEL(512, float)
+{
+        const MAT_KERNEL(float) *const zero = (const MAT_KERNEL(float) *)&SN_ZERO_KERNEL;
+        const MAT_KERNEL(float) *k[4][4];
+
+        for (int kc = 0; kc < 4; ++kc) {
+                for (int kr = 0; kr < 4; ++kr) {
+                        const int col = col_index + kc;
+                        const int row = row_index + kr;
+                        k[kc][kr] = (col < col_length && row < row_length) ?
+                                            (matrix + (col * row_length + row)) :
+                                            zero;
+                }
+        }
+
+        /* Each panel column corresponds to a kernel column (0..3) within a block col (0..3). */
+        for (int pc = 0; pc < 16; ++pc) {
+                const int blk_col = pc >> 2;     /* which kernel column block (0..3) */
+                const int kc_col = pc & 3;       /* column inside a 4x4 kernel */
+
+                const __m128 r0 = k[blk_col][0]->xmm[kc_col];
+                const __m128 r1 = k[blk_col][1]->xmm[kc_col];
+                const __m128 r2 = k[blk_col][2]->xmm[kc_col];
+                const __m128 r3 = k[blk_col][3]->xmm[kc_col];
+
+                const __m256 top = _mm256_set_m128(r1, r0);
+                const __m256 bot = _mm256_set_m128(r3, r2);
+
+                __m512 col = _mm512_castps256_ps512(top);
+                col = _mm512_insertf32x8(col, bot, 1);
+
+                PANEL_COL_VEC((*panel), pc, 0) = col;
+        }
+}
+
 DECLARE_STORE_PANEL(avx, float)
 {
         const int Q = 8; /* fixed for AVX float panel */
@@ -96,3 +141,49 @@ DECLARE_STORE_PANEL(avx, float)
                 k_lb->xmm[kc_col] = bot128;
         }
 }
+
+DECLARE_STORE_PANEL(sse, float)
+{
+        MAT_KERNEL(float) *dst = matrix + (col_index * row_length + row_index);
+
+        for (int pc = 0; pc < 4; ++pc) {
+                const __m128 col = PANEL_COL_VEC((*panel), pc, 0);
+                dst->xmm[pc] = col;
+        }
+}
+
+DECLARE_STORE_PANEL(512, float)
+{
+        MAT_KERNEL(float) zero;
+        MAT_KERNEL(float) *k[4][4];
+
+        for (int kc = 0; kc < 4; ++kc) {
+                for (int kr = 0; kr < 4; ++kr) {
+                        const int col = col_index + kc;
+                        const int row = row_index + kr;
+                        k[kc][kr] = (col < col_length && row < row_length) ?
+                                            (matrix + (col * row_length + row)) :
+                                            &zero;
+                }
+        }
+
+        for (int pc = 0; pc < 16; ++pc) {
+                const int blk_col = pc >> 2;     /* which kernel column block (0..3) */
+                const int kc_col = pc & 3;       /* column inside a 4x4 kernel */
+
+                const __m512 col = PANEL_COL_VEC((*panel), pc, 0);
+
+                const __m256 top = _mm512_castps512_ps256(col);
+                const __m256 bot = _mm512_extractf32x8_ps(col, 1);
+
+                const __m128 r0 = _mm256_castps256_ps128(top);
+                const __m128 r1 = _mm256_extractf128_ps(top, 1);
+                const __m128 r2 = _mm256_castps256_ps128(bot);
+                const __m128 r3 = _mm256_extractf128_ps(bot, 1);
+
+                k[blk_col][0]->xmm[kc_col] = r0;
+                k[blk_col][1]->xmm[kc_col] = r1;
+                k[blk_col][2]->xmm[kc_col] = r2;
+                k[blk_col][3]->xmm[kc_col] = r3;
+        }
+}

StgSharp.Native/sn-src/load_intrinsic.c

@@ -65,8 +65,10 @@ void load_intrinsic_sse(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
+        context->f32_clear_panel = PNL_PROC(float, sse, clear);
+        context->f32_build_panel = BUILD_PANEL(sse, float);
         context->f32_buffer_add = BUF_PROC(float, sse, add);
-        context->f32_panel_fma = KER_PROC(float, sse, fma);
+        context->f32_panel_fma = PNL_PROC(float, sse, fma);
         context->f32_buffer_fill = BUF_PROC(float, sse, fill);
         context->variant = BUF_PROC(float, sse, scalar_mul);
         context->f32_buffer_sub = BUF_PROC(float, sse, sub);
@@ -80,8 +82,10 @@ void load_intrinsic_avx(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
+        context->f32_clear_panel = PNL_PROC(float, avx, clear);
+        context->f32_build_panel = BUILD_PANEL(avx, float);
         context->f32_buffer_add = BUF_PROC(float, avx, add);
-        context->f32_panel_fma = KER_PROC(float, sse, fma);
+        context->f32_panel_fma = PNL_PROC(float, avx, fma);
         context->f32_buffer_fill = BUF_PROC(float, avx, fill);
         context->variant = BUF_PROC(float, avx, scalar_mul);
         context->f32_buffer_sub = BUF_PROC(float, avx, sub);
@@ -95,8 +99,10 @@ void load_intrinsic_512(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
+        context->f32_clear_panel = PNL_PROC(float, sse, clear);
+        context->f32_build_panel = BUILD_PANEL(512, float);
         context->f32_buffer_add = BUF_PROC(float, 512, add);
-        context->f32_panel_fma = KER_PROC(float, sse, fma);
+        context->f32_panel_fma = PNL_PROC(float, 512, fma);
         context->f32_buffer_fill = BUF_PROC(float, 512, fill);
         context->variant = BUF_PROC(float, 512, scalar_mul);
         context->f32_buffer_sub = BUF_PROC(float, 512, sub);