Better FMA

StgSharp.Benchmark/HlsfAllocatorBench.cs

@@ -2,8 +2,8 @@
 // -----------------------------------------------------------------------
 // file="HlsfAllocatorBench"
 // Project: StgSharp
-// AuthorGroup: Nitload Space
+// AuthorGroup: Nitload
-// Copyright (c) Nitload Space. All rights reserved.
+// 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
@@ -27,6 +27,8 @@
 // -----------------------------------------------------------------------
 using BenchmarkDotNet.Attributes;
 using StgSharp.Mathematics.Memory;
+using System;
+using System.Collections.Generic;
 using System.Runtime.InteropServices;
 
 namespace StgSharp.Benchmark
@@ -40,7 +42,15 @@ namespace StgSharp.Benchmark
 
         private const int AllocationCount = 2048;
 
-        private static readonly int[] AllocSizes = [64, 96, 150, 200, 400, 1024, 4096];
+        private static readonly int[] AllocSizes = [
+            64,
+            96,
+            150,
+            200,
+            400,
+            1024,
+            4096
+        ];
         private HybridLayerSegregatedFitAllocator _hlsfAllocator = null!;
         private readonly List<HybridLayerSegregatedFitAllocationHandle> _hlsfLocalList = new();
         private readonly List<nint> _libcLocalList = new();
@@ -96,7 +106,9 @@ namespace StgSharp.Benchmark
             _hlsfAllocator = new HybridLayerSegregatedFitAllocator(GetArenaBytes(), align:(int)Alignment);
         }
 
-        private static uint AlignSize(uint size)
+        private static uint AlignSize(
+                            uint size
+        )
         {
             uint mask = Alignment - 1;
             return (size + mask) & ~mask;
@@ -109,7 +121,9 @@ namespace StgSharp.Benchmark
             int mb = defaultMb;
 
             string? env = Environment.GetEnvironmentVariable("HLSF_BENCH_ARENA_MB");
-            if (!string.IsNullOrWhiteSpace(env) && int.TryParse(env, out int parsed) && parsed >= minMb) {
+            if (!string.IsNullOrWhiteSpace(env) &&
+                int.TryParse(env, out int parsed) &&
+                parsed >= minMb) {
                 mb = parsed;
             }
 
@@ -118,7 +132,10 @@ namespace StgSharp.Benchmark
 
         private uint GetNextRandomSize() => (uint)AllocSizes[_random.Next(AllocSizes.Length)];
 
-        private void RunHlsf(HybridLayerSegregatedFitAllocator allocator, FreePolicy policy)
+        private void RunHlsf(
+                     HybridLayerSegregatedFitAllocator allocator,
+                     FreePolicy policy
+        )
         {
             _hlsfLocalList.Clear();
 

StgSharp.Mathematics/Geometries/PlainGeometry/Define/PlainGeometry.cs

@@ -25,7 +25,7 @@
 //     
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
-using StgSharp.Graphics;
+using StgSharp;
 using StgSharp.Mathematics;
 using StgSharp.Mathematics.Graphics;
 
@@ -43,12 +43,16 @@ namespace StgSharp.Geometries
         protected Vec4[] vertexMat;
         internal CoordinationBase coordinate;
 
-        internal PlainGeometry(CoordinationBase coordination)
+        internal PlainGeometry(
+                 CoordinationBase coordination
+        )
         {
             coordinate = coordination;
         }
 
-        public Point this[int index]
+        public Point this[
+                     int index
+        ]
         {
             get => new Point(vertexMat[index]);
             set => vertexMat[index].vec = value.coordVec;

StgSharp.Mathematics/Internal/MatrixPanel.cs

@@ -93,7 +93,7 @@ namespace StgSharp.Mathematics
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         internal static int KernelSideCount<T>() where T: unmanaged, INumber<T>
         {
-            SIMDID id = MathematicsModule.IntrinsicMask;
+            SIMDID id = NumericalModule.IntrinsicMask;
             switch (id.MaskByte[0] & 0b_00001111)
             {
                 case 1:
@@ -148,7 +148,7 @@ namespace StgSharp.Mathematics
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         private static int PanelSideLength<T>() where T: unmanaged, INumber<T>
         {
-            SIMDID id = MathematicsModule.IntrinsicMask;
+            SIMDID id = NumericalModule.IntrinsicMask;
             switch (id.MaskByte[0] & 0b_00001111)
             {
                 case 1:

StgSharp.Mathematics/Internal/MatrixParallelTask.cs

@@ -39,9 +39,9 @@ namespace StgSharp.Mathematics.Internal
         [FieldOffset(0)]public MatrixKernel* Mat1;    // 8
         [FieldOffset(8)] public MatrixKernel* Mat2;    // 8
         [FieldOffset(16)] public MatrixKernel* Mat3;    // 16
-        [FieldOffset(32)] public int Common_i;
+        [FieldOffset(32)] public int CommonX;
-        [FieldOffset(36)] public int Common_j;
         [FieldOffset(40)] public int Common_k;
+        [FieldOffset(36)] public int CommonY;
         [FieldOffset(44)] public int i;
         [FieldOffset(48)] public int j;
         [FieldOffset(52)] public int k;

StgSharp.Mathematics/Internal/UnsafeCompute.cs

@@ -35,19 +35,19 @@ namespace StgSharp.Mathematics.Internal
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static unsafe int CityHashSimplify(char* span, int length)
         {
-            return MathematicsModule.GlobalContext.city_hash_simplify(span, length);
+            return NumericalModule.GlobalContext.city_hash_simplify(span, length);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static unsafe int CityHashSimplify(char* span, int begin, int length)
         {
-            return MathematicsModule.GlobalContext.city_hash_simplify(span + begin, length);
+            return NumericalModule.GlobalContext.city_hash_simplify(span + begin, length);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static unsafe int IndexOfCharPair(char* source, int pairValue, int length)
         {
-            return MathematicsModule.GlobalContext.index_pair(source, pairValue, length);
+            return NumericalModule.GlobalContext.index_pair(source, pairValue, length);
         }
 
     }

StgSharp.Mathematics/Moduel.cs

@@ -1,80 +0,0 @@
-//-----------------------------------------------------------------------
-// -----------------------------------------------------------------------
-// file="Moduel"
-// Project: StgSharp
-// AuthorGroup: Nitload Space
-// Copyright (c) Nitload Space. 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.Mathematics.Internal;
-using System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Runtime.CompilerServices;
-using System.Text;
-using System.Threading.Tasks;
-
-namespace StgSharp.Mathematics
-{
-    public class MathematicsModule : IStaticModule
-    {
-
-        public static unsafe int DefaultSIMDAlignment
-        {
-            get
-            {
-                SIMDID id = GlobalIntrinsicMask;
-                switch (id.MaskByte[0] & 0b_00001111)
-                {
-                    case 1:
-                        return id.MaskByte[1] switch
-                        {
-                            1 => 16,
-                            2 => 32,
-                            3 => 64,
-                            4 => 64,
-                            5 => 64,
-                            _ => 8,
-                        };
-                    default:
-                        return 8;
-                }
-            }
-        }
-
-        public string ModuleName => "Mathematics";
-
-        internal static IntrinsicContext GlobalContext
-        {
-            [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get;
-        } = new();
-
-        internal static SIMDID GlobalIntrinsicMask { get; set; }
-
-        public void InitializeModule()
-        {
-            throw new NotImplementedException();
-        }
-
-    }
-}

StgSharp.Mathematics/Numeric/Graphics/Matrix44.cs

@@ -189,7 +189,7 @@ namespace StgSharp.Mathematics.Graphics
         {
             const int f32 = 0;
             fixed (MatrixKernel<float>* source = &this.kernel, target = &transpose.kernel) {
-                MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Transpose,
+                NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Transpose,
                                                                null, (MatrixKernel*)source,
                                                                (MatrixKernel*)target, 0);
             }
@@ -203,7 +203,7 @@ namespace StgSharp.Mathematics.Graphics
         {
             Matrix44 ret = new();
             int f32 = MatrixElementType.F32.IntrinsicNode;
-            MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Sub,
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Sub,
                                                            (MatrixKernel*)&left,
                                                            (MatrixKernel*)&right,
                                                            (MatrixKernel*)&ret, 1);
@@ -237,7 +237,7 @@ namespace StgSharp.Mathematics.Graphics
         {
             Matrix44 ret = new();
             int f32 = MatrixElementType.F32.IntrinsicNode;
-            MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Fma,
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Fma,
                                                            (MatrixKernel*)&left,
                                                            (MatrixKernel*)&right,
                                                            (MatrixKernel*)&ret, 0);
@@ -252,7 +252,7 @@ namespace StgSharp.Mathematics.Graphics
         {
             Matrix44 ret = new();
             int f32 = MatrixElementType.F32.IntrinsicNode;
-            MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Add,
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Add,
                                                            (MatrixKernel*)&left,
                                                            (MatrixKernel*)&right,
                                                            (MatrixKernel*)&ret, 1);

StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs

@@ -77,7 +77,7 @@ namespace StgSharp.Mathematics.Numeric
                 ulong* pivots = stackalloc ulong[4];
                 for (int kBlock = 0; kBlock < n; kBlock += 4)
                 {
-                    MathematicsModule.GlobalContext
+                    NumericalModule.GlobalContext
                                    .mat_sk[(int)MatrixElementType.F32].pivot(matrix.Buffer,
                                                                             (ulong)(kBlock >> 2),
                                                                             (ulong)(matrix.KernelRowLength - (kBlock >> 2)),

StgSharp.Mathematics/Numeric/MatrixCompute.Float.cs

@@ -49,7 +49,7 @@ namespace StgSharp.Mathematics.Numeric
                                   Matrix<float> ans
         )
         {
-            int f32 = (int)MatrixElementType.F32;
+            MatrixElementType f32 = MatrixElementType.F32;
             if (left.ColumnLength != right.ColumnLength ||
                 left.RowLength != right.RowLength ||
                 left.ColumnLength != ans.ColumnLength ||
@@ -62,7 +62,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
                                                                      count);
 
                 return;
@@ -73,21 +73,21 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
                                                                      count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
             }
             MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
-            package->Left = (nint)left.Buffer;
+            package->Mat1 = (MatrixKernel*)left.Buffer;
-            package->Right = (nint)right.Buffer;
+            package->Mat2 = (MatrixKernel*)right.Buffer;
-            package->Answer = (nint)ans.Buffer;
+            package->Mat3 = (MatrixKernel*)ans.Buffer;
             package->ElementSize = sizeof(float);
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f32].buffer_add;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_add;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -100,18 +100,18 @@ namespace StgSharp.Mathematics.Numeric
                                   float value
         )
         {
-            int f32 = (int)MatrixElementType.F32;
+            MatrixElementType f32 = MatrixElementType.F32;
             long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
             ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
             scalar->Data<float>(0) = value;
             if (count <= 256)
             {
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
                 return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
             }
@@ -121,7 +121,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f32].buffer_fill;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_fill;
             package->Scalar = scalar;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -135,7 +135,7 @@ namespace StgSharp.Mathematics.Numeric
                                   Matrix<float> ans
         )
         {
-            int f32 = (int)MatrixElementType.F32;
+            MatrixElementType f32 = MatrixElementType.F32;
 
             // For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
             // ColumnLength = number of columns (N), RowLength = number of rows (M)
@@ -168,7 +168,7 @@ namespace StgSharp.Mathematics.Numeric
                         enumeration->Member<int>(0) = i;
                         enumeration->Member<int>(1) = j;
 
-                        MathematicsModule.GlobalContext
+                        NumericalModule.GlobalContext
                                            .mat_sk[f32].kernel_tile_fma((MatrixKernel*)left.Buffer,
                                                                         (MatrixKernel*)right.Buffer,
                                                                         (MatrixKernel*)ans.Buffer,
@@ -205,7 +205,7 @@ namespace StgSharp.Mathematics.Numeric
                                   Matrix<float> ans
         )
         {
-            int f32 = (int)MatrixElementType.F32;
+            MatrixElementType f32 = MatrixElementType.F32;
             if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
                 throw new IndexOutOfRangeException("Matrix dimensions must agree.");
             }
@@ -220,7 +220,7 @@ namespace StgSharp.Mathematics.Numeric
             {
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
                                                                             count);
 
                 return;
@@ -230,7 +230,7 @@ namespace StgSharp.Mathematics.Numeric
 
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
                                                                             count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
@@ -244,7 +244,7 @@ namespace StgSharp.Mathematics.Numeric
             package->SecCount = left.KernelRowLength;
             package->Scalar = scalar;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f32].buffer_scalar_mul;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -259,7 +259,7 @@ namespace StgSharp.Mathematics.Numeric
                                   Matrix<float> ans
         )
         {
-            int f32 = (int)MatrixElementType.F32;
+            MatrixElementType f32 = MatrixElementType.F32;
             if (left.ColumnLength != right.ColumnLength ||
                 left.RowLength != right.RowLength ||
                 left.ColumnLength != ans.ColumnLength ||
@@ -277,7 +277,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
                                                                      count);
 
                 return;
@@ -288,7 +288,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
                                                                      count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
@@ -302,7 +302,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = left.KernelColumnLength;
             package->SecCount = left.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f32].buffer_sub;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_sub;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);

StgSharp.Mathematics/Numeric/MatrixCompute.double.cs

@@ -58,7 +58,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* rightPtr = right.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
                                                                     count);
 
                 return;
@@ -69,7 +69,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* rightPtr = right.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
                                                                     count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
@@ -83,7 +83,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f64].buffer_add;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_add;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -102,12 +102,12 @@ namespace StgSharp.Mathematics.Numeric
             scalar->Data<double>(0) = value;
             if (count <= 256)
             {
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
                 return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
             }
@@ -117,7 +117,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f64].buffer_fill;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_fill;
             package->Scalar = scalar;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -158,7 +158,7 @@ namespace StgSharp.Mathematics.Numeric
                         enumeration->Member<int>(0) = i;
                         enumeration->Member<int>(1) = j;
 
-                        MathematicsModule.GlobalContext
+                        NumericalModule.GlobalContext
                                            .mat_sk[f64].kernel_tile_fma((MatrixKernel*)left.Buffer,
                                                                      (MatrixKernel*)right.Buffer,
                                                                      (MatrixKernel*)ans.Buffer,
@@ -209,7 +209,7 @@ namespace StgSharp.Mathematics.Numeric
             {
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
                                                                            count);
 
                 return;
@@ -219,7 +219,7 @@ namespace StgSharp.Mathematics.Numeric
 
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
                                                                            count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
@@ -233,7 +233,7 @@ namespace StgSharp.Mathematics.Numeric
             package->SecCount = left.KernelRowLength;
             package->Scalar = scalar;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f64].buffer_scalar_mul;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);
@@ -266,7 +266,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* rightPtr = right.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
                                                                     count);
 
                 return;
@@ -277,7 +277,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<double>* leftPtr = left.Buffer;
                 MatrixKernel<double>* rightPtr = right.Buffer;
                 MatrixKernel<double>* ansPtr = ans.Buffer;
-                MathematicsModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
                                                                     count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
@@ -291,7 +291,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = left.KernelColumnLength;
             package->SecCount = left.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)MathematicsModule.GlobalContext.mat_sk[f64].buffer_sub;
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_sub;
             MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
 
             taskGroup.LaunchParallel(SleepMode.DeepSleep);

