rewrite of f32 fma, and other optmization of simd kernel ops

StgSharp.Common/Collections/CapacityFixedStack.cs

@@ -91,7 +91,7 @@ namespace StgSharp.Collections
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public ReadOnlySpan<T> PopRange(int count)
+        public T[] PopRange(int count)
         {
             int i = int.Min(count, _index) - 1;
             return _values[0..i];

StgSharp.Common/HighPerformance/IRegisterType.cs

@@ -34,12 +34,65 @@ using System.Threading.Tasks;
 
 namespace StgSharp.HighPerformance
 {
-    public interface IRegisterType
+    /// <summary>
+    ///   Abstraction for a wide SIMD register or equivalent packed numeric storage. Provides raw
+    ///   reference access, broadcast initialization and element level access. Typical
+    ///   implementations map to 128/256/512-bit hardware registers or software emulated blocks.
+    /// </summary>
+
+    public interface IRegisterPresentation
     {
 
+        /// <summary>
+        ///   Gets a reference to the first element of the register reinterpreted as type <typeparamref
+        ///   name="T" />.
+        /// </summary>
+        /// <typeparam name="T">
+        ///   Unmanaged numeric element type implementing <see cref="INumber{T}" />.
+        /// </typeparam>
+        /// <returns>
+        ///   Writable reference to element 0 (lowest lane / lowest bits).
+        /// </returns>
+        /// <remarks>
+        ///   Intended for treating the register as a contiguous block of elements. The returned reference
+        ///   must remain valid for the lifetime of the register and be identical to <see
+        ///   cref="Member{T}(int)" /> with index 0. Caller must choose a type layout compatible with the
+        ///   actual register representation.
+        /// </remarks>
         ref T AsRef<T>() where T: unmanaged, INumber<T>;
 
+        /// <summary>
+        ///   Fills all elements (lanes) of the register with the provided scalar <paramref name="value" />.
+        /// </summary>
+        /// <typeparam name="T">
+        ///   Unmanaged numeric element type implementing <see cref="INumber{T}" />.
+        /// </typeparam>
+        /// <param name="value">
+        ///   Scalar value to replicate to every element.
+        /// </param>
+        /// <remarks>
+        ///   Should map to a single broadcast/shuffle instruction where possible. After completion all
+        ///   lanes read back the same value.
+        /// </remarks>
+        void BroadCastFrom<T>(T value) where T: unmanaged, INumber<T>;
+
+        /// <summary>
+        ///   Gets a writable reference to the element at the specified <paramref name="index" />.
+        /// </summary>
+        /// <typeparam name="T">
+        ///   Unmanaged numeric element type implementing <see cref="INumber{T}" />.
+        /// </typeparam>
+        /// <param name="index">
+        ///   Zero-based lane index (range defined by implementation).
+        /// </param>
+        /// <returns>
+        ///   Writable reference to the lane.
+        /// </returns>
+        /// <remarks>
+        ///   Enables scalar lane edits without unpacking the entire register. Implementations may omit
+        ///   bounds checks for performance. Caller is responsible for passing a valid index.
+        /// </remarks>
         ref T Member<T>(int index) where T: unmanaged, INumber<T>;
 
     }
-}
+}

StgSharp.Common/HighPerformance/M128.cs

@@ -68,6 +68,11 @@ namespace StgSharp.HighPerformance
             return Vec.AsVector4();
         }
 
