youmuliu/stgsharp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

tile fma added

Browse Source
This commit is contained in:
Youmu Liu
2026-01-09 18:13:46 +08:00
parent 4ff1d856a7
commit 6827259891
12 changed files with 736 additions and 265 deletions
Download Patch File Download Diff File Expand all files Collapse all files
StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs
+1
View File
@@ -46,6 +46,7 @@ namespace StgSharp.Internal.Intrinsic
public delegate* unmanaged[Cdecl]<MatrixPanel*, MatrixKernel*, int, int, int, int, void> build_panel; //f32_build_panel
public delegate* unmanaged[Cdecl]<MatrixPanel*, void> clear_panel; //f32_clear_panel
public delegate* unmanaged[Cdecl]<void*, void*, void*, void> kernel_fma; //f32_panel_fma
public delegate* unmanaged[Cdecl]<void*, void*, void*, M128*, void> kernel_tile_fma; //f32_panel_fma
public delegate* unmanaged[Cdecl]<void*, ulong, ulong, ulong, void> pivot; //f32_panel_mul
public delegate* unmanaged[Cdecl]<void*, void*, int, int, int, int, void> store_panel; //f32_store_panel
StgSharp.Mathematics/Numeric/MatrixCompute.cs
+2 -4
View File
@@ -95,9 +95,6 @@ namespace StgSharp.Mathematics.Numeric
{
int i = (int)(c / ansWidth);
int j = (int)(c % ansWidth);
MatrixKernel<double>* ansPanel = GetKernelAddressUnsafe((MatrixKernel<double>*)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<double>*)rightBuffer,
commonK, j, k),
(MatrixPanel*)ansPanel);
(MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel<double>*)ansBuffer,
ansHeight, i, j));
}
}
}
StgSharp.Mathematics/Numeric/MatrixCompute.double.cs
+6 -51
View File
@@ -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<double>* 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<double>* leftPanel = MatrixPanel.Create<double>();
MatrixPanel<double>* rightPanel = MatrixPanel.Create<double>();
MatrixPanel<double>* ansPanel = MatrixPanel.Create<double>();
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!);
StgSharp.Native/StgSharp.Native.vcxproj
+1
View File
@@ -160,6 +160,7 @@
<ClCompile Include="sn-src\dllmain.c" />
<ClCompile Include="sn-src\m.bufop.x86.c" />
<ClCompile Include="sn-src\m.pivot.x86.c" />
<ClCompile Include="sn-src\m.fma_tile.x86.c" />
<ClCompile Include="sn-src\s.factorial.x86.c" />
<ClCompile Include="sn-src\m.fma.x86.c" />
<ClCompile Include="sn-src\s.hash.x86.c" />
StgSharp.Native/StgSharp.Native.vcxproj.filters
+6 -3
View File
@@ -63,6 +63,9 @@
<ClCompile Include="sn-src\cpuid.x86.c">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="sn-src\m.fma_tile.x86.c">
<Filter>源文件</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="sn-head\gl.h">
@@ -80,9 +83,6 @@
<ClInclude Include="sn-head\wgl.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="sn-head\sn_intrinsic.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="sn-head\glfw_function.h">
<Filter>头文件</Filter>
</ClInclude>
@@ -104,6 +104,9 @@
<ClInclude Include="sn-head\sn_intrincontext.x86.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="sn-head\sn_intrinsic.h">
<Filter>头文件</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<Library Include="lib\glfw\glfw3.lib">
StgSharp.Native/sn-head/sn_internal.h
+3
View File
@@ -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 <windows.h>
#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
StgSharp.Native/sn-head/sn_intrincontext.x86.h
+10
View File
@@ -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
StgSharp.Native/sn-head/sn_intrinsic.h
+34 -1
View File
@@ -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;
StgSharp.Native/sn-src/load_intrinsic.x86.c
+16 -8
View File
@@ -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);
StgSharp.Native/sn-src/m.fma.x86.c
+68 -192
View File
@@ -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
StgSharp.Native/sn-src/m.fma_tile.x86.c
+583
View File
@@ -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<T>* GetKernelAddressUnsafe<T>(
MatrixKernel<T>* buffer,
int colLength,
int x,
int y
) where T : unmanaged, INumber<T>
{
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<double>*)leftBuffer,
ansHeight, k, j),
(MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel<double>*)rightBuffer,
commonK, j, k),
(MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel<double>*)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
StgSharpDebug/Program.cs
+6 -6
View File
@@ -56,14 +56,14 @@ namespace StgSharp
unchecked {
hlsf = new((nuint)(64 * (ulong)10000 * 10000));
}
Matrix<float> m = Matrix<float>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
Matrix<double> m = Matrix<double>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(m, 1.0f);
Matrix<float> n = Matrix<float>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
Matrix<double> n = Matrix<double>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(n, 1.0f);
Matrix<float> p = Matrix<float>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
Matrix<double> p = Matrix<double>.Create(4000, 4000, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(p, 0.0f);
/**/