StgSharp.Mathematics/Numeric/NumericalModule.cs

@@ -1,4 +1,4 @@
-//-----------------------------------------------------------------------
+//-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
 // file="NumericalModule"
 // Project: StgSharp
@@ -25,9 +25,11 @@
 //     
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
+using StgSharp.Mathematics.Internal;
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Runtime.CompilerServices;
 using System.Text;
 using System.Threading.Tasks;
 
@@ -36,12 +38,53 @@ namespace StgSharp.Mathematics.Numeric
     public class NumericalModule : IStaticModule
     {
 
+        public static unsafe int DefaultSIMDAlignment
+        {
+            get
+            {
+                SIMDID id = GlobalIntrinsicMask;
+                return (id.MaskByte[0] & 0b_00001111) switch
+                {
+                    1 => id.MaskByte[1] switch
+                    {
+                        1 => 16,
+                        2 => 32,
+                        3 => 64,
+                        4 => 64,
+                        5 => 64,
+                        _ => 8,
+                    },
+                    _ => 8,
+                };
+            }
+        }
+
         public string ModuleName => "Math.Numerical";
 
+        internal static IntrinsicContext GlobalContext
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get;
+        } = new();
+
+        internal static SIMDID GlobalIntrinsicMask { get; set; }
+
         public void InitializeModule()
+        {
+            throw new NotImplementedException();
+        }
+
+        public void InitializeModule(
+                    IModuleInitializeProfile profile
+        )
         {
             MatrixParallel.Init();
         }
 
+        public void UninitializeModule()
+        {
+            throw new NotImplementedException();
+        }
+
     }
 }

StgSharp.Mathematics/Numeric/Runtime/MatrixParallelQueue.cs

@@ -33,11 +33,12 @@ using System.Collections.Generic;
 using System.Linq;
 using System.Runtime.CompilerServices;
 using System.Text;
+using System.Threading;
 using System.Threading.Tasks;
 
 namespace StgSharp.Mathematics.Numeric
 {
-    public unsafe class MatrixParallelQueue
+    internal unsafe class MatrixParallelQueue
     {
 
         private static ConcurrentStack<MatrixParallelQueue> _recycle = [];
@@ -57,7 +58,7 @@ namespace StgSharp.Mathematics.Numeric
                 queue = new MatrixParallelQueue();
             }
             queue._package = pack;
-            queue._maxCount = (ulong)pack->PrimCount * (ulong)pack->SecCount;
+            queue._maxCount = (ulong)pack->CommonX * (ulong)pack->CommonY;
             queue._cursor = 0;
             return queue;
         }

StgSharp.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs

@@ -65,7 +65,7 @@ namespace StgSharp.Mathematics.Numeric
                 int j = (int)(c % ansHeight);
                 enumeration->Member<int>(0) = i;
                 enumeration->Member<int>(1) = j;