+        public void BroadCastFrom<T>(T value) where T: unmanaged, INumber<T>
+        {
+            Vec = Vector128.CreateScalarUnsafe(value).As<T, float>();
+        }
+
         public override readonly bool Equals([NotNullWhen(true)] object obj)
         {
             return obj is M128 other && this == other;

StgSharp.Common/HighPerformance/M256.cs

@@ -31,6 +31,7 @@ using System.Linq;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
 using System.Text;
 using System.Threading.Tasks;
 
@@ -41,6 +42,7 @@ namespace StgSharp.HighPerformance
     {
 
         [FieldOffset(0)] internal fixed byte Buffer[32];
+        [FieldOffset(0)] internal Vector256<float> Vec;
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public ref T AsRef<T>() where T: unmanaged, INumber<T>
@@ -48,6 +50,11 @@ namespace StgSharp.HighPerformance
             return ref Unsafe.As<byte, T>(ref Buffer[0]);
         }
 
+        public void BroadCastFrom<T>(T value) where T: unmanaged, INumber<T>
+        {
+            Vec = Vector256.CreateScalarUnsafe(value).As<T, float>();
+        }
+
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public ref T Member<T>(int index) where T: unmanaged, INumber<T>
         {

StgSharp.Common/HighPerformance/M512.cs

@@ -31,6 +31,7 @@ using System.Linq;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
 using System.Text;
 using System.Threading.Tasks;
 
@@ -41,12 +42,18 @@ namespace StgSharp.HighPerformance
     {
 
         [FieldOffset(0)]private fixed byte _buffer[64];
+        [FieldOffset(0)]private Vector512<float> Vec;
 
         public ref T AsRef<T>() where T: unmanaged, INumber<T>
         {
             return ref Unsafe.As<byte, T>(ref _buffer[0]);
         }
 
+        public void BroadCastFrom<T>(T value) where T: unmanaged, INumber<T>
+        {
+            Vec = Vector512.CreateScalarUnsafe(value).As<T, float>();
+        }
+
         public ref T Member<T>(int index) where T: unmanaged, INumber<T>
         {
             return ref Unsafe.As<byte, T>(ref _buffer[0]);

StgSharp.Common/HighPerformance/ScalarPacket.cs

@@ -1,6 +1,6 @@
 //-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
-// file="IRegisterPresentation"
+// file="ScalarPacket"
 // Project: StgSharp
 // AuthorGroup: Nitload
 // Copyright (c) Nitload. All rights reserved.
@@ -29,17 +29,36 @@ using System;
 using System.Collections.Generic;
 using System.Linq;
 using System.Numerics;
+using System.Runtime.CompilerServices;
 using System.Text;
 using System.Threading.Tasks;
 
 namespace StgSharp.HighPerformance
 {
-    public interface IRegisterPresentation
+    public unsafe struct ScalarPacket
     {
 
-        public ref T AsRef<T>() where T: unmanaged, INumber<T>;
+        private fixed ulong data[8];
 
-        public ref T Member<T>(int index) where T: unmanaged, INumber<T>;
+        public void BroadCast<TScalar, TRegister>(int slotIndex, TScalar value)
+            where TScalar: unmanaged, INumber<TScalar>
+            where TRegister: unmanaged, IRegisterPresentation
+        {
+            if (slotIndex is < 0 or >= 8) {
+                throw new IndexOutOfRangeException();
+            }
+            ref TRegister register = ref Unsafe.As<ulong, TRegister>(ref data[slotIndex]);
+            register.BroadCastFrom(value);
+        }
+
+        public ref T Data<T>(int index) where T: unmanaged, INumber<T>
+        {
+            if (index is < 0 or >= 8) {
+                throw new IndexOutOfRangeException();
+            }
+
+            return ref Unsafe.As<ulong, T>(ref data[index]);
+        }
 
     }
 }

StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs

@@ -26,6 +26,7 @@
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
 using StgSharp.HighPerformance;
+using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
@@ -38,8 +39,9 @@ namespace StgSharp.Internal.Intrinsic
 
         public delegate* unmanaged[Cdecl]<char*, int, int> city_hash_simplify;
         public delegate* unmanaged[Cdecl]<void*, void*, void*, void> f32_add;
+        public delegate* unmanaged[Cdecl]<void*, ScalarPacket*, void> f32_fill;
         public delegate* unmanaged[Cdecl]<void*, void*, void*, void> f32_fma;
-        public delegate* unmanaged[Cdecl]<void*, float, void*, void> f32_scalar_mul;
+        public delegate* unmanaged[Cdecl]<void*, void*, ScalarPacket*, void> f32_scalar_mul;
         public delegate* unmanaged[Cdecl]<void*, void*, void*, void> f32_sub;
         public delegate* unmanaged[Cdecl]<void*, void*, void> f32_transpose;
         public delegate* unmanaged[Cdecl]<int, ulong> factorial_simd;
@@ -52,4 +54,4 @@ namespace StgSharp.Internal.Intrinsic
         }
 
     }
-}
+}

StgSharp.Mathematics/Numeric/MatrixCompute.Basic.cs

@@ -31,6 +31,7 @@ using StgSharp.Internal;
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Numerics;
 using System.Security.AccessControl;
 using System.Text;
 using System.Threading.Tasks;
@@ -59,23 +60,21 @@ namespace StgSharp.Mathematics.Numeric
             }
             if (count <= 1024)
             {
-                ReadOnlySpan<MatrixParallelThread> pool = MatrixParallel.LeadParallel(1);
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
-                if (pool.Length == 0)
+
-                {
-                    /*
-                     * wait
-                     */
-                }
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
                 for (int i = 0; i < count; i++) {
                     NativeIntrinsic.Intrinsic.f32_add(leftPtr + i, rightPtr + i, ansPtr + i);
                 }
+                MatrixParallel.ReturnParallelResources(ref pool);
                 return ans;
             }
