seems tile_fma kernel has index error and indexer of matrix has error as well

This commit is contained in:
Youmu Liu
2026-01-11 16:45:38 +08:00
parent cd8210d9a1
commit b8bf1e9f08
8 changed files with 56 additions and 29 deletions
+2 -1
View File
@@ -1,4 +1,4 @@
//-----------------------------------------------------------------------
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="Global"
// Project: StgSharp
@@ -33,6 +33,7 @@ using System.Diagnostics.CodeAnalysis;
using System.Runtime.CompilerServices;
[assembly: InternalsVisibleTo("StgSharp.Mathematics")]
[assembly: InternalsVisibleTo("StgSharpDebug")]
[assembly: SuppressMessage(
"Design",
+9 -9
View File
@@ -158,12 +158,12 @@ namespace StgSharp.Mathematics.Numeric
)
{
// Caller guarantees: indices in range, layout is dense square.
int k_c = y >> 2;
int k_r = x >> 2;
int c = y & 3;
int r = x & 3;
int k_c = y / 4;
int k_r = x / 4;
int c = y % 4;
int r = x % 4;
ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in _size)];
ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_r, k_c, in _size)];
return ref kernel[c, r];
}
@@ -183,10 +183,10 @@ namespace StgSharp.Mathematics.Numeric
)
{
// Caller guarantees: indices in range, layout is dense square.
int k_c = y >> 2;
int k_r = x >> 2;
int c = y & 3;
int r = x & 3;
int k_c = y / 4;
int k_r = x / 4;
int c = y % 4;
int r = x % 4;
return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in _size)];
}
@@ -121,7 +121,7 @@ namespace StgSharp.Mathematics.Numeric
long y = yKernel;
long col = enumerator.KernelColumnLength;
Console.Write($"Attempt to get offset at ({x},{y}), offset is {x * col + y} ");
// Console.Write($"Attempt to get offset at ({x},{y}), offset is {x * col + y} ");
return x * col + y;
}
@@ -155,7 +155,7 @@ namespace StgSharp.Mathematics.Numeric
if (count > 256) {
pool = MatrixParallel.LeadParallel(1);
}
M128* enumeration = stackalloc M128[4];
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++)
@@ -162,9 +162,11 @@ namespace StgSharp.Mathematics.Numeric
(MatrixKernel*)right.Buffer,
(MatrixKernel*)ans.Buffer,
enumeration);
Console.WriteLine($"{i},{j},{ans[4*i,4*j]:0.0E00}\t");
// Console.WriteLine($"{i},{j},{ans[4*i,4*j]:0.0E00}\t");
}
Console.WriteLine();
// Console.WriteLine();
}
if (count > 256) {
MatrixParallel.ReturnParallelResources(ref pool!);
@@ -45,7 +45,11 @@ namespace StgSharp.Mathematics.Numeric
private int _remainThreadCount;
private readonly MatrixParallelHandle _handle;
internal MatrixParallelWrap(MatrixParallelHandle handle, int begin, int count)
internal MatrixParallelWrap(
MatrixParallelHandle handle,
int begin,
int count
)
{
WrapTaskCapacity = count;
_handle = handle;
@@ -121,7 +125,9 @@ namespace StgSharp.Mathematics.Numeric
_ = Interlocked.Exchange(ref _remainThreadCount, WrapTaskCapacity);
}
private void PublishBufferTaskToWorkers(MatrixParallelTaskPackage* package)
private void PublishBufferTaskToWorkers(
MatrixParallelTaskPackage* package
)
{
ReadOnlySpan<MatrixParallelThread> pool = Threads;
if (pool.Length == 1)
@@ -157,7 +163,7 @@ namespace StgSharp.Mathematics.Numeric
""";
/**/
current = MatrixParallelFactory.FromExistPackage(current);
offset = length;
offset += length;
}
// Console.WriteLine(msg);
+19 -11
View File
@@ -55,8 +55,9 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, sse, , fma)
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, i, k); // GetKernelAddressUnsafe
const *right_cur = right + GET_INDEX(commonK, i, k), // GetKernelAddressUnsafe
const *left_cur =
left + GET_INDEX(ansHeight, k, j); // 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]);
@@ -182,9 +183,9 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, avx, _fma, fma)
for (size_t k = 0; k < commonK; k++) {
MAT_KERNEL(float)
const *left_cur = left + GET_INDEX(commonK, i, k), // GetKernelAddressUnsafe
const *left_cur = left + GET_INDEX(ansHeight, k, j), // GetKernelAddressUnsafe
const *right_cur =
right + GET_INDEX(ansHeight, k, j); // GetKernelAddressUnsafe
right + GET_INDEX(commonK, i, 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]);
@@ -239,8 +240,13 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma)
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);
MAT_KERNEL(double) *const left_cur = left + GET_INDEX(ansHeight, k, j),
const *right_cur =
right + GET_INDEX(commonK, i,
k); // GetKernelAddressUnsafe
register __m512 a = _mm512_load_ps((float *)left_cur);
register __m512 b = _mm512_load_ps((float *)right_cur);
b = _mm512_permutexvar_ps(b_shift, b);
register __m512 a_perm = _mm512_shuffle_f32x4(a, a, _MM_SHUFFLE(3, 2, 1, 0));
@@ -260,7 +266,7 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(float, 512, , fma)
c = _mm512_fmadd_ps(a_perm, b_perm, c);
}
_mm512_store_ps((float *)ans, c);
_mm512_store_ps((float *)(ans + GET_INDEX(ansHeight, i, j)), c);
}
// double
@@ -358,8 +364,9 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, , fma)
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, i, k); // GetKernelAddressUnsafe
const *right_cur = right + GET_INDEX(commonK, i, k), // GetKernelAddressUnsafe
const *left_cur =
left + GET_INDEX(ansHeight, k, j); // GetKernelAddressUnsafe
const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]);
const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]);
@@ -435,8 +442,9 @@ DECLARE_KERTILE_PROC_LEFT_RIGHT_ANS(double, avx, _fma, fma)
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, i, k); // GetKernelAddressUnsafe
const *right_cur = right + GET_INDEX(commonK, i, k), // GetKernelAddressUnsafe
const *left_cur =
left + GET_INDEX(ansHeight, k, j); // GetKernelAddressUnsafe
const __m256d a0 = _mm256_load_pd((double *)&left_cur->f64_y[0]);
const __m256d a1 = _mm256_load_pd((double *)&left_cur->f64_y[1]);
+11 -1
View File
@@ -83,13 +83,23 @@ namespace StgSharp
}
for (int j = 0; j < 40; j++)
{
string extra = j % 4 == 3 ? " \t" : ",\t";
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 ---