-                MathematicsModule.GlobalContext.mat_sk[f32].kernel_tile_fma(leftBuffer, rightBuffer,
+                NumericalModule.GlobalContext.mat_sk[f32].kernel_tile_fma(leftBuffer, rightBuffer,
                                                                          ansBuffer, enumeration);
             }
         }
@@ -97,7 +97,7 @@ namespace StgSharp.Mathematics.Numeric
                 int j = (int)(c % ansHeight);
                 enumeration->Member<int>(0) = i;
                 enumeration->Member<int>(1) = j;
-                MathematicsModule.GlobalContext.mat_sk[f64].kernel_tile_fma(leftBuffer, rightBuffer,
+                NumericalModule.GlobalContext.mat_sk[f64].kernel_tile_fma(leftBuffer, rightBuffer,
                                                                          ansBuffer, enumeration);
             }
         }

StgSharp.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs

@@ -63,8 +63,6 @@ namespace StgSharp.Mathematics.Numeric
             _count = threadArray.Length;
         }
 
-        public MatrixParallelQueue TaskQueue { get; set; }
-
         public ReadOnlySpan<MatrixParallelThread> Threads
         {
             get
@@ -77,6 +75,8 @@ namespace StgSharp.Mathematics.Numeric
 
         public ThreadPriority AfterWork { get; private set; }
 
+        internal MatrixParallelQueue TaskQueue { get; set; }
+
         private bool IsMultithread { get; init; }
 
         public void Dispose()

StgSharp.Mathematics/Numeric/Scalar/Factorial.cs

@@ -54,7 +54,7 @@ namespace StgSharp.Mathematics
                 return factorialResult[n];
             }
             if (n < 30) {
-                return MathematicsModule.GlobalContext.factorial_simd(n);
+                return NumericalModule.GlobalContext.factorial_simd(n);
             }
             double 
                 dn = n,

StgSharp.Native/StgSharp.Native.vcxproj

@@ -176,9 +176,9 @@
     <ClCompile Include="sn-src\dllmain.c" />
     <ClCompile Include="sn-src\m.bufop.x86.c" />
     <ClCompile Include="sn-src\m.solve.x86.c" />
-    <ClCompile Include="sn-src\m.fma_tile.x86.c" />
+    <ClCompile Include="sn-src\mk.fma.x86.c" />
     <ClCompile Include="sn-src\s.factorial.x86.c" />
-    <ClCompile Include="sn-src\m.fma.x86.c" />
+    <ClCompile Include="sn-src\sk.fma.x86.c" />
     <ClCompile Include="sn-src\s.hash.x86.c" />
     <ClCompile Include="sn-src\load_intrinsic.x86.c" />
     <ClCompile Include="sn-src\load_glfw.c" />

StgSharp.Native/StgSharp.Native.vcxproj.filters

@@ -42,7 +42,7 @@
     <ClCompile Include="sn-src\m.basic.x86.c">
       <Filter>源文件</Filter>
     </ClCompile>
-    <ClCompile Include="sn-src\m.fma.x86.c">
+    <ClCompile Include="sn-src\sk.fma.x86.c">
       <Filter>源文件</Filter>
     </ClCompile>
     <ClCompile Include="sn-src\load_intrinsic.x86.c">
@@ -57,7 +57,7 @@
     <ClCompile Include="sn-src\simdid.x86.c">
       <Filter>源文件</Filter>
     </ClCompile>
-    <ClCompile Include="sn-src\m.fma_tile.x86.c">
+    <ClCompile Include="sn-src\mk.fma.x86.c">
       <Filter>源文件</Filter>
     </ClCompile>
     <ClCompile Include="sn-src\t.affinity.c">

StgSharp.Native/sn-head/sn_intrinsic.std.h

@@ -14,12 +14,12 @@
         T_type##_multi_kernel_proc_##T_op##_##T_arch##_##T_feature
 #define SN_MK_PROC_DECL(T_type, P_name) INTERNAL void SN_DECL P_name(MAT_TASK(T_type) const *task)
 #define SN_MK_PROC_DECL_STD(T_type, T_arch, T_feature, T_op)                          \
