From c6d53d97eafad7ae4b15167f79b495512be2eabf Mon Sep 17 00:00:00 2001 From: Youmu Liu <105042599+Nitload-NSI@users.noreply.github.com> Date: Fri, 26 Dec 2025 17:16:12 +0800 Subject: [PATCH] add double compute kernel --- .../HighPerformance/MatrixPanel.cs | 49 +-- .../Internal/Intrinsic/IntrinsicContext.cs | 33 +- .../Mathematics/Graphic/Matrix44.cs | 6 +- StgSharp.Common/Mathematics/Linear/Matrix.cs | 52 +++- .../Numeric/MatrixCompute.Additional.cs | 90 ++++-- .../Numeric/MatrixCompute.Factorize.cs | 6 +- StgSharp.Mathematics/Numeric/MatrixCompute.cs | 108 ++++--- StgSharp.Native/sn-head/sn_intrinsic.h | 57 +++- StgSharp.Native/sn-head/sn_matkernel.h | 64 ++-- StgSharp.Native/sn-src/k_basicop.c | 53 +++- StgSharp.Native/sn-src/k_bufop.c | 281 +++++++++++++++--- StgSharp.Native/sn-src/k_panel.c | 173 +++++++++-- StgSharp.Native/sn-src/k_pivot.c | 4 +- StgSharp.Native/sn-src/load_intrinsic.c | 105 +++++-- StgSharp.Native/sn-src/vector.c | 8 +- StgSharp.Native/sn-src/vector_dot.c | 32 +- 16 files changed, 821 insertions(+), 300 deletions(-) diff --git a/StgSharp.Common/HighPerformance/MatrixPanel.cs b/StgSharp.Common/HighPerformance/MatrixPanel.cs index 7f3ae6d0..5a99f17d 100644 --- a/StgSharp.Common/HighPerformance/MatrixPanel.cs +++ b/StgSharp.Common/HighPerformance/MatrixPanel.cs @@ -61,17 +61,9 @@ namespace StgSharp.HighPerformance [MethodImpl(MethodImplOptions.AggressiveInlining)] public static int ElementSideLength() where T : unmanaged, INumber { - return Math.Max(8, AlignmentBytes / sizeof(T)); + return PanelSideLength(); } - // View an existing pointer; caller guarantees lifetime/alignment. - public static unsafe MatrixPanel* FromPointer(void* ptr) - where T : unmanaged, INumber - { - return (MatrixPanel*)Align64(ptr); - } - - [MethodImpl(MethodImplOptions.AggressiveInlining)] public static int MaxElementCount() where T : unmanaged, INumber { int side = ElementSideLength(); @@ -115,28 +107,43 @@ namespace StgSharp.HighPerformance [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int KernelSideCountAvx2() where T : unmanaged, INumber { - const int minQ = 8; - int lanes = 32 / sizeof(T); - int side = Math.Max(lanes, minQ); - return side / 4; + return PanelSideFromBits(256) / 4; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int KernelSideCountAvx512() where T : unmanaged, INumber { - const int minQ = 8; - int lanes = 64 / sizeof(T); - int side = Math.Max(lanes, minQ); - return side / 4; + return PanelSideFromBits(512) / 4; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static int KernelSideCountSse() where T : unmanaged, INumber { - const int minQ = 4; - int lanes = 16 / sizeof(T); - int side = Math.Max(lanes, minQ); - return side / 4; + return PanelSideFromBits(128) / 4; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int PanelSideFromBits(int vectorBits) where T : unmanaged, INumber + { + int elemBits = sizeof(T) * 8; + + // Q = max(4, vectorBits/64, vectorBits/elemBits) + int by64 = vectorBits / 64; + int byElem = vectorBits / elemBits; + int side = Math.Max(4, Math.Max(by64, byElem)); + return side; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static int PanelSideLength() where T : unmanaged, INumber + { + return NativeIntrinsic.IntrinsicLevel switch + { + IntrinsicLevel.SSE => PanelSideFromBits(128), + IntrinsicLevel.AVX2 => PanelSideFromBits(256), + IntrinsicLevel.AVX512 => PanelSideFromBits(512), + _ => 0 + }; } } diff --git a/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs b/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs index 53475398..a1fea492 100644 --- a/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs +++ b/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs @@ -39,20 +39,29 @@ namespace StgSharp.Internal.Intrinsic { public delegate* unmanaged[Cdecl] city_hash_simplify; //city_hash_simplify - public delegate* unmanaged[Cdecl] f32_buffer_add; //f32_buffer_add - public delegate* unmanaged[Cdecl] f32_buffer_fill; //f32_buffer_fill - public delegate* unmanaged[Cdecl] f32_buffer_scalar_mul; //f32_buffer_s_m - public delegate* unmanaged[Cdecl] f32_buffer_sub; //f32_kernel_sub - public delegate* unmanaged[Cdecl] f32_buffer_transpose; //f32_kernel_transpose - public delegate* unmanaged[Cdecl] f32_build_panel; //f32_build_panel - public delegate* unmanaged[Cdecl] f32_clear_panel; - public delegate* unmanaged[Cdecl] f32_normalize_3; //f32_normalize_3 - public delegate* unmanaged[Cdecl] f32_panel_fma; //f32_panel_fma - public delegate* unmanaged[Cdecl] f32_pivot; //f32_panel_mul - public delegate* unmanaged[Cdecl] f32_store_panel; //f32_store_panel public delegate* unmanaged[Cdecl] factorial_simd; //factorial_simd public delegate* unmanaged[Cdecl] index_pair; //index_pair - public delegate* unmanaged[Cdecl] variant ; //variant + public delegate* unmanaged[Cdecl] mat_f32_buffer_add; //f32_buffer_add + public delegate* unmanaged[Cdecl] mat_f32_buffer_fill; //f32_buffer_fill + public delegate* unmanaged[Cdecl] mat_f32_buffer_scalar_mul; //f32_buffer_s_m + public delegate* unmanaged[Cdecl] mat_f32_buffer_sub; //f32_kernel_sub + public delegate* unmanaged[Cdecl] mat_f32_buffer_transpose; //f32_kernel_transpose + public delegate* unmanaged[Cdecl] mat_f32_build_panel; //f32_build_panel + public delegate* unmanaged[Cdecl] mat_f32_clear_panel; //f32_clear_panel + public delegate* unmanaged[Cdecl] mat_f32_panel_fma; //f32_panel_fma + public delegate* unmanaged[Cdecl] mat_f32_pivot; //f32_panel_mul + public delegate* unmanaged[Cdecl] mat_f32_store_panel; //f32_store_panel + public delegate* unmanaged[Cdecl] mat_f64_buffer_add; //f32_buffer_add + public delegate* unmanaged[Cdecl] mat_f64_buffer_fill; //f32_buffer_fill + public delegate* unmanaged[Cdecl] mat_f64_buffer_scalar_mul; //f32_buffer_s_m + public delegate* unmanaged[Cdecl] mat_f64_buffer_sub; //f32_kernel_sub + public delegate* unmanaged[Cdecl] mat_f64_buffer_transpose; //f32_kernel_transpose + public delegate* unmanaged[Cdecl] mat_f64_build_panel; //f32_build_panel + public delegate* unmanaged[Cdecl] mat_f64_clear_panel; //f32_clear_panel + public delegate* unmanaged[Cdecl] mat_f64_panel_fma; //f32_panel_fma + public delegate* unmanaged[Cdecl] mat_f64_pivot; //f32_panel_mul + public delegate* unmanaged[Cdecl] mat_f64_store_panel; //f32_store_panel + public delegate* unmanaged[Cdecl] vec_f32_normalize_3; //f32_normalize_3 public IntrinsicContext() { diff --git a/StgSharp.Common/Mathematics/Graphic/Matrix44.cs b/StgSharp.Common/Mathematics/Graphic/Matrix44.cs index 3d60ecb5..90a55590 100644 --- a/StgSharp.Common/Mathematics/Graphic/Matrix44.cs +++ b/StgSharp.Common/Mathematics/Graphic/Matrix44.cs @@ -169,7 +169,7 @@ namespace StgSharp.Mathematics.Graphic internal unsafe void InternalTranspose(ref Matrix44 transpose) { fixed (MatrixKernel* source = &this.kernel, target = &transpose.kernel) { - NativeIntrinsic.Context.f32_buffer_transpose(source, target); + NativeIntrinsic.Context.mat_f32_buffer_transpose(source, target); } } @@ -177,7 +177,7 @@ namespace StgSharp.Mathematics.Graphic public static unsafe Matrix44 operator -(Matrix44 left, Matrix44 right) { Matrix44 ret = new(); - NativeIntrinsic.Context.f32_buffer_sub(&left, &right, &ret, 1); + NativeIntrinsic.Context.mat_f32_buffer_sub(&left, &right, &ret, 1); return ret; } @@ -219,7 +219,7 @@ namespace StgSharp.Mathematics.Graphic public static unsafe Matrix44 operator +(Matrix44 left, Matrix44 right) { Matrix44 ret = new(); - NativeIntrinsic.Context.f32_buffer_add(&left, &right, &ret, 1); + NativeIntrinsic.Context.mat_f32_buffer_add(&left, &right, &ret, 1); return ret; } diff --git a/StgSharp.Common/Mathematics/Linear/Matrix.cs b/StgSharp.Common/Mathematics/Linear/Matrix.cs index 26bfd118..c2527bec 100644 --- a/StgSharp.Common/Mathematics/Linear/Matrix.cs +++ b/StgSharp.Common/Mathematics/Linear/Matrix.cs @@ -50,14 +50,17 @@ namespace StgSharp.Mathematics.Numeric MatrixSize size, MatrixStorage buffer - ) + ) { _boarder = boarder; _size = size; _buffer = buffer; } - public T this[int x, int y] + public T this[ + int x, + int y + ] { get { @@ -108,7 +111,8 @@ namespace StgSharp.Mathematics.Numeric /// public unsafe Matrix this[ ReadOnlySpan column, - ReadOnlySpan row] + ReadOnlySpan row + ] { get { throw new NotImplementedException(); } set { throw new NotImplementedException(); } @@ -137,7 +141,10 @@ namespace StgSharp.Mathematics.Numeric protected internal int KernelRowLength => _size.KernelRowLength; [MethodImpl(MethodImplOptions.AggressiveInlining)] - public unsafe ref MatrixKernel GetMatrixKernel(int x, int y) + public unsafe ref MatrixKernel GetMatrixKernel( + int x, + int y + ) { MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength, x, y); @@ -147,7 +154,10 @@ namespace StgSharp.Mathematics.Numeric } [MethodImpl(MethodImplOptions.AggressiveInlining)] - internal unsafe ref T GetElementUnsafe(int x, int y) + internal unsafe ref T GetElementUnsafe( + int x, + int y + ) { // Caller guarantees: indices in range, layout is dense square. int k_c = y >> 2; @@ -161,5 +171,37 @@ namespace StgSharp.Mathematics.Numeric return ref kernel[c, r]; } + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal unsafe MatrixKernel* GetKernelAddressUnsafe( + int x, + int y + ) + { + // Caller guarantees: indices in range, layout is dense square. + int k_c = y >> 2; + int k_r = x >> 2; + + MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength, + k_c, k_r); + return Buffer + _boarder.GetIndexOffset(k_c, k_r, in enumerator); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal unsafe ref MatrixKernel GetKernelUnsafe( + int x, + int y + ) + { + // Caller guarantees: indices in range, layout is dense square. + int k_c = y >> 2; + int k_r = x >> 2; + int c = y & 3; + int r = x & 3; + + MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength, + k_c, k_r); + return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in enumerator)]; + } + } } diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs index e7f26ad0..0a46703c 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs @@ -30,6 +30,7 @@ using StgSharp.Internal; using System; using System.Collections.Generic; using System.Linq; +using System.Numerics; using System.Runtime.CompilerServices; using System.Text; using System.Threading.Tasks; @@ -39,45 +40,84 @@ namespace StgSharp.Mathematics.Numeric public partial class MatrixCompute { - internal static unsafe void ComputeFmaPanelF32(MatrixParallelTaskPackage* p) + internal static unsafe void ComputeFmaPanelF32( + MatrixParallelTaskPackage* p + ) { int size = p->ElementSize; MatrixKernel* leftBuffer = (MatrixKernel*)p->Left; MatrixKernel* rightBuffer = (MatrixKernel*)p->Right; MatrixKernel* ansBuffer = (MatrixKernel*)p->Answer; - MatrixPanel* leftPanel = MatrixPanel.Create(); - MatrixPanel* rightPanel = MatrixPanel.Create(); - MatrixPanel* ansPanel = MatrixPanel.Create(); ref long offsetRef = ref Unsafe.As(ref p->PrimTileOffset); ref long lengthRef = ref Unsafe.As(ref p->PrimCount); - int ansWidth = p-> AnsPrimOffset; + int ansWidth = p->AnsPrimOffset; int ansHeight = p->AnsSecOffset; int commonK = p->LeftPrimOffset; int blocksPerRow = ansWidth / size; - - for (long c = offsetRef; c < offsetRef + lengthRef; c++) + if (NativeIntrinsic.Context.mat_f32_build_panel == null) { - int i = (int)(c / blocksPerRow) * size; - int j = (int)(c % blocksPerRow) * size; - NativeIntrinsic.Context.f32_clear_panel((MatrixPanel*)ansPanel); - for (int k = 0; k < commonK; k += size) + // in this occasion size must be 1 + for (long c = offsetRef; c < offsetRef + lengthRef; c++) { - NativeIntrinsic.Context - .f32_build_panel((MatrixPanel*)leftPanel, leftBuffer, commonK, - ansHeight, k, i); - NativeIntrinsic.Context - .f32_build_panel((MatrixPanel*)rightPanel, rightBuffer, ansWidth, - commonK, j, k); - NativeIntrinsic.Context - .f32_panel_fma((MatrixPanel*)leftPanel, (MatrixPanel*)rightPanel, (MatrixPanel*)ansPanel); + int i = (int)(c / blocksPerRow); + int j = (int)(c % blocksPerRow); + MatrixKernel* ansPanel = GetKernelAddressUnsafe((MatrixKernel*)ansBuffer, + ansHeight, i, j); + NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel); + for (int k = 0; k < commonK; k += size) { + NativeIntrinsic.Context + .mat_f32_panel_fma( + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)leftBuffer, + ansHeight, k, j), + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)rightBuffer, + commonK, j, k), + (MatrixPanel*)ansPanel); + } } - NativeIntrinsic.Context - .f32_store_panel((MatrixPanel*)ansPanel, ansBuffer, ansWidth, - ansHeight, j, i); + } else + { + MatrixPanel* leftPanel = MatrixPanel.Create(); + MatrixPanel* rightPanel = MatrixPanel.Create(); + MatrixPanel* ansPanel = MatrixPanel.Create(); + for (long c = offsetRef; c < offsetRef + lengthRef; c++) + { + int i = (int)(c / blocksPerRow) * size; + int j = (int)(c % blocksPerRow) * size; + NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel); + for (int k = 0; k < commonK; k += size) + { + NativeIntrinsic.Context + .mat_f32_build_panel((MatrixPanel*)leftPanel, leftBuffer, + commonK, ansHeight, k, i); + NativeIntrinsic.Context + .mat_f32_build_panel((MatrixPanel*)rightPanel, rightBuffer, + ansWidth, commonK, j, k); + NativeIntrinsic.Context + .mat_f32_panel_fma((MatrixPanel*)leftPanel, (MatrixPanel*)rightPanel, (MatrixPanel*)ansPanel); + } + NativeIntrinsic.Context + .mat_f32_store_panel((MatrixPanel*)ansPanel, ansBuffer, ansWidth, + ansHeight, j, i); + } + MatrixPanel.Destroy(leftPanel); + MatrixPanel.Destroy(rightPanel); + MatrixPanel.Destroy(ansPanel); } - MatrixPanel.Destroy(leftPanel); - MatrixPanel.Destroy(rightPanel); - MatrixPanel.Destroy(ansPanel); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static unsafe MatrixKernel* GetKernelAddressUnsafe( + MatrixKernel* buffer, + int colLength, + int x, + int y + ) where T : unmanaged, INumber + { + // Caller guarantees: indices in range, layout is dense square. + int k_c = y >> 2; + int k_r = x >> 2; + + return buffer + colLength * x + y; } } diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs index f4ebcf7b..0614187a 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs @@ -77,9 +77,9 @@ namespace StgSharp.Mathematics.Numeric for (int kBlock = 0; kBlock < n; kBlock += 4) { NativeIntrinsic.Context - .f32_pivot(matrix.Buffer, (ulong)(kBlock >> 2), - (ulong)(matrix.KernelRowLength - (kBlock >> 2)), - (ulong)pivots); + .mat_f32_pivot(matrix.Buffer, (ulong)(kBlock >> 2), + (ulong)(matrix.KernelRowLength - (kBlock >> 2)), + (ulong)pivots); int colsInBlock = Math.Min(4, n - kBlock); for (int offset = 0; offset < colsInBlock; offset++) { diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.cs index e5caaaa1..80c0082a 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.cs @@ -56,7 +56,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_add(leftPtr, rightPtr, ansPtr, count); + NativeIntrinsic.Context.mat_f32_buffer_add(leftPtr, rightPtr, ansPtr, count); return; } else if (count <= 1024) @@ -66,7 +66,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_add(leftPtr, rightPtr, ansPtr, count); + NativeIntrinsic.Context.mat_f32_buffer_add(leftPtr, rightPtr, ansPtr, count); MatrixParallel.ReturnParallelResources(ref pool); return; @@ -79,7 +79,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)NativeIntrinsic.Context.f32_buffer_add; + package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_add; MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package); MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep); @@ -94,12 +94,12 @@ namespace StgSharp.Mathematics.Numeric scalar->Data(0) = value; if (count <= 256) { - NativeIntrinsic.Context.f32_buffer_fill(ans.Buffer, scalar, count); + NativeIntrinsic.Context.mat_f32_buffer_fill(ans.Buffer, scalar, count); return; } else if (count <= 1024) { MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1); - NativeIntrinsic.Context.f32_buffer_fill(ans.Buffer, scalar, count); + NativeIntrinsic.Context.mat_f32_buffer_fill(ans.Buffer, scalar, count); MatrixParallel.ReturnParallelResources(ref pool); return; } @@ -109,7 +109,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)NativeIntrinsic.Context.f32_buffer_fill; + package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_fill; package->Scalar = scalar; MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package); MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep); @@ -135,43 +135,65 @@ namespace StgSharp.Mathematics.Numeric if (count > 256) { pool = MatrixParallel.LeadParallel(1); } - MatrixPanel* leftPanel = MatrixPanel.Create(); - MatrixPanel* rightPanel = MatrixPanel.Create(); - MatrixPanel* ansPanel = MatrixPanel.Create(); - for (int i = 0; i < ans.KernelRowLength; i += size) + if (NativeIntrinsic.Context.mat_f32_build_panel == null) { - for (int j = 0; j < ans.KernelColumnLength; j += size) + for (int i = 0; i < ans.KernelRowLength; i += size) { - NativeIntrinsic.Context.f32_clear_panel((MatrixPanel*)ansPanel); - - // Console.WriteLine($"{i}{j}"); - for (int k = 0; k < left.KernelColumnLength; k += size) + for (int j = 0; j < ans.KernelColumnLength; j += size) { - NativeIntrinsic.Context - .f32_build_panel((MatrixPanel*)leftPanel, - (MatrixKernel*)left.Buffer, - left.KernelColumnLength, - left.KernelRowLength, i, k); - NativeIntrinsic.Context - .f32_build_panel((MatrixPanel*)rightPanel, - (MatrixKernel*)right.Buffer, - right.KernelColumnLength, - right.KernelRowLength, k, j); - NativeIntrinsic.Context - .f32_panel_fma((MatrixPanel*)leftPanel, - (MatrixPanel*)rightPanel, - (MatrixPanel*)ansPanel); + MatrixKernel* ansPanel = ans.GetKernelAddressUnsafe(i, j); + NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel); + + // Console.WriteLine($"{i}{j}"); + for (int k = 0; k < left.KernelColumnLength; k += size) + { + NativeIntrinsic.Context + .mat_f32_panel_fma((MatrixPanel*)left.GetKernelAddressUnsafe(i, k), + (MatrixPanel*)right.GetKernelAddressUnsafe(k, j), + (MatrixPanel*)ansPanel); + } } - NativeIntrinsic.Context - .f32_store_panel((MatrixPanel*)ansPanel, - (MatrixKernel*)ans.Buffer, - ans.KernelColumnLength, ans.KernelRowLength, - i, j); } + } else + { + MatrixPanel* leftPanel = MatrixPanel.Create(); + MatrixPanel* rightPanel = MatrixPanel.Create(); + MatrixPanel* ansPanel = MatrixPanel.Create(); + for (int i = 0; i < ans.KernelRowLength; i += size) + { + for (int j = 0; j < ans.KernelColumnLength; j += size) + { + NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel); + + // Console.WriteLine($"{i}{j}"); + for (int k = 0; k < left.KernelColumnLength; k += size) + { + NativeIntrinsic.Context + .mat_f32_build_panel((MatrixPanel*)leftPanel, + (MatrixKernel*)left.Buffer, + left.KernelColumnLength, + left.KernelRowLength, i, k); + NativeIntrinsic.Context + .mat_f32_build_panel((MatrixPanel*)rightPanel, + (MatrixKernel*)right.Buffer, + right.KernelColumnLength, + right.KernelRowLength, k, j); + NativeIntrinsic.Context + .mat_f32_panel_fma((MatrixPanel*)leftPanel, + (MatrixPanel*)rightPanel, + (MatrixPanel*)ansPanel); + } + NativeIntrinsic.Context + .mat_f32_store_panel((MatrixPanel*)ansPanel, + (MatrixKernel*)ans.Buffer, + ans.KernelColumnLength, + ans.KernelRowLength, i, j); + } + } + MatrixPanel.Destroy(leftPanel); + MatrixPanel.Destroy(rightPanel); + MatrixPanel.Destroy(ansPanel); } - MatrixPanel.Destroy(leftPanel); - MatrixPanel.Destroy(rightPanel); - MatrixPanel.Destroy(ansPanel); if (count > 256) { MatrixParallel.ReturnParallelResources(ref pool!); } @@ -211,7 +233,7 @@ namespace StgSharp.Mathematics.Numeric { MatrixKernel* leftPtr = left.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count); + NativeIntrinsic.Context.mat_f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count); return; } else if (count <= 1024) @@ -220,7 +242,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count); + NativeIntrinsic.Context.mat_f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count); MatrixParallel.ReturnParallelResources(ref pool); return; @@ -233,7 +255,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)NativeIntrinsic.Context.f32_buffer_scalar_mul; + package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_scalar_mul; MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package); MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep); @@ -261,7 +283,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count); + NativeIntrinsic.Context.mat_f32_buffer_sub(leftPtr, rightPtr, ansPtr, count); return; } else if (count <= 1024) @@ -271,7 +293,7 @@ namespace StgSharp.Mathematics.Numeric MatrixKernel* leftPtr = left.Buffer; MatrixKernel* rightPtr = right.Buffer; MatrixKernel* ansPtr = ans.Buffer; - NativeIntrinsic.Context.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count); + NativeIntrinsic.Context.mat_f32_buffer_sub(leftPtr, rightPtr, ansPtr, count); MatrixParallel.ReturnParallelResources(ref pool); return; @@ -284,7 +306,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)NativeIntrinsic.Context.f32_buffer_sub; + package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_sub; MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package); MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep); diff --git a/StgSharp.Native/sn-head/sn_intrinsic.h b/StgSharp.Native/sn-head/sn_intrinsic.h index 1b05a915..ebd87154 100644 --- a/StgSharp.Native/sn-head/sn_intrinsic.h +++ b/StgSharp.Native/sn-head/sn_intrinsic.h @@ -215,10 +215,18 @@ DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, add); DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, add); DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, add); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, add); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, add); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, add); + DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, sub); DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, sub); DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, sub); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, sub); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, sub); +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, sub); + DECLARE_KER_PROC_RIGHT_ANS(float, sse, transpose); DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose); DECLARE_KER_PROC_RIGHT_ANS(float, 512, transpose); @@ -227,10 +235,18 @@ DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul); DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul); DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul); +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, sse, scalar_mul); +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, avx, scalar_mul); +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, 512, scalar_mul); + DECLARE_BUF_PROC_ANS_SCALAR(float, sse, fill); DECLARE_BUF_PROC_ANS_SCALAR(float, avx, fill); DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill); +DECLARE_BUF_PROC_ANS_SCALAR(double, sse, fill); +DECLARE_BUF_PROC_ANS_SCALAR(double, avx, fill); +DECLARE_BUF_PROC_ANS_SCALAR(double, 512, fill); + #pragma endregion #pragma region fma @@ -239,6 +255,10 @@ DECLARE_PNL_PROC_ANS(float, sse, clear); DECLARE_PNL_PROC_ANS(float, avx, clear); DECLARE_PNL_PROC_ANS(float, 512, clear); +DECLARE_PNL_PROC_ANS(double, sse, clear); +DECLARE_PNL_PROC_ANS(double, avx, clear); +DECLARE_PNL_PROC_ANS(double, 512, clear); + DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, sse, fma); DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx, fma); DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx_fma, fma); @@ -278,20 +298,33 @@ INTERNAL int SN_DECL index_pair_sse(short const *str, uint32_t target, int lengt // Mirrors C# IntrinsicContext field order so managed/unmanaged layouts stay in sync. typedef struct sn_intrinsic { HASH city_hash_simplify; - BUF_PROC_LRA f32_buffer_add; - BUF_PROC_AS f32_buffer_fill; - BUF_PROC_RAS f32_buffer_scalar_mul; - BUF_PROC_LRA f32_buffer_sub; - KER_PROC_RA f32_buffer_transpose; - PANEL_PROC f32_build_panel; - PANEL_PROC f32_clear_panel; - VECTORNORMALIZEPROC f32_normalize_3; - PNL_PROC_LRA f32_panel_fma; - PIVOTPROC f32_pivot; - PANEL_PROC f32_store_panel; FACTORIALROC factorial_simd; INDEX_PAIR index_pair; - BUF_PROC_RAS variant; + + BUF_PROC_LRA mat_f32_buffer_add; + BUF_PROC_AS mat_f32_buffer_fill; + BUF_PROC_RAS mat_f32_buffer_scalar_mul; + BUF_PROC_LRA mat_f32_buffer_sub; + KER_PROC_RA mat_f32_buffer_transpose; + PANEL_PROC mat_f32_build_panel; + PANEL_PROC mat_f32_clear_panel; + PNL_PROC_LRA mat_f32_panel_fma; + PIVOTPROC mat_f32_pivot; + PANEL_PROC mat_f32_store_panel; + + // f64 slots (currently unimplemented; filled with NULL in loader) + BUF_PROC_LRA mat_f64_buffer_add; + BUF_PROC_AS mat_f64_buffer_fill; + BUF_PROC_RAS mat_f64_buffer_scalar_mul; + BUF_PROC_LRA mat_f64_buffer_sub; + KER_PROC_RA mat_f64_buffer_transpose; + PANEL_PROC mat_f64_build_panel; + PANEL_PROC mat_f64_clear_panel; + PNL_PROC_LRA mat_f64_panel_fma; + PIVOTPROC mat_f64_pivot; + PANEL_PROC mat_f64_store_panel; + + VECTORNORMALIZEPROC vec_f32_normalize_3; } sn_intrinsic; typedef enum most_advanced_instruction { diff --git a/StgSharp.Native/sn-head/sn_matkernel.h b/StgSharp.Native/sn-head/sn_matkernel.h index 52b58d11..939ff308 100644 --- a/StgSharp.Native/sn-head/sn_matkernel.h +++ b/StgSharp.Native/sn-head/sn_matkernel.h @@ -12,20 +12,28 @@ #define YMMCOUNT(T) (sizeof(T) / 2) #define ZMMCOUNT(T) (sizeof(T) / 4) -#define MAT_KERNEL(T) \ - union { \ - __m128 xmm[XMMCOUNT(T) / 4]; \ - __m256 ymm[YMMCOUNT(T) / 4]; \ - __m512 zmm[ZMMCOUNT(T) / 4]; \ - T m[4][4]; \ +#define MAT_KERNEL(T) \ + union /* MAT_KERNEL_##T## /**/ { \ + __m128 f32_x[XMMCOUNT(T)]; \ + __m256 f32_y[YMMCOUNT(T)]; \ + __m512 f32_z[ZMMCOUNT(T)]; \ + \ + __m128d f64_x[XMMCOUNT(T)]; \ + __m256d f64_y[YMMCOUNT(T)]; \ + __m512d f64_z[ZMMCOUNT(T)]; \ + T m[4][4]; \ } -#define VEC(T) \ - union { \ - __m128 xmm[XMMCOUNT(T) / 4]; \ - __m256 ymm[YMMCOUNT(T) / 4]; \ - __m512 zmm[ZMMCOUNT(T) / 4]; \ - T v[4]; \ +#define VEC(T) \ + union { \ + __m128 f32_x[XMMCOUNT(T) / 4]; \ + __m256 f32_y[YMMCOUNT(T) / 4]; \ + __m512 f32_z[ZMMCOUNT(T) / 4]; \ + \ + __m128d f64_x[XMMCOUNT(T) / 4]; \ + __m256d f64_y[YMMCOUNT(T) / 4]; \ + __m512d f64_z[ZMMCOUNT(T) / 4]; \ + T v[4]; \ } /* Provide a named alias for the zero kernel type to ensure consistent linkage in C++. */ @@ -35,6 +43,8 @@ typedef MAT_KERNEL(uint64_t) sn_zero_kernel_u64; #define SN_LANES(T_VEC, T) (sizeof(T_VEC) / sizeof(T)) #define SN_MAX(a, b) (((a)) > ((b)) ? (a) : (b)) #define SN_CEIL_DIV(a, b) (((a) + (b) - 1) / (b)) +#define SN_BITS(x) ((int)(sizeof(x) * 8)) +#define PANEL_MINQ(T_VEC, T) SN_MAX(4, SN_MAX(SN_BITS(T_VEC) / 64, SN_BITS(T_VEC) / SN_BITS(T))) /* Panel geometry helpers (all compile-time): */ #define PANEL_Q(T, T_VEC, MINQ) SN_MAX(SN_LANES(T_VEC, T), (MINQ)) @@ -60,26 +70,25 @@ typedef MAT_KERNEL(uint64_t) sn_zero_kernel_u64; /* ISA presets with minimum square side: - SSE : MINQ=4 (avoid 2x2 for double) - - AVX2 : MINQ=8 - - AVX-512: MINQ=8 (16 for float/i32 naturally, so this is a no-op there) + - AVX2 : MINQ computed as max(4, vector_bits/64, vector_bits/element_bits) + - AVX-512: same rule (float/i32 naturally expand to 16) */ /* SSE panel layout: column-major buffer holding 4 logical columns with MINQ=4. Columns: each column packs up to four scalars so ALIGN equals sizeof(__m128). Alignment: every column boundary lands on a 16-byte multiple for XMM loads. */ -#define MAT_PANEL_sse(T) SN_ALIGNAS(16) MAT_PANEL_SQ_MIN_NOALIGN(T, __m128, 4) -/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ=8. - Columns: each column stores eight elements via two 256-bit chunks when T is 32-bit. +#define MAT_PANEL_sse(T) SN_ALIGNAS(16) MAT_PANEL_SQ_MIN_NOALIGN(T, __m128, PANEL_MINQ(__m128, T)) +/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ per PANEL_MINQ. + Columns: each column stores up to lanes elements via contiguous chunks. Alignment: column starts are 32-byte aligned so AVX2 loads/stores remain contiguous. */ -#define MAT_PANEL_avx(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, 8) -/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ=8. - Columns: each column stores eight elements via two 256-bit chunks when T is 32-bit. - Alignment: column starts are 32-byte aligned so AVX2 loads/stores remain contiguous. */ -#define MAT_PANEL_avx_fma(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, 8) +#define MAT_PANEL_avx(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, PANEL_MINQ(__m256, T)) +/* AVX2/FMA panel layout shares the same geometry as AVX2. */ +#define MAT_PANEL_avx_fma(T) \ + SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, PANEL_MINQ(__m256, T)) /* AVX-512 panel layout: column-major buffer holding 8 or more columns depending on T. - Columns: each column expands to ceil(Q/lanes) ZMM chunks; default MINQ=8 keeps doublse-wide tiles. + Columns: each column expands to ceil(Q/lanes) ZMM chunks. Alignment: column starts sit on 64-byte boundaries to match AVX-512 streaming access. */ -#define MAT_PANEL_512(T) SN_ALIGNAS(64) MAT_PANEL_SQ_MIN_NOALIGN(T, __m512, 8) +#define MAT_PANEL_512(T) SN_ALIGNAS(64) MAT_PANEL_SQ_MIN_NOALIGN(T, __m512, PANEL_MINQ(__m512, T)) /* Global zero kernel declaration (defined in dllmain.c), aligned to 64 bytes */ #ifdef __cplusplus @@ -117,11 +126,14 @@ extern __declspec(align(64)) sn_zero_kernel_u64 SN_ZERO_KERNEL; int col_index, int row_index) /* Declarations for AVX float variant */ -DECLARE_BUILD_PANEL(sse, float); DECLARE_BUILD_PANEL(avx, float); DECLARE_BUILD_PANEL(512, float); -DECLARE_STORE_PANEL(sse, float); + DECLARE_STORE_PANEL(avx, float); DECLARE_STORE_PANEL(512, float); +DECLARE_BUILD_PANEL(512, double); + +DECLARE_STORE_PANEL(512, double); + #endif /* SN_MATKERNEL */ diff --git a/StgSharp.Native/sn-src/k_basicop.c b/StgSharp.Native/sn-src/k_basicop.c index 07d2ef06..c7682591 100644 --- a/StgSharp.Native/sn-src/k_basicop.c +++ b/StgSharp.Native/sn-src/k_basicop.c @@ -25,28 +25,53 @@ DECLARE_PNL_PROC_ANS(float, 512, clear) } } +DECLARE_PNL_PROC_ANS(double, sse, clear) +{ + register __m128 zero = _mm_setzero_ps(); + PANEL_COL_VEC((*ans), 0, 0) = zero; + PANEL_COL_VEC((*ans), 1, 0) = zero; + PANEL_COL_VEC((*ans), 2, 0) = zero; + PANEL_COL_VEC((*ans), 3, 0) = zero; +} + +DECLARE_PNL_PROC_ANS(double, avx, clear) +{ + register __m256 zero = _mm256_setzero_ps(); + for (size_t i = 0; i < 8; i++) { + PANEL_COL_VEC((*ans), i, 0) = zero; + } +} + +DECLARE_PNL_PROC_ANS(double, 512, clear) +{ + register __m512 zero = _mm512_setzero_ps(); + for (size_t i = 0; i < 16; i++) { + PANEL_COL_VEC((*ans), i, 0) = zero; + } +} + // Matrix transpose implementations DECLARE_KER_PROC_RIGHT_ANS(float, sse, transpose) { - register __m128 col0 = right->xmm[0]; - register __m128 col1 = right->xmm[1]; - register __m128 col2 = right->xmm[2]; - register __m128 col3 = right->xmm[3]; + register __m128 col0 = right->f32_x[0]; + register __m128 col1 = right->f32_x[1]; + register __m128 col2 = right->f32_x[2]; + register __m128 col3 = right->f32_x[3]; register __m128 tmp0 = _mm_unpacklo_ps(col0, col1); register __m128 tmp1 = _mm_unpackhi_ps(col0, col1); register __m128 tmp2 = _mm_unpacklo_ps(col2, col3); register __m128 tmp3 = _mm_unpackhi_ps(col2, col3); - ans->xmm[0] = _mm_movelh_ps(tmp0, tmp2); - ans->xmm[1] = _mm_movehl_ps(tmp2, tmp0); - ans->xmm[2] = _mm_movelh_ps(tmp1, tmp3); - ans->xmm[3] = _mm_movehl_ps(tmp3, tmp1); + ans->f32_x[0] = _mm_movelh_ps(tmp0, tmp2); + ans->f32_x[1] = _mm_movehl_ps(tmp2, tmp0); + ans->f32_x[2] = _mm_movelh_ps(tmp1, tmp3); + ans->f32_x[3] = _mm_movehl_ps(tmp3, tmp1); } DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose) { - register __m256 col_low = right->ymm[0]; - register __m256 col_high = right->ymm[1]; + register __m256 col_low = right->f32_y[0]; + register __m256 col_high = right->f32_y[1]; register __m256 unpack_low = _mm256_unpacklo_ps(col_low, col_high); register __m256 unpack_high = _mm256_unpackhi_ps(col_low, col_high); @@ -54,16 +79,16 @@ DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose) col_low = _mm256_permute2f128_ps(unpack_low, unpack_high, 0x20); col_high = _mm256_permute2f128_ps(unpack_low, unpack_high, 0x31); - ans->ymm[0] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(2, 0, 2, 0)); - ans->ymm[1] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(3, 1, 3, 1)); + ans->f32_y[0] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(2, 0, 2, 0)); + ans->f32_y[1] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(3, 1, 3, 1)); } DECLARE_KER_PROC_RIGHT_ANS(float, 512, transpose) { - register __m512 matrix = right->zmm[0]; + register __m512 matrix = right->f32_z[0]; const __m512i transpose_indices = _mm512_setr_epi32(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15); - ans->zmm[0] = _mm512_permutexvar_ps(transpose_indices, matrix); + ans->f32_z[0] = _mm512_permutexvar_ps(transpose_indices, matrix); } diff --git a/StgSharp.Native/sn-src/k_bufop.c b/StgSharp.Native/sn-src/k_bufop.c index ce3d29ad..0bfa9ca7 100644 --- a/StgSharp.Native/sn-src/k_bufop.c +++ b/StgSharp.Native/sn-src/k_bufop.c @@ -4,10 +4,10 @@ DECLARE_BUF_PROC_ANS_SCALAR(float, sse, fill) { register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]); for (size_t i = 0; i < count; i++) { - ans[i].xmm[0] = vec; - ans[i].xmm[1] = vec; - ans[i].xmm[2] = vec; - ans[i].xmm[3] = vec; + ans[i].f32_x[0] = vec; + ans[i].f32_x[1] = vec; + ans[i].f32_x[2] = vec; + ans[i].f32_x[3] = vec; } } @@ -15,8 +15,8 @@ DECLARE_BUF_PROC_ANS_SCALAR(float, avx, fill) { register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data))); for (size_t i = 0; i < count; i++) { - ans[i].ymm[0] = vec; - ans[i].ymm[1] = vec; + ans[i].f32_y[0] = vec; + ans[i].f32_y[1] = vec; } } @@ -24,22 +24,57 @@ DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill) { register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]); for (size_t i = 0; i < count; i++) { - ans[i].zmm[0] = vec; + ans[i].f32_z[0] = vec; + } +} + +DECLARE_BUF_PROC_ANS_SCALAR(double, sse, fill) +{ + register __m128d vec = _mm_set1_pd(((double *)(scalar->data))[0]); + for (size_t i = 0; i < count; i++) { + ans[i].f64_x[0] = vec; + ans[i].f64_x[1] = vec; + ans[i].f64_x[2] = vec; + ans[i].f64_x[3] = vec; + ans[i].f64_x[4] = vec; + ans[i].f64_x[5] = vec; + ans[i].f64_x[6] = vec; + ans[i].f64_x[7] = vec; + } +} + +DECLARE_BUF_PROC_ANS_SCALAR(double, avx, fill) +{ + register __m256d vec = _mm256_broadcast_sd(((double *)(scalar->data))); + for (size_t i = 0; i < count; i++) { + ans[i].f64_y[0] = vec; + ans[i].f64_y[1] = vec; + ans[i].f64_y[2] = vec; + ans[i].f64_y[3] = vec; + } +} + +DECLARE_BUF_PROC_ANS_SCALAR(double, 512, fill) +{ + register __m512d vec = _mm512_broadcastsd_pd(_mm_loadu_pd(scalar)); + for (size_t i = 0; i < count; i++) { + ans[i].f64_z[0] = vec; + ans[i].f64_z[1] = vec; } } DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, add) { for (size_t i = 0; i < count; i++) { - register __m128 c0 = _mm_add_ps(left[i].xmm[0], right[i].xmm[0]); - register __m128 c1 = _mm_add_ps(left[i].xmm[1], right[i].xmm[1]); - register __m128 c2 = _mm_add_ps(left[i].xmm[2], right[i].xmm[2]); - register __m128 c3 = _mm_add_ps(left[i].xmm[3], right[i].xmm[3]); + register __m128 c0 = _mm_add_ps(left[i].f32_x[0], right[i].f32_x[0]); + register __m128 c1 = _mm_add_ps(left[i].f32_x[1], right[i].f32_x[1]); + register __m128 c2 = _mm_add_ps(left[i].f32_x[2], right[i].f32_x[2]); + register __m128 c3 = _mm_add_ps(left[i].f32_x[3], right[i].f32_x[3]); - ans[i].xmm[0] = c0; - ans[i].xmm[1] = c1; - ans[i].xmm[2] = c2; - ans[i].xmm[3] = c3; + ans[i].f32_x[0] = c0; + ans[i].f32_x[1] = c1; + ans[i].f32_x[2] = c2; + ans[i].f32_x[3] = c3; } } @@ -47,51 +82,147 @@ DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, add) { for (size_t i = 0; i < count; i++) { size_t s = sizeof(MAT_KERNEL(float)); - register __m256 c0 = _mm256_add_ps(left[i].ymm[0], right[i].ymm[0]); - register __m256 c1 = _mm256_add_ps(left[i].ymm[1], right[i].ymm[1]); - ans[i].ymm[0] = c0; - ans[i].ymm[1] = c1; + register __m256 c0 = _mm256_add_ps(left[i].f32_y[0], right[i].f32_y[0]); + register __m256 c1 = _mm256_add_ps(left[i].f32_y[1], right[i].f32_y[1]); + ans[i].f32_y[0] = c0; + ans[i].f32_y[1] = c1; } } DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, add) { for (size_t i = 0; i < count; i++) { - register __m512 c0 = _mm512_add_ps(left[i].zmm[0], right[i].zmm[0]); - ans[i].zmm[0] = c0; + register __m512 c0 = _mm512_add_ps(left[i].f32_z[0], right[i].f32_z[0]); + ans[i].f32_z[0] = c0; + } +} +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, add) +{ + for (size_t i = 0; i < count; i++) { + register __m128d c0 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[0]), + _mm_castps_pd(right[i].f32_x[0])); + register __m128d c1 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[1]), + _mm_castps_pd(right[i].f32_x[1])); + register __m128d c2 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[2]), + _mm_castps_pd(right[i].f32_x[2])); + register __m128d c3 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[3]), + _mm_castps_pd(right[i].f32_x[3])); + register __m128d c4 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[4]), + _mm_castps_pd(right[i].f32_x[4])); + register __m128d c5 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[5]), + _mm_castps_pd(right[i].f32_x[5])); + register __m128d c6 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[6]), + _mm_castps_pd(right[i].f32_x[6])); + register __m128d c7 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[7]), + _mm_castps_pd(right[i].f32_x[7])); + + ans[i].f32_x[0] = _mm_castpd_ps(c0); + ans[i].f32_x[1] = _mm_castpd_ps(c1); + ans[i].f32_x[2] = _mm_castpd_ps(c2); + ans[i].f32_x[3] = _mm_castpd_ps(c3); + ans[i].f32_x[4] = _mm_castpd_ps(c4); + ans[i].f32_x[5] = _mm_castpd_ps(c5); + ans[i].f32_x[6] = _mm_castpd_ps(c6); + ans[i].f32_x[7] = _mm_castpd_ps(c7); + } +} + +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, add) +{ + for (size_t i = 0; i < count; i++) { + size_t s = sizeof(MAT_KERNEL(float)); + register __m256d c0 = _mm256_add_pd(left[i].f64_y[0], right[i].f64_y[0]); + register __m256d c1 = _mm256_add_pd(left[i].f64_y[1], right[i].f64_y[1]); + register __m256d c2 = _mm256_add_pd(left[i].f64_y[0], right[i].f64_y[0]); + register __m256d c3 = _mm256_add_pd(left[i].f64_y[1], right[i].f64_y[1]); + ans[i].f64_y[0] = c0; + ans[i].f64_y[1] = c1; + ans[i].f64_y[2] = c2; + ans[i].f64_y[3] = c3; + } +} + +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, add) +{ + for (size_t i = 0; i < count; i++) { + register __m512d c0 = _mm512_add_pd(left[i].f64_z[0], right[i].f64_z[0]); + register __m512d c1 = _mm512_add_pd(left[i].f64_z[1], right[i].f64_z[1]); + ans[i].f64_z[0] = c0; + ans[i].f64_z[1] = c1; } } DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, sub) { for (size_t i = 0; i < count; i++) { - register __m128 c0 = _mm_sub_ps(left[i].xmm[0], right[i].xmm[0]); - register __m128 c1 = _mm_sub_ps(left[i].xmm[1], right[i].xmm[1]); - register __m128 c2 = _mm_sub_ps(left[i].xmm[2], right[i].xmm[2]); - register __m128 c3 = _mm_sub_ps(left[i].xmm[3], right[i].xmm[3]); + register __m128 c0 = _mm_sub_ps(left[i].f32_x[0], right[i].f32_x[0]); + register __m128 c1 = _mm_sub_ps(left[i].f32_x[1], right[i].f32_x[1]); + register __m128 c2 = _mm_sub_ps(left[i].f32_x[2], right[i].f32_x[2]); + register __m128 c3 = _mm_sub_ps(left[i].f32_x[3], right[i].f32_x[3]); - ans[i].xmm[0] = c0; - ans[i].xmm[1] = c1; - ans[i].xmm[2] = c2; - ans[i].xmm[3] = c3; + ans[i].f32_x[0] = c0; + ans[i].f32_x[1] = c1; + ans[i].f32_x[2] = c2; + ans[i].f32_x[3] = c3; } } DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, sub) { for (size_t i = 0; i < count; i++) { - register __m256 c0 = _mm256_sub_ps(left[i].ymm[0], right[i].ymm[0]); - register __m256 c1 = _mm256_sub_ps(left[i].ymm[1], right[i].ymm[1]); - ans[i].ymm[0] = c0; - ans[i].ymm[1] = c1; + register __m256 c0 = _mm256_sub_ps(left[i].f32_y[0], right[i].f32_y[0]); + register __m256 c1 = _mm256_sub_ps(left[i].f32_y[1], right[i].f32_y[1]); + ans[i].f32_y[0] = c0; + ans[i].f32_y[1] = c1; } } DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, sub) { for (size_t i = 0; i < count; i++) { - register __m512 c0 = _mm512_sub_ps(left[i].zmm[0], right[i].zmm[0]); - ans[i].zmm[0] = c0; + register __m512 c0 = _mm512_sub_ps(left[i].f32_z[0], right[i].f32_z[0]); + ans[i].f32_z[0] = c0; + } +} + +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, sub) +{ + for (size_t i = 0; i < count; i++) { + register __m128d c0 = _mm_sub_pd(left[i].f64_x[0], right[i].f64_x[0]); + register __m128d c1 = _mm_sub_pd(left[i].f64_x[1], right[i].f64_x[1]); + register __m128d c2 = _mm_sub_pd(left[i].f64_x[2], right[i].f64_x[2]); + register __m128d c3 = _mm_sub_pd(left[i].f64_x[3], right[i].f64_x[3]); + register __m128d c4 = _mm_sub_pd(left[i].f64_x[4], right[i].f64_x[4]); + register __m128d c5 = _mm_sub_pd(left[i].f64_x[5], right[i].f64_x[5]); + register __m128d c6 = _mm_sub_pd(left[i].f64_x[6], right[i].f64_x[6]); + register __m128d c7 = _mm_sub_pd(left[i].f64_x[7], right[i].f64_x[7]); + + ans[i].f64_x[0] = c0; + ans[i].f64_x[1] = c1; + ans[i].f64_x[2] = c2; + ans[i].f64_x[3] = c3; + ans[i].f64_x[4] = c4; + ans[i].f64_x[5] = c5; + ans[i].f64_x[6] = c6; + ans[i].f64_x[7] = c7; + } +} + +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, sub) +{ + for (size_t i = 0; i < count; i++) { + register __m256 c0 = _mm256_sub_ps(left[i].f32_y[0], right[i].f32_y[0]); + register __m256 c1 = _mm256_sub_ps(left[i].f32_y[1], right[i].f32_y[1]); + ans[i].f32_y[0] = c0; + ans[i].f32_y[1] = c1; + } +} + +DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, sub) +{ + for (size_t i = 0; i < count; i++) { + register __m512 c0 = _mm512_sub_ps(left[i].f32_z[0], right[i].f32_z[0]); + ans[i].f32_z[0] = c0; } } @@ -99,15 +230,15 @@ DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul) { register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]); for (size_t i = 0; i < count; i++) { - register __m128 c0 = _mm_mul_ps(right[i].xmm[0], vec); - register __m128 c1 = _mm_mul_ps(right[i].xmm[1], vec); - register __m128 c2 = _mm_mul_ps(right[i].xmm[2], vec); - register __m128 c3 = _mm_mul_ps(right[i].xmm[3], vec); + register __m128 c0 = _mm_mul_ps(right[i].f32_x[0], vec); + register __m128 c1 = _mm_mul_ps(right[i].f32_x[1], vec); + register __m128 c2 = _mm_mul_ps(right[i].f32_x[2], vec); + register __m128 c3 = _mm_mul_ps(right[i].f32_x[3], vec); - ans[i].xmm[0] = c0; - ans[i].xmm[1] = c1; - ans[i].xmm[2] = c2; - ans[i].xmm[3] = c3; + ans[i].f32_x[0] = c0; + ans[i].f32_x[1] = c1; + ans[i].f32_x[2] = c2; + ans[i].f32_x[3] = c3; } } @@ -115,11 +246,11 @@ DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul) { register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data))); for (size_t i = 0; i < count; i++) { - register __m256 c0 = _mm256_mul_ps(right[i].ymm[0], vec); - register __m256 c1 = _mm256_mul_ps(right[i].ymm[1], vec); + register __m256 c0 = _mm256_mul_ps(right[i].f32_y[0], vec); + register __m256 c1 = _mm256_mul_ps(right[i].f32_y[1], vec); - ans[i].ymm[0] = c0; - ans[i].ymm[1] = c1; + ans[i].f32_y[0] = c0; + ans[i].f32_y[1] = c1; } } @@ -127,8 +258,60 @@ DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul) { register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]); for (size_t i = 0; i < count; i++) { - register __m512 c0 = _mm512_mul_ps(right[i].zmm[0], vec); + register __m512 c0 = _mm512_mul_ps(right[i].f32_z[0], vec); - ans[i].zmm[0] = c0; + ans[i].f32_z[0] = c0; + } +} + +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, sse, scalar_mul) +{ + register __m128d vec = _mm_set1_pd(((double *)(scalar->data))[0]); + for (size_t i = 0; i < count; i++) { + register __m128d c0 = _mm_mul_pd(right[i].f64_x[0], vec); + register __m128d c1 = _mm_mul_pd(right[i].f64_x[1], vec); + register __m128d c2 = _mm_mul_pd(right[i].f64_x[2], vec); + register __m128d c3 = _mm_mul_pd(right[i].f64_x[3], vec); + register __m128d c4 = _mm_mul_pd(right[i].f64_x[4], vec); + register __m128d c5 = _mm_mul_pd(right[i].f64_x[5], vec); + register __m128d c6 = _mm_mul_pd(right[i].f64_x[6], vec); + register __m128d c7 = _mm_mul_pd(right[i].f64_x[7], vec); + + ans[i].f64_x[0] = c0; + ans[i].f64_x[1] = c1; + ans[i].f64_x[2] = c2; + ans[i].f64_x[3] = c3; + ans[i].f64_x[4] = c4; + ans[i].f64_x[5] = c5; + ans[i].f64_x[6] = c6; + ans[i].f64_x[7] = c7; + } +} + +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, avx, scalar_mul) +{ + register __m256d vec = _mm256_broadcast_sd(((double *)(scalar->data))); + for (size_t i = 0; i < count; i++) { + register __m256d c0 = _mm256_mul_pd(right[i].f64_y[0], vec); + register __m256d c1 = _mm256_mul_pd(right[i].f64_y[1], vec); + register __m256d c2 = _mm256_mul_pd(right[i].f64_y[2], vec); + register __m256d c3 = _mm256_mul_pd(right[i].f64_y[3], vec); + + ans[i].f64_y[0] = c0; + ans[i].f64_y[1] = c1; + ans[i].f64_y[2] = c2; + ans[i].f64_y[3] = c3; + } +} + +DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, 512, scalar_mul) +{ + register __m512d vec = _mm512_broadcastsd_pd(_mm_load_pd(scalar)); + for (size_t i = 0; i < count; i++) { + register __m512d c0 = _mm512_mul_pd(right[i].f64_z[0], vec); + register __m512d c1 = _mm512_mul_pd(right[i].f64_z[1], vec); + + ans[i].f64_z[0] = c0; + ans[i].f64_z[1] = c1; } } diff --git a/StgSharp.Native/sn-src/k_panel.c b/StgSharp.Native/sn-src/k_panel.c index 916a6f23..aae9ff26 100644 --- a/StgSharp.Native/sn-src/k_panel.c +++ b/StgSharp.Native/sn-src/k_panel.c @@ -31,8 +31,8 @@ DECLARE_BUILD_PANEL(avx, float) const MAT_KERNEL(float) *src_bot = is_right_half ? k_rb : k_lb; /* Vector loads: each kernel column is a compact __m128 */ - top128 = src_top->xmm[kc_col]; - bot128 = src_bot->xmm[kc_col]; + top128 = src_top->f32_x[kc_col]; + bot128 = src_bot->f32_x[kc_col]; /* Combine into an AVX column: low 128 = top, high 128 = bottom */ col256 = _mm256_set_m128(bot128, top128); @@ -41,15 +41,6 @@ DECLARE_BUILD_PANEL(avx, float) } } -DECLARE_BUILD_PANEL(sse, float) -{ - /* Single kernel -> 4 columns, each is one __m128 column vector. */ - const MAT_KERNEL(float) *src = matrix + (row_index * col_length + col_index); - - for (int pc = 0; pc < 4; ++pc) { - PANEL_COL_VEC((*panel), pc, 0) = src->xmm[pc]; - } -} DECLARE_BUILD_PANEL(512, float) { @@ -71,10 +62,10 @@ DECLARE_BUILD_PANEL(512, float) const int blk_col = pc >> 2; /* which kernel column block (0..3) */ const int kc_col = pc & 3; /* column inside a 4x4 kernel */ - const __m128 r0 = k[blk_col][0]->xmm[kc_col]; - const __m128 r1 = k[blk_col][1]->xmm[kc_col]; - const __m128 r2 = k[blk_col][2]->xmm[kc_col]; - const __m128 r3 = k[blk_col][3]->xmm[kc_col]; + const __m128 r0 = k[blk_col][0]->f32_x[kc_col]; + const __m128 r1 = k[blk_col][1]->f32_x[kc_col]; + const __m128 r2 = k[blk_col][2]->f32_x[kc_col]; + const __m128 r3 = k[blk_col][3]->f32_x[kc_col]; const __m256 top = _mm256_set_m128(r1, r0); const __m256 bot = _mm256_set_m128(r3, r2); @@ -86,6 +77,43 @@ DECLARE_BUILD_PANEL(512, float) } } + +DECLARE_BUILD_PANEL(512, double) +{ + const int Q = 8; /* max(4, lanes) for double with ZMM -> 8 columns */ + + const int has_right = (col_index + 1) < col_length; + const int has_bottom = (row_index + 1) < row_length; + + const MAT_KERNEL(double) *const k_lt = matrix + (row_index * col_length + col_index); + const MAT_KERNEL(double) *const k_rt = + has_right ? (matrix + (row_index * col_length + (col_index + 1))) : + (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL; + const MAT_KERNEL(double) *const k_lb = + has_bottom ? (matrix + ((row_index + 1) * col_length + col_index)) : + (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL; + const MAT_KERNEL(double) *const k_rb = + (has_right && has_bottom) ? + (matrix + ((row_index + 1) * col_length + (col_index + 1))) : + (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL; + + for (int pc = 0; pc < Q; ++pc) { + const int is_right_half = (pc >= 4); + const int kc_col = pc & 3; /* column index inside 4x4 kernel */ + + const MAT_KERNEL(double) *src_top = is_right_half ? k_rt : k_lt; + const MAT_KERNEL(double) *src_bot = is_right_half ? k_rb : k_lb; + + const __m256d top256 = src_top->f64_y[kc_col]; + const __m256d bot256 = src_bot->f64_y[kc_col]; + + __m512d col = _mm512_castpd256_pd512(top256); + col = _mm512_insertf64x4(col, bot256, 1); + + PANEL_COL_VEC((*panel), pc, 0) = _mm512_castpd_ps(col); + } +} + DECLARE_STORE_PANEL(avx, float) { const int Q = 8; /* fixed for AVX float panel */ @@ -124,8 +152,8 @@ DECLARE_STORE_PANEL(avx, float) bot128 = _mm256_extractf128_ps(col256, 1); /* high 128 */ /* Store back into kernels directly (out-of-bounds writes go to dummy) */ - k_rt->xmm[kc_col] = top128; - k_rb->xmm[kc_col] = bot128; + k_rt->f32_x[kc_col] = top128; + k_rb->f32_x[kc_col] = bot128; } for (; i < Q; i++) { const int kc_col = i - 4; /* column index inside 4x4 kernel */ @@ -137,20 +165,11 @@ DECLARE_STORE_PANEL(avx, float) bot128 = _mm256_extractf128_ps(col256, 1); /* high 128 */ /* Store back into kernels directly (out-of-bounds writes go to dummy) */ - k_lt->xmm[kc_col] = top128; - k_lb->xmm[kc_col] = bot128; + k_lt->f32_x[kc_col] = top128; + k_lb->f32_x[kc_col] = bot128; } } -DECLARE_STORE_PANEL(sse, float) -{ - MAT_KERNEL(float) *dst = matrix + (row_index * col_length + col_index); - - for (int pc = 0; pc < 4; ++pc) { - const __m128 col = PANEL_COL_VEC((*panel), pc, 0); - dst->xmm[pc] = col; - } -} DECLARE_STORE_PANEL(512, float) { @@ -181,9 +200,99 @@ DECLARE_STORE_PANEL(512, float) const __m128 r2 = _mm256_castps256_ps128(bot); const __m128 r3 = _mm256_extractf128_ps(bot, 1); - k[blk_col][0]->xmm[kc_col] = r0; - k[blk_col][1]->xmm[kc_col] = r1; - k[blk_col][2]->xmm[kc_col] = r2; - k[blk_col][3]->xmm[kc_col] = r3; + k[blk_col][0]->f32_x[kc_col] = r0; + k[blk_col][1]->f32_x[kc_col] = r1; + k[blk_col][2]->f32_x[kc_col] = r2; + k[blk_col][3]->f32_x[kc_col] = r3; + } +} + +DECLARE_STORE_PANEL(sse, double) +{ + MAT_KERNEL(double) *dst = matrix + (row_index * col_length + col_index); + + for (int pc = 0; pc < 4; ++pc) { + const int base = pc * 2; + dst->f64_x[base] = _mm_castps_pd(PANEL_COL_VEC((*panel), pc, 0)); + dst->f64_x[base + 1] = _mm_castps_pd(PANEL_COL_VEC((*panel), pc, 1)); + } +} + +DECLARE_STORE_PANEL(avx, double) +{ + const int Q = 8; /* fixed for AVX double panel */ + + MAT_KERNEL(double) zero; + + const int has_right = (col_index + 1) < col_length; + const int has_bottom = (row_index + 1) < row_length; + + MAT_KERNEL(double) *const k_lt = (col_index < col_length && row_index < row_length) ? + (matrix + (row_index * col_length + col_index)) : + &zero; + MAT_KERNEL(double) *const k_rt = + has_right && (row_index < row_length) ? + (matrix + (row_index * col_length + (col_index + 1))) : + &zero; + MAT_KERNEL(double) *const k_lb = + has_bottom && (col_index < col_length) ? + (matrix + ((row_index + 1) * col_length + col_index)) : + &zero; + MAT_KERNEL(double) *const k_rb = + (has_right && has_bottom) ? + (matrix + ((row_index + 1) * col_length + (col_index + 1))) : + &zero; + + for (int i = 0; i < Q; ++i) { + const int kc_col = (i < 4) ? i : (i - 4); + const __m256 top_chunk = PANEL_COL_VEC((*panel), i, 0); + const __m256 bot_chunk = PANEL_COL_VEC((*panel), i, 1); + + MAT_KERNEL(double) *const top_dst = (i < 4) ? k_rt : k_lt; + MAT_KERNEL(double) *const bot_dst = (i < 4) ? k_rb : k_lb; + + top_dst->f64_y[kc_col] = _mm256_castps_pd(top_chunk); + bot_dst->f64_y[kc_col] = _mm256_castps_pd(bot_chunk); + } +} + +DECLARE_STORE_PANEL(512, double) +{ + const int Q = 8; /* max(4, lanes) -> 8 columns, one ZMM per column */ + + MAT_KERNEL(double) zero; + + const int has_right = (col_index + 1) < col_length; + const int has_bottom = (row_index + 1) < row_length; + + MAT_KERNEL(double) *const k_lt = (col_index < col_length && row_index < row_length) ? + (matrix + (row_index * col_length + col_index)) : + &zero; + MAT_KERNEL(double) *const k_rt = + has_right && (row_index < row_length) ? + (matrix + (row_index * col_length + (col_index + 1))) : + &zero; + MAT_KERNEL(double) *const k_lb = + has_bottom && (col_index < col_length) ? + (matrix + ((row_index + 1) * col_length + col_index)) : + &zero; + MAT_KERNEL(double) *const k_rb = + (has_right && has_bottom) ? + (matrix + ((row_index + 1) * col_length + (col_index + 1))) : + &zero; + + for (int i = 0; i < Q; ++i) { + const int kc_col = (i < 4) ? i : (i - 4); + const __m512 col = PANEL_COL_VEC((*panel), i, 0); + const __m512d col_pd = _mm512_castps_pd(col); + + const __m256d top256 = _mm512_castpd512_pd256(col_pd); + const __m256d bot256 = _mm512_extractf64x4_pd(col_pd, 1); + + MAT_KERNEL(double) *const top_dst = (i < 4) ? k_rt : k_lt; + MAT_KERNEL(double) *const bot_dst = (i < 4) ? k_rb : k_lb; + + top_dst->f64_y[kc_col] = top256; + bot_dst->f64_y[kc_col] = bot256; } } diff --git a/StgSharp.Native/sn-src/k_pivot.c b/StgSharp.Native/sn-src/k_pivot.c index 3abf4c66..cc6a1733 100644 --- a/StgSharp.Native/sn-src/k_pivot.c +++ b/StgSharp.Native/sn-src/k_pivot.c @@ -2,8 +2,8 @@ // Hybrid implementation: SIMD for Abs, Scalar for Max // Assumes Column-Major layout for both MAT_KERNEL content and Kernel sequence -static void pivot_float(MAT_KERNEL(float) const *source, size_t col_begin, size_t col_length, - size_t *restrict ans) +INTERNAL void SN_DECL pivot_float(MAT_KERNEL(float) const *source, size_t col_begin, + size_t col_length, size_t *restrict ans) { float max_v[4] = { -1.0f, -1.0f, -1.0f, -1.0f }; diff --git a/StgSharp.Native/sn-src/load_intrinsic.c b/StgSharp.Native/sn-src/load_intrinsic.c index 22b93f11..949bc51c 100644 --- a/StgSharp.Native/sn-src/load_intrinsic.c +++ b/StgSharp.Native/sn-src/load_intrinsic.c @@ -65,17 +65,30 @@ void load_intrinsic_sse(sn_intrinsic *context) context->city_hash_simplify = city_hash_simplify_sse; context->factorial_simd = factorial_simd_sse; context->index_pair = index_pair_sse; - context->f32_store_panel = STORE_PANEL(sse, float); - context->f32_clear_panel = PNL_PROC(float, sse, clear); - context->f32_build_panel = BUILD_PANEL(sse, float); - context->f32_buffer_add = BUF_PROC(float, sse, add); - context->f32_panel_fma = PNL_PROC(float, sse, fma); - context->f32_buffer_fill = BUF_PROC(float, sse, fill); - context->variant = BUF_PROC(float, sse, scalar_mul); - context->f32_buffer_sub = BUF_PROC(float, sse, sub); - context->f32_buffer_transpose = KER_PROC(float, sse, transpose); - context->f32_buffer_scalar_mul = BUF_PROC(float, sse, scalar_mul); - context->f32_normalize_3 = f32_normalize; + + context->mat_f32_buffer_add = BUF_PROC(float, sse, add); + context->mat_f32_buffer_fill = BUF_PROC(float, sse, fill); + context->mat_f32_buffer_scalar_mul = BUF_PROC(float, sse, scalar_mul); + context->mat_f32_buffer_sub = BUF_PROC(float, sse, sub); + context->mat_f32_buffer_transpose = KER_PROC(float, sse, transpose); + context->mat_f32_build_panel = NULL; + context->mat_f32_clear_panel = PNL_PROC(float, sse, clear); + context->mat_f32_panel_fma = PNL_PROC(float, sse, fma); + context->mat_f32_pivot = pivot_float; + context->mat_f32_store_panel = NULL; + + context->mat_f64_buffer_add = NULL; + context->mat_f64_buffer_fill = NULL; + context->mat_f64_buffer_scalar_mul = NULL; + context->mat_f64_buffer_sub = NULL; + context->mat_f64_buffer_transpose = NULL; + context->mat_f64_build_panel = NULL; + context->mat_f64_clear_panel = NULL; + context->mat_f64_panel_fma = NULL; + context->mat_f64_pivot = NULL; + context->mat_f64_store_panel = NULL; + + context->vec_f32_normalize_3 = f32_normalize; } void load_intrinsic_avx(sn_intrinsic *context) @@ -83,17 +96,30 @@ void load_intrinsic_avx(sn_intrinsic *context) context->city_hash_simplify = city_hash_simplify_sse; context->factorial_simd = factorial_simd_sse; context->index_pair = index_pair_sse; - context->f32_store_panel = STORE_PANEL(avx, float); - context->f32_clear_panel = PNL_PROC(float, avx, clear); - context->f32_build_panel = BUILD_PANEL(avx, float); - context->f32_buffer_add = BUF_PROC(float, avx, add); - context->f32_panel_fma = PNL_PROC(float, avx, fma); - context->f32_buffer_fill = BUF_PROC(float, avx, fill); - context->variant = BUF_PROC(float, avx, scalar_mul); - context->f32_buffer_sub = BUF_PROC(float, avx, sub); - context->f32_buffer_transpose = KER_PROC(float, avx, transpose); - context->f32_buffer_scalar_mul = BUF_PROC(float, avx, scalar_mul); - context->f32_normalize_3 = f32_normalize; + + context->mat_f32_buffer_add = BUF_PROC(float, avx, add); + context->mat_f32_buffer_fill = BUF_PROC(float, avx, fill); + context->mat_f32_buffer_scalar_mul = BUF_PROC(float, avx, scalar_mul); + context->mat_f32_buffer_sub = BUF_PROC(float, avx, sub); + context->mat_f32_buffer_transpose = KER_PROC(float, avx, transpose); + context->mat_f32_build_panel = BUILD_PANEL(avx, float); + context->mat_f32_clear_panel = PNL_PROC(float, avx, clear); + context->mat_f32_panel_fma = PNL_PROC(float, avx, fma); + context->mat_f32_pivot = pivot_float; + context->mat_f32_store_panel = STORE_PANEL(avx, float); + + context->mat_f64_buffer_add = BUF_PROC(double, avx, add); + context->mat_f64_buffer_fill = BUF_PROC(double, avx, fill); + context->mat_f64_buffer_scalar_mul = BUF_PROC(double, avx, scalar_mul); + context->mat_f64_buffer_sub = BUF_PROC(double, avx, sub); + context->mat_f64_buffer_transpose = NULL; + context->mat_f64_build_panel = NULL; + context->mat_f64_clear_panel = PNL_PROC(double, avx, clear); + context->mat_f64_panel_fma = NULL; + context->mat_f64_pivot = NULL; + context->mat_f64_store_panel = NULL; + + context->vec_f32_normalize_3 = f32_normalize; } void load_intrinsic_512(sn_intrinsic *context) @@ -101,15 +127,28 @@ void load_intrinsic_512(sn_intrinsic *context) context->city_hash_simplify = city_hash_simplify_sse; context->factorial_simd = factorial_simd_sse; context->index_pair = index_pair_sse; - context->f32_store_panel = STORE_PANEL(512, float); - context->f32_clear_panel = PNL_PROC(float, 512, clear); - context->f32_build_panel = BUILD_PANEL(512, float); - context->f32_buffer_add = BUF_PROC(float, 512, add); - context->f32_panel_fma = PNL_PROC(float, 512, fma); - context->f32_buffer_fill = BUF_PROC(float, 512, fill); - context->variant = BUF_PROC(float, 512, scalar_mul); - context->f32_buffer_sub = BUF_PROC(float, 512, sub); - context->f32_buffer_transpose = KER_PROC(float, 512, transpose); - context->f32_buffer_scalar_mul = BUF_PROC(float, 512, scalar_mul); - context->f32_normalize_3 = f32_normalize; + + context->mat_f32_buffer_add = BUF_PROC(float, 512, add); + context->mat_f32_buffer_fill = BUF_PROC(float, 512, fill); + context->mat_f32_buffer_scalar_mul = BUF_PROC(float, 512, scalar_mul); + context->mat_f32_buffer_sub = BUF_PROC(float, 512, sub); + context->mat_f32_buffer_transpose = KER_PROC(float, 512, transpose); + context->mat_f32_build_panel = BUILD_PANEL(512, float); + context->mat_f32_clear_panel = PNL_PROC(float, 512, clear); + context->mat_f32_panel_fma = PNL_PROC(float, 512, fma); + context->mat_f32_pivot = pivot_float; + context->mat_f32_store_panel = STORE_PANEL(512, float); + + context->mat_f64_buffer_add = NULL; + context->mat_f64_buffer_fill = NULL; + context->mat_f64_buffer_scalar_mul = NULL; + context->mat_f64_buffer_sub = NULL; + context->mat_f64_buffer_transpose = NULL; + context->mat_f64_build_panel = NULL; + context->mat_f64_clear_panel = NULL; + context->mat_f64_panel_fma = NULL; + context->mat_f64_pivot = NULL; + context->mat_f64_store_panel = NULL; + + context->vec_f32_normalize_3 = f32_normalize; } diff --git a/StgSharp.Native/sn-src/vector.c b/StgSharp.Native/sn-src/vector.c index 3cacaaca..60bae8b3 100644 --- a/StgSharp.Native/sn-src/vector.c +++ b/StgSharp.Native/sn-src/vector.c @@ -3,7 +3,7 @@ INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *restrict target) { - register __m128 s = _mm_load_ps(&source->xmm); + register __m128 s = _mm_load_ps(&source->f32_x); register __m128 a_length_sq = _mm_mul_ps(s, s); a_length_sq = _mm_add_ps(a_length_sq, _mm_shuffle_ps(a_length_sq, a_length_sq, _MM_SHUFFLE(2, 3, 0, 1))); @@ -13,13 +13,13 @@ INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *res __m128 a_inv_length = _mm_rsqrt_ps(a_length_sq); s = _mm_mul_ps(s, a_inv_length); - __m128 dot_product = _mm_mul_ps(target->xmm[0], s); + __m128 dot_product = _mm_mul_ps(target->f32_x[0], s); dot_product = _mm_add_ps(dot_product, _mm_shuffle_ps(dot_product, dot_product, _MM_SHUFFLE(2, 3, 0, 1))); dot_product = _mm_add_ss(dot_product, _mm_shuffle_ps(dot_product, dot_product, _MM_SHUFFLE(0, 0, 3, 2))); - __m128 b_orthogonal = _mm_sub_ps(target->xmm[0], _mm_mul_ps(s, dot_product)); + __m128 b_orthogonal = _mm_sub_ps(target->f32_x[0], _mm_mul_ps(s, dot_product)); __m128 b_length_sq = _mm_mul_ps(b_orthogonal, b_orthogonal); b_length_sq = _mm_add_ps(b_length_sq, @@ -28,5 +28,5 @@ INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *res _mm_shuffle_ps(b_length_sq, b_length_sq, _MM_SHUFFLE(0, 0, 3, 2))); __m128 b_inv_length = _mm_rsqrt_ps(b_length_sq); - target->xmm[0] = _mm_mul_ps(b_orthogonal, b_inv_length); + target->f32_x[0] = _mm_mul_ps(b_orthogonal, b_inv_length); } diff --git a/StgSharp.Native/sn-src/vector_dot.c b/StgSharp.Native/sn-src/vector_dot.c index 7b18cf27..6b7ca295 100644 --- a/StgSharp.Native/sn-src/vector_dot.c +++ b/StgSharp.Native/sn-src/vector_dot.c @@ -6,12 +6,12 @@ static const int avxdot_indices[] = { 0, 4, 1, 5, 2, 3, 6, 7 }; INTERNAL void FORCEINLINE SN_DECL dot_41_sse(MAT_KERNEL(float) const *transpose, __m128 const *vector, __m128 *ans) { - *ans = _mm_dp_ps(transpose->xmm[0], *vector, 0B11110001); - *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], *vector, 0B11110001), + *ans = _mm_dp_ps(transpose->f32_x[0], *vector, 0B11110001); + *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], *vector, 0B11110001), SSE_INSERT_SIMPLE(1, 1)); - *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], *vector, 0B11110001), + *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], *vector, 0B11110001), SSE_INSERT_SIMPLE(1, 3)); - *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], *vector, 0B11110001), + *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], *vector, 0B11110001), SSE_INSERT_SIMPLE(1, 4)); } @@ -19,8 +19,8 @@ INTERNAL void SN_DECL dot_41_avx(MAT_KERNEL(float) const *transpose, __m128 cons __m128 *ans) { __m256 _double = _mm256_set_m128(vector[0], vector[0]); - __m256 tmp = _mm256_dp_ps(transpose->ymm[0], _double, 0B11110001); - tmp = _mm256_blend_ps(tmp, _mm256_dp_ps(transpose->ymm[0], _double, 0B11110010), + __m256 tmp = _mm256_dp_ps(transpose->f32_y[0], _double, 0B11110001); + tmp = _mm256_blend_ps(tmp, _mm256_dp_ps(transpose->f32_y[0], _double, 0B11110010), 0B00100010); *ans = _mm256_castps256_ps128( _mm256_permutevar8x32_ps(tmp, _mm256_loadu_epi32(avxdot_indices))); @@ -29,34 +29,34 @@ INTERNAL void SN_DECL dot_41_avx(MAT_KERNEL(float) const *transpose, __m128 cons INTERNAL void SN_DECL dot_42_sse(MAT_KERNEL(float) const *transpose, __m128 const *vector, __m128 *ans) { - ans[0] = _mm_dp_ps(transpose->xmm[0], vector[0], 0B11110001); - ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], vector[0], 0B11110001), + ans[0] = _mm_dp_ps(transpose->f32_x[0], vector[0], 0B11110001); + ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], vector[0], 0B11110001), SSE_INSERT_SIMPLE(1, 2)); - ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], vector[0], 0B11110001), + ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], vector[0], 0B11110001), SSE_INSERT_SIMPLE(1, 3)); - ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], vector[0], 0B11110001), + ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], vector[0], 0B11110001), SSE_INSERT_SIMPLE(1, 4)); - ans[1] = _mm_dp_ps(transpose->xmm[0], vector[1], 0B11110001); - ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], vector[1], 0B11110001), + ans[1] = _mm_dp_ps(transpose->f32_x[0], vector[1], 0B11110001); + ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], vector[1], 0B11110001), SSE_INSERT_SIMPLE(1, 2)); - ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], vector[1], 0B11110001), + ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], vector[1], 0B11110001), SSE_INSERT_SIMPLE(1, 3)); - ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], vector[1], 0B11110001), + ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], vector[1], 0B11110001), SSE_INSERT_SIMPLE(1, 4)); } INTERNAL void SN_DECL dot_42_avx(MAT_KERNEL(float) const *transpose, __m128 const *vector, __m128 *ans) { - __m256 temp = transpose->ymm[0]; + __m256 temp = transpose->f32_y[0]; __m256 vector_stream = _mm256_set_m128(vector[1], vector[0]); __m256 result = _mm256_dp_ps(temp, vector_stream, 0B11110001); temp = _mm256_permute2f128_ps(temp, temp, 0B00000001); result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11110010), 0B00100010); - temp = transpose->ymm[1]; + temp = transpose->f32_y[1]; result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11110100), 0B01000100); temp = _mm256_permute2f128_ps(temp, temp, 0B00000001); result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11111000), 0B10001000);