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

tiled fma is ok, mat mul is now available again

Browse Source
This commit is contained in:
Youmu Liu
2026-01-12 15:09:39 +08:00
parent b8bf1e9f08
commit ad3b2212fa
6 changed files with 172 additions and 119 deletions
Download Patch File Download Diff File Expand all files Collapse all files
StgSharp.Common/Mathematics/Linear/Matrix.cs
+27 -17
View File
@@ -67,30 +67,22 @@ namespace StgSharp.Mathematics.Numeric
if (y < 0 || y > _size.ElementColumnLength || x < 0 || x > _size.ElementRowLength) {
throw new ArgumentOutOfRangeException($"Index [{x},{y}] is out of matrix index range.");
}
int
k_c = y / 4,
c = y % 4,
k_r = x / 4,
r = x % 4;
ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_r, k_c);
return Unsafe.IsNullRef(ref kernel) ? T.Zero : kernel[r, c];
long offset = GetElementOffset(x, y);
if (offset < 0) {
return T.Zero;
}
return ((T*)Buffer)[offset];
}
set
{
if (y < 0 || y > _size.ElementColumnLength || x < 0 || x > _size.ElementRowLength) {
throw new ArgumentOutOfRangeException($"Index [{x},{y}] is out of matrix index range.");
}
int
k_c = y / 4,
c = y % 4,
k_r = x / 4,
r = x % 4;
ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_r, k_c);
if (!Unsafe.IsNullRef(ref kernel)) {
kernel[r, c] = value;
long offset = GetElementOffset(x, y);
if (offset < 0) {
return;
}
((T*)Buffer)[offset] = value;
}
}
@@ -151,6 +143,24 @@ namespace StgSharp.Mathematics.Numeric
ref Unsafe.NullRef<MatrixKernel<T>>();
}
internal long GetElementOffset(
int x,
int y
)
{
int
k_c = y / 4,
c = y % 4,
k_r = x / 4,
r = x % 4;
if (!_boarder.IsInBuffer(k_r, k_c, in _size)) {
return -1;
}
long kernelOffset = _boarder.GetIndexOffset(k_r, k_c, in _size);
return kernelOffset * 16 + 4 * r + c;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal unsafe ref T GetElementUnsafe(
int x,
StgSharp.Mathematics/Numeric/MatrixCompute.cs
+9 -13
View File
@@ -54,23 +54,19 @@ namespace StgSharp.Mathematics.Numeric
int ansHeight = p->AnsSecOffset;
int commonK = p->LeftPrimOffset;
// int blocksPerRow = ansWidth / size;
M128* enumeration = stackalloc M128[1];
enumeration->Member<int>(2) = commonK;
enumeration->Member<int>(3) = ansHeight;
// in this occasion size must be 1
for (long c = offsetRef; c < offsetRef + lengthRef; c++)
{
int i = (int)(c / ansHeight);
int j = (int)(c % ansHeight);
for (int k = 0; k < commonK; k++) {
NativeIntrinsic.Context
.mat[f32].kernel_fma(
GetKernelAddressUnsafe((MatrixKernel<float>*)leftBuffer,
ansHeight, k, j),
GetKernelAddressUnsafe((MatrixKernel<float>*)rightBuffer,
commonK, j, k),
GetKernelAddressUnsafe((MatrixKernel<float>*)ansBuffer,
ansHeight, i, j));
}
enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j;
NativeIntrinsic.Context.mat[f32].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer,
enumeration);
}
}
@@ -91,6 +87,8 @@ namespace StgSharp.Mathematics.Numeric
int commonK = p->LeftPrimOffset;
M128* enumeration = stackalloc M128[1];
enumeration->Member<int>(2) = commonK;
enumeration->Member<int>(3) = ansHeight;
// in this occasion size must be 1
for (long c = offsetRef; c < offsetRef + lengthRef; c++)
@@ -99,8 +97,6 @@ namespace StgSharp.Mathematics.Numeric
int j = (int)(c % ansHeight);
enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j;
enumeration->Member<int>(2) = commonK;
enumeration->Member<int>(3) = ansHeight;
NativeIntrinsic.Context.mat[f64].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer,
enumeration);
}
StgSharp.Mathematics/Numeric/MatrixParallelWrap.cs
+1 -1
View File
@@ -163,7 +163,7 @@ namespace StgSharp.Mathematics.Numeric
""";
/**/
current = MatrixParallelFactory.FromExistPackage(current);
offset += length;
offset = length;
}
// Console.WriteLine(msg);
StgSharp.Native/sn-src/m.fma_tile.x86.c
+31 -37
View File
@@ -175,8 +175,8 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, _fma, fma)
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 __m256 c01 = _mm256_load_ps((float *)&ans_cur->f32_y[0]);
register __m256 c23 = _mm256_load_ps((float *)&ans_cur->f32_y[1]);
register __m256i b01_shift = _mm256_load_si256(&b01_shift_source);
register __m256i b23_shift = _mm256_load_si256(&b23_shift_source);
@@ -195,9 +195,6 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, _fma, fma)
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));
@@ -270,35 +267,35 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma)
}
// 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)); \
#define DOUBLE_KERNEL_FMA_CYCLE(i) \
/* k = 0 */ \
register __m128d b_lo = _mm_load_pd(right_cur + (2 * i + 0)); \
register __m128d b_hi = _mm_load_pd(right_cur + (2 * i + 1)); \
register __m128d a_lo = _mm_load_pd(left_cur + (0 * 2 + 0)); \
register __m128d a_hi = _mm_load_pd(left_cur + (0 * 2 + 1)); \
register __m128d b = _mm_shuffle_pd(b_lo, b_lo, 0x0); \
c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b)); \
c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b)); \
\
/* group = 1 */ \
a_lo = _mm_load_pd(left_cur + (1 * 2 + 0)); \
a_hi = _mm_load_pd(left_cur + (1 * 2 + 1)); \
b = _mm_shuffle_pd(b_lo, b_lo, 0x3); \
c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b)); \
c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b)); \
\
/* groupp = 2 */ \
a_lo = _mm_load_pd(left_cur + (2 * 2 + 0)); \
a_hi = _mm_load_pd(left_cur + (2 * 2 + 1)); \
b = _mm_shuffle_pd(b_hi, b_hi, 0x0); \
c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b)); \
c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b)); \
\
/* froup = 3 */ \
a_lo = _mm_load_pd(left_cur + (3 * 2 + 0)); \
a_hi = _mm_load_pd(left_cur + (3 * 2 + 1)); \
b = _mm_shuffle_pd(b_hi, b_hi, 0x3); \
c##i##_lo = _mm_add_pd(c##i##_lo, _mm_mul_pd(a_lo, b)); \
c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));
DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma)
@@ -519,9 +516,6 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, 512, , fma)
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));
StgSharpDebug/MatrixTest.cs
+103
View File
@@ -0,0 +1,103 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="MatrixTest"
// Project: StgSharp
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using Microsoft.CodeAnalysis.CSharp.Syntax;
using StgSharp.HighPerformance.Memory;
using StgSharp.Mathematics.Numeric;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
using T = double;
namespace StgSharp
{
internal static partial class Program
{
private static int sideLength = 400;
internal static unsafe void TestMatrix(
)
{
// TestRender();
/**/
HybridLayerSegregatedFitAllocator hlsf;
unchecked {
hlsf = new((nuint)(4 * (ulong)1024 * 1024 * 1024));
}
Matrix<T> m = Matrix<T>.Create(sideLength, sideLength, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(m, 1.0f);
Matrix<T> n = Matrix<T>.Create(sideLength, sideLength, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(n, 1.0f);
Matrix<T> p = Matrix<T>.Create(sideLength, sideLength, MatrixLayout.DenseRectangle,
hlsf);
MatrixCompute.Fill(p, 0.0f);
/*
for (int j = 0; j < 20; j++)
{
for (int k = 0; k < 10; k++) {
Console.Write($"{m[j, k]:0.0E+00}\t");
}
Console.WriteLine();
}
/**/
Console.WriteLine("===================================================");
MatrixCompute.Mul(n, m, p);
/**/
for (int i = 0; i < 400; i++)
{
if (i % 4 == 0) {
Console.WriteLine();
}
for (int j = 0; j < 30; j++)
{
string extra = j % 4 == 3 ? " " : " ";
T value = p[j, i];
// value = T.Abs(value) < 1e-6 ? 0 : value;
// Console.Write($"{p.GetElementOffset(j,i)}{extra}");
Console.Write($"{value:000}{extra}");
}
Console.WriteLine();
}
Console.WriteLine("===================================================");
/**/
}
}
}
StgSharpDebug/Program.cs
+1 -51
View File
@@ -49,58 +49,8 @@ namespace StgSharp
.UseModule<NumericalModule>()
.UseModule<GraphicModule>()
.InitializeModule();
TestMatrix();
// TestRender();
/**/
HybridLayerSegregatedFitAllocator hlsf;
unchecked {
hlsf = new((nuint)(4 * (ulong)1024 * 1024 * 1024));
}
Matrix<double> m = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf);
MatrixCompute.Fill(m, 1.0f);
Matrix<double> n = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf);
MatrixCompute.Fill(n, 1.0f);
Matrix<double> p = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf);
MatrixCompute.Fill(p, 0.0f);
/*
for (int j = 0; j < 20; j++)
{
for (int k = 0; k < 10; k++) {
Console.Write($"{m[j, k]:0.0E+00}\t");
}
Console.WriteLine();
}
/**/
Console.WriteLine("===================================================");
MatrixCompute.Mul(n, m, p);
/**/
for (int i = 0; i < 40; i++)
{
if (i % 4 == 0) {
Console.WriteLine();
}
for (int j = 0; j < 40; j++)
{
string extra = j % 4 == 3 ? " " : ",";
double value = p[i, j];
value = double.Abs(value) < 1e-6 ? 0 : value;
Console.Write($"{value:0E+0}{extra}");
}
Console.WriteLine();
}
Console.WriteLine("===================================================");
double *buffer = (double*)p.Buffer;
for (int i = 0; i < 40; i++)
{
for (int j = 0; j < 40; j++) {
Console.Write($"{buffer[i * 40 + j]:0.0E0} ");
}
Console.WriteLine();
}
/**/
// --- Logarithmic time sequence test: print seconds at each emitted span ---
}