-        INTERNAL void SN_DECL MK_PROC_NAME_STD(T_type, T_arch, T_feature, \
+        INTERNAL FORCEINLINE void SN_DECL MK_PROC_NAME_STD(T_type, T_arch, T_feature, \
                                                            T_op)(MAT_TASK(T_type) const *task)
 #define SK_PROC_NAME_STD(T_type, T_arch, T_feature, T_op) \
         T_type##_single_kernel_proc_##T_op##_##T_arch##_##T_feature
 #define SN_SK_PROC_DECL_STD(T_type, T_arch, T_feature, T_op)                                 \
-        INTERNAL void SN_DECL SK_PROC_NAME_STD(T_type, T_arch, T_feature, T_op)( \
+        INTERNAL FORCEINLINE void SN_DECL SK_PROC_NAME_STD(T_type, T_arch, T_feature, T_op)( \
                 MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right,             \
                 MAT_KERNEL(T_type) *restrict ans, uint64_t scalar)
 

StgSharp.Native/sn-src/m.fma.x86.c

@@ -1,367 +0,0 @@
-#include "sn_target.h"
-
-#if SN_IS_ARCH(SN_ARCH_X86_64)
-
-#include "sn_intrinsic.h"
-#include "sn_intrinsic.context.x86.h"
-#include "sn_intrinsic.std.h"
-
-#define _FMA_SELECT(i) _MM_SHUFFLE(i, i, i, i)
-
-// float
-
-SN_SK_PROC_DECL_STD(float, sse, , fma)
-{
-        const __m128 a0 = _mm_load_ps((float *)&left->f32_x[0]);
-        const __m128 a1 = _mm_load_ps((float *)&left->f32_x[1]);
-        const __m128 a2 = _mm_load_ps((float *)&left->f32_x[2]);
-        const __m128 a3 = _mm_load_ps((float *)&left->f32_x[3]);
-
-        __m128 c0 = _mm_load_ps((float *)&ans->f32_x[0]);
-        __m128 c1 = _mm_load_ps((float *)&ans->f32_x[1]);
-        __m128 c2 = _mm_load_ps((float *)&ans->f32_x[2]);
-        __m128 c3 = _mm_load_ps((float *)&ans->f32_x[3]);
-
-        __m128 bcol;
-
-        /* Column 0 */
-        bcol = _mm_load_ps((float *)&right->f32_x[0]);
-        c0 = _mm_add_ps(c0, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
-        c0 = _mm_add_ps(c0, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
-        c0 = _mm_add_ps(c0, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
-        c0 = _mm_add_ps(c0, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
-
-        /* Column 1 */
-        bcol = _mm_load_ps((float *)&right->f32_x[1]);
-        c1 = _mm_add_ps(c1, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
-        c1 = _mm_add_ps(c1, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
-        c1 = _mm_add_ps(c1, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
-        c1 = _mm_add_ps(c1, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
-
-        /* Column 2 */
-        bcol = _mm_load_ps((float *)&right->f32_x[2]);
-        c2 = _mm_add_ps(c2, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
-        c2 = _mm_add_ps(c2, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
-        c2 = _mm_add_ps(c2, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
-        c2 = _mm_add_ps(c2, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
-
-        /* Column 3 */
-        bcol = _mm_load_ps((float *)&right->f32_x[3]);
-        c3 = _mm_add_ps(c3, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
-        c3 = _mm_add_ps(c3, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
-        c3 = _mm_add_ps(c3, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
-        c3 = _mm_add_ps(c3, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
-
-        _mm_store_ps((float *)&ans->f32_x[0], c0);
-        _mm_store_ps((float *)&ans->f32_x[1], c1);
-        _mm_store_ps((float *)&ans->f32_x[2], c2);
-        _mm_store_ps((float *)&ans->f32_x[3], c3);
-}
-
-SN_SK_PROC_DECL_STD(float, avx, , fma)
-{
-        static const __m256i b01_shift = { .m256i_i32 = { 0, 1, 2, 3, 5, 6, 7, 4 } };
-        static const __m256i b23_shift = { .m256i_i32 = { 2, 3, 0, 1, 7, 4, 5, 6 } };
-
-        register __m256 a01 = _mm256_load_ps((float *)&left->f32_y[0]);
-        register __m256 a23 = _mm256_load_ps((float *)&left->f32_y[1]);
-
-        register __m256 b01 = _mm256_load_ps((float *)&right->f32_y[0]);
-        register __m256 b23 = _mm256_load_ps((float *)&right->f32_y[1]);
-        b01 = _mm256_permutevar_ps(b01, b01_shift);
-        b23 = _mm256_permutevar_ps(b23, b23_shift);
-
-        register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]);
-        register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]);
-
-        //group 0
-        register __m256 b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(0));
-        register __m256 b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(0));
-        c01 = _mm256_add_ps(_mm256_mul_ps(a01, b01_perm), c01);
-        c23 = _mm256_add_ps(_mm256_mul_ps(a23, b23_perm), c23);
-
-        //group 2
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(2));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(2));
-        c01 = _mm256_add_ps(_mm256_mul_ps(a23, b01_perm), c01);
-        c23 = _mm256_add_ps(_mm256_mul_ps(a01, b23_perm), c23);
-
-        //group 1
-        register __m256 a12 = _mm256_permute2f128_ps(a01, a23, 0x21);
-        register __m256 a30 = _mm256_permute2f128_ps(a01, a23, 0x30);
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(1));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(1));
-        c01 = _mm256_add_ps(_mm256_mul_ps(a12, b01_perm), c01);
-        c23 = _mm256_add_ps(_mm256_mul_ps(a30, b23_perm), c23);
-
-        //group 3
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(3));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(3));
-        c01 = _mm256_add_ps(_mm256_mul_ps(a30, b01_perm), c01);
-        c23 = _mm256_add_ps(_mm256_mul_ps(a12, b23_perm), c23);
-
-        ans->f32_y[0] = c01;
-        ans->f32_y[1] = c23;
-}
-
-SN_SK_PROC_DECL_STD(float, avx, _fma, fma)
-{
-        static const __m256i b01_shift = { .m256i_i32 = { 0, 1, 2, 3, 5, 6, 7, 4 } };
-        static const __m256i b23_shift = { .m256i_i32 = { 2, 3, 0, 1, 7, 4, 5, 6 } };
-
-        register __m256 a01 = _mm256_load_ps((float *)&left->f32_y[0]);
-        register __m256 a23 = _mm256_load_ps((float *)&left->f32_y[1]);
-
-        register __m256 b01 = _mm256_load_ps((float *)&right->f32_y[0]);
-        register __m256 b23 = _mm256_load_ps((float *)&right->f32_y[1]);
-        b01 = _mm256_permutevar_ps(b01, b01_shift);
-        b23 = _mm256_permutevar_ps(b23, b23_shift);
-
-        register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]);
-        register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]);
-
-        //group 0
-        register __m256 b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(0));
-        register __m256 b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(0));
-        c01 = _mm256_fmadd_ps(a01, b01_perm, c01);
-        c23 = _mm256_fmadd_ps(a23, b23_perm, c23);
-
-        //group 2
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(2));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(2));
-        c01 = _mm256_fmadd_ps(a23, b01_perm, c01);
-        c23 = _mm256_fmadd_ps(a01, b23_perm, c23);
-
-        //group 1
-        register __m256 a12 = _mm256_permute2f128_ps(a01, a23, 0x21);
-        register __m256 a30 = _mm256_permute2f128_ps(a01, a23, 0x30);
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(1));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(1));
-        c01 = _mm256_fmadd_ps(a12, b01_perm, c01);
-        c23 = _mm256_fmadd_ps(a30, b23_perm, c23);
-
-        //group 3
-        b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(3));
-        b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(3));
-        c01 = _mm256_fmadd_ps(a30, b01_perm, c01);
-        c23 = _mm256_fmadd_ps(a12, b23_perm, c23);
-
-        ans->f32_y[0] = c01;
-        ans->f32_y[1] = c23;
-}
-
-SN_SK_PROC_DECL_STD(float, 512, , fma)
-{
-        static const __m512i b_shift = { .m512i_i32 = {
-                                                 // shuff order
-                                                 0, 1, 2, 3, // col0
-                                                 5, 6, 7, 4, //col1
-                                                 10, 11, 8, 9, //col2
-                                                 15, 12, 13, 14 //col3
-                                         } };
-
-        register __m512 a = _mm512_load_ps((float *)left);
-        register __m512 b = _mm512_load_ps((float *)right);
-        register __m512 c = _mm512_load_ps((float *)ans);
-
-        b = _mm512_permutexvar_ps(b_shift, b);
-        register __m512 a_perm0 = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(3, 2, 1, 0));
-        register __m512 b_perm0 = _mm512_permute_ps(b, _FMA_SELECT(0));
-        c = _mm512_fmadd_ps(a_perm0, b_perm0, c);
-
-        register __m512 a_perm1 = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(2, 1, 0, 3));
-        register __m512 b_perm1 = _mm512_permute_ps(b, _FMA_SELECT(1));
-        c = _mm512_fmadd_ps(a_perm1, b_perm1, c);
-
-        register __m512 a_perm2 = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(1, 0, 3, 2));
-        register __m512 b_perm2 = _mm512_permute_ps(b, _FMA_SELECT(2));
-        c = _mm512_fmadd_ps(a_perm2, b_perm2, c);
-
-        register __m512 a_perm3 = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(0, 3, 2, 1));
-        register __m512 b_perm3 = _mm512_permute_ps(b, _FMA_SELECT(3));
-        c = _mm512_fmadd_ps(a_perm3, b_perm3, c);
-
-        _mm512_store_ps((float *)ans, c);
-}
-
-// double
-
-SN_SK_PROC_DECL_STD(double, sse, , fma)
-{
-        for (int col = 0; col < 4; ++col) {
-                /* current answer column (two chunks: rows 0-1 and 2-3) */
-                register __m128d c_lo = _mm_load_pd(ans->f64_x + (col * 2 + 0));
-                register __m128d c_hi = _mm_load_pd(ans->f64_x + (col * 2 + 1));
-
-                /* right matrix column split into two chunks */
-                register __m128d b_lo = _mm_load_pd(right + (2 * col + 0));
-                register __m128d b_hi = _mm_load_pd(right + (2 * col + 1));
-
-                /* k = 0 */
-                register __m128d a_lo = _mm_load_pd(left + (0 * 2 + 0));
-                register __m128d a_hi = _mm_load_pd(left + (0 * 2 + 1));
-                register __m128d b = _mm_shuffle_pd(b_lo, b_lo, 0x0);
-                c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b));
-                c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b));
-
-                /* k = 1 */
-                a_lo = _mm_load_pd(left + (1 * 2 + 0));
-                a_hi = _mm_load_pd(left + (1 * 2 + 1));
-                b = _mm_shuffle_pd(b_lo, b_lo, 0x3);
-                c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b));
-                c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b));
-
-                /* k = 2 */
-                a_lo = _mm_load_pd(left + (1 * 2 + 0));
-                a_hi = _mm_load_pd(left + (1 * 2 + 1));
-                b = _mm_shuffle_pd(b_hi, b_hi, 0x0);
-                c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b));
-                c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b));
-
-                /* k = 3 */
-                a_lo = _mm_load_pd(left + (1 * 2 + 0));
-                a_hi = _mm_load_pd(left + (1 * 2 + 1));
-                b = _mm_shuffle_pd(b_hi, b_hi, 0x3);
-                c_lo = _mm_add_pd(c_lo, _mm_mul_pd(a_lo, b));
-                c_hi = _mm_add_pd(c_hi, _mm_mul_pd(a_hi, b));
-
-                ans->f64_x[col * 2 + 0] = c_lo;
-                ans->f64_x[col * 2 + 1] = c_hi;
-        }
-}
-
-SN_SK_PROC_DECL_STD(double, avx, , fma)
-{
-        const __m256d a0 = _mm256_load_pd((double *)&left->f64_y[0]);
-        const __m256d a1 = _mm256_load_pd((double *)&left->f64_y[1]);
-        const __m256d a2 = _mm256_load_pd((double *)&left->f64_y[2]);
-        const __m256d a3 = _mm256_load_pd((double *)&left->f64_y[3]);
-
-        __m256d b0 = _mm256_load_pd((double *)&right->f64_y[0]);
-        __m256d b1 = _mm256_load_pd((double *)&right->f64_y[1]);
-        __m256d b2 = _mm256_load_pd((double *)&right->f64_y[2]);
-        __m256d b3 = _mm256_load_pd((double *)&right->f64_y[3]);
-
-        register __m256d ccol;
-
-        /* Column 0 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[0]);
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a0, _mm256_permute4x64_pd(b0, _FMA_SELECT(0))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a1, _mm256_permute4x64_pd(b0, _FMA_SELECT(1))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a2, _mm256_permute4x64_pd(b0, _FMA_SELECT(2))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a3, _mm256_permute4x64_pd(b0, _FMA_SELECT(3))));
-        _mm256_store_pd((double *)&ans->f64_y[0], ccol);
-        /* Column 1 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[1]);
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a0, _mm256_permute4x64_pd(b1, _FMA_SELECT(0))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a1, _mm256_permute4x64_pd(b1, _FMA_SELECT(1))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a2, _mm256_permute4x64_pd(b1, _FMA_SELECT(2))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a3, _mm256_permute4x64_pd(b1, _FMA_SELECT(3))));
-        _mm256_store_pd((double *)&ans->f64_y[1], ccol);
-        /* Column 2 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[2]);
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a0, _mm256_permute4x64_pd(b2, _FMA_SELECT(0))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a1, _mm256_permute4x64_pd(b2, _FMA_SELECT(1))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a2, _mm256_permute4x64_pd(b2, _FMA_SELECT(2))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a3, _mm256_permute4x64_pd(b2, _FMA_SELECT(3))));
-        _mm256_store_pd((double *)&ans->f64_y[2], ccol);
-        /* Column 3 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[3]);
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a0, _mm256_permute4x64_pd(b3, _FMA_SELECT(0))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a1, _mm256_permute4x64_pd(b3, _FMA_SELECT(1))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a2, _mm256_permute4x64_pd(b3, _FMA_SELECT(2))));
-        ccol = _mm256_add_pd(ccol, _mm256_mul_pd(a3, _mm256_permute4x64_pd(b3, _FMA_SELECT(3))));
-        _mm256_store_pd((double *)&ans->f64_y[3], ccol);
-}
-
-SN_SK_PROC_DECL_STD(double, avx, _fma, fma)
-{
-        const __m256d a0 = _mm256_load_pd((double *)&left->f64_y[0]);
-        const __m256d a1 = _mm256_load_pd((double *)&left->f64_y[1]);
-        const __m256d a2 = _mm256_load_pd((double *)&left->f64_y[2]);
-        const __m256d a3 = _mm256_load_pd((double *)&left->f64_y[3]);
-
-        __m256d b0 = _mm256_load_pd((double *)&right->f64_y[0]);
-        __m256d b1 = _mm256_load_pd((double *)&right->f64_y[1]);
-        __m256d b2 = _mm256_load_pd((double *)&right->f64_y[2]);
-        __m256d b3 = _mm256_load_pd((double *)&right->f64_y[3]);
-
-        register __m256d ccol;
-
-        /* Column 0 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[0]);
-        ccol = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(b0, _FMA_SELECT(0)), ccol);
-        ccol = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(b0, _FMA_SELECT(1)), ccol);
-        ccol = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(b0, _FMA_SELECT(2)), ccol);
-        ccol = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(b0, _FMA_SELECT(3)), ccol);
-        _mm256_store_pd((double *)&ans->f64_y[0], ccol);
-        /* Column 1 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[1]);
-        ccol = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(b1, _FMA_SELECT(0)), ccol);
-        ccol = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(b1, _FMA_SELECT(1)), ccol);
-        ccol = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(b1, _FMA_SELECT(2)), ccol);
-        ccol = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(b1, _FMA_SELECT(3)), ccol);
-        _mm256_store_pd((double *)&ans->f64_y[1], ccol);
-        /* Column 2 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[2]);
-        ccol = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(b2, _FMA_SELECT(0)), ccol);
-        ccol = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(b2, _FMA_SELECT(1)), ccol);
-        ccol = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(b2, _FMA_SELECT(2)), ccol);
-        ccol = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(b2, _FMA_SELECT(3)), ccol);
-        _mm256_store_pd((double *)&ans->f64_y[2], ccol);
-        /* Column 3 */
-        ccol = _mm256_load_pd((double *)&ans->f64_y[3]);
-        ccol = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(b3, _FMA_SELECT(0)), ccol);
-        ccol = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(b3, _FMA_SELECT(1)), ccol);
-        ccol = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(b3, _FMA_SELECT(2)), ccol);
-        ccol = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(b3, _FMA_SELECT(3)), ccol);
-        _mm256_store_pd((double *)&ans->f64_y[3], ccol);
-}
-
-SN_SK_PROC_DECL_STD(double, 512, , fma)
-{
-        static const __m512i b01_shift = { .m512i_i64 = { 0, 1, 2, 3, 5, 6, 7, 4 } };
-        static const __m512i b23_shift = { .m512i_i64 = { 2, 3, 0, 1, 7, 4, 5, 6 } };
-
-        register __m512d a01 = _mm512_load_pd((double *)&left->f64_z[0]);
-        register __m512d a23 = _mm512_load_pd((double *)&left->f64_z[1]);
-
-        register __m512d b01 = _mm512_load_pd((double *)&right->f64_z[0]);
-        register __m512d b23 = _mm512_load_pd((double *)&right->f64_z[1]);
-        b01 = _mm512_permutexvar_pd(b01_shift, b01);
-        b23 = _mm512_permutexvar_pd(b23_shift, b23);
-
-        register __m512d c01 = _mm512_load_pd((float *)&ans->f64_z[0]);
-        register __m512d c23 = _mm512_load_pd((float *)&ans->f64_z[1]);
-
-        //group 0
-        register __m512d b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(0));
-        register __m512d b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(0));
-        c01 = _mm512_fmadd_pd(a01, b01_perm, c01);
-        c23 = _mm512_fmadd_pd(a23, b23_perm, c23);
-
-        //group 2
-        b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(2));
-        b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(2));
-        c01 = _mm512_fmadd_pd(a23, b01_perm, c01);
-        c23 = _mm512_fmadd_pd(a01, b23_perm, c23);
-
-        //group 1
-        register __m512d a12 = _mm512_shuffle_f64x2(a01, a23, _MM_SHUFFLE(1, 0, 3, 2));
-        register __m512d a30 = _mm512_shuffle_f64x2(a23, a01, _MM_SHUFFLE(1, 0, 3, 2));
-        b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(1));
-        b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(1));
-        c01 = _mm512_fmadd_pd(a12, b01_perm, c01);
-        c23 = _mm512_fmadd_pd(a30, b23_perm, c23);
-
-        //group 3
-        b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(3));
-        b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(3));
-        c01 = _mm512_fmadd_pd(a30, b01_perm, c01);
-        c23 = _mm512_fmadd_pd(a12, b23_perm, c23);
-
-        _mm512_store_pd((double *)&ans->f64_z[0], c01);
-        _mm512_store_pd((double *)&ans->f64_z[1], c23);
-}
-
-#endif

