From 6827259891b7c80b8a97f4442002d150764109b7 Mon Sep 17 00:00:00 2001 From: Youmu Liu <105042599+Nitload-NSI@users.noreply.github.com> Date: Fri, 9 Jan 2026 18:13:46 +0800 Subject: [PATCH] tile fma added --- .../Internal/Intrinsic/IntrinsicContext.cs | 1 + StgSharp.Mathematics/Numeric/MatrixCompute.cs | 6 +- .../Numeric/MatrixCompute.double.cs | 57 +- StgSharp.Native/StgSharp.Native.vcxproj | 1 + .../StgSharp.Native.vcxproj.filters | 9 +- StgSharp.Native/sn-head/sn_internal.h | 3 + .../sn-head/sn_intrincontext.x86.h | 10 + StgSharp.Native/sn-head/sn_intrinsic.h | 35 +- StgSharp.Native/sn-src/load_intrinsic.x86.c | 24 +- StgSharp.Native/sn-src/m.fma.x86.c | 260 ++------ StgSharp.Native/sn-src/m.fma_tile.x86.c | 583 ++++++++++++++++++ StgSharpDebug/Program.cs | 12 +- 12 files changed, 736 insertions(+), 265 deletions(-) create mode 100644 StgSharp.Native/sn-src/m.fma_tile.x86.c diff --git a/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs b/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs index 93f5020c..c11f0fc2 100644 --- a/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs +++ b/StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs @@ -46,6 +46,7 @@ namespace StgSharp.Internal.Intrinsic public delegate* unmanaged[Cdecl] build_panel; //f32_build_panel public delegate* unmanaged[Cdecl] clear_panel; //f32_clear_panel public delegate* unmanaged[Cdecl] kernel_fma; //f32_panel_fma + public delegate* unmanaged[Cdecl] kernel_tile_fma; //f32_panel_fma public delegate* unmanaged[Cdecl] pivot; //f32_panel_mul public delegate* unmanaged[Cdecl] store_panel; //f32_store_panel diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.cs index 902046e8..19a13abd 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.cs @@ -95,9 +95,6 @@ namespace StgSharp.Mathematics.Numeric { int i = (int)(c / ansWidth); int j = (int)(c % ansWidth); - MatrixKernel* ansPanel = GetKernelAddressUnsafe((MatrixKernel*)ansBuffer, - ansHeight, i, j); - NativeIntrinsic.Context.mat[f64].clear_panel((MatrixPanel*)ansPanel); for (int k = 0; k < commonK; k += 1) { NativeIntrinsic.Context .mat[f64].kernel_fma( @@ -105,7 +102,8 @@ namespace StgSharp.Mathematics.Numeric ansHeight, k, j), (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)rightBuffer, commonK, j, k), - (MatrixPanel*)ansPanel); + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)ansBuffer, + ansHeight, i, j)); } } } diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs index b782278d..e26fe7e7 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.double.cs @@ -146,64 +146,19 @@ namespace StgSharp.Mathematics.Numeric if (count > 256) { pool = MatrixParallel.LeadParallel(1); } - if (NativeIntrinsic.Context.mat[f64].build_panel == null) + for (int i = 0; i < ans.KernelRowLength; i += size) { - for (int i = 0; i < ans.KernelRowLength; i += size) + for (int j = 0; j < ans.KernelColumnLength; j += size) { - for (int j = 0; j < ans.KernelColumnLength; j += size) + // Console.WriteLine($"{i}{j}"); + for (int k = 0; k < left.KernelColumnLength; k += size) { - MatrixKernel* ansPanel = ans.GetKernelAddressUnsafe(i, j); - NativeIntrinsic.Context.mat[f64].clear_panel((MatrixPanel*)ansPanel); - - // Console.WriteLine($"{i}{j}"); - for (int k = 0; k < left.KernelColumnLength; k += size) - { - NativeIntrinsic.Context + NativeIntrinsic.Context .mat[f64].kernel_fma((MatrixPanel*)left.GetKernelAddressUnsafe(i, k), (MatrixPanel*)right.GetKernelAddressUnsafe(k, j), - (MatrixPanel*)ansPanel); - } + (MatrixPanel*)ans.GetKernelAddressUnsafe(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[f64].clear_panel((MatrixPanel*)ansPanel); - - // Console.WriteLine($"{i}{j}"); - for (int k = 0; k < left.KernelColumnLength; k += size) - { - NativeIntrinsic.Context - .mat[f64].build_panel((MatrixPanel*)leftPanel, - (MatrixKernel*)left.Buffer, - left.KernelColumnLength, - left.KernelRowLength, i, k); - NativeIntrinsic.Context - .mat[f64].build_panel((MatrixPanel*)rightPanel, - (MatrixKernel*)right.Buffer, - right.KernelColumnLength, - right.KernelRowLength, k, j); - NativeIntrinsic.Context - .mat[f64].kernel_fma((MatrixPanel*)leftPanel, - (MatrixPanel*)rightPanel, - (MatrixPanel*)ansPanel); - } - NativeIntrinsic.Context - .mat[f64].store_panel((MatrixPanel*)ansPanel, - (MatrixKernel*)ans.Buffer, - ans.KernelColumnLength, - ans.KernelRowLength, i, j); - } - } - MatrixPanel.Destroy(leftPanel); - MatrixPanel.Destroy(rightPanel); - MatrixPanel.Destroy(ansPanel); } if (count > 256) { MatrixParallel.ReturnParallelResources(ref pool!); diff --git a/StgSharp.Native/StgSharp.Native.vcxproj b/StgSharp.Native/StgSharp.Native.vcxproj index a267628a..69f7554e 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj +++ b/StgSharp.Native/StgSharp.Native.vcxproj @@ -160,6 +160,7 @@ + diff --git a/StgSharp.Native/StgSharp.Native.vcxproj.filters b/StgSharp.Native/StgSharp.Native.vcxproj.filters index f5e7ef25..fc090e7e 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj.filters +++ b/StgSharp.Native/StgSharp.Native.vcxproj.filters @@ -63,6 +63,9 @@ 源文件 + + 源文件 + @@ -80,9 +83,6 @@ 头文件 - - 头文件 - 头文件 @@ -104,6 +104,9 @@ 头文件 + + 头文件 + diff --git a/StgSharp.Native/sn-head/sn_internal.h b/StgSharp.Native/sn-head/sn_internal.h index c2d751f0..12e80b68 100644 --- a/StgSharp.Native/sn-head/sn_internal.h +++ b/StgSharp.Native/sn-head/sn_internal.h @@ -8,16 +8,19 @@ #define _CRT_SECURE_NO_WARNINGS 1 #define WIN32_LEAN_AND_MEAN #define SN_DECL __cdecl +#define SN_MAT_DECL __cdecl __forceinline #define SN_API __declspec(dllexport) #define FORCEINLINE __forceinline #include #elif defined(__clang__) || defined(__GNUC__) #define SN_DECL __attribute__((cdecl)) +#define SN_MAT_DECL __attribute__((cdecl)) __attribute__((always_inline)) #define SN_API __attribute__((dllexport)) #define FORCEINLINE __attribute__((always_inline)) #else #define SN_API #define SN_DECL +#define SN_MAT_DECL #define FORCEINLINE #endif diff --git a/StgSharp.Native/sn-head/sn_intrincontext.x86.h b/StgSharp.Native/sn-head/sn_intrincontext.x86.h index b68e0e63..83d1bb99 100644 --- a/StgSharp.Native/sn-head/sn_intrincontext.x86.h +++ b/StgSharp.Native/sn-head/sn_intrincontext.x86.h @@ -159,6 +159,16 @@ DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, , fma); DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma); DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, 512, , fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, sse, , fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, , fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, _fma, fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma); + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, , fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma); +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, 512, , fma); + #pragma endregion #pragma region scaler_simd diff --git a/StgSharp.Native/sn-head/sn_intrinsic.h b/StgSharp.Native/sn-head/sn_intrinsic.h index 36a6b2f7..380a1992 100644 --- a/StgSharp.Native/sn-head/sn_intrinsic.h +++ b/StgSharp.Native/sn-head/sn_intrinsic.h @@ -188,6 +188,38 @@ typedef void(SN_DECL *KER_PROC_LAS)(mat_kernel *matrix, mat_kernel *ans, scalar_ #pragma endregion +#pragma region matrix kernel tile function type + +typedef void(SN_DECL *KERTILE_PROC_AS)(mat_kernel *ans, scalar_pack *scalar, __m128 *enumeration); +typedef void(SN_DECL *KERTILE_PROC_RA)(mat_kernel *left, mat_kernel *ans, __m128 *enumeration); +typedef void(SN_DECL *KERTILE_PROC_LRA)(mat_kernel *left, mat_kernel *right, mat_kernel *ans, + __m128 *enumeration); +typedef void(SN_DECL *KERTILE_PROC_LAS)(mat_kernel *matrix, mat_kernel *ans, scalar_pack *scalar, + __m128 *enumeration); + +#define KERTILE_PROC(T_type, T_arch, T_feature, T_op_name) \ + T_type##_kernel_tile_##T_op_name##_##T_arch##T_feature + +#define DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS_SCALAR(T_type, T_arch, T_feature, T_op_name) \ + INTERNAL void SN_DECL KERTILE_PROC(T_type, T_arch, T_feature, T_op_name)( \ + MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \ + MAT_KERNEL(T_type) *restrict ans, scalar_pack const *scalar, \ + __m128 const *enumeration) +#define DECLARE_KERTILE_PROC_RIGHT_ANS_SCALAR(T_type, T_arch, T_feature, T_op_name) \ + INTERNAL void SN_DECL KERTILE_PROC(T_type, T_arch, T_feature, T_op_name)( \ + MAT_KERNEL(T_type) const *right, MAT_KERNEL(T_type) *restrict ans, \ + scalar_pack const *scalar, __m128 const *enumeration) +#define DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(T_type, T_arch, T_feature, T_op_name) \ + INTERNAL void SN_DECL KERTILE_PROC(T_type, T_arch, T_feature, T_op_name)( \ + MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \ + MAT_KERNEL(T_type) *restrict ans, __m128 const *enumeration) +#define DECLARE_KERTILE_PROC_RIGHT_ANS(T_type, T_arch, T_feature, T_op_name) \ + INTERNAL void SN_DECL KERTILE_PROC(T_type, T_arch, T_feature, T_op_name)( \ + MAT_KERNEL(T_type) const *right, MAT_KERNEL(T_type) *restrict ans, \ + __m128 const *enumeration) + +#pragma endregion + #pragma region special matrix kernel operation typedef void(SN_DECL *PIVOTPROC)(mat_kernel const *source, size_t col_begin, size_t col_length, @@ -216,7 +248,8 @@ typedef struct mat_intrinsic { KER_PROC_RA buffer_transpose; PANEL_PROC build_panel; PANEL_PROC clear_panel; - PNL_PROC_LRA panel_fma; + PNL_PROC_LRA ker_fma; + KERTILE_PROC_LRA ker_tile_fma; PIVOTPROC pivot; PANEL_PROC store_panel; diff --git a/StgSharp.Native/sn-src/load_intrinsic.x86.c b/StgSharp.Native/sn-src/load_intrinsic.x86.c index b17ec566..a1d70eb0 100644 --- a/StgSharp.Native/sn-src/load_intrinsic.x86.c +++ b/StgSharp.Native/sn-src/load_intrinsic.x86.c @@ -52,7 +52,8 @@ void load_intrinsic_sse(sn_intrinsic *context, SIMDID id) 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 = KER_PROC(float, sse, , fma); + context->mat[F32].ker_fma = KER_PROC(float, sse, , fma); + context->mat[F32].ker_tile_fma = KERTILE_PROC(float, sse, , fma); context->mat[F32].pivot = pivot_float; context->mat[F32].store_panel = NULL; @@ -63,7 +64,8 @@ void load_intrinsic_sse(sn_intrinsic *context, SIMDID id) context->mat[F64].buffer_transpose = NULL; context->mat[F64].build_panel = NULL; context->mat[F64].clear_panel = PNL_PROC(double, sse, , clear); - context->mat[F64].panel_fma = KER_PROC(double, sse, , fma); + context->mat[F64].ker_fma = KER_PROC(double, sse, , fma); + context->mat[F64].ker_tile_fma = KERTILE_PROC(double, sse, , fma); context->mat[F64].pivot = NULL; context->mat[F64].store_panel = NULL; @@ -87,8 +89,10 @@ void load_intrinsic_avx(sn_intrinsic *context, SIMDID id) 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 = is_fma ? KER_PROC(float, avx, , fma) : - KER_PROC(float, avx, _fma, fma); + context->mat[F32].ker_fma = is_fma ? KER_PROC(float, avx, , fma) : + KER_PROC(float, avx, _fma, fma); + context->mat[F32].ker_tile_fma = is_fma ? KERTILE_PROC(float, avx, , fma) : + KERTILE_PROC(float, avx, _fma, fma); context->mat[F32].pivot = pivot_float; context->mat[F32].store_panel = STORE_PANEL(avx, float); @@ -99,8 +103,10 @@ void load_intrinsic_avx(sn_intrinsic *context, SIMDID id) 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 = is_fma ? KER_PROC(double, avx, , fma) : - KER_PROC(double, avx, _fma, fma); + context->mat[F64].ker_fma = is_fma ? KER_PROC(double, avx, , fma) : + KER_PROC(double, avx, _fma, fma); + context->mat[F64].ker_tile_fma = is_fma ? KERTILE_PROC(double, avx, , fma) : + KERTILE_PROC(double, avx, _fma, fma); context->mat[F64].pivot = NULL; context->mat[F64].store_panel = NULL; @@ -120,7 +126,8 @@ void load_intrinsic_512(sn_intrinsic *context, SIMDID id) 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 = KER_PROC(float, 512, , fma); + context->mat[F32].ker_fma = KER_PROC(float, 512, , fma); + context->mat[F32].ker_tile_fma = KERTILE_PROC(float, 512, , fma); context->mat[F32].pivot = pivot_float; context->mat[F32].store_panel = STORE_PANEL(512, float); @@ -131,7 +138,8 @@ void load_intrinsic_512(sn_intrinsic *context, SIMDID id) context->mat[F64].buffer_transpose = NULL; context->mat[F64].build_panel = BUILD_PANEL(512, double); context->mat[F64].clear_panel = PNL_PROC(double, 512, , clear); - context->mat[F64].panel_fma = KER_PROC(double, 512, , fma); + context->mat[F64].ker_fma = KER_PROC(double, 512, , fma); + context->mat[F64].ker_tile_fma = KERTILE_PROC(double, 512, , fma); context->mat[F64].pivot = NULL; context->mat[F64].store_panel = STORE_PANEL(512, double); diff --git a/StgSharp.Native/sn-src/m.fma.x86.c b/StgSharp.Native/sn-src/m.fma.x86.c index 4a931a8e..4f93dcc9 100644 --- a/StgSharp.Native/sn-src/m.fma.x86.c +++ b/StgSharp.Native/sn-src/m.fma.x86.c @@ -180,6 +180,50 @@ DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, 512, , fma) // double +DECLARE_KER_PROC_LEFT_RIGHT_ANS(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; + } +} + DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, , fma) { const __m256d a0 = _mm256_load_pd((double *)&left->f64_y[0]); @@ -201,27 +245,27 @@ DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, , fma) 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)))); ans->f64_y[0] = ccol; - /* Column 0 */ + /* Column 1 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; - /* Column 0 */ + ans->f64_y[1] = ccol; + /* Column 2 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; - /* Column 0 */ + ans->f64_y[2] = ccol; + /* Column 3 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; + ans->f64_y[3] = ccol; } DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma) @@ -245,27 +289,27 @@ DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma) 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); ans->f64_y[0] = ccol; - /* Column 0 */ + /* Column 1 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; - /* Column 0 */ + 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); + ans->f64_y[1] = ccol; + /* Column 2 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; - /* Column 0 */ + 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); + ans->f64_y[2] = ccol; + /* Column 3 */ ccol = _mm256_load_pd((double *)&ans->f32_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)))); - ans->f64_y[0] = ccol; + 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); + ans->f64_y[3] = ccol; } DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, 512, , fma) @@ -314,172 +358,4 @@ DECLARE_KER_PROC_LEFT_RIGHT_ANS(double, 512, , fma) ans->f64_z[1] = c23; } -DECLARE_KER_PROC_LEFT_RIGHT_ANS(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; - } -} - -DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, avx, , fma) -{ - register __m256d a0 = _mm256_load_pd((double *)&left->vec[0][0]); - register __m256d a1 = _mm256_load_pd((double *)&left->vec[1][0]); - register __m256d a2 = _mm256_load_pd((double *)&left->vec[2][0]); - register __m256d a3 = _mm256_load_pd((double *)&left->vec[3][0]); - - register __m256d c0 = _mm256_load_pd((double *)&ans->vec[0][0]); - register __m256d c1 = _mm256_load_pd((double *)&ans->vec[1][0]); - register __m256d c2 = _mm256_load_pd((double *)&ans->vec[2][0]); - register __m256d c3 = _mm256_load_pd((double *)&ans->vec[3][0]); - - register __m256d bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c0 = _mm256_add_pd(c0, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); - c0 = _mm256_add_pd(c0, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); - c0 = _mm256_add_pd(c0, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); - c0 = _mm256_add_pd(c0, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); - _mm256_store_pd((double *)&ans->vec[0][0], c0); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c1 = _mm256_add_pd(c1, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); - c1 = _mm256_add_pd(c1, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); - c1 = _mm256_add_pd(c1, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); - c1 = _mm256_add_pd(c1, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); - _mm256_store_pd((double *)&ans->vec[1][0], c1); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c2 = _mm256_add_pd(c2, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); - c2 = _mm256_add_pd(c2, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); - c2 = _mm256_add_pd(c2, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); - c2 = _mm256_add_pd(c2, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); - _mm256_store_pd((double *)&ans->vec[2][0], c2); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c3 = _mm256_add_pd(c3, _mm256_mul_pd(a0, _mm256_permute_pd(bcol, 0x0))); - c3 = _mm256_add_pd(c3, _mm256_mul_pd(a1, _mm256_permute_pd(bcol, 0x3))); - c3 = _mm256_add_pd(c3, _mm256_mul_pd(a2, _mm256_permute_pd(bcol, 0x0))); - c3 = _mm256_add_pd(c3, _mm256_mul_pd(a3, _mm256_permute_pd(bcol, 0x3))); - _mm256_store_pd((double *)&ans->vec[3][0], c3); -} - -DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma) -{ - register __m256d a0 = _mm256_load_pd((double *)&left->vec[0][0]); - register __m256d a1 = _mm256_load_pd((double *)&left->vec[1][0]); - register __m256d a2 = _mm256_load_pd((double *)&left->vec[2][0]); - register __m256d a3 = _mm256_load_pd((double *)&left->vec[3][0]); - - register __m256d c0 = _mm256_load_pd((double *)&ans->vec[0][0]); - register __m256d c1 = _mm256_load_pd((double *)&ans->vec[1][0]); - register __m256d c2 = _mm256_load_pd((double *)&ans->vec[2][0]); - register __m256d c3 = _mm256_load_pd((double *)&ans->vec[3][0]); - - register __m256d bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c0 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c0); - c0 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c0); - c0 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c0); - c0 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c0); - _mm256_store_pd((double *)&ans->vec[0][0], c0); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c1 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c1); - c1 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c1); - c1 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c1); - c1 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c1); - _mm256_store_pd((double *)&ans->vec[1][0], c1); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c2 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c2); - c2 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c2); - c2 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c2); - c2 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c2); - _mm256_store_pd((double *)&ans->vec[2][0], c2); - - bcol = _mm256_load_pd((double *)&right->vec[0][0]); - c3 = _mm256_fmadd_pd(a0, _mm256_permute_pd(bcol, 0x0), c3); - c3 = _mm256_fmadd_pd(a1, _mm256_permute_pd(bcol, 0x3), c3); - c3 = _mm256_fmadd_pd(a2, _mm256_permute_pd(bcol, 0x0), c3); - c3 = _mm256_fmadd_pd(a3, _mm256_permute_pd(bcol, 0x3), c3); - _mm256_store_pd((double *)&ans->vec[3][0], c3); -} - -// ans = left * right + ans -// right is not transposed -DECLARE_PNL_PROC_LEFT_RIGHT_ANS(double, 512, , fma) -{ - const __m512i ones = _mm512_set1_epi32(1); - - register __m512d a0 = _mm512_load_pd((double *)&left->vec[0][0]); - register __m512d a1 = _mm512_load_pd((double *)&left->vec[1][0]); - register __m512d a2 = _mm512_load_pd((double *)&left->vec[2][0]); - register __m512d a3 = _mm512_load_pd((double *)&left->vec[3][0]); - register __m512d a4 = _mm512_load_pd((double *)&left->vec[4][0]); - register __m512d a5 = _mm512_load_pd((double *)&left->vec[5][0]); - register __m512d a6 = _mm512_load_pd((double *)&left->vec[6][0]); - register __m512d a7 = _mm512_load_pd((double *)&left->vec[7][0]); - - for (int col = 0; col < 8; ++col) { - register __m512d bcol = _mm512_load_pd((double *)&right->vec[col][0]); - register __m512d ccol = _mm512_load_pd((double *)&ans->vec[col][0]); - register __m512i idx = _mm512_setzero_si512(); - register __m512i increament = _mm512_load_si512(&ones); - - ccol = _mm512_fmadd_pd(a0, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a1, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a2, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a3, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a4, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a5, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a6, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - ccol = _mm512_fmadd_pd(a7, _mm512_permutexvar_pd(idx, bcol), ccol); - idx = _mm512_add_epi64(idx, increament); - - _mm512_store_pd((double *)&ans->vec[col][0], ccol); - } -} - #endif diff --git a/StgSharp.Native/sn-src/m.fma_tile.x86.c b/StgSharp.Native/sn-src/m.fma_tile.x86.c new file mode 100644 index 00000000..ffdb2c0d --- /dev/null +++ b/StgSharp.Native/sn-src/m.fma_tile.x86.c @@ -0,0 +1,583 @@ +#include "sn_target.h" + +#if SN_IS_ARCH(SN_ARCH_X86_64) + +#include "sn_intrinsic.h" +#include "sn_intrincontext.x86.h" + +#define _FMA_SELECT(i) _MM_SHUFFLE(i, i, i, i) + +/* +[MethodImpl(MethodImplOptions.AggressiveInlining)] +internal static unsafe MatrixKernel* GetKernelAddressUnsafe( + MatrixKernel* buffer, + int colLength, + int x, + int y +) where T : unmanaged, INumber +{ + return buffer + (colLength * x + y); +} +*/ + +#define GET_INDEX(col_len, x, y) col_len *x + y + +/* +for (int k = 0; k < commonK; k += 1) { + NativeIntrinsic.Context + .mat[f64].kernel_fma( + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)leftBuffer, + ansHeight, k, j), + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)rightBuffer, + commonK, j, k), + (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel*)ansBuffer, + ansHeight, i, j)); +} +*/ + +// float + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, sse, , fma) +{ + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(float) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(float) + const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe + + 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); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(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 } }; + + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(float) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(float) + const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *right_cur = + right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe + + 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); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(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 } }; + + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(float) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(float) + const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *right_cur = + right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe + + 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); + + 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); + } + _mm256_store_ps((float *)&ans_cur->f32_y[0], c01); + _mm256_store_ps((float *)&ans_cur->f32_y[1], c23); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma) +{ + static const __m512i b_shift_source = { .m512i_i32 = { 0, 1, 2, 3, 5, 6, 7, 4, 10, 11, 8, 9, + 15, 12, 13, 14 } }; + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + register __m512i b_shift = _mm512_load_si512(&b_shift_source); + register __m512 c = _mm512_load_ps((float *)(ans + GET_INDEX(ansHeight, i, j))); + for (size_t k = 0; k < commonK; k++) { + register __m512 a = _mm512_load_ps((float *)left); + register __m512 b = _mm512_load_ps((float *)right); + + 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, c); +} + +// double +#define DOUBLE_KERNEL_FMA_CYCLE(i) \ + /* 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_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)); + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma) +{ + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + 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)); + + for (size_t k = 0; k < commonK; k++) { + MAT_KERNEL(double) + const *left_cur = left + GET_INDEX(ansHeight, k, j), + const *right_cur = right + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe + + 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); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, , fma) +{ + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(double) + const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe + + const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_z[0]); + const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_z[1]); + const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_z[2]); + const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_z[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); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma) +{ + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(double) + const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe + + const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_z[0]); + const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_z[1]); + const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_z[2]); + const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_z[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); +} + +DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(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 } }; + + int i = enumeration->m128_i32[0]; + int j = enumeration->m128_i32[1]; + int commonK = enumeration->m128_i32[2]; + int ansHeight = enumeration->m128_i32[3]; + + MAT_KERNEL(double) + const *ans_cur = ans + GET_INDEX(ansHeight, i, j); // 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 k = 0; k < commonK; k++) { + MAT_KERNEL(double) + const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe + const *right_cur = + right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe + + 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); + + register __m512d c01 = _mm512_load_pd((double *)&ans->f64_z[0]); + register __m512d c23 = _mm512_load_pd((double *)&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_permute2f128_pd(a01, a23, 0x21); + register __m512d a30 = _mm512_permute2f128_pd(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); + /* + 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); + + ans->f64_z[0] = c01; + ans->f64_z[1] = c23; + /**/ +} + +#endif diff --git a/StgSharpDebug/Program.cs b/StgSharpDebug/Program.cs index 8878d99a..22b9bc50 100644 --- a/StgSharpDebug/Program.cs +++ b/StgSharpDebug/Program.cs @@ -56,14 +56,14 @@ namespace StgSharp unchecked { hlsf = new((nuint)(64 * (ulong)10000 * 10000)); } - Matrix m = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, - hlsf); + Matrix m = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, + hlsf); MatrixCompute.Fill(m, 1.0f); - Matrix n = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, - hlsf); + Matrix n = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, + hlsf); MatrixCompute.Fill(n, 1.0f); - Matrix p = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, - hlsf); + Matrix p = Matrix.Create(4000, 4000, MatrixLayout.DenseRectangle, + hlsf); MatrixCompute.Fill(p, 0.0f); /**/