-            MatrixParallelHandle taskGroup = MatrixParallel.PublicTask(left.GetColumnEnumeration(), right.GetColumnEnumeration(), ans.GetColumnEnumeration(),
+            MatrixParallelHandle taskGroup = MatrixParallel.PublicTask(left.GetColumnEnumeration(),
-                (delegate*<MatrixKernel<float>*, MatrixKernel<float>*, MatrixKernel<float>*, void>)NativeIntrinsic.Intrinsic.f32_add);
+                                                                       right.GetColumnEnumeration(),
+                                                                       ans.GetColumnEnumeration(),
+                                                                       (delegate*<MatrixKernel<float>*, MatrixKernel<float>*, MatrixKernel<float>*, void>)NativeIntrinsic.Intrinsic.f32_add);
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
             /*
@@ -85,5 +84,37 @@ namespace StgSharp.Mathematics.Numeric
             return ans;
         }
 
+        public static unsafe void Fill<T>(Matrix<float> left, T value) where T: unmanaged, INumber<T>
+        {
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            ScalarPacket* scaler = MatrixParallelFactory.CreateScalarPacket();
+            scaler->BroadCast<T, M512>(0, value);
+            if (count <= 256)
+            {
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                for (int i = 0; i < count; i++) {
+                    NativeIntrinsic.Intrinsic.f32_fill(leftPtr + i, scaler);
+                }
+            }
+            if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                for (int i = 0; i < count; i++) {
+                    NativeIntrinsic.Intrinsic.f32_fill(leftPtr + i, scaler);
+                }
+                MatrixParallel.ReturnParallelResources(ref pool);
+            }
+            MatrixParallelHandle taskGroup = MatrixParallel.PublicTask(left.GetColumnEnumeration(),
+                                                                       (delegate*<MatrixKernel<float>*, MatrixKernel<float>*, MatrixKernel<float>*, void>)NativeIntrinsic.Intrinsic.f32_add);
+
+            MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
+            /*
+             * run parallel
+             */
+            taskGroup.Dispose();
+        }
+
     }
 }

StgSharp.Mathematics/Numeric/MatrixParallel.Await.cs

@@ -2,8 +2,8 @@
 // -----------------------------------------------------------------------
 // file="MatrixParallel.Await"
 // 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
@@ -56,7 +56,16 @@ namespace StgSharp.Mathematics.Numeric
                     min,
                     recommendCount
                 );