StgSharp.Native/sn-src/mk.fma.x86.c

@@ -82,46 +82,25 @@ SN_MK_PROC_DECL_STD(float, avx, , fma)
         register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]);
         register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]);
 
-        register __m256i b01_shift = _mm256_load_si256(&b01_shift_source);
-        register __m256i b23_shift = _mm256_load_si256(&b23_shift_source);
-
         for (size_t z = 0; z < common_z; z++) {
                 SN_MAT_KERNEL_LOCAL const *right_cur = right + GET_INDEX(common_z, x, z),
                                           const *left_cur = left + GET_INDEX(common_y, z, y);
 
-                register __m256 a01 = _mm256_load_ps((float *)&left_cur->f32_y[0]);
-                register __m256 a23 = _mm256_load_ps((float *)&left_cur->f32_y[1]);
-
                 register __m256 b01 = _mm256_load_ps((float *)&right_cur->f32_y[0]);
                 register __m256 b23 = _mm256_load_ps((float *)&right_cur->f32_y[1]);
-                b01 = _mm256_permutevar_ps(b01, b01_shift);
-                b23 = _mm256_permutevar_ps(b23, b23_shift);
 
-                //group 0
+                register __m256 a0 = _mm256_broadcast_ps(&left_cur->f32_x[0]);
-                register __m256 b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(0));
+                c01 = _mm256_add_ps(_mm256_mul_ps(a0, b01), c01);
-                register __m256 b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(0));
+                c23 = _mm256_add_ps(_mm256_mul_ps(a0, b23), c23);
-                c01 = _mm256_add_ps(_mm256_mul_ps(a01, b01_perm), c01);
+                register __m256 a1 = _mm256_broadcast_ps(&left_cur->f32_x[1]);
-                c23 = _mm256_add_ps(_mm256_mul_ps(a23, b23_perm), c23);
+                c01 = _mm256_add_ps(_mm256_mul_ps(a1, b01), c01);
-
+                c23 = _mm256_add_ps(_mm256_mul_ps(a1, b23), c23);
-                //group 2
+                register __m256 a2 = _mm256_broadcast_ps(&left_cur->f32_x[2]);
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(2));
+                c01 = _mm256_add_ps(_mm256_mul_ps(a2, b01), c01);
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(2));
+                c23 = _mm256_add_ps(_mm256_mul_ps(a2, b23), c23);
-                c01 = _mm256_add_ps(_mm256_mul_ps(a23, b01_perm), c01);
+                register __m256 a3 = _mm256_broadcast_ps(&left_cur->f32_x[3]);
-                c23 = _mm256_add_ps(_mm256_mul_ps(a01, b23_perm), c23);
+                c01 = _mm256_add_ps(_mm256_mul_ps(a3, b01), c01);
-
+                c23 = _mm256_add_ps(_mm256_mul_ps(a3, b23), c23);
-                //group 1
-                register __m256 a12 = _mm256_permute2f128_ps(a01, a23, 0x21);
-                register __m256 a30 = _mm256_permute2f128_ps(a01, a23, 0x30);
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(1));
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(1));
-                c01 = _mm256_add_ps(_mm256_mul_ps(a12, b01_perm), c01);
-                c23 = _mm256_add_ps(_mm256_mul_ps(a30, b23_perm), c23);
-
-                //group 3
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(3));
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(3));
-                c01 = _mm256_add_ps(_mm256_mul_ps(a30, b01_perm), c01);
-                c23 = _mm256_add_ps(_mm256_mul_ps(a12, b23_perm), c23);
         }
         _mm256_store_ps((float *)&ans_cur->f32_y[0], c01);
         _mm256_store_ps((float *)&ans_cur->f32_y[1], c23);
