fixing tile fma

This commit is contained in:
Youmu Liu
2026-01-11 15:28:30 +08:00
parent 4be5c1c7c9
commit cd8210d9a1
8 changed files with 198 additions and 149 deletions
@@ -53,17 +53,22 @@ namespace StgSharp.HighPerformance
Unsafe.SkipInit(out this); Unsafe.SkipInit(out this);
} }
public MatrixKernel(bool clear) public MatrixKernel(
bool clear
)
{ {
if (!clear) { if (!clear) {
Unsafe.SkipInit(out this); Unsafe.SkipInit(out this);
} }
} }
public ref T this[int row, int column] public ref T this[
int x,
int y
]
{ {
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
get => ref Unsafe.Add(ref c0r0, (row * 4) + column); get => ref Unsafe.Add(ref c0r0, (x * 4) + y);
} }
} }
+11 -21
View File
@@ -70,11 +70,11 @@ namespace StgSharp.Mathematics.Numeric
int int
k_c = y / 4, k_c = y / 4,
k_r = x / 4,
c = y % 4, c = y % 4,
k_r = x / 4,
r = x % 4; r = x % 4;
ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_c, k_r); ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_r, k_c);
return Unsafe.IsNullRef(ref kernel) ? T.Zero : kernel[c, r]; return Unsafe.IsNullRef(ref kernel) ? T.Zero : kernel[r, c];
} }
set set
{ {
@@ -84,12 +84,12 @@ namespace StgSharp.Mathematics.Numeric
int int
k_c = y / 4, k_c = y / 4,
k_r = x / 4,
c = y % 4, c = y % 4,
k_r = x / 4,
r = x % 4; r = x % 4;
ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_c, k_r); ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_r, k_c);
if (!Unsafe.IsNullRef(ref kernel)) { if (!Unsafe.IsNullRef(ref kernel)) {
kernel[c, r] = value; kernel[r, c] = value;
} }
} }
} }
@@ -146,10 +146,8 @@ namespace StgSharp.Mathematics.Numeric
int y int y
) )
{ {
MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength, return ref _boarder.IsInBuffer(x, y, in _size) ?
x, y); ref _buffer[_boarder.GetIndexOffset(x, y, in _size)] :
return ref _boarder.IsInBuffer(x, y, in enumerator) ?
ref _buffer[_boarder.GetIndexOffset(x, y, in enumerator)] :
ref Unsafe.NullRef<MatrixKernel<T>>(); ref Unsafe.NullRef<MatrixKernel<T>>();
} }
@@ -165,9 +163,7 @@ namespace StgSharp.Mathematics.Numeric
int c = y & 3; int c = y & 3;
int r = x & 3; int r = x & 3;
MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength, ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in _size)];
k_c, k_r);
ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in enumerator)];
return ref kernel[c, r]; return ref kernel[c, r];
} }
@@ -177,10 +173,7 @@ namespace StgSharp.Mathematics.Numeric
int y int y
) )
{ {
// Caller guarantees: indices in range, layout is dense square. return Buffer + _boarder.GetIndexOffset(x, y, in _size);
MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength,
x, y);
return Buffer + _boarder.GetIndexOffset(x, y, in enumerator);
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
@@ -194,10 +187,7 @@ namespace StgSharp.Mathematics.Numeric
int k_r = x >> 2; int k_r = x >> 2;
int c = y & 3; int c = y & 3;
int r = x & 3; int r = x & 3;
return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in _size)];
MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength,
k_c, k_r);
return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in enumerator)];
} }
} }
@@ -48,7 +48,9 @@ namespace StgSharp.Mathematics.Numeric
private static readonly MatrixKernelBoarder UpperTriangleInstance = new UpperTriangleKernelBoarder(); private static readonly MatrixKernelBoarder UpperTriangleInstance = new UpperTriangleKernelBoarder();
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public static MatrixKernelBoarder Create(MatrixLayout layout) public static MatrixKernelBoarder Create(
MatrixLayout layout
)
{ {
return layout switch return layout switch
{ {
@@ -63,25 +65,30 @@ namespace StgSharp.Mathematics.Numeric
public abstract long GetIndexOffset( public abstract long GetIndexOffset(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator); in MatrixSize size
);
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public abstract bool IsInBuffer( public abstract bool IsInBuffer(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator); in MatrixSize enumerator
);
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
protected static long GetTriangularOrder(in MatrixKernelEnumerator enumerator) protected static long GetTriangularOrder(
in MatrixSize enumerator
)
{ {
return Math.Min(enumerator.KernelColumnLength, enumerator.KernelRowLength); return Math.Min(enumerator.KernelColumnLength, enumerator.KernelRowLength);
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
protected static bool IsWithinKernelDomain( protected static bool IsWithinKernelDomain(
in MatrixKernelEnumerator enumerator, in MatrixSize enumerator,
int xKernel, int xKernel,
int yKernel) int yKernel
)
{ {
return xKernel >= 0 && return xKernel >= 0 &&
yKernel >= 0 && yKernel >= 0 &&
@@ -91,9 +98,10 @@ namespace StgSharp.Mathematics.Numeric
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
protected static bool IsWithinTriangularDomain( protected static bool IsWithinTriangularDomain(
in MatrixKernelEnumerator enumerator, in MatrixSize enumerator,
int xKernel, int xKernel,
int yKernel) int yKernel
)
{ {
long order = GetTriangularOrder(enumerator); long order = GetTriangularOrder(enumerator);
return xKernel >= 0 && yKernel >= 0 && xKernel < order && yKernel < order; return xKernel >= 0 && yKernel >= 0 && xKernel < order && yKernel < order;
@@ -106,19 +114,23 @@ namespace StgSharp.Mathematics.Numeric
public override long GetIndexOffset( public override long GetIndexOffset(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
long x = xKernel; long x = xKernel;
long y = yKernel; long y = yKernel;
long rows = enumerator.KernelRowLength; long col = enumerator.KernelColumnLength;
return x * rows + y;
Console.Write($"Attempt to get offset at ({x},{y}), offset is {x * col + y} ");
return x * col + y;
} }
[MethodImpl(MethodImplOptions.AggressiveInlining)] [MethodImpl(MethodImplOptions.AggressiveInlining)]
public override bool IsInBuffer( public override bool IsInBuffer(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
return IsWithinKernelDomain(enumerator, xKernel, yKernel); return IsWithinKernelDomain(enumerator, xKernel, yKernel);
} }
@@ -132,7 +144,8 @@ namespace StgSharp.Mathematics.Numeric
public override long GetIndexOffset( public override long GetIndexOffset(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
long order = GetTriangularOrder(enumerator); long order = GetTriangularOrder(enumerator);
long x = xKernel; long x = xKernel;
@@ -146,7 +159,8 @@ namespace StgSharp.Mathematics.Numeric
public override bool IsInBuffer( public override bool IsInBuffer(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
return IsWithinTriangularDomain(enumerator, xKernel, yKernel) && yKernel >= xKernel; return IsWithinTriangularDomain(enumerator, xKernel, yKernel) && yKernel >= xKernel;
} }
@@ -160,7 +174,8 @@ namespace StgSharp.Mathematics.Numeric
public override long GetIndexOffset( public override long GetIndexOffset(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
long x = xKernel; long x = xKernel;
long y = yKernel; long y = yKernel;
@@ -173,7 +188,8 @@ namespace StgSharp.Mathematics.Numeric
public override bool IsInBuffer( public override bool IsInBuffer(
int xKernel, int xKernel,
int yKernel, int yKernel,
in MatrixKernelEnumerator enumerator) in MatrixSize enumerator
)
{ {
return IsWithinTriangularDomain(enumerator, xKernel, yKernel) && yKernel <= xKernel; return IsWithinTriangularDomain(enumerator, xKernel, yKernel) && yKernel <= xKernel;
} }
@@ -197,7 +213,8 @@ namespace StgSharp.Mathematics.Numeric
int elementColumnLength, int elementColumnLength,
int elementRowLength, int elementRowLength,
int kernelColumnLength, int kernelColumnLength,
int kernelRowLength) int kernelRowLength
)
{ {
ElementColumnLength = elementColumnLength; ElementColumnLength = elementColumnLength;
ElementRowLength = elementRowLength; ElementRowLength = elementRowLength;
@@ -222,7 +239,8 @@ namespace StgSharp.Mathematics.Numeric
int kernelColumnLength, int kernelColumnLength,
int kernelRowLength, int kernelRowLength,
int currentColumn = 0, int currentColumn = 0,
int currentRow = 0) int currentRow = 0
)
{ {
CurrentColumn = currentColumn; CurrentColumn = currentColumn;
CurrentRow = currentRow; CurrentRow = currentRow;
@@ -1,9 +1,9 @@
//----------------------------------------------------------------------- //-----------------------------------------------------------------------
// ----------------------------------------------------------------------- // -----------------------------------------------------------------------
// file="MatrixCompute.Float" // file="MatrixCompute.Float"
// Project: StgSharp // Project: StgSharp
// AuthorGroup: Nitload Space // AuthorGroup: Nitload
// Copyright (c) Nitload Space. All rights reserved. // Copyright (c) Nitload. All rights reserved.
// //
// Permission is hereby granted, free of charge, to any person obtaining a copy // Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal // of this software and associated documentation files (the "Software"), to deal
@@ -42,7 +42,11 @@ namespace StgSharp.Mathematics.Numeric
public static partial class MatrixCompute public static partial class MatrixCompute
{ {
public static unsafe void Add(Matrix<float> left, Matrix<float> right, Matrix<float> ans) public static unsafe void Add(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
if (left.ColumnLength != right.ColumnLength || if (left.ColumnLength != right.ColumnLength ||
@@ -88,7 +92,10 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Fill(Matrix<float> ans, float value) public static unsafe void Fill(
Matrix<float> ans,
float value
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
long count = (long)ans.KernelColumnLength * ans.KernelRowLength; long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
@@ -119,7 +126,11 @@ namespace StgSharp.Mathematics.Numeric
taskGroup.Dispose(); taskGroup.Dispose();
} }
public static unsafe void Mul(Matrix<float> left, Matrix<float> right, Matrix<float> ans) public static unsafe void Mul(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
@@ -145,31 +156,20 @@ namespace StgSharp.Mathematics.Numeric
pool = MatrixParallel.LeadParallel(1); pool = MatrixParallel.LeadParallel(1);
} }
M128* enumeration = stackalloc M128[4]; M128* enumeration = stackalloc M128[4];
enumeration->Member<int>(2) = left.KernelColumnLength;
enumeration->Member<int>(3) = ans.KernelColumnLength;
for (int i = 0; i < ans.KernelRowLength; i++) for (int i = 0; i < ans.KernelRowLength; i++)
{ {
for (int j = 0; j < ans.KernelColumnLength; j++) for (int j = 0; j < ans.KernelColumnLength; j++)
{ {
enumeration->Member<int>(0) = i; enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j; enumeration->Member<int>(1) = j;
enumeration->Member<int>(2) = ans.KernelRowLength;
enumeration->Member<int>(3) = left.KernelColumnLength;
NativeIntrinsic.Context NativeIntrinsic.Context
.mat[f32].kernel_tile_fma((MatrixKernel*)left.Buffer, .mat[f32].kernel_tile_fma((MatrixKernel*)left.Buffer,
(MatrixKernel*)right.Buffer, (MatrixKernel*)right.Buffer,
(MatrixKernel*)ans.Buffer, enumeration); (MatrixKernel*)ans.Buffer,
/* enumeration);
for (int k = 0; k < left.KernelColumnLength; k++) {
NativeIntrinsic.Context
.mat[f32].kernel_fma(GetKernelAddressUnsafe(left.Buffer,
left.KernelColumnLength, i,
k),
GetKernelAddressUnsafe(right.Buffer,
right.KernelColumnLength, k, j),
GetKernelAddressUnsafe(ans.Buffer,
ans.KernelColumnLength, i, j));
}
/**/
} }
} }
@@ -196,13 +196,18 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Mul(Matrix<float> left, float right, Matrix<float> ans) public static unsafe void Mul(
Matrix<float> left,
float right,
Matrix<float> ans
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) { if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree."); throw new IndexOutOfRangeException("Matrix dimensions must agree.");
} }
if (left.Layout != MatrixLayout.DenseRectangle || ans.Layout != MatrixLayout.DenseRectangle) { if (left.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported."); throw new NotSupportedException("Only dense rectangle matrices are supported.");
} }
long count = (long)left.KernelColumnLength * left.KernelRowLength; long count = (long)left.KernelColumnLength * left.KernelRowLength;
@@ -243,7 +248,11 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Sub(Matrix<float> left, Matrix<float> right, Matrix<float> ans) public static unsafe void Sub(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
if (left.ColumnLength != right.ColumnLength || if (left.ColumnLength != right.ColumnLength ||
+23 -24
View File
@@ -1,9 +1,9 @@
//----------------------------------------------------------------------- //-----------------------------------------------------------------------
// ----------------------------------------------------------------------- // -----------------------------------------------------------------------
// file="MatrixCompute" // file="MatrixCompute"
// Project: StgSharp // Project: StgSharp
// AuthorGroup: Nitload Space // AuthorGroup: Nitload
// Copyright (c) Nitload Space. All rights reserved. // Copyright (c) Nitload. All rights reserved.
// //
// Permission is hereby granted, free of charge, to any person obtaining a copy // Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal // of this software and associated documentation files (the "Software"), to deal
@@ -40,7 +40,9 @@ namespace StgSharp.Mathematics.Numeric
public partial class MatrixCompute public partial class MatrixCompute
{ {
internal static unsafe void ComputeFmaPanelF32(MatrixParallelTaskPackage* p) internal static unsafe void ComputeFmaPanelF32(
MatrixParallelTaskPackage* p
)
{ {
int f32 = (int)MatrixElementType.F32; int f32 = (int)MatrixElementType.F32;
MatrixKernel* leftBuffer = (MatrixKernel*)p->Left; MatrixKernel* leftBuffer = (MatrixKernel*)p->Left;
@@ -57,19 +59,24 @@ namespace StgSharp.Mathematics.Numeric
// in this occasion size must be 1 // in this occasion size must be 1
for (long c = offsetRef; c < offsetRef + lengthRef; c++) for (long c = offsetRef; c < offsetRef + lengthRef; c++)
{ {
int i = (int)(c / ansWidth); int i = (int)(c / ansHeight);
int j = (int)(c % ansWidth); int j = (int)(c % ansHeight);
for (int k = 0; k < commonK; k++) { for (int k = 0; k < commonK; k++) {
NativeIntrinsic.Context NativeIntrinsic.Context
.mat[f32].kernel_fma( .mat[f32].kernel_fma(
GetKernelAddressUnsafe((MatrixKernel<float>*)leftBuffer, ansHeight, k, j), GetKernelAddressUnsafe((MatrixKernel<float>*)leftBuffer,
GetKernelAddressUnsafe((MatrixKernel<float>*)rightBuffer, commonK, j, k), ansHeight, k, j),
GetKernelAddressUnsafe((MatrixKernel<float>*)ansBuffer, ansHeight, i, j)); GetKernelAddressUnsafe((MatrixKernel<float>*)rightBuffer,
commonK, j, k),
GetKernelAddressUnsafe((MatrixKernel<float>*)ansBuffer,
ansHeight, i, j));
} }
} }
} }
internal static unsafe void ComputeFmaPanelF64(MatrixParallelTaskPackage* p) internal static unsafe void ComputeFmaPanelF64(
MatrixParallelTaskPackage* p
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
@@ -88,22 +95,14 @@ namespace StgSharp.Mathematics.Numeric
// in this occasion size must be 1 // in this occasion size must be 1
for (long c = offsetRef; c < offsetRef + lengthRef; c++) for (long c = offsetRef; c < offsetRef + lengthRef; c++)
{ {
int i = (int)(c / ansWidth); int i = (int)(c / ansHeight);
int j = (int)(c % ansWidth); int j = (int)(c % ansHeight);
enumeration->Member<int>(0) = i; enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j; enumeration->Member<int>(1) = j;
enumeration->Member<int>(2) = commonK; enumeration->Member<int>(2) = commonK;
enumeration->Member<int>(3) = ansHeight; enumeration->Member<int>(3) = ansHeight;
NativeIntrinsic.Context.mat[f64].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer, enumeration); NativeIntrinsic.Context.mat[f64].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer,
/* enumeration);
for (int k = 0; k < commonK; k += 1) {
NativeIntrinsic.Context
.mat[f64].kernel_fma(
GetKernelAddressUnsafe((MatrixKernel<double>*)leftBuffer, ansHeight, k, j),
GetKernelAddressUnsafe((MatrixKernel<double>*)rightBuffer, commonK, j, k),
GetKernelAddressUnsafe((MatrixKernel<double>*)ansBuffer, ansHeight, i, j));
}
/**/
} }
} }
@@ -112,8 +111,8 @@ namespace StgSharp.Mathematics.Numeric
MatrixKernel<T>* buffer, MatrixKernel<T>* buffer,
int colLength, int colLength,
int x, int x,
int y) int y
where T: unmanaged, INumber<T> ) where T : unmanaged, INumber<T>
{ {
return (MatrixKernel*)(buffer + (colLength * x + y)); return (MatrixKernel*)(buffer + (colLength * x + y));
} }
@@ -1,9 +1,9 @@
//----------------------------------------------------------------------- //-----------------------------------------------------------------------
// ----------------------------------------------------------------------- // -----------------------------------------------------------------------
// file="MatrixCompute.double" // file="MatrixCompute.double"
// Project: StgSharp // Project: StgSharp
// AuthorGroup: Nitload Space // AuthorGroup: Nitload
// Copyright (c) Nitload Space. All rights reserved. // Copyright (c) Nitload. All rights reserved.
// //
// Permission is hereby granted, free of charge, to any person obtaining a copy // Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal // of this software and associated documentation files (the "Software"), to deal
@@ -38,7 +38,11 @@ namespace StgSharp.Mathematics.Numeric
public partial class MatrixCompute public partial class MatrixCompute
{ {
public static unsafe void Add(Matrix<double> left, Matrix<double> right, Matrix<double> ans) public static unsafe void Add(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
if (left.ColumnLength != right.ColumnLength || if (left.ColumnLength != right.ColumnLength ||
@@ -84,7 +88,10 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Fill(Matrix<double> ans, double value) public static unsafe void Fill(
Matrix<double> ans,
double value
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
long count = (long)ans.KernelColumnLength * ans.KernelRowLength; long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
@@ -115,7 +122,11 @@ namespace StgSharp.Mathematics.Numeric
taskGroup.Dispose(); taskGroup.Dispose();
} }
public static unsafe void Mul(Matrix<double> left, Matrix<double> right, Matrix<double> ans) public static unsafe void Mul(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
@@ -137,34 +148,23 @@ namespace StgSharp.Mathematics.Numeric
pool = MatrixParallel.LeadParallel(1); pool = MatrixParallel.LeadParallel(1);
} }
M128* enumeration = stackalloc M128[1]; M128* enumeration = stackalloc M128[1];
enumeration->Member<int>(2) = left.KernelColumnLength;
enumeration->Member<int>(3) = ans.KernelColumnLength;
for (int i = 0; i < ans.KernelRowLength; i++) for (int i = 0; i < ans.KernelRowLength; i++)
{ {
for (int j = 0; j < ans.KernelColumnLength; j += 1) for (int j = 0; j < ans.KernelColumnLength; j++)
{ {
enumeration->Member<int>(0) = i; enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j; enumeration->Member<int>(1) = j;
enumeration->Member<int>(2) = ans.KernelRowLength;
enumeration->Member<int>(3) = left.KernelColumnLength;
NativeIntrinsic.Context NativeIntrinsic.Context
.mat[f64].kernel_tile_fma((MatrixKernel*)left.Buffer, .mat[f64].kernel_tile_fma((MatrixKernel*)left.Buffer,
(MatrixKernel*)right.Buffer, (MatrixKernel*)right.Buffer,
(MatrixKernel*)ans.Buffer, enumeration); (MatrixKernel*)ans.Buffer,
/* enumeration);
// Console.WriteLine($"{i}{j}"); Console.WriteLine($"{i},{j},{ans[4*i,4*j]:0.0E00}\t");
for (int k = 0; k < left.KernelColumnLength; k++)
{
NativeIntrinsic.Context
.mat[f64].kernel_fma(GetKernelAddressUnsafe(left.Buffer,
left.KernelColumnLength, i,
k),
GetKernelAddressUnsafe(right.Buffer,
right.KernelColumnLength, k, j),
GetKernelAddressUnsafe(ans.Buffer,
ans.KernelColumnLength, i, j));
}
/**/
} }
Console.WriteLine();
} }
if (count > 256) { if (count > 256) {
MatrixParallel.ReturnParallelResources(ref pool!); MatrixParallel.ReturnParallelResources(ref pool!);
@@ -189,13 +189,18 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Mul(Matrix<double> left, double right, Matrix<double> ans) public static unsafe void Mul(
Matrix<double> left,
double right,
Matrix<double> ans
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) { if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree."); throw new IndexOutOfRangeException("Matrix dimensions must agree.");
} }
if (left.Layout != MatrixLayout.DenseRectangle || ans.Layout != MatrixLayout.DenseRectangle) { if (left.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported."); throw new NotSupportedException("Only dense rectangle matrices are supported.");
} }
long count = (long)left.KernelColumnLength * left.KernelRowLength; long count = (long)left.KernelColumnLength * left.KernelRowLength;
@@ -236,7 +241,11 @@ namespace StgSharp.Mathematics.Numeric
return; return;
} }
public static unsafe void Sub(Matrix<double> left, Matrix<double> right, Matrix<double> ans) public static unsafe void Sub(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{ {
int f64 = (int)MatrixElementType.F64; int f64 = (int)MatrixElementType.F64;
if (left.ColumnLength != right.ColumnLength || if (left.ColumnLength != right.ColumnLength ||
+28 -22
View File
@@ -1,4 +1,4 @@
#include "sn_target.h" #include "sn_target.h"
#if SN_IS_ARCH(SN_ARCH_X86_64) #if SN_IS_ARCH(SN_ARCH_X86_64)
@@ -20,7 +20,7 @@ internal static unsafe MatrixKernel<T>* GetKernelAddressUnsafe<T>(
} }
*/ */
#define GET_INDEX(col_len, x, y) col_len *x + y #define GET_INDEX(col_len, x, y) (col_len * x + y)
/* /*
for (int k = 0; k < commonK; k += 1) { for (int k = 0; k < commonK; k += 1) {
@@ -56,7 +56,7 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, sse, , fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(float) MAT_KERNEL(float)
const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe
const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(commonK, i, k); // GetKernelAddressUnsafe
const __m128 a0 = _mm_load_ps((float *)&left_cur->f32_x[0]); 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 a1 = _mm_load_ps((float *)&left_cur->f32_x[1]);
@@ -119,7 +119,7 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, , fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(float) MAT_KERNEL(float)
const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(ansHeight, i, k), // GetKernelAddressUnsafe
const *right_cur = const *right_cur =
right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe
@@ -182,9 +182,9 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, _fma, fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(float) MAT_KERNEL(float)
const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(commonK, i, k), // GetKernelAddressUnsafe
const *right_cur = const *right_cur =
right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe right + GET_INDEX(ansHeight, k, j); // GetKernelAddressUnsafe
register __m256 a01 = _mm256_load_ps((float *)&left_cur->f32_y[0]); 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 a23 = _mm256_load_ps((float *)&left_cur->f32_y[1]);
@@ -293,7 +293,7 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma)
a_hi = _mm_load_pd(left_cur + (3 * 2 + 1)); \ a_hi = _mm_load_pd(left_cur + (3 * 2 + 1)); \
b = _mm_shuffle_pd(b_hi, b_hi, 0x3); \ 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##_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));\ c##i##_hi = _mm_add_pd(c##i##_hi, _mm_mul_pd(a_hi, b));
DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma) DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma)
{ {
@@ -313,12 +313,18 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, sse, , fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(double) MAT_KERNEL(double)
const *left_cur = left + GET_INDEX(ansHeight, k, j), const *left_cur = left + GET_INDEX(ansHeight, i, k),
const *right_cur = right + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe const *right_cur = right + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe
{DOUBLE_KERNEL_FMA_CYCLE(0)} {
{DOUBLE_KERNEL_FMA_CYCLE(1)} DOUBLE_KERNEL_FMA_CYCLE(0)
{DOUBLE_KERNEL_FMA_CYCLE(2)} }
{
DOUBLE_KERNEL_FMA_CYCLE(1)
}
{
DOUBLE_KERNEL_FMA_CYCLE(2)
}
{ {
DOUBLE_KERNEL_FMA_CYCLE(3) DOUBLE_KERNEL_FMA_CYCLE(3)
} }
@@ -353,12 +359,12 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, , fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(double) MAT_KERNEL(double)
const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe
const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(commonK, i, k); // GetKernelAddressUnsafe
const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_z[0]); const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]);
const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_z[1]); const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]);
const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_z[2]); const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]);
const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_z[3]); const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]);
/* Column 0 */ /* Column 0 */
bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]); bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]);
@@ -430,12 +436,12 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(double) MAT_KERNEL(double)
const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *right_cur = right + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe
const *left_cur = left + GET_INDEX(commonK, j, k); // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(commonK, i, k); // GetKernelAddressUnsafe
const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_z[0]); const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]);
const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_z[1]); const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]);
const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_z[2]); const __m256d a2 = _mm256_load_pd((double *)&left_cur->f64_y[2]);
const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_z[3]); const __m256d a3 = _mm256_load_pd((double *)&left_cur->f64_y[3]);
/* Column 0 */ /* Column 0 */
bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]); bcol = _mm256_load_pd((double *)&right_cur->f64_y[0]);
@@ -493,7 +499,7 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, 512, , fma)
for (size_t k = 0; k < commonK; k++) { for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(double) MAT_KERNEL(double)
const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe const *left_cur = left + GET_INDEX(ansHeight, i, k), // GetKernelAddressUnsafe
const *right_cur = const *right_cur =
right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe right + GET_INDEX(ansHeight, j, k); // GetKernelAddressUnsafe
+23 -10
View File
@@ -1,9 +1,9 @@
//----------------------------------------------------------------------- //-----------------------------------------------------------------------
// ----------------------------------------------------------------------- // -----------------------------------------------------------------------
// file="Program" // file="Program"
// Project: StgSharp // Project: StgSharp
// AuthorGroup: Nitload Space // AuthorGroup: Nitload
// Copyright (c) Nitload Space. All rights reserved. // Copyright (c) Nitload. All rights reserved.
// //
// Permission is hereby granted, free of charge, to any person obtaining a copy // Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal // of this software and associated documentation files (the "Software"), to deal
@@ -39,17 +39,22 @@ namespace StgSharp
private static string testExpPath = @"D:\Programmer\stepcode\data\ap203\ap203.exp"; private static string testExpPath = @"D:\Programmer\stepcode\data\ap203\ap203.exp";
internal static unsafe void Main(string[] args) internal static unsafe void Main(
string[] args
)
{ {
try try
{ {
World.Initialize().UseModule<NumericalModule>().UseModule<GraphicModule>().InitializeModule(); World.Initialize()
.UseModule<NumericalModule>()
.UseModule<GraphicModule>()
.InitializeModule();
// TestRender(); // TestRender();
/**/ /**/
HybridLayerSegregatedFitAllocator hlsf; HybridLayerSegregatedFitAllocator hlsf;
unchecked { unchecked {
hlsf = new((nuint)(64 * (ulong)10000 * 10000)); hlsf = new((nuint)(4 * (ulong)1024 * 1024 * 1024));
} }
Matrix<double> m = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf); Matrix<double> m = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf);
MatrixCompute.Fill(m, 1.0f); MatrixCompute.Fill(m, 1.0f);
@@ -58,7 +63,7 @@ namespace StgSharp
Matrix<double> p = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf); Matrix<double> p = Matrix<double>.Create(40, 40, MatrixLayout.DenseRectangle, hlsf);
MatrixCompute.Fill(p, 0.0f); MatrixCompute.Fill(p, 0.0f);
/**/ /*
for (int j = 0; j < 20; j++) for (int j = 0; j < 20; j++)
{ {
for (int k = 0; k < 10; k++) { for (int k = 0; k < 10; k++) {
@@ -70,10 +75,18 @@ namespace StgSharp
/**/ /**/
Console.WriteLine("==================================================="); Console.WriteLine("===================================================");
MatrixCompute.Mul(n, m, p); MatrixCompute.Mul(n, m, p);
for (int i = 0; i < 20; i++) /**/
for (int i = 0; i < 40; i++)
{ {
for (int j = 0; j < 10; j++) { if (i % 4 == 0) {
Console.Write($"{p[i, j]:0.0E+00}\t"); Console.WriteLine();
}
for (int j = 0; j < 40; j++)
{
string extra = j % 4 == 3 ? " \t" : ",\t";
double value = p[i, j];
value = double.Abs(value) < 1e-6 ? 0 : value;
Console.Write($"{value:0E+0}{extra}");
} }
Console.WriteLine(); Console.WriteLine();
} }