-                _awaitingLeaderHandle.Enqueue(handle);
+                do
+                {
+                    if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                    {
+                        _awaitingLeaderHandle.Enqueue(handle);
+                        _ = Interlocked.Exchange(ref _sortState, 0);
+                        Interlocked.Increment(ref _createGeneration);
+                        break;
+                    }
+                } while (true);
                 int cycleGeneration = _createGeneration;
                 do
                 {
@@ -66,7 +75,7 @@ namespace StgSharp.Mathematics.Numeric
                     {
                         if (Interlocked.Exchange(ref _sortState, 1) == 0)
                         {
-                            if (_awaitingLeaderHandle.TryPopIfRefEqual(cycleGeneration, ref handle))
+                            if (_awaitingLeaderHandle.TryPopIfRefEqual(cycleGeneration, in handle))
                             {
                                 // Sort is implemented in TryPopIfRefEqual
                                 ret = _threads.PopRange(handle.BestCount);
@@ -101,7 +110,16 @@ namespace StgSharp.Mathematics.Numeric
 
         public static void ReturnParallelResources(ref MatrixParallelThread[] pool)
         {
-            _threads.PushRange(pool);
+            while (true)
+            {
+                if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                {
+                    _threads.PushRange(pool);
+                    Interlocked.Exchange(ref _sortState, 0);
+                    break;
+                }
+            }
+
             pool = null!;
         }
 

StgSharp.Mathematics/Numeric/MatrixParallel.AwaitingQueue.cs

@@ -2,8 +2,8 @@
 // -----------------------------------------------------------------------
 // file="MatrixParallel.AwaitingQueue"
 // 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
@@ -47,7 +47,6 @@ namespace StgSharp.Mathematics.Numeric
 
             private Element _quickElements = new();
             private volatile int _count;
-            private volatile int _lastGeneration;
             private readonly Queue<WaitHandle> _fallback = [];
 
             public AwaitingQueue() { }
@@ -68,22 +67,16 @@ namespace StgSharp.Mathematics.Numeric
             {
                 if (_count == 0)
                 {
-                    _lastGeneration = generation;
                     handle = null!;
                     return false;
                 }
-                if (_lastGeneration == generation)
+                for (int i = 0; i < _count; i++) {
-                {
+                    _quickElements[i].UpdatePriority(generation);
-                    for (int i = 0; i < _count; i++) {
-                        _quickElements[i].UpdatePriority(generation);
-                    }
-                    Element.Sort(_quickElements);
-                } else
-                {
-                    _lastGeneration = generation;
                 }
+                Element.Sort(ref _quickElements);
                 handle = _quickElements[0];
                 _count--;
+                _quickElements[0] = default!;
                 if (_count > 7) {
                     _ = _fallback.TryDequeue(out _quickElements[0]!);
                 }
@@ -94,40 +87,31 @@ namespace StgSharp.Mathematics.Numeric
             {
                 if (_count == 0)
                 {
-                    _lastGeneration = generation;
                     handle = null!;
                     return false;
                 }
                 for (int i = 0; i < _count; i++) {
                     _quickElements[i].UpdatePriority(generation);
                 }
-                Element.Sort(_quickElements);
+
-                _lastGeneration = generation;
+                Element.Sort(ref _quickElements);
                 handle = _quickElements[0];
-                _count--;
-                if (_count > 7) {
-                    _ = _fallback.TryDequeue(out _quickElements[0]!);
-                }
                 return true;
             }
 
-            public bool TryPopIfRefEqual(int generation, ref WaitHandle handle)
+            public bool TryPopIfRefEqual(int generation, in WaitHandle handle)
             {
-                if (_count == 0)
+                if (_count == 0) {
-                {
-                    handle = null!;
                     return false;
                 }
                 for (int i = 0; i < _count; i++) {
                     _quickElements[i].UpdatePriority(generation);
                 }
                 Element.Sort(ref _quickElements);
-                _lastGeneration = generation;
-                handle = _quickElements[0];
-                _count--;
                 if (!ReferenceEquals(_quickElements[0], handle)) {
                     return false;
                 }
+                _count--;
                 _quickElements[0] = null!;
                 if (_count > 7) {
                     _ = _fallback.TryDequeue(out _quickElements[0]!);
@@ -147,7 +131,7 @@ namespace StgSharp.Mathematics.Numeric
             {
                 int maxIdx = 0;
                 WaitHandle first = e[0];
-                int maxPr = first?.Priority ?? int.MaxValue;
+                int maxPr = first?.Priority ?? int.MinValue;
 
                 for (int i = 1; i < 8; i++)
                 {

StgSharp.Mathematics/Numeric/MatrixParallel.cs

@@ -2,8 +2,8 @@
 // -----------------------------------------------------------------------
 // file="MatrixParallel"
 // 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
@@ -40,7 +40,7 @@ namespace StgSharp.Mathematics.Numeric
 
         private static bool _meaningful;
 
-        private static ConcurrentCapacityFixedStack<MatrixParallelThread> _threads;
+        private static CapacityFixedStack<MatrixParallelThread> _threads;
 
         private static volatile int RemainThreadCount = 0;
         private static readonly object _lock = new();

StgSharp.Mathematics/Numeric/MatrixParallelHandle.cs

@@ -1,6 +1,6 @@
 //-----------------------------------------------------------------------
 // -----------------------------------------------------------------------
-// file="MatrixParallelTask"
+// file="MatrixParallelHandle"
 // Project: StgSharp
 // AuthorGroup: Nitload
 // Copyright (c) Nitload. All rights reserved.
@@ -34,22 +34,6 @@ using System.Runtime.InteropServices;
 
 namespace StgSharp.Mathematics.Numeric
 {
-    public unsafe struct ScalarPacket
-    {
-
-        private fixed ulong data[8];
-
-        public ref T Data<T>(int index) where T: unmanaged,INumber<T>
-        {
-            if (index is < 0 or >= 8) {
-                throw new IndexOutOfRangeException();
-            }
-
-            return ref Unsafe.As<ulong, T>(ref data[index]);
-        }
-
-    }
-
     [StructLayout(LayoutKind.Sequential)]
     internal unsafe struct MatrixParallelTaskPackageNonGeneric
     {
@@ -59,11 +43,8 @@ namespace StgSharp.Mathematics.Numeric
                                   MatrixParallelTaskPackageNonGeneric* source,
                                   MatrixParallelTaskPackageNonGeneric* target)
         {
-            Buffer.MemoryCopy(
+            Buffer.MemoryCopy(source, target, sizeof(MatrixParallelTaskPackageNonGeneric),
-                source,
+                              sizeof(MatrixParallelTaskPackageNonGeneric));
-                target,
-                sizeof(MatrixParallelTaskPackageNonGeneric),
-                sizeof(MatrixParallelTaskPackageNonGeneric));
         }
 
         #region source
@@ -148,11 +129,8 @@ namespace StgSharp.Mathematics.Numeric
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static unsafe void Copy(MatrixParallelTaskPackage<T>* source, MatrixParallelTaskPackage<T>* target)
         {
-            Buffer.MemoryCopy(
+            Buffer.MemoryCopy(source, target, sizeof(MatrixParallelTaskPackage<T>),
-                source,
+                              sizeof(MatrixParallelTaskPackage<T>));
-                target,
-                sizeof(MatrixParallelTaskPackage<T>),
-                sizeof(MatrixParallelTaskPackage<T>));
         }
 
         #region source

StgSharp.Native/StgSharp.Native.vcxproj

@@ -109,7 +109,9 @@
     <ClCompile Include="sn-src\str_factorial.c" />
     <ClCompile Include="sn-src\k_fma.c" />
     <ClCompile Include="sn-src\str_hash.c" />
-    <ClCompile Include="sn-src\load_intrinsic.c" />
+    <ClCompile Include="sn-src\load_intrinsic.c">
+      <FileType>CppCode</FileType>
+    </ClCompile>
     <ClCompile Include="sn-src\load_glfw.c" />
     <ClCompile Include="sn-src\stb_image.c" />
     <ClCompile Include="sn-src\sn_native.c" />

StgSharp.Native/StgSharp.Native.vcxproj.filters

@@ -24,9 +24,6 @@
     <ClCompile Include="sn-src\str_hash.c">
       <Filter>源文件</Filter>
     </ClCompile>
-    <ClCompile Include="sn-src\load_intrinsic.c">
-      <Filter>源文件</Filter>
-    </ClCompile>
     <ClCompile Include="sn-src\stb_image.c">
       <Filter>源文件</Filter>
     </ClCompile>
@@ -51,6 +48,9 @@
     <ClCompile Include="sn-src\k_fma.c">
       <Filter>源文件</Filter>
     </ClCompile>
+    <ClCompile Include="sn-src\load_intrinsic.c">
+      <Filter>源文件</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="sn-head\gl.h">

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -31,12 +31,7 @@ typedef union uint64_u {
         uint32_t u32[2];
 } uint64_u;
 
-typedef union mat_kernel {
+typedef struct mat_kernel mat_kernel;
-        __m128 xmm[4];
-        __m256 ymm[2];
-        __m512 zmm[1];
-        float m[4][4];
-} mat_kernel;
 
 typedef union scalar_pack {
         uint64_t data[8];
@@ -44,13 +39,10 @@ typedef union scalar_pack {
 
 typedef void(SN_DECL *VECTORNORMALIZEPROC)(__m128 *source, __m128 *target);
 typedef void(SN_DECL *DOTPROC)(void *transpose, __m128 *vector, __m128 *ans);
-typedef void(SN_DECL *MATRIXTRANSPOSEPROC)(void *source, void *target);
+typedef void(SN_DECL *MAT_PROC_AS)(mat_kernel *ans, scalar_pack *scalar);
-typedef void(SN_DECL *MATRIXDETPROC)(void *mat, void *transpose);
+typedef void(SN_DECL *MAT_PROC_LA)(mat_kernel *left, mat_kernel *ans);
-typedef void(SN_DECL *MATARITHMETIC)(void *left, void *right, void *ans);
+typedef void(SN_DECL *MAT_PROC_LRA)(mat_kernel *left, mat_kernel *right, mat_kernel *ans);
-typedef void(SN_DECL *SCALERMATARITHMETIC)(void *left, void *right, void *ans);
+typedef void(SN_DECL *MAT_PROC_LSA)(mat_kernel *matrix, scalar_pack *scalar, mat_kernel *ans);
-typedef void(SN_DECL *MATFMAPROC)(mat_kernel *left, mat_kernel *right, mat_kernel *ans);
-typedef void(SN_DECL *SCALARMATARITHMETIC)(mat_kernel *matrix, float scalar, mat_kernel *ans);
-typedef void(SN_DECL *MATTRANSPOSEPROC)(mat_kernel *source, mat_kernel *target);
 typedef uint64_t(SN_DECL *FACTORIALROC)(int n);
 
 #pragma region matix function
@@ -75,12 +67,18 @@ INTERNAL void SN_DECL f32_transpose_avx(MAT_KERNEL(float) const *source,
                                         MAT_KERNEL(float) *restrict target);
 INTERNAL void SN_DECL f32_transpose_512(MAT_KERNEL(float) const *source,
                                         MAT_KERNEL(float) *restrict target);
-INTERNAL void SN_DECL f32_scalar_mul_sse(MAT_KERNEL(float) const *matrix, float const scalar,
+INTERNAL void SN_DECL f32_scalar_mul_sse(MAT_KERNEL(float) const *matrix,
-                                         MAT_KERNEL(float) *restrict ans);
+                                         MAT_KERNEL(float) *restrict ans,
-INTERNAL void SN_DECL f32_scalar_mul_avx(MAT_KERNEL(float) const *matrix, float const scalar,
+                                         scalar_pack const *scalar);
-                                         MAT_KERNEL(float) *restrict ans);
+INTERNAL void SN_DECL f32_scalar_mul_avx(MAT_KERNEL(float) const *matrix,
-INTERNAL void SN_DECL f32_scalar_mul_512(MAT_KERNEL(float) const *matrix, float const scalar,
+                                         MAT_KERNEL(float) *restrict ans,
-                                         MAT_KERNEL(float) *restrict ans);
+                                         scalar_pack const *scalar);
+INTERNAL void SN_DECL f32_scalar_mul_512(MAT_KERNEL(float) const *matrix,
+                                         MAT_KERNEL(float) *restrict ans,
+                                         scalar_pack const *scalar);
+INTERNAL void SN_DECL f32_fill_sse(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler);
+INTERNAL void SN_DECL f32_fill_avx(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler);
+INTERNAL void SN_DECL f32_fill_512(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler);
 #pragma endregion
 
 #pragma region fma
@@ -88,6 +86,12 @@ INTERNAL void SN_DECL f32_scalar_mul_512(MAT_KERNEL(float) const *matrix, float
 INTERNAL void SN_DECL f32_fma_sse(MAT_KERNEL(float) const *restrict left,
                                   MAT_KERNEL(float) const *restrict right,
                                   MAT_KERNEL(float) *restrict ans);
+INTERNAL void SN_DECL f32_fma_avx(MAT_KERNEL(float) const *restrict left,
+                                  MAT_KERNEL(float) const *restrict right,
+                                  MAT_KERNEL(float) *restrict ans);
+INTERNAL void SN_DECL f32_fma_512(MAT_KERNEL(float) const *restrict left,
+                                  MAT_KERNEL(float) const *restrict right,
+                                  MAT_KERNEL(float) *restrict ans);
 
 #pragma endregion
 
@@ -119,11 +123,12 @@ typedef struct sn_intrinsic {
         int (*city_hash_simplify)(byte *str, int length);
         FACTORIALROC factorial_simd;
         int (*index_pair)(short const *str, uint32_t target, int length);
-        MATARITHMETIC f32_add;
+        MAT_PROC_LRA f32_add;
-        void (*f32_fma)(void *left, void *right, void *ans);
+        MAT_PROC_AS f32_fill;
-        SCALARMATARITHMETIC f32_scalar_mul;
+        MAT_PROC_LRA f32_fma;
-        MATARITHMETIC f32_sub;
+        MAT_PROC_LRA f32_scalar_mul;
-        MATTRANSPOSEPROC f32_transpose;
+        MAT_PROC_LRA f32_sub;
+        MAT_PROC_LA f32_transpose;
         VECTORNORMALIZEPROC f32_normalize_3;
 } sn_intrinsic;
 

StgSharp.Native/sn-src/k_basicop.c

@@ -1,5 +1,27 @@
 #include "sn_intrinsic.h"
 
+INTERNAL void SN_DECL f32_fill_sse(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler)
+{
+        register __m128 vec = _mm_set_ps1(((float *)(scaler->data))[0]);
+        dest->xmm[0] = vec;
+        dest->xmm[1] = vec;
+        dest->xmm[2] = vec;
+        dest->xmm[3] = vec;
+}
+
+INTERNAL void SN_DECL f32_fill_avx(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler)
+{
+        register __m256 vec = _mm256_broadcast_ss(((float *)(scaler->data)));
+        dest->ymm[0] = vec;
+        dest->ymm[1] = vec;
+}
+
+INTERNAL void SN_DECL f32_fill_512(MAT_KERNEL(float) *restrict dest, scalar_pack *scaler)
+{
+        register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scaler->data))[0]);
+        dest->zmm[0] = vec;
+}
+
 INTERNAL void SN_DECL f32_add_sse(MAT_KERNEL(float) const *left, MAT_KERNEL(float) const *right,
                                   MAT_KERNEL(float) *restrict ans)
 {
@@ -106,11 +128,10 @@ INTERNAL void SN_DECL f32_transpose_512(MAT_KERNEL(float) const *source,
         target->zmm[0] = _mm512_permutexvar_ps(transpose_indices, matrix);
 }
 
-INTERNAL void SN_DECL f32_scalar_mul_sse(MAT_KERNEL(float) const *matrix, float const scalar,
+INTERNAL void SN_DECL f32_scalar_mul_sse(MAT_KERNEL(float) const *matrix,
-                                         MAT_KERNEL(float) *restrict ans)
+                                         MAT_KERNEL(float) *restrict ans, scalar_pack const *scalar)
 {
-        register __m128 temp = _mm_load_ss(&scalar);
+        register __m128 scalar_vec = ((__m128 *)(scalar->data))[0];
-        register __m128 scalar_vec = _mm_shuffle_ps(temp, temp, _MM_SHUFFLE(0, 0, 0, 0));
 
         register __m128 c0 = _mm_mul_ps(matrix->xmm[0], scalar_vec);
         register __m128 c1 = _mm_mul_ps(matrix->xmm[1], scalar_vec);
@@ -123,10 +144,10 @@ INTERNAL void SN_DECL f32_scalar_mul_sse(MAT_KERNEL(float) const *matrix, float
         ans->xmm[3] = c3;
 }
 
-INTERNAL void SN_DECL f32_scalar_mul_avx(MAT_KERNEL(float) const *matrix, float const scalar,
+INTERNAL void SN_DECL f32_scalar_mul_avx(MAT_KERNEL(float) const *matrix,
-                                         MAT_KERNEL(float) *restrict ans)
+                                         MAT_KERNEL(float) *restrict ans, scalar_pack const *scalar)
 {
-        register __m256 scalar_vec = _mm256_broadcast_ss(&scalar);
+        register __m256 scalar_vec = ((__m256 *)(scalar->data))[0];
 
         register __m256 c0 = _mm256_mul_ps(matrix->ymm[0], scalar_vec);
         register __m256 c1 = _mm256_mul_ps(matrix->ymm[1], scalar_vec);
@@ -135,10 +156,10 @@ INTERNAL void SN_DECL f32_scalar_mul_avx(MAT_KERNEL(float) const *matrix, float
         ans->ymm[1] = c1;
 }
 
-INTERNAL void SN_DECL f32_scalar_mul_512(MAT_KERNEL(float) const *matrix, float const scalar,
+INTERNAL void SN_DECL f32_scalar_mul_512(MAT_KERNEL(float) const *matrix,
-                                         MAT_KERNEL(float) *restrict ans)
+                                         MAT_KERNEL(float) *restrict ans, scalar_pack const *scalar)
 {
-        __m512 scalar_vec = _mm512_broadcastss_ps(_mm_load_ss(&scalar));
+        register __m512 scalar_vec = ((__m512 *)(scalar->data))[0];
 
         register __m512 c0 = _mm512_mul_ps(matrix->zmm[0], scalar_vec);
 

StgSharp.Native/sn-src/k_fma.c

@@ -1,6 +1,7 @@
 #include "sn_intrinsic.h"
 
 // ans = left * right + ans
+// right is not transposed
 INTERNAL void SN_DECL f32_fma_sse(MAT_KERNEL(float) const *restrict left,
                                   MAT_KERNEL(float) const *restrict right,
                                   MAT_KERNEL(float) *restrict ans)
@@ -10,34 +11,47 @@ INTERNAL void SN_DECL f32_fma_sse(MAT_KERNEL(float) const *restrict left,
         register __m128 source2 = _mm_load_ps(&left->xmm[2]);
         register __m128 source3 = _mm_load_ps(&left->xmm[3]);
 
+        register __m128 ans0 = _mm_load_ps(&ans->xmm[0]);
+        register __m128 ans1 = _mm_load_ps(&ans->xmm[1]);
+        register __m128 ans2 = _mm_load_ps(&ans->xmm[2]);
+        register __m128 ans3 = _mm_load_ps(&ans->xmm[3]);
+
         register __m128 broadcast;
-
+        for (size_t i = 0; i < 4; i++) {
-        register __m128 ans0 = _mm_setzero_ps();
-        register __m128 ans1 = _mm_setzero_ps();
-        register __m128 ans2 = _mm_setzero_ps();
-        register __m128 ans3 = _mm_setzero_ps();
-        for (size_t i = 0; i < 3; i++) {
                 broadcast = _mm_load_ps(&right->xmm[i]);
-                _mm_add_ps(ans0, _mm_mul_ps(source0, broadcast));
-                _mm_add_ps(ans1, _mm_mul_ps(source1, broadcast));
-                _mm_add_ps(ans2, _mm_mul_ps(source2, broadcast));
-                _mm_add_ps(ans3, _mm_mul_ps(source3, broadcast));
-                for (size_t t = 0; t < 2; t++) {
-                        broadcast = _mm_shuffle_ps(broadcast, broadcast, _MM_SHUFFLE(2, 1, 0, 3));
 
-                        _mm_add_ps(ans0, _mm_mul_ps(source0, broadcast));
+                ans0 = _mm_add_ps(ans0,
-                        _mm_add_ps(ans1, _mm_mul_ps(source1, broadcast));
+                                  _mm_mul_ps(source0, _mm_shuffle_ps(broadcast, broadcast, 0x00)));
-                        _mm_add_ps(ans2, _mm_mul_ps(source2, broadcast));
+                ans1 = _mm_add_ps(ans1,
-                        _mm_add_ps(ans3, _mm_mul_ps(source3, broadcast));
+                                  _mm_mul_ps(source1, _mm_shuffle_ps(broadcast, broadcast, 0x55)));
-                }
+                ans2 = _mm_add_ps(ans2,
+                                  _mm_mul_ps(source2, _mm_shuffle_ps(broadcast, broadcast, 0xAA)));
+                ans3 = _mm_add_ps(ans3,
+                                  _mm_mul_ps(source3, _mm_shuffle_ps(broadcast, broadcast, 0xFF)));
         }
-        source0 = _mm_load_ps(&ans->xmm[0]);
-        source1 = _mm_load_ps(&ans->xmm[1]);
-        source2 = _mm_load_ps(&ans->xmm[2]);
-        source3 = _mm_load_ps(&ans->xmm[3]);
 
-        _mm_store_ps(&ans->xmm[0], _mm_add_ps(source0, ans0));
+        _mm_store_ps(&ans->xmm[0], ans0);
-        _mm_store_ps(&ans->xmm[1], _mm_add_ps(source1, ans1));
+        _mm_store_ps(&ans->xmm[1], ans1);
-        _mm_store_ps(&ans->xmm[2], _mm_add_ps(source2, ans2));
+        _mm_store_ps(&ans->xmm[2], ans2);
-        _mm_store_ps(&ans->xmm[3], _mm_add_ps(source3, ans3));
+        _mm_store_ps(&ans->xmm[3], ans3);
+}
+
+// ans = left * right + ans
+INTERNAL void SN_DECL f32_fma_avx(MAT_KERNEL(float) const *restrict left,
+                                  MAT_KERNEL(float) const *restrict right,
+                                  MAT_KERNEL(float) *restrict ans)
+{
+        /*
+        * to be implemented
+        */
+}
+
+// ans = left * right + ans
+INTERNAL void SN_DECL f32_fma_512(MAT_KERNEL(float) const *restrict left,
+                                  MAT_KERNEL(float) const *restrict right,
+                                  MAT_KERNEL(float) *restrict ans)
+{
+        /*
+        * to be implemented
+        */
 }

StgSharp.Native/sn-src/load_intrinsic.c

@@ -65,11 +65,11 @@ 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_add = (MATARITHMETIC)f32_add_sse;
+        context->f32_add = (MAT_PROC_LRA)f32_add_sse;
-        context->f32_fma = (void (*)(void *, void *, void *))f32_fma_sse;
+        context->f32_fma = (MAT_PROC_LRA)f32_fma_sse;
         context->f32_scalar_mul = f32_scalar_mul_sse;
-        context->f32_sub = (MATARITHMETIC)f32_sub_sse;
+        context->f32_sub = (MAT_PROC_LRA)f32_sub_sse;
-        context->f32_transpose = (MATTRANSPOSEPROC)f32_transpose_sse;
+        context->f32_transpose = (MAT_PROC_LA)f32_transpose_sse;
         context->f32_normalize_3 = f32_normalize;
 }
 
@@ -78,11 +78,11 @@ 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_add = (MATARITHMETIC)f32_add_avx;
+        context->f32_add = (MAT_PROC_LRA)f32_add_avx;
-        context->f32_fma = (void (*)(void *, void *, void *))f32_fma_sse;
+        context->f32_fma = (MAT_PROC_LRA)f32_fma_avx;
         context->f32_scalar_mul = f32_scalar_mul_avx;
-        context->f32_sub = (MATARITHMETIC)f32_sub_avx;
+        context->f32_sub = (MAT_PROC_LRA)f32_sub_avx;
-        context->f32_transpose = (MATTRANSPOSEPROC)f32_transpose_avx;
+        context->f32_transpose = (MAT_PROC_LA)f32_transpose_avx;
         context->f32_normalize_3 = f32_normalize;
 }
 
@@ -91,10 +91,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_add = (MATARITHMETIC)f32_add_512;
+        context->f32_add = (MAT_PROC_LRA)f32_add_512;
-        context->f32_fma = (void (*)(void *, void *, void *))f32_fma_sse;
+        context->f32_fma = (MAT_PROC_LRA)f32_fma_512;
         context->f32_scalar_mul = f32_scalar_mul_512;
-        context->f32_sub = (MATARITHMETIC)f32_sub_512;
+        context->f32_sub = (MAT_PROC_LRA)f32_sub_512;
-        context->f32_transpose = (MATTRANSPOSEPROC)f32_transpose_512;
+        context->f32_transpose = (MAT_PROC_LA)f32_transpose_512;
         context->f32_normalize_3 = f32_normalize;
 }