@@ -129,9 +108,6 @@ SN_MK_PROC_DECL_STD(float, avx, , fma)
 
 SN_MK_PROC_DECL_STD(float, avx, _fma, fma)
 {
-        static const __m256i b01_shift_source = { .m256i_i32 = { 0, 1, 2, 3, 5, 6, 7, 4 } };
-        static const __m256i b23_shift_source = { .m256i_i32 = { 2, 3, 0, 1, 7, 4, 5, 6 } };
-
         __mk_param_std(TILE, LEFT_RIGHT_ANS, float);
 
         MAT_KERNEL(float)
@@ -140,46 +116,24 @@ SN_MK_PROC_DECL_STD(float, avx, _fma, fma)
         register __m256 c01 = _mm256_load_ps((float *)&ans_cur->f32_y[0]);
         register __m256 c23 = _mm256_load_ps((float *)&ans_cur->f32_y[1]);
 
-        register __m256i b01_shift = _mm256_load_si256(&b01_shift_source);
-        register __m256i b23_shift = _mm256_load_si256(&b23_shift_source);
-
         for (size_t z = 0; z < common_z; z++) {
                 SN_MAT_KERNEL_LOCAL const *right_cur = right + GET_INDEX(common_z, x, z),
                                           const *left_cur = left + GET_INDEX(common_y, z, y);
-
-                register __m256 a01 = _mm256_load_ps((float *)&left_cur->f32_y[0]);
-                register __m256 a23 = _mm256_load_ps((float *)&left_cur->f32_y[1]);
-
                 register __m256 b01 = _mm256_load_ps((float *)&right_cur->f32_y[0]);
                 register __m256 b23 = _mm256_load_ps((float *)&right_cur->f32_y[1]);
-                b01 = _mm256_permutevar_ps(b01, b01_shift);
-                b23 = _mm256_permutevar_ps(b23, b23_shift);
 
-                //group 0
+                register __m256 a0 = _mm256_broadcast_ps((float *)&left_cur->f32_x[0]);
-                register __m256 b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(0));
+                c01 = _mm256_fmadd_ps(a0, b01, c01);
-                register __m256 b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(0));
+                c23 = _mm256_fmadd_ps(a0, b23, c23);
-                c01 = _mm256_fmadd_ps(a01, b01_perm, c01);
+                register __m256 a1 = _mm256_broadcast_ps((float *)&left_cur->f32_x[1]);
-                c23 = _mm256_fmadd_ps(a23, b23_perm, c23);
+                c01 = _mm256_fmadd_ps(a1, b01, c01);
-
+                c23 = _mm256_fmadd_ps(a1, b23, c23);
-                //group 2
+                register __m256 a2 = _mm256_broadcast_ps((float *)&left_cur->f32_x[2]);
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(2));
+                c01 = _mm256_fmadd_ps(a2, b01, c01);
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(2));
+                c23 = _mm256_fmadd_ps(a2, b23, c23);
-                c01 = _mm256_fmadd_ps(a23, b01_perm, c01);
+                register __m256 a3 = _mm256_broadcast_ps((float *)&left_cur->f32_x[3]);
-                c23 = _mm256_fmadd_ps(a01, b23_perm, c23);
+                c01 = _mm256_fmadd_ps(a3, b01, c01);
-
+                c23 = _mm256_fmadd_ps(a3, b23, c23);
-                //group 1
-                register __m256 a12 = _mm256_permute2f128_ps(a01, a23, 0x21);
-                register __m256 a30 = _mm256_permute2f128_ps(a01, a23, 0x30);
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(1));
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(1));
-                c01 = _mm256_fmadd_ps(a12, b01_perm, c01);
-                c23 = _mm256_fmadd_ps(a30, b23_perm, c23);
-
-                //group 3
-                b01_perm = _mm256_shuffle_ps(b01, b01, _FMA_SELECT(3));
-                b23_perm = _mm256_shuffle_ps(b23, b23, _FMA_SELECT(3));
-                c01 = _mm256_fmadd_ps(a30, b01_perm, c01);
-                c23 = _mm256_fmadd_ps(a12, b23_perm, c23);
         }
         _mm256_store_ps((float *)&ans_cur->f32_y[0], c01);
         _mm256_store_ps((float *)&ans_cur->f32_y[1], c23);
@@ -187,42 +141,34 @@ SN_MK_PROC_DECL_STD(float, avx, _fma, fma)
 
 SN_MK_PROC_DECL_STD(float, 512, , fma)
 {
-        static const __m512i b_shift_source = { .m512i_i32 = {
-                                                        // shuff order
-                                                        0, 1, 2, 3, // col0
-                                                        5, 6, 7, 4, //col1
-                                                        10, 11, 8, 9, //col2
-                                                        15, 12, 13, 14 //col3
-                                                } };
-        register __m512i b_shift = _mm512_load_si512(&b_shift_source);
         __mk_param_std(TILE, LEFT_RIGHT_ANS, float);
         SN_MAT_KERNEL_LOCAL const *ans_cur = ans + GET_INDEX(common_y, x, y);
-        register __m512 c = _mm512_load_ps(ans_cur);
+        register __m512 c = _mm512_load_ps((float *)ans_cur);
-        for (size_t z = 0; z < common_z; z++) {
-                SN_MAT_KERNEL_LOCAL const *right_cur = right + GET_INDEX(common_z, x, z),
-                                          const *left_cur = left + GET_INDEX(common_y, z, y);
 
-                register __m512 a = _mm512_load_ps((float *)left_cur);
+        for (size_t z = 0; z < common_z; z++) {
+                SN_MAT_KERNEL_LOCAL const *right_cur = right + GET_INDEX(common_z, x, z);
+                SN_MAT_KERNEL_LOCAL const *left_cur = left + GET_INDEX(common_y, z, y);
+
+                // full right load：[kb_col0 | kb_col1 | kb_col2 | kb_col3]
                 register __m512 b = _mm512_load_ps((float *)right_cur);
 
-                b = _mm512_permutexvar_ps(b_shift, b);
+                // k=0: broadcast ka_col0 到4个lane，permute_ps 取 kb(0,j)
-                register __m512 a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(3, 2, 1, 0));
+                register __m512 ak;
-                register __m512 b_perm = _mm512_permute_ps(b, _FMA_SELECT(0));
+                ak = _mm512_broadcast_f32x4(left_cur->f32_x[0]);
-                c = _mm512_fmadd_ps(a_perm, b_perm, c);
+                c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(0)), c);
 
-                a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(2, 1, 0, 3));
+                // k=1
-                b_perm = _mm512_permute_ps(b, _FMA_SELECT(1));
+                ak = _mm512_broadcast_f32x4(left_cur->f32_x[1]);
-                c = _mm512_fmadd_ps(a_perm, b_perm, c);
+                c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(1)), c);
 
-                a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(1, 0, 3, 2));
+                // k=2
-                b_perm = _mm512_permute_ps(b, _FMA_SELECT(2));
+                ak = _mm512_broadcast_f32x4(left_cur->f32_x[2]);
-                c = _mm512_fmadd_ps(a_perm, b_perm, c);
+                c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(2)), c);
 
-                a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(0, 3, 2, 1));
+                // k=3
-                b_perm = _mm512_permute_ps(b, _FMA_SELECT(3));
+                ak = _mm512_broadcast_f32x4(left_cur->f32_x[3]);
-                c = _mm512_fmadd_ps(a_perm, b_perm, c);
+                c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(3)), c);
         }
-
         _mm512_store_ps((float *)ans_cur, c);
 }
 
@@ -431,20 +377,14 @@ SN_MK_PROC_DECL_STD(double, avx, _fma, fma)
 
 SN_MK_PROC_DECL_STD(double, 512, , fma)
 {
-        static const __m512i b01_shift_source = { .m512i_i64 = { 0, 1, 2, 3, 5, 6, 7, 4 } };
-        static const __m512i b23_shift_source = { .m512i_i64 = { 2, 3, 0, 1, 7, 4, 5, 6 } };
-
         __mk_param_std(TILE, LEFT_RIGHT_ANS, double);
 
         MAT_KERNEL(double)
-        const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe
+        const *ans_cur = ans + GET_INDEX(common_y, x, y);
 
         register __m512d c01 = _mm512_load_pd((double *)&ans->f64_z[0]);
         register __m512d c23 = _mm512_load_pd((double *)&ans->f64_z[1]);
 
-        register __m512i b01_shift = _mm512_load_si512(&b01_shift_source);
-        register __m512i b23_shift = _mm512_load_si512(&b23_shift_source);
-
         for (size_t z = 0; z < common_z; z++) {
                 SN_MAT_KERNEL_LOCAL const *right_cur = right + GET_INDEX(common_z, x, z),
                                           const *left_cur = left + GET_INDEX(common_y, z, y);
@@ -452,39 +392,23 @@ SN_MK_PROC_DECL_STD(double, 512, , fma)
                 register __m512d a01 = _mm512_load_pd((double *)&left_cur->f64_z[0]);
                 register __m512d a23 = _mm512_load_pd((double *)&left_cur->f64_z[1]);
 
-                register __m512d b01 = _mm512_load_pd((double *)&right_cur->f64_z[0]);
+                register __m512d b0 = _mm512_broadcast_f64x4(right_cur->f64_y[0]);
-                register __m512d b23 = _mm512_load_pd((double *)&right_cur->f64_z[1]);
+                c01 = _mm512_fmadd_pd(a01, b0, c01);
-                b01 = _mm512_permutexvar_pd(b01_shift, b01);
+                c23 = _mm512_fmadd_pd(a23, b0, c23);
-                b23 = _mm512_permutexvar_pd(b23_shift, b23);
 
-                //group 0
+                register __m512d b1 = _mm512_broadcast_f64x4(right_cur->f64_y[1]);
-                register __m512d b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(0));
+                c01 = _mm512_fmadd_pd(a01, b1, c01);
-                register __m512d b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(0));
+                c23 = _mm512_fmadd_pd(a23, b1, c23);
-                c01 = _mm512_fmadd_pd(a01, b01_perm, c01);
-                c23 = _mm512_fmadd_pd(a23, b23_perm, c23);
 
