diff --git a/StgSharp.Common/Mathematics/Linear/Matrix.cs b/StgSharp.Common/Mathematics/Linear/Matrix.cs index e30103ba..0271b9b8 100644 --- a/StgSharp.Common/Mathematics/Linear/Matrix.cs +++ b/StgSharp.Common/Mathematics/Linear/Matrix.cs @@ -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 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 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>(); } + 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, diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.cs index d2ba217f..d378d0df 100644 --- a/StgSharp.Mathematics/Numeric/MatrixCompute.cs +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.cs @@ -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(2) = commonK; + enumeration->Member(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*)leftBuffer, - ansHeight, k, j), - GetKernelAddressUnsafe((MatrixKernel*)rightBuffer, - commonK, j, k), - GetKernelAddressUnsafe((MatrixKernel*)ansBuffer, - ansHeight, i, j)); - } + enumeration->Member(0) = i; + enumeration->Member(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(2) = commonK; + enumeration->Member(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(0) = i; enumeration->Member(1) = j; - enumeration->Member(2) = commonK; - enumeration->Member(3) = ansHeight; NativeIntrinsic.Context.mat[f64].kernel_tile_fma(leftBuffer, rightBuffer, ansBuffer, enumeration); } diff --git a/StgSharp.Mathematics/Numeric/MatrixParallelWrap.cs b/StgSharp.Mathematics/Numeric/MatrixParallelWrap.cs index 8e1a2ddd..59d4efa0 100644 --- a/StgSharp.Mathematics/Numeric/MatrixParallelWrap.cs +++ b/StgSharp.Mathematics/Numeric/MatrixParallelWrap.cs @@ -163,7 +163,7 @@ namespace StgSharp.Mathematics.Numeric """; /**/ current = MatrixParallelFactory.FromExistPackage(current); - offset += length; + offset = length; } // Console.WriteLine(msg); diff --git a/StgSharp.Native/sn-src/m.fma_tile.x86.c b/StgSharp.Native/sn-src/m.fma_tile.x86.c index 7e877368..445cf9f6 100644 --- a/StgSharp.Native/sn-src/m.fma_tile.x86.c +++ b/StgSharp.Native/sn-src/m.fma_tile.x86.c @@ -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)); diff --git a/StgSharpDebug/MatrixTest.cs b/StgSharpDebug/MatrixTest.cs new file mode 100644 index 00000000..a33c3bba --- /dev/null +++ b/StgSharpDebug/MatrixTest.cs @@ -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 m = Matrix.Create(sideLength, sideLength, MatrixLayout.DenseRectangle, + hlsf); + MatrixCompute.Fill(m, 1.0f); + Matrix n = Matrix.Create(sideLength, sideLength, MatrixLayout.DenseRectangle, + hlsf); + MatrixCompute.Fill(n, 1.0f); + Matrix p = Matrix.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("==================================================="); + + /**/ + } + + } +} \ No newline at end of file diff --git a/StgSharpDebug/Program.cs b/StgSharpDebug/Program.cs index 515fb0da..a7327df0 100644 --- a/StgSharpDebug/Program.cs +++ b/StgSharpDebug/Program.cs @@ -49,58 +49,8 @@ namespace StgSharp .UseModule() .UseModule() .InitializeModule(); + TestMatrix(); - // TestRender(); - /**/ - HybridLayerSegregatedFitAllocator hlsf; - unchecked { - hlsf = new((nuint)(4 * (ulong)1024 * 1024 * 1024)); - } - Matrix m = Matrix.Create(40, 40, MatrixLayout.DenseRectangle, hlsf); - MatrixCompute.Fill(m, 1.0f); - Matrix n = Matrix.Create(40, 40, MatrixLayout.DenseRectangle, hlsf); - MatrixCompute.Fill(n, 1.0f); - Matrix p = Matrix.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 --- }