From 83997f64d00249cbe23f7fd63533588a08960774 Mon Sep 17 00:00:00 2001 From: Youmu Liu <105042599+Nitload-NSI@users.noreply.github.com> Date: Thu, 5 Mar 2026 18:04:34 +0800 Subject: [PATCH] Better FMA --- StgSharp.Benchmark/HlsfAllocatorBench.cs | 29 +- .../PlainGeometry/Define/PlainGeometry.cs | 10 +- StgSharp.Mathematics/Internal/MatrixPanel.cs | 4 +- .../Internal/MatrixParallelTask.cs | 4 +- .../Internal/UnsafeCompute.cs | 6 +- StgSharp.Mathematics/Moduel.cs | 80 -- .../Numeric/Graphics/Matrix44.cs | 30 +- .../Numeric/MatrixCompute.Factorize.cs | 2 +- .../Numeric/MatrixCompute.Float.cs | 60 +- .../Numeric/MatrixCompute.double.cs | 26 +- .../Numeric/NumericalModule.cs | 45 +- .../Numeric/Runtime/MatrixParallelQueue.cs | 5 +- .../Runtime/MatrixParallelThread.Fma.cs | 4 +- .../Numeric/Runtime/MatrixParallelWrap.cs | 4 +- .../Numeric/Scalar/Factorial.cs | 2 +- StgSharp.Native/StgSharp.Native.vcxproj | 4 +- .../StgSharp.Native.vcxproj.filters | 4 +- StgSharp.Native/sn-head/sn_intrinsic.std.h | 12 +- StgSharp.Native/sn-src/m.fma.x86.c | 367 ------- .../sn-src/{m.fma_tile.x86.c => mk.fma.x86.c} | 904 ++++++++---------- StgSharp.Native/sn-src/sk.fma.x86.c | 321 +++++++ 21 files changed, 893 insertions(+), 1030 deletions(-) delete mode 100644 StgSharp.Mathematics/Moduel.cs delete mode 100644 StgSharp.Native/sn-src/m.fma.x86.c rename StgSharp.Native/sn-src/{m.fma_tile.x86.c => mk.fma.x86.c} (69%) create mode 100644 StgSharp.Native/sn-src/sk.fma.x86.c diff --git a/StgSharp.Benchmark/HlsfAllocatorBench.cs b/StgSharp.Benchmark/HlsfAllocatorBench.cs index 0d23d002..2d38e93b 100644 --- a/StgSharp.Benchmark/HlsfAllocatorBench.cs +++ b/StgSharp.Benchmark/HlsfAllocatorBench.cs @@ -2,8 +2,8 @@ // ----------------------------------------------------------------------- // file="HlsfAllocatorBench" // Project: StgSharp -// AuthorGroup: Nitload Space -// Copyright (c) Nitload Space. All rights reserved. +// AuthorGroup: Nitload +// Copyright (c) Nitload. All rights reserved. // // Permission is hereby granted, free of charge, to any person obtaining a copy // of this software and associated documentation files (the "Software"), to deal @@ -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 _hlsfLocalList = new(); private readonly List _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(); diff --git a/StgSharp.Mathematics/Geometries/PlainGeometry/Define/PlainGeometry.cs b/StgSharp.Mathematics/Geometries/PlainGeometry/Define/PlainGeometry.cs index d7d57787..270341c4 100644 --- a/StgSharp.Mathematics/Geometries/PlainGeometry/Define/PlainGeometry.cs +++ b/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; diff --git a/StgSharp.Mathematics/Internal/MatrixPanel.cs b/StgSharp.Mathematics/Internal/MatrixPanel.cs index 24ad172e..87a94ecb 100644 --- a/StgSharp.Mathematics/Internal/MatrixPanel.cs +++ b/StgSharp.Mathematics/Internal/MatrixPanel.cs @@ -93,7 +93,7 @@ namespace StgSharp.Mathematics [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static int KernelSideCount() where T: unmanaged, INumber { - 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() where T: unmanaged, INumber { - SIMDID id = MathematicsModule.IntrinsicMask; + SIMDID id = NumericalModule.IntrinsicMask; switch (id.MaskByte[0] & 0b_00001111) { case 1: diff --git a/StgSharp.Mathematics/Internal/MatrixParallelTask.cs b/StgSharp.Mathematics/Internal/MatrixParallelTask.cs index 735c6102..b1ea87d3 100644 --- a/StgSharp.Mathematics/Internal/MatrixParallelTask.cs +++ b/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(36)] public int Common_j; + [FieldOffset(32)] public int CommonX; [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; diff --git a/StgSharp.Mathematics/Internal/UnsafeCompute.cs b/StgSharp.Mathematics/Internal/UnsafeCompute.cs index 6fccfaa8..95829078 100644 --- a/StgSharp.Mathematics/Internal/UnsafeCompute.cs +++ b/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); } } diff --git a/StgSharp.Mathematics/Moduel.cs b/StgSharp.Mathematics/Moduel.cs deleted file mode 100644 index bf62a3ad..00000000 --- a/StgSharp.Mathematics/Moduel.cs +++ /dev/null @@ -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(); - } - - } -} diff --git a/StgSharp.Mathematics/Numeric/Graphics/Matrix44.cs b/StgSharp.Mathematics/Numeric/Graphics/Matrix44.cs index 7ea83e07..e7f5be49 100644 --- a/StgSharp.Mathematics/Numeric/Graphics/Matrix44.cs +++ b/StgSharp.Mathematics/Numeric/Graphics/Matrix44.cs @@ -189,9 +189,9 @@ namespace StgSharp.Mathematics.Graphics { const int f32 = 0; fixed (MatrixKernel* source = &this.kernel, target = &transpose.kernel) { - MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Transpose, - null, (MatrixKernel*)source, - (MatrixKernel*)target, 0); + NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Transpose, + null, (MatrixKernel*)source, + (MatrixKernel*)target, 0); } } @@ -203,10 +203,10 @@ namespace StgSharp.Mathematics.Graphics { Matrix44 ret = new(); int f32 = MatrixElementType.F32.IntrinsicNode; - MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Sub, - (MatrixKernel*)&left, - (MatrixKernel*)&right, - (MatrixKernel*)&ret, 1); + NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Sub, + (MatrixKernel*)&left, + (MatrixKernel*)&right, + (MatrixKernel*)&ret, 1); return ret; } @@ -237,10 +237,10 @@ namespace StgSharp.Mathematics.Graphics { Matrix44 ret = new(); int f32 = MatrixElementType.F32.IntrinsicNode; - MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Fma, - (MatrixKernel*)&left, - (MatrixKernel*)&right, - (MatrixKernel*)&ret, 0); + NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Fma, + (MatrixKernel*)&left, + (MatrixKernel*)&right, + (MatrixKernel*)&ret, 0); return ret; } @@ -252,10 +252,10 @@ namespace StgSharp.Mathematics.Graphics { Matrix44 ret = new(); int f32 = MatrixElementType.F32.IntrinsicNode; - MathematicsModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Add, - (MatrixKernel*)&left, - (MatrixKernel*)&right, - (MatrixKernel*)&ret, 1); + NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Add, + (MatrixKernel*)&left, + (MatrixKernel*)&right, + (MatrixKernel*)&ret, 1); return ret; } diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs index b5d28d39..0b4d1d49 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs +++ b/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)), diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.Float.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.Float.cs index 1aff6019..68af94b1 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.Float.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.Float.cs @@ -49,7 +49,7 @@ namespace StgSharp.Mathematics.Numeric Matrix ans ) { - int f32 = (int)MatrixElementType.F32; + MatrixElementType f32 = MatrixElementType.F32; if (left.ColumnLength != right.ColumnLength || left.RowLength != right.RowLength || left.ColumnLength != ans.ColumnLength || @@ -62,8 +62,8 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_add(leftPtr, rightPtr, ansPtr, - count); + NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr, + count); return; } else if (count <= 1024) @@ -73,21 +73,21 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_add(leftPtr, rightPtr, ansPtr, - count); + 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->Right = (nint)right.Buffer; - package->Answer = (nint)ans.Buffer; + package->Mat1 = (MatrixKernel*)left.Buffer; + package->Mat2 = (MatrixKernel*)right.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(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 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,11 +168,11 @@ namespace StgSharp.Mathematics.Numeric enumeration->Member(0) = i; enumeration->Member(1) = j; - MathematicsModule.GlobalContext + NumericalModule.GlobalContext .mat_sk[f32].kernel_tile_fma((MatrixKernel*)left.Buffer, - (MatrixKernel*)right.Buffer, - (MatrixKernel*)ans.Buffer, - enumeration); + (MatrixKernel*)right.Buffer, + (MatrixKernel*)ans.Buffer, + enumeration); } } @@ -205,7 +205,7 @@ namespace StgSharp.Mathematics.Numeric Matrix 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,8 +220,8 @@ namespace StgSharp.Mathematics.Numeric { MatrixKernel* leftPtr = left.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar, - count); + NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar, + count); return; } else if (count <= 1024) @@ -230,8 +230,8 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar, - count); + NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar, + count); MatrixParallel.ReturnParallelResources(ref pool); return; @@ -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 ans ) { - int f32 = (int)MatrixElementType.F32; + MatrixElementType f32 = MatrixElementType.F32; if (left.ColumnLength != right.ColumnLength || left.RowLength != right.RowLength || left.ColumnLength != ans.ColumnLength || @@ -277,8 +277,8 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr, - count); + NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr, + count); return; } else if (count <= 1024) @@ -288,8 +288,8 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - MathematicsModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr, - count); + NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr, + count); MatrixParallel.ReturnParallelResources(ref pool); return; @@ -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); diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs index e2aea435..7c1d6a4a 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs @@ -58,7 +58,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* 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* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* 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(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(0) = i; enumeration->Member(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* leftPtr = left.Buffer; MatrixKernel* 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* leftPtr = left.Buffer; MatrixKernel* 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* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* 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* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* 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); diff --git a/StgSharp.Mathematics/Numeric/NumericalModule.cs b/StgSharp.Mathematics/Numeric/NumericalModule.cs index 0b1d6318..dff8ab4e 100644 --- a/StgSharp.Mathematics/Numeric/NumericalModule.cs +++ b/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(); + } + } } diff --git a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelQueue.cs b/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelQueue.cs index c4510594..d6526909 100644 --- a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelQueue.cs +++ b/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 _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; } diff --git a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs b/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs index e72329d1..986512f3 100644 --- a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs +++ b/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs @@ -65,7 +65,7 @@ namespace StgSharp.Mathematics.Numeric int j = (int)(c % ansHeight); enumeration->Member(0) = i; enumeration->Member(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(0) = i; enumeration->Member(1) = j; - MathematicsModule.GlobalContext.mat_sk[f64].kernel_tile_fma(leftBuffer, rightBuffer, + NumericalModule.GlobalContext.mat_sk[f64].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer, enumeration); } } diff --git a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs b/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs index 322251eb..354d0fa8 100644 --- a/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs +++ b/StgSharp.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs @@ -63,8 +63,6 @@ namespace StgSharp.Mathematics.Numeric _count = threadArray.Length; } - public MatrixParallelQueue TaskQueue { get; set; } - public ReadOnlySpan 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() diff --git a/StgSharp.Mathematics/Numeric/Scalar/Factorial.cs b/StgSharp.Mathematics/Numeric/Scalar/Factorial.cs index 7920c761..24b3bfbb 100644 --- a/StgSharp.Mathematics/Numeric/Scalar/Factorial.cs +++ b/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, diff --git a/StgSharp.Native/StgSharp.Native.vcxproj b/StgSharp.Native/StgSharp.Native.vcxproj index be1f116a..7e899e45 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj +++ b/StgSharp.Native/StgSharp.Native.vcxproj @@ -176,9 +176,9 @@ - + - + diff --git a/StgSharp.Native/StgSharp.Native.vcxproj.filters b/StgSharp.Native/StgSharp.Native.vcxproj.filters index 8c79c37c..43e5a309 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj.filters +++ b/StgSharp.Native/StgSharp.Native.vcxproj.filters @@ -42,7 +42,7 @@ 源文件 - + 源文件 @@ -57,7 +57,7 @@ 源文件 - + 源文件 diff --git a/StgSharp.Native/sn-head/sn_intrinsic.std.h b/StgSharp.Native/sn-head/sn_intrinsic.std.h index 0f981626..23c4d6bd 100644 --- a/StgSharp.Native/sn-head/sn_intrinsic.std.h +++ b/StgSharp.Native/sn-head/sn_intrinsic.std.h @@ -13,14 +13,14 @@ #define MK_PROC_NAME_STD(T_type, T_arch, T_feature, T_op) \ 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, \ - T_op)(MAT_TASK(T_type) const *task) +#define SN_MK_PROC_DECL_STD(T_type, T_arch, T_feature, T_op) \ + 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)( \ - MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \ +#define SN_SK_PROC_DECL_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) #define _SN_STD__PROBE_0 0 diff --git a/StgSharp.Native/sn-src/m.fma.x86.c b/StgSharp.Native/sn-src/m.fma.x86.c deleted file mode 100644 index a357ef71..00000000 --- a/StgSharp.Native/sn-src/m.fma.x86.c +++ /dev/null @@ -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 diff --git a/StgSharp.Native/sn-src/m.fma_tile.x86.c b/StgSharp.Native/sn-src/mk.fma.x86.c similarity index 69% rename from StgSharp.Native/sn-src/m.fma_tile.x86.c rename to StgSharp.Native/sn-src/mk.fma.x86.c index 41b418b4..12251ad1 100644 --- a/StgSharp.Native/sn-src/m.fma_tile.x86.c +++ b/StgSharp.Native/sn-src/mk.fma.x86.c @@ -1,490 +1,414 @@ -#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_MK_PROC_DECL_STD(float, sse, , fma) -{ - __mk_param_std(TILE, LEFT_RIGHT_ANS, float); - - MAT_KERNEL(float) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - register __m128 c0 = _mm_load_ps((float *)&ans_cur->f32_x[0]); - register __m128 c1 = _mm_load_ps((float *)&ans_cur->f32_x[1]); - register __m128 c2 = _mm_load_ps((float *)&ans_cur->f32_x[2]); - register __m128 c3 = _mm_load_ps((float *)&ans_cur->f32_x[3]); - register __m128 bcol; - - 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); - - const __m128 a0 = _mm_load_ps((float *)&left_cur->f32_x[0]); - const __m128 a1 = _mm_load_ps((float *)&left_cur->f32_x[1]); - const __m128 a2 = _mm_load_ps((float *)&left_cur->f32_x[2]); - const __m128 a3 = _mm_load_ps((float *)&left_cur->f32_x[3]); - - /* Column 0 */ - bcol = _mm_load_ps((float *)&right_cur->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_cur->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_cur->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_cur->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_cur->f32_x[0], c0); - _mm_store_ps((float *)&ans_cur->f32_x[1], c1); - _mm_store_ps((float *)&ans_cur->f32_x[2], c2); - _mm_store_ps((float *)&ans_cur->f32_x[3], c3); -} - -SN_MK_PROC_DECL_STD(float, avx, , 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) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - 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 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); - } - _mm256_store_ps((float *)&ans_cur->f32_y[0], c01); - _mm256_store_ps((float *)&ans_cur->f32_y[1], c23); -} - -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) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - 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 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); - } - _mm256_store_ps((float *)&ans_cur->f32_y[0], c01); - _mm256_store_ps((float *)&ans_cur->f32_y[1], c23); -} - -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); - 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); - register __m512 b = _mm512_load_ps((float *)right_cur); - - b = _mm512_permutexvar_ps(b_shift, b); - register __m512 a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(3, 2, 1, 0)); - register __m512 b_perm = _mm512_permute_ps(b, _FMA_SELECT(0)); - c = _mm512_fmadd_ps(a_perm, b_perm, c); - - a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(2, 1, 0, 3)); - b_perm = _mm512_permute_ps(b, _FMA_SELECT(1)); - c = _mm512_fmadd_ps(a_perm, b_perm, c); - - a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(1, 0, 3, 2)); - b_perm = _mm512_permute_ps(b, _FMA_SELECT(2)); - c = _mm512_fmadd_ps(a_perm, b_perm, c); - - a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(0, 3, 2, 1)); - b_perm = _mm512_permute_ps(b, _FMA_SELECT(3)); - c = _mm512_fmadd_ps(a_perm, b_perm, c); - } - - _mm512_store_ps((float *)ans_cur, c); -} - -// double -#define DOUBLE_KERNEL_FMA_CYCLE(i) \ - do { \ - /* k = 0 */ \ - register __m128d b_lo = _mm_load_pd(right_cur + (2 * i + 0)); \ - register __m128d b_hi = _mm_load_pd(right_cur + (2 * i + 1)); \ - register __m128d a_lo = _mm_load_pd(left_cur + (0 * 2 + 0)); \ - register __m128d a_hi = _mm_load_pd(left_cur + (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_cur + (1 * 2 + 0)); \ - a_hi = _mm_load_pd(left_cur + (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_cur + (2 * 2 + 0)); \ - a_hi = _mm_load_pd(left_cur + (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_cur + (3 * 2 + 0)); \ - a_hi = _mm_load_pd(left_cur + (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) - -SN_MK_PROC_DECL_STD(double, sse, , fma) -{ - __mk_param_std(TILE, LEFT_RIGHT_ANS, double); - - MAT_KERNEL(double) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - register __m128d c0_lo = _mm_load_pd(ans_cur->f64_x + (0 * 2 + 0)); - register __m128d c0_hi = _mm_load_pd(ans_cur->f64_x + (0 * 2 + 1)); - register __m128d c1_lo = _mm_load_pd(ans_cur->f64_x + (1 * 2 + 0)); - register __m128d c1_hi = _mm_load_pd(ans_cur->f64_x + (1 * 2 + 1)); - register __m128d c2_lo = _mm_load_pd(ans_cur->f64_x + (2 * 2 + 0)); - register __m128d c2_hi = _mm_load_pd(ans_cur->f64_x + (2 * 2 + 1)); - register __m128d c3_lo = _mm_load_pd(ans_cur->f64_x + (3 * 2 + 0)); - register __m128d c3_hi = _mm_load_pd(ans_cur->f64_x + (3 * 2 + 1)); - - 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); - - DOUBLE_KERNEL_FMA_CYCLE(0); - - DOUBLE_KERNEL_FMA_CYCLE(1); - - DOUBLE_KERNEL_FMA_CYCLE(2); - - DOUBLE_KERNEL_FMA_CYCLE(3); - } - - _mm_store_pd(ans_cur->f64_x + (0 * 2 + 0), c0_lo); - _mm_store_pd(ans_cur->f64_x + (0 * 2 + 1), c0_hi); - _mm_store_pd(ans_cur->f64_x + (1 * 2 + 0), c1_lo); - _mm_store_pd(ans_cur->f64_x + (1 * 2 + 1), c1_hi); - _mm_store_pd(ans_cur->f64_x + (2 * 2 + 0), c2_lo); - _mm_store_pd(ans_cur->f64_x + (2 * 2 + 1), c2_hi); - _mm_store_pd(ans_cur->f64_x + (3 * 2 + 0), c3_lo); - _mm_store_pd(ans_cur->f64_x + (3 * 2 + 1), c3_hi); -} - -SN_MK_PROC_DECL_STD(double, avx, , fma) -{ - __mk_param_std(TILE, LEFT_RIGHT_ANS, double); - - MAT_KERNEL(double) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - register __m256d c0 = _mm256_load_pd((double *)&ans_cur->f64_y[0]); - register __m256d c1 = _mm256_load_pd((double *)&ans_cur->f64_y[1]); - register __m256d c2 = _mm256_load_pd((double *)&ans_cur->f64_y[2]); - register __m256d c3 = _mm256_load_pd((double *)&ans_cur->f64_y[3]); - register __m256d bcol; - - 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); - - const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]); - const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]); - const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]); - const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]); - - /* Column 0 */ - bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]); - c0 = _mm256_add_pd(c0, - _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); - c0 = _mm256_add_pd(c0, - _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); - c0 = _mm256_add_pd(c0, - _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); - c0 = _mm256_add_pd(c0, - _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); - - /* Column 1 */ - bcol = _mm256_load_pd((double *)&right_cur->f64_y[1]); - c1 = _mm256_add_pd(c1, - _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); - c1 = _mm256_add_pd(c1, - _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); - c1 = _mm256_add_pd(c1, - _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); - c1 = _mm256_add_pd(c1, - _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); - - /* Column 2 */ - bcol = _mm256_load_pd((double *)&right_cur->f64_y[2]); - c2 = _mm256_add_pd(c2, - _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); - c2 = _mm256_add_pd(c2, - _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); - c2 = _mm256_add_pd(c2, - _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); - c2 = _mm256_add_pd(c2, - _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); - - /* Column 3 */ - bcol = _mm256_load_pd((double *)&right_cur->f64_y[3]); - c3 = _mm256_add_pd(c3, - _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); - c3 = _mm256_add_pd(c3, - _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); - c3 = _mm256_add_pd(c3, - _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); - c3 = _mm256_add_pd(c3, - _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); - } - - _mm256_store_pd((double *)&ans_cur->f64_y[0], c0); - _mm256_store_pd((double *)&ans_cur->f64_y[1], c1); - _mm256_store_pd((double *)&ans_cur->f64_y[2], c2); - _mm256_store_pd((double *)&ans_cur->f64_y[3], c3); -} - -SN_MK_PROC_DECL_STD(double, avx, _fma, fma) -{ - __mk_param_std(TILE, LEFT_RIGHT_ANS, double); - - MAT_KERNEL(double) - const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe - - register __m256d c0 = _mm256_load_pd((double *)&ans_cur->f64_y[0]); - register __m256d c1 = _mm256_load_pd((double *)&ans_cur->f64_y[1]); - register __m256d c2 = _mm256_load_pd((double *)&ans_cur->f64_y[2]); - register __m256d c3 = _mm256_load_pd((double *)&ans_cur->f64_y[3]); - register __m256d bcol; - - 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); - - const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]); - const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]); - const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]); - const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]); - - /* Column 0 */ - bcol = _mm256_load_pd((double *)&right_cur->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_cur->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_cur->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_cur->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_cur->f64_y[0], c0); - _mm256_store_pd((double *)&ans_cur->f64_y[1], c1); - _mm256_store_pd((double *)&ans_cur->f64_y[2], c2); - _mm256_store_pd((double *)&ans_cur->f64_y[3], c3); -} - -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 - - 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); - - 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 b23 = _mm512_load_pd((double *)&right_cur->f64_z[1]); - b01 = _mm512_permutexvar_pd(b01_shift, b01); - b23 = _mm512_permutexvar_pd(b23_shift, b23); - - //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, 0x21); - register __m512d a30 = _mm512_shuffle_f64x2(a01, a23, 0x30); - 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 +#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_MK_PROC_DECL_STD(float, sse, , fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, float); + + MAT_KERNEL(float) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m128 c0 = _mm_load_ps((float *)&ans_cur->f32_x[0]); + register __m128 c1 = _mm_load_ps((float *)&ans_cur->f32_x[1]); + register __m128 c2 = _mm_load_ps((float *)&ans_cur->f32_x[2]); + register __m128 c3 = _mm_load_ps((float *)&ans_cur->f32_x[3]); + register __m128 bcol; + + 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); + + const __m128 a0 = _mm_load_ps((float *)&left_cur->f32_x[0]); + const __m128 a1 = _mm_load_ps((float *)&left_cur->f32_x[1]); + const __m128 a2 = _mm_load_ps((float *)&left_cur->f32_x[2]); + const __m128 a3 = _mm_load_ps((float *)&left_cur->f32_x[3]); + + /* Column 0 */ + bcol = _mm_load_ps((float *)&right_cur->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_cur->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_cur->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_cur->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_cur->f32_x[0], c0); + _mm_store_ps((float *)&ans_cur->f32_x[1], c1); + _mm_store_ps((float *)&ans_cur->f32_x[2], c2); + _mm_store_ps((float *)&ans_cur->f32_x[3], c3); +} + +SN_MK_PROC_DECL_STD(float, avx, , 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) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m256 c01 = _mm256_load_ps((float *)&ans->f32_y[0]); + register __m256 c23 = _mm256_load_ps((float *)&ans->f32_y[1]); + + 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 b01 = _mm256_load_ps((float *)&right_cur->f32_y[0]); + register __m256 b23 = _mm256_load_ps((float *)&right_cur->f32_y[1]); + + register __m256 a0 = _mm256_broadcast_ps(&left_cur->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_cur->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_cur->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_cur->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_cur->f32_y[0], c01); + _mm256_store_ps((float *)&ans_cur->f32_y[1], c23); +} + +SN_MK_PROC_DECL_STD(float, avx, _fma, fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, float); + + MAT_KERNEL(float) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m256 c01 = _mm256_load_ps((float *)&ans_cur->f32_y[0]); + register __m256 c23 = _mm256_load_ps((float *)&ans_cur->f32_y[1]); + + 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 b01 = _mm256_load_ps((float *)&right_cur->f32_y[0]); + register __m256 b23 = _mm256_load_ps((float *)&right_cur->f32_y[1]); + + register __m256 a0 = _mm256_broadcast_ps((float *)&left_cur->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_cur->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_cur->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_cur->f32_x[3]); + c01 = _mm256_fmadd_ps(a3, b01, c01); + c23 = _mm256_fmadd_ps(a3, b23, c23); + } + _mm256_store_ps((float *)&ans_cur->f32_y[0], c01); + _mm256_store_ps((float *)&ans_cur->f32_y[1], c23); +} + +SN_MK_PROC_DECL_STD(float, 512, , fma) +{ + __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((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); + 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); + + // k=0: broadcast ka_col0 到4个lane，permute_ps 取 kb(0,j) + register __m512 ak; + ak = _mm512_broadcast_f32x4(left_cur->f32_x[0]); + c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(0)), c); + + // k=1 + ak = _mm512_broadcast_f32x4(left_cur->f32_x[1]); + c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(1)), c); + + // k=2 + ak = _mm512_broadcast_f32x4(left_cur->f32_x[2]); + c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(2)), c); + + // k=3 + ak = _mm512_broadcast_f32x4(left_cur->f32_x[3]); + c = _mm512_fmadd_ps(ak, _mm512_permute_ps(b, _FMA_SELECT(3)), c); + } + _mm512_store_ps((float *)ans_cur, c); +} + +// double +#define DOUBLE_KERNEL_FMA_CYCLE(i) \ + do { \ + /* k = 0 */ \ + register __m128d b_lo = _mm_load_pd(right_cur + (2 * i + 0)); \ + register __m128d b_hi = _mm_load_pd(right_cur + (2 * i + 1)); \ + register __m128d a_lo = _mm_load_pd(left_cur + (0 * 2 + 0)); \ + register __m128d a_hi = _mm_load_pd(left_cur + (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_cur + (1 * 2 + 0)); \ + a_hi = _mm_load_pd(left_cur + (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_cur + (2 * 2 + 0)); \ + a_hi = _mm_load_pd(left_cur + (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_cur + (3 * 2 + 0)); \ + a_hi = _mm_load_pd(left_cur + (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) + +SN_MK_PROC_DECL_STD(double, sse, , fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, double); + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m128d c0_lo = _mm_load_pd(ans_cur->f64_x + (0 * 2 + 0)); + register __m128d c0_hi = _mm_load_pd(ans_cur->f64_x + (0 * 2 + 1)); + register __m128d c1_lo = _mm_load_pd(ans_cur->f64_x + (1 * 2 + 0)); + register __m128d c1_hi = _mm_load_pd(ans_cur->f64_x + (1 * 2 + 1)); + register __m128d c2_lo = _mm_load_pd(ans_cur->f64_x + (2 * 2 + 0)); + register __m128d c2_hi = _mm_load_pd(ans_cur->f64_x + (2 * 2 + 1)); + register __m128d c3_lo = _mm_load_pd(ans_cur->f64_x + (3 * 2 + 0)); + register __m128d c3_hi = _mm_load_pd(ans_cur->f64_x + (3 * 2 + 1)); + + 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); + + DOUBLE_KERNEL_FMA_CYCLE(0); + + DOUBLE_KERNEL_FMA_CYCLE(1); + + DOUBLE_KERNEL_FMA_CYCLE(2); + + DOUBLE_KERNEL_FMA_CYCLE(3); + } + + _mm_store_pd(ans_cur->f64_x + (0 * 2 + 0), c0_lo); + _mm_store_pd(ans_cur->f64_x + (0 * 2 + 1), c0_hi); + _mm_store_pd(ans_cur->f64_x + (1 * 2 + 0), c1_lo); + _mm_store_pd(ans_cur->f64_x + (1 * 2 + 1), c1_hi); + _mm_store_pd(ans_cur->f64_x + (2 * 2 + 0), c2_lo); + _mm_store_pd(ans_cur->f64_x + (2 * 2 + 1), c2_hi); + _mm_store_pd(ans_cur->f64_x + (3 * 2 + 0), c3_lo); + _mm_store_pd(ans_cur->f64_x + (3 * 2 + 1), c3_hi); +} + +SN_MK_PROC_DECL_STD(double, avx, , fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, double); + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m256d c0 = _mm256_load_pd((double *)&ans_cur->f64_y[0]); + register __m256d c1 = _mm256_load_pd((double *)&ans_cur->f64_y[1]); + register __m256d c2 = _mm256_load_pd((double *)&ans_cur->f64_y[2]); + register __m256d c3 = _mm256_load_pd((double *)&ans_cur->f64_y[3]); + register __m256d bcol; + + 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); + + const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]); + const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]); + const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]); + const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]); + + /* Column 0 */ + bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]); + c0 = _mm256_add_pd(c0, + _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); + c0 = _mm256_add_pd(c0, + _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); + c0 = _mm256_add_pd(c0, + _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); + c0 = _mm256_add_pd(c0, + _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); + + /* Column 1 */ + bcol = _mm256_load_pd((double *)&right_cur->f64_y[1]); + c1 = _mm256_add_pd(c1, + _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); + c1 = _mm256_add_pd(c1, + _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); + c1 = _mm256_add_pd(c1, + _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); + c1 = _mm256_add_pd(c1, + _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); + + /* Column 2 */ + bcol = _mm256_load_pd((double *)&right_cur->f64_y[2]); + c2 = _mm256_add_pd(c2, + _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); + c2 = _mm256_add_pd(c2, + _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); + c2 = _mm256_add_pd(c2, + _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); + c2 = _mm256_add_pd(c2, + _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); + + /* Column 3 */ + bcol = _mm256_load_pd((double *)&right_cur->f64_y[3]); + c3 = _mm256_add_pd(c3, + _mm256_mul_pd(a0, _mm256_permute4x64_pd(bcol, _FMA_SELECT(0)))); + c3 = _mm256_add_pd(c3, + _mm256_mul_pd(a1, _mm256_permute4x64_pd(bcol, _FMA_SELECT(1)))); + c3 = _mm256_add_pd(c3, + _mm256_mul_pd(a2, _mm256_permute4x64_pd(bcol, _FMA_SELECT(2)))); + c3 = _mm256_add_pd(c3, + _mm256_mul_pd(a3, _mm256_permute4x64_pd(bcol, _FMA_SELECT(3)))); + } + + _mm256_store_pd((double *)&ans_cur->f64_y[0], c0); + _mm256_store_pd((double *)&ans_cur->f64_y[1], c1); + _mm256_store_pd((double *)&ans_cur->f64_y[2], c2); + _mm256_store_pd((double *)&ans_cur->f64_y[3], c3); +} + +SN_MK_PROC_DECL_STD(double, avx, _fma, fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, double); + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(common_y, x, y); // GetKernelAddressUnsafe + + register __m256d c0 = _mm256_load_pd((double *)&ans_cur->f64_y[0]); + register __m256d c1 = _mm256_load_pd((double *)&ans_cur->f64_y[1]); + register __m256d c2 = _mm256_load_pd((double *)&ans_cur->f64_y[2]); + register __m256d c3 = _mm256_load_pd((double *)&ans_cur->f64_y[3]); + register __m256d bcol; + + 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); + + const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]); + const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]); + const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]); + const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]); + + /* Column 0 */ + bcol = _mm256_load_pd((double *)&right_cur->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_cur->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_cur->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_cur->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_cur->f64_y[0], c0); + _mm256_store_pd((double *)&ans_cur->f64_y[1], c1); + _mm256_store_pd((double *)&ans_cur->f64_y[2], c2); + _mm256_store_pd((double *)&ans_cur->f64_y[3], c3); +} + +SN_MK_PROC_DECL_STD(double, 512, , fma) +{ + __mk_param_std(TILE, LEFT_RIGHT_ANS, double); + + MAT_KERNEL(double) + 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]); + + 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 __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 b0 = _mm512_broadcast_f64x4(right_cur->f64_y[0]); + c01 = _mm512_fmadd_pd(a01, b0, c01); + c23 = _mm512_fmadd_pd(a23, b0, c23); + + register __m512d b1 = _mm512_broadcast_f64x4(right_cur->f64_y[1]); + c01 = _mm512_fmadd_pd(a01, b1, c01); + c23 = _mm512_fmadd_pd(a23, b1, c23); + + register __m512d b2 = _mm512_broadcast_f64x4(right_cur->f64_y[2]); + c01 = _mm512_fmadd_pd(a01, b2, c01); + c23 = _mm512_fmadd_pd(a23, b2, c23); + + register __m512d b3 = _mm512_broadcast_f64x4(right_cur->f64_y[3]); + c01 = _mm512_fmadd_pd(a01, b3, c01); + c23 = _mm512_fmadd_pd(a23, b3, c23); + } + _mm512_store_pd((double *)&ans_cur->f64_z[0], c01); + _mm512_store_pd((double *)&ans_cur->f64_z[1], c23); +} +#endif diff --git a/StgSharp.Native/sn-src/sk.fma.x86.c b/StgSharp.Native/sn-src/sk.fma.x86.c new file mode 100644 index 00000000..52433927 --- /dev/null +++ b/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