-                //group 2
+                register __m512d b2 = _mm512_broadcast_f64x4(right_cur->f64_y[2]);
-                b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(2));
+                c01 = _mm512_fmadd_pd(a01, b2, c01);
-                b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(2));
+                c23 = _mm512_fmadd_pd(a23, b2, c23);
-                c01 = _mm512_fmadd_pd(a23, b01_perm, c01);
-                c23 = _mm512_fmadd_pd(a01, b23_perm, c23);
 
-                //group 1
+                register __m512d b3 = _mm512_broadcast_f64x4(right_cur->f64_y[3]);
-                register __m512d a12 = _mm512_shuffle_f64x2(a01, a23, 0x21);
+                c01 = _mm512_fmadd_pd(a01, b3, c01);
-                register __m512d a30 = _mm512_shuffle_f64x2(a01, a23, 0x30);
+                c23 = _mm512_fmadd_pd(a23, b3, c23);
-                b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(1));
-                b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(1));
-                c01 = _mm512_fmadd_pd(a12, b01_perm, c01);
-                c23 = _mm512_fmadd_pd(a30, b23_perm, c23);
-
-                //group 3
-                b01_perm = _mm512_shuffle_pd(b01, b01, _FMA_SELECT(3));
-                b23_perm = _mm512_shuffle_pd(b23, b23, _FMA_SELECT(3));
-                c01 = _mm512_fmadd_pd(a30, b01_perm, c01);
-                c23 = _mm512_fmadd_pd(a12, b23_perm, c23);
         }
         _mm512_store_pd((double *)&ans_cur->f64_z[0], c01);
         _mm512_store_pd((double *)&ans_cur->f64_z[1], c23);
 }
-
 #endif

StgSharp.Native/sn-src/sk.fma.x86.c

