From 688a420574ddcdb0b74c17dc1576123acfd92a6f Mon Sep 17 00:00:00 2001 From: Youmu Liu <105042599+Nitload-NSI@users.noreply.github.com> Date: Wed, 17 Dec 2025 16:56:13 +0800 Subject: [PATCH] Working on Factorize --- StgSharp.Common/Mathematics/Linear/Matrix.cs | 14 +++- .../Mathematics/Linear/Matrix_T.HLSF.cs | 21 +++-- .../Numeric/MatrixCompute.Factorize.cs | 56 +++++++++++++ StgSharp.Native/StgSharp.Native.vcxproj | 1 + .../StgSharp.Native.vcxproj.filters | 3 + StgSharp.Native/sn-head/sn_intrinsic.h | 15 ++++ StgSharp.Native/sn-src/k_pivot.c | 79 +++++++++++++++++++ StgSharpDebug/Program.cs | 6 +- 8 files changed, 181 insertions(+), 14 deletions(-) create mode 100644 StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs create mode 100644 StgSharp.Native/sn-src/k_pivot.c diff --git a/StgSharp.Common/Mathematics/Linear/Matrix.cs b/StgSharp.Common/Mathematics/Linear/Matrix.cs index c75d2bf3..5738ecca 100644 --- a/StgSharp.Common/Mathematics/Linear/Matrix.cs +++ b/StgSharp.Common/Mathematics/Linear/Matrix.cs @@ -1,4 +1,4 @@ -//----------------------------------------------------------------------- +//----------------------------------------------------------------------- // ----------------------------------------------------------------------- // file="Matrix" // Project: StgSharp @@ -52,7 +52,8 @@ namespace StgSharp.Mathematics.Numeric public static MatrixAllocation DefaultAllocation { get; set; } = MatrixAllocation.GC; - public static Matrix FromDefault(int columnLength, int rowLength) where T : unmanaged, INumber + public static Matrix FromDefault(int columnLength, int rowLength) + where T : unmanaged, INumber { return DefaultAllocation switch { @@ -62,15 +63,20 @@ namespace StgSharp.Mathematics.Numeric }; } - public static Matrix FromGC(int columnLength, int rowLength) where T : unmanaged, INumber + public static Matrix FromGC(int columnLength, int rowLength) + where T : unmanaged, INumber { return new GCMatrix(columnLength, rowLength); } [MethodImpl(MethodImplOptions.AggressiveInlining)] - public static Matrix FromHlsf(int columnLength, int rowLength, hlsfAllocator allocator) + public static Matrix FromHlsf( + int columnLength, + int rowLength, + hlsfAllocator allocator) where T : unmanaged, INumber { + ArgumentNullException.ThrowIfNull(allocator, nameof(allocator)); return new HlsfMatrix(columnLength, rowLength, allocator); } diff --git a/StgSharp.Common/Mathematics/Linear/Matrix_T.HLSF.cs b/StgSharp.Common/Mathematics/Linear/Matrix_T.HLSF.cs index e1f1f530..c0b28d22 100644 --- a/StgSharp.Common/Mathematics/Linear/Matrix_T.HLSF.cs +++ b/StgSharp.Common/Mathematics/Linear/Matrix_T.HLSF.cs @@ -1,4 +1,4 @@ -//----------------------------------------------------------------------- +//----------------------------------------------------------------------- // ----------------------------------------------------------------------- // file="Matrix_T.HLSF" // Project: StgSharp @@ -52,7 +52,11 @@ namespace StgSharp.Mathematics.Numeric KernelRowLength = (row + 3) / t_count; ColumnLength = column; RowLength = row; - this.handle = allocator.Alloc((uint)(KernelColumnLength * KernelRowLength * sizeof(MatrixKernel))); + long size = (long)KernelColumnLength * (long)KernelRowLength * (long)sizeof(MatrixKernel); + this.handle = allocator.Alloc((uint)(size)); + if (handle.AllocSize < size) { + throw new InsufficientMemoryException(); + } k_buffer = (MatrixKernel*)handle.Pointer; } @@ -68,7 +72,9 @@ namespace StgSharp.Mathematics.Numeric } } - public override Matrix this[ReadOnlySpan column, ReadOnlySpan row] + public override Matrix this[ + ReadOnlySpan column, + ReadOnlySpan row] { get => throw new NotImplementedException(); set => throw new NotImplementedException(); @@ -106,14 +112,15 @@ namespace StgSharp.Mathematics.Numeric public override unsafe MatrixSegmentEnumeration ColumnEnumeration { [MethodImpl(MethodImplOptions.AggressiveInlining)] - get => new MatrixSegmentEnumeration(k_buffer, KernelColumnLength, KernelColumnLength * KernelRowLength, - 1, KernelColumnLength); + get => new MatrixSegmentEnumeration(k_buffer, KernelColumnLength, + KernelColumnLength * KernelRowLength, 1, + KernelColumnLength); } public override unsafe MatrixSegmentEnumeration GetRowEnumeration() { - return new MatrixSegmentEnumeration(k_buffer, 1, KernelColumnLength, KernelColumnLength, - KernelRowLength); + return new MatrixSegmentEnumeration(k_buffer, 1, KernelColumnLength, + KernelColumnLength, KernelRowLength); } public override unsafe MatrixSegmentEnumeration SequentialEnumeration diff --git a/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs b/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs new file mode 100644 index 00000000..ebc58bdd --- /dev/null +++ b/StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs @@ -0,0 +1,56 @@ +//----------------------------------------------------------------------- +// ----------------------------------------------------------------------- +// file="MatrixCompute.Factorize" +// 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 System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace StgSharp.Mathematics.Numeric +{ + public static partial class MatrixCompute + { + + public static void FactorizePLU( + Matrix matrix, + Matrix p, + Matrix lower, + Matrix upper) + { + if (matrix.ColumnLength != matrix.RowLength) { + throw new InvalidOperationException("Matrix must be square."); + } + int size = matrix.KernelColumnLength * matrix.KernelRowLength; + if (size < 256) + { + for (int i = 0; i < matrix.KernelColumnLength; i++) { } + } else if (size < 1024) { } else { } + } + + } +} diff --git a/StgSharp.Native/StgSharp.Native.vcxproj b/StgSharp.Native/StgSharp.Native.vcxproj index 81129ad1..2858b025 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj +++ b/StgSharp.Native/StgSharp.Native.vcxproj @@ -152,6 +152,7 @@ + diff --git a/StgSharp.Native/StgSharp.Native.vcxproj.filters b/StgSharp.Native/StgSharp.Native.vcxproj.filters index 92bd11f3..8161ef1c 100644 --- a/StgSharp.Native/StgSharp.Native.vcxproj.filters +++ b/StgSharp.Native/StgSharp.Native.vcxproj.filters @@ -54,6 +54,9 @@ 源文件 + + 源文件 + diff --git a/StgSharp.Native/sn-head/sn_intrinsic.h b/StgSharp.Native/sn-head/sn_intrinsic.h index 3a34a6f1..8ba80fc4 100644 --- a/StgSharp.Native/sn-head/sn_intrinsic.h +++ b/StgSharp.Native/sn-head/sn_intrinsic.h @@ -7,6 +7,7 @@ #include "sn_matkernel.h" #include #include +#include #ifdef _MSC_VER #include @@ -193,6 +194,20 @@ typedef void(SN_DECL *KER_PROC_LAS)(mat_kernel *matrix, mat_kernel *ans, scalar_ #pragma endregion +#pragma region special matrix kernel operation + +typedef void(SN_DECL *POVITPROC)(mat_kernel const *source, mat_kernel *ans); +#define MAT_POVIT_PROC(T_type, T_arch) povit_##T_type##_##T_arch +#define DECLARE_MAT_POVIT_PROC(T_type, T_arch) \ + INTERNAL void SN_DECL MAT_POVIT_PROC(T_type, T_arch)(MAT_KERNEL(T_type) const *source, \ + size_t col_begin, size_t col_length, \ + size_t *restrict ans) +DECLARE_MAT_POVIT_PROC(float, sse); +DECLARE_MAT_POVIT_PROC(float, avx); +DECLARE_MAT_POVIT_PROC(float, 512); + +#pragma endregion + typedef void(SN_DECL *VECTORNORMALIZEPROC)(__m128 *source, __m128 *target); typedef void(SN_DECL *DOTPROC)(void *transpose, __m128 *vector, __m128 *ans); typedef uint64_t(SN_DECL *FACTORIALROC)(int n); diff --git a/StgSharp.Native/sn-src/k_pivot.c b/StgSharp.Native/sn-src/k_pivot.c new file mode 100644 index 00000000..6617c14b --- /dev/null +++ b/StgSharp.Native/sn-src/k_pivot.c @@ -0,0 +1,79 @@ +#include "sn_intrinsic.h" + +// Hybrid implementation: SIMD for Abs, Scalar for Max +// Assumes Column-Major layout for both MAT_KERNEL content and Kernel sequence +static void povit_float_simd_abs(MAT_KERNEL(float) const *source, size_t col_begin, + size_t col_length, size_t *restrict ans) +{ + float max_v[4] = { -1.0f, -1.0f, -1.0f, -1.0f }; + + // Initialize ans with the first row index of the search range + size_t start_row = col_begin * 4; + ans[0] = start_row; + ans[1] = start_row; + ans[2] = start_row; + ans[3] = start_row; + + const MAT_KERNEL(float) *k_ptr = source + col_begin; + + // Mask for absolute value: 0x7FFFFFFF (clears the sign bit) + __m128i abs_mask = _mm_set1_epi32(0x7FFFFFFF); + + for (size_t i = 0; i < col_length; ++i) { + const float *mat_ptr = (const float *)(k_ptr + i); + size_t base_row = (col_begin + i) * 4; + + // Iterate over 4 columns of the kernel + for (int c = 0; c < 4; ++c) { + // Load column c: 4 floats (Column-Major: contiguous in memory) + __m128 col_vec = _mm_load_ps(mat_ptr + c * 4); + + // Convert to integer vector, AND with mask, convert back + __m128i col_i = _mm_castps_si128(col_vec); + col_i = _mm_and_si128(col_i, abs_mask); + col_vec = _mm_castsi128_ps(col_i); + + VEC(float) vals; + _mm_store_ps(vals.v, col_vec); + + // Local reduction to find the max within this column chunk (4 elements) + // This breaks the dependency chain and reduces branches/cmovs on the global max_v + float v0 = vals.v[0]; + float v1 = vals.v[1]; + float v2 = vals.v[2]; + float v3 = vals.v[3]; + + // Compare pairs (0 vs 1) and (2 vs 3) - these can execute in parallel + int i01 = (v1 > v0) ? 1 : 0; + float m01 = (v1 > v0) ? v1 : v0; + + int i23 = (v3 > v2) ? 3 : 2; + float m23 = (v3 > v2) ? v3 : v2; + + // Compare winners + int i_local = (m23 > m01) ? i23 : i01; + float m_local = (m23 > m01) ? m23 : m01; + + // Single check against global maximum + if (m_local > max_v[c]) { + max_v[c] = m_local; + ans[c] = base_row + i_local; + } + } + } +} + +DECLARE_MAT_POVIT_PROC(float, sse) +{ + povit_float_simd_abs(source, col_begin, col_length, ans); +} + +DECLARE_MAT_POVIT_PROC(float, avx) +{ + povit_float_simd_abs(source, col_begin, col_length, ans); +} + +DECLARE_MAT_POVIT_PROC(float, 512) +{ + povit_float_simd_abs(source, col_begin, col_length, ans); +} diff --git a/StgSharpDebug/Program.cs b/StgSharpDebug/Program.cs index 11299973..bbae7ad6 100644 --- a/StgSharpDebug/Program.cs +++ b/StgSharpDebug/Program.cs @@ -58,11 +58,11 @@ namespace StgSharp unchecked { hlsf = new((nuint)(64 * (ulong)10000 * 10000)); } - Matrix m = Matrix.FromHlsf(400, 400, hlsf); + Matrix m = Matrix.FromHlsf(4000, 4000, hlsf); MatrixCompute.Fill(m, 1.0f); - Matrix n = Matrix.FromHlsf(400, 400, hlsf); + Matrix n = Matrix.FromHlsf(4000, 4000, hlsf); MatrixCompute.Fill(n, 1.0f); - Matrix p = Matrix.FromHlsf(400, 400, hlsf); + Matrix p = Matrix.FromHlsf(4000, 4000, hlsf); // MatrixCompute.Fill(p, 0.0f); /*