@@ -0,0 +1,321 @@
+#include "sn_target.h"
+
+#if SN_IS_ARCH(SN_ARCH_X86_64)
+
+#include "sn_intrinsic.h"
+#include "sn_intrinsic.std.h"
+#include "sn_intrinsic.context.x86.h"
+
+#define _FMA_SELECT(i) _MM_SHUFFLE(i, i, i, i)
+
+#define GET_INDEX(col_len, x, y) (col_len * x + y)
+
+// float
+
+SN_SK_PROC_DECL_STD(float, sse, , fma)
+{
+        register __m128 c0 = _mm_load_ps((float *)&ans->f32_x[0]);
+        register __m128 c1 = _mm_load_ps((float *)&ans->f32_x[1]);
+        register __m128 c2 = _mm_load_ps((float *)&ans->f32_x[2]);
+        register __m128 c3 = _mm_load_ps((float *)&ans->f32_x[3]);
+        register __m128 bcol;
+
+        const __m128 a0 = _mm_load_ps((float *)&left->f32_x[0]);
+        const __m128 a1 = _mm_load_ps((float *)&left->f32_x[1]);
+        const __m128 a2 = _mm_load_ps((float *)&left->f32_x[2]);
+        const __m128 a3 = _mm_load_ps((float *)&left->f32_x[3]);
+
+        /* Column 0 */
+        bcol = _mm_load_ps((float *)&right->f32_x[0]);
+        c0 = _mm_add_ps(c0, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
+        c0 = _mm_add_ps(c0, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
+        c0 = _mm_add_ps(c0, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
+        c0 = _mm_add_ps(c0, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
+
+        /* Column 1 */
+        bcol = _mm_load_ps((float *)&right->f32_x[1]);
+        c1 = _mm_add_ps(c1, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
+        c1 = _mm_add_ps(c1, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
+        c1 = _mm_add_ps(c1, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
+        c1 = _mm_add_ps(c1, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
+
+        /* Column 2 */
+        bcol = _mm_load_ps((float *)&right->f32_x[2]);
+        c2 = _mm_add_ps(c2, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
+        c2 = _mm_add_ps(c2, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
+        c2 = _mm_add_ps(c2, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
+        c2 = _mm_add_ps(c2, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
+
+        /* Column 3 */
+        bcol = _mm_load_ps((float *)&right->f32_x[3]);
+        c3 = _mm_add_ps(c3, _mm_mul_ps(a0, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(0))));
+        c3 = _mm_add_ps(c3, _mm_mul_ps(a1, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(1))));
+        c3 = _mm_add_ps(c3, _mm_mul_ps(a2, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(2))));
+        c3 = _mm_add_ps(c3, _mm_mul_ps(a3, _mm_shuffle_ps(bcol, bcol, _FMA_SELECT(3))));
+
+        _mm_store_ps((float *)&ans->f32_x[0], c0);
+        _mm_store_ps((float *)&ans->f32_x[1], c1);
+        _mm_store_ps((float *)&ans->f32_x[2], c2);
+        _mm_store_ps((float *)&ans->f32_x[3], c3);
+}
+
+SN_SK_PROC_DECL_STD(float, avx, , fma)
+{
+        register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]);
+        register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]);
+
+        register __m256 b01 = _mm256_load_ps((float *)&right->f32_y[0]);
+        register __m256 b23 = _mm256_load_ps((float *)&right->f32_y[1]);
+
+        register __m256 a0 = _mm256_broadcast_ps(&left->f32_x[0]);
+        c01 = _mm256_add_ps(_mm256_mul_ps(a0, b01), c01);
+        c23 = _mm256_add_ps(_mm256_mul_ps(a0, b23), c23);
+        register __m256 a1 = _mm256_broadcast_ps(&left->f32_x[1]);
+        c01 = _mm256_add_ps(_mm256_mul_ps(a1, b01), c01);
+        c23 = _mm256_add_ps(_mm256_mul_ps(a1, b23), c23);
+        register __m256 a2 = _mm256_broadcast_ps(&left->f32_x[2]);
+        c01 = _mm256_add_ps(_mm256_mul_ps(a2, b01), c01);
+        c23 = _mm256_add_ps(_mm256_mul_ps(a2, b23), c23);
+        register __m256 a3 = _mm256_broadcast_ps(&left->f32_x[3]);
+        c01 = _mm256_add_ps(_mm256_mul_ps(a3, b01), c01);
+        c23 = _mm256_add_ps(_mm256_mul_ps(a3, b23), c23);
+
+        _mm256_store_ps((float *)&ans->f32_y[0], c01);
+        _mm256_store_ps((float *)&ans->f32_y[1], c23);
+}
+
+SN_SK_PROC_DECL_STD(float, avx, _fma, fma)
+{
+        register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]);
+        register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]);
+
+        register __m256 b01 = _mm256_load_ps((float *)&right->f32_y[0]);
+        register __m256 b23 = _mm256_load_ps((float *)&right->f32_y[1]);
+
+        register __m256 a0 = _mm256_broadcast_ps((float *)&left->f32_x[0]);
+        c01 = _mm256_fmadd_ps(a0, b01, c01);
+        c23 = _mm256_fmadd_ps(a0, b23, c23);
+        register __m256 a1 = _mm256_broadcast_ps((float *)&left->f32_x[1]);
+        c01 = _mm256_fmadd_ps(a1, b01, c01);
+        c23 = _mm256_fmadd_ps(a1, b23, c23);
+        register __m256 a2 = _mm256_broadcast_ps((float *)&left->f32_x[2]);
+        c01 = _mm256_fmadd_ps(a2, b01, c01);
+        c23 = _mm256_fmadd_ps(a2, b23, c23);
+        register __m256 a3 = _mm256_broadcast_ps((float *)&left->f32_x[3]);
+        c01 = _mm256_fmadd_ps(a3, b01, c01);
+        c23 = _mm256_fmadd_ps(a3, b23, c23);
+
+        _mm256_store_ps((float *)&ans->f32_y[0], c01);
+        _mm256_store_ps((float *)&ans->f32_y[1], c23);
+}
+
+SN_SK_PROC_DECL_STD(float, 512, , fma)
+{
+        register __m512 c = _mm512_load_ps((float *)ans);
+        // full right load：[kb_col0 | kb_col1 | kb_col2 | kb_col3]
+        register __m512 b = _mm512_load_ps((float *)right);
+        // k=0: broadcast ka_col0 到4个lane，permute_ps 取 kb(0,j)
+        register __m512 ak;
+        ak = _mm512_broadcast_f32x4(left->f32_x[0]);
+        c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(0)), c);
+        // k=1
+        ak = _mm512_broadcast_f32x4(left->f32_x[1]);
+        c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(1)), c);
+        // k=2
+        ak = _mm512_broadcast_f32x4(left->f32_x[2]);
+        c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(2)), c);
+        // k=3
+        ak = _mm512_broadcast_f32x4(left->f32_x[3]);
+        c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(3)), c);
+        _mm512_store_ps((float *)ans, c);
+}
+
+// double
+
+// clang-format off
+
+#define DOUBLE_KERNEL_FMA_CYCLE(i)                                            \
+        do {                                                                  \
+                /* k = 0 */                                                   \
+                register __m128d b_lo = _mm_load_pd(right + (2 * i + 0)); \
+                register __m128d b_hi = _mm_load_pd(right + (2 * i + 1)); \
+                register __m128d a_lo = _mm_load_pd(left + (0 * 2 + 0));  \
+                register __m128d a_hi = _mm_load_pd(left + (0 * 2 + 1));  \
+                register __m128d b = _mm_shuffle_pd(b_lo, b_lo, 0x0);         \
+                c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b));       \
+                c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));       \
+                                                                              \
+                /* group = 1 */                                               \
+                a_lo = _mm_load_pd(left + (1 * 2 + 0));                   \
+                a_hi = _mm_load_pd(left + (1 * 2 + 1));                   \
+                b = _mm_shuffle_pd(b_lo, b_lo, 0x3);                          \
+                c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b));       \
+                c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));       \
+                                                                              \
+                /* groupp = 2 */                                              \
+                a_lo = _mm_load_pd(left + (2 * 2 + 0));                   \
+                a_hi = _mm_load_pd(left + (2 * 2 + 1));                   \
+                b = _mm_shuffle_pd(b_hi, b_hi, 0x0);                          \
+                c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b));       \
+                c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));       \
+                                                                              \
+                /* froup = 3 */                                               \
+                a_lo = _mm_load_pd(left + (3 * 2 + 0));                   \
+                a_hi = _mm_load_pd(left + (3 * 2 + 1));                   \
+                b = _mm_shuffle_pd(b_hi, b_hi, 0x3);                          \
+                c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b));       \
+                c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));       \
+        } while (0)
+
+// clang-format on
+
+SN_SK_PROC_DECL_STD(double, sse, , fma)
+{
+        register __m128d c0_lo = _mm_load_pd(ans->f64_x + (0 * 2 + 0));
+        register __m128d c0_hi = _mm_load_pd(ans->f64_x + (0 * 2 + 1));
+        register __m128d c1_lo = _mm_load_pd(ans->f64_x + (1 * 2 + 0));
+        register __m128d c1_hi = _mm_load_pd(ans->f64_x + (1 * 2 + 1));
+        register __m128d c2_lo = _mm_load_pd(ans->f64_x + (2 * 2 + 0));
+        register __m128d c2_hi = _mm_load_pd(ans->f64_x + (2 * 2 + 1));
+        register __m128d c3_lo = _mm_load_pd(ans->f64_x + (3 * 2 + 0));
+        register __m128d c3_hi = _mm_load_pd(ans->f64_x + (3 * 2 + 1));
+
+        DOUBLE_KERNEL_FMA_CYCLE(0);
+        DOUBLE_KERNEL_FMA_CYCLE(1);
+        DOUBLE_KERNEL_FMA_CYCLE(2);
+        DOUBLE_KERNEL_FMA_CYCLE(3);
+
+        _mm_store_pd(ans->f64_x + (0 * 2 + 0), c0_lo);
+        _mm_store_pd(ans->f64_x + (0 * 2 + 1), c0_hi);
+        _mm_store_pd(ans->f64_x + (1 * 2 + 0), c1_lo);
+        _mm_store_pd(ans->f64_x + (1 * 2 + 1), c1_hi);
+        _mm_store_pd(ans->f64_x + (2 * 2 + 0), c2_lo);
+        _mm_store_pd(ans->f64_x + (2 * 2 + 1), c2_hi);
+        _mm_store_pd(ans->f64_x + (3 * 2 + 0), c3_lo);
+        _mm_store_pd(ans->f64_x + (3 * 2 + 1), c3_hi);
+}
+
+// clang-format off
+SN_SK_PROC_DECL_STD(double, avx, , fma)
+{
+        register __m256d c0 = _mm256_load_pd((double *)&ans->f64_y[0]);
+        register __m256d c1 = _mm256_load_pd((double *)&ans->f64_y[1]);
+        register __m256d c2 = _mm256_load_pd((double *)&ans->f64_y[2]);
+        register __m256d c3 = _mm256_load_pd((double *)&ans->f64_y[3]);
+        register __m256d bcol;
+
+        const __m256d a0 = _mm256_load_pd((double *)&left->f64_y[0]);
+        const __m256d a1 = _mm256_load_pd((double *)&left->f64_y[1]);
+        const __m256d a2 = _mm256_load_pd((double *)&left->f64_y[2]);
+        const __m256d a3 = _mm256_load_pd((double *)&left->f64_y[3]);
+
+        /* Column 0 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[0]);
+        c0 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c0);
+        c0 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c0);
+        c0 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c0);
+        c0 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c0);
+
+        /* Column 1 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[1]);
+        c1 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c1);
+        c1 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c1);
+        c1 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c1);
+        c1 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c1);
+
+        /* Column 2 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[2]);
+        c2 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c2);
+        c2 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c2);
+        c2 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c2);
+        c2 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c2);
+
+        /* Column 3 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[3]);
+        c3 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c3);
+        c3 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c3);
+        c3 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c3);
+        c3 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c3);
+
+        _mm256_store_pd((double *)&ans->f64_y[0], c0);
+        _mm256_store_pd((double *)&ans->f64_y[1], c1);
+        _mm256_store_pd((double *)&ans->f64_y[2], c2);
+        _mm256_store_pd((double *)&ans->f64_y[3], c3);
+}
+// clang-format on
+
+SN_SK_PROC_DECL_STD(double, avx, _fma, fma)
+{
+        register __m256d c0 = _mm256_load_pd((double *)&ans->f64_y[0]);
+        register __m256d c1 = _mm256_load_pd((double *)&ans->f64_y[1]);
+        register __m256d c2 = _mm256_load_pd((double *)&ans->f64_y[2]);
+        register __m256d c3 = _mm256_load_pd((double *)&ans->f64_y[3]);
+        register __m256d bcol;
+
+        const __m256d a0 = _mm256_load_pd((double *)&left->f64_y[0]);
+        const __m256d a1 = _mm256_load_pd((double *)&left->f64_y[1]);
+        const __m256d a2 = _mm256_load_pd((double *)&left->f64_y[2]);
+        const __m256d a3 = _mm256_load_pd((double *)&left->f64_y[3]);
+
+        /* Column 0 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[0]);
+        c0 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c0);
+        c0 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c0);
+        c0 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c0);
+        c0 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c0);
+
+        /* Column 1 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[1]);
+        c1 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c1);
+        c1 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c1);
+        c1 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c1);
+        c1 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c1);
+
+        /* Column 2 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[2]);
+        c2 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c2);
+        c2 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c2);
+        c2 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c2);
+        c2 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c2);
+
+        /* Column 3 */
+        bcol = _mm256_load_pd((double *)&right->f64_y[3]);
+        c3 = _mm256_fmadd_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)), c3);
+        c3 = _mm256_fmadd_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)), c3);
+        c3 = _mm256_fmadd_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)), c3);
+        c3 = _mm256_fmadd_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)), c3);
+
+        _mm256_store_pd((double *)&ans->f64_y[0], c0);
+        _mm256_store_pd((double *)&ans->f64_y[1], c1);
+        _mm256_store_pd((double *)&ans->f64_y[2], c2);
+        _mm256_store_pd((double *)&ans->f64_y[3], c3);
+}
+
+SN_SK_PROC_DECL_STD(double, 512, , fma)
+{
+        register __m512d c01 = _mm512_load_pd((double *)&ans->f64_z[0]);
+        register __m512d c23 = _mm512_load_pd((double *)&ans->f64_z[1]);
+
+        register __m512d a01 = _mm512_load_pd((double *)&left->f64_z[0]);
+        register __m512d a23 = _mm512_load_pd((double *)&left->f64_z[1]);
+
+        register __m512d b0 = _mm512_broadcast_f64x4(right->f64_y[0]);
+        c01 = _mm512_fmadd_pd(a01, b0, c01);
+        c23 = _mm512_fmadd_pd(a23, b0, c23);
+
+        register __m512d b1 = _mm512_broadcast_f64x4(right->f64_y[1]);
+        c01 = _mm512_fmadd_pd(a01, b1, c01);
+        c23 = _mm512_fmadd_pd(a23, b1, c23);
+
+        register __m512d b2 = _mm512_broadcast_f64x4(right->f64_y[2]);
+        c01 = _mm512_fmadd_pd(a01, b2, c01);
+        c23 = _mm512_fmadd_pd(a23, b2, c23);
+
+        register __m512d b3 = _mm512_broadcast_f64x4(right->f64_y[3]);
+        c01 = _mm512_fmadd_pd(a01, b3, c01);
+        c23 = _mm512_fmadd_pd(a23, b3, c23);
+
+        _mm512_store_pd((double *)&ans->f64_z[0], c01);
+        _mm512_store_pd((double *)&ans->f64_z[1], c23);
+}
+#endif