Merge branch 'giga' of ssh://ssh.github.com:443/Nitload-NSI/CCD_Github into local

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<T> FromDefault<T>(int columnLength, int rowLength) where T : unmanaged, INumber<T>
+        public static Matrix<T> FromDefault<T>(int columnLength, int rowLength)
+            where T : unmanaged, INumber<T>
         {
             return DefaultAllocation switch
             {
@@ -62,15 +63,20 @@ namespace StgSharp.Mathematics.Numeric
             };
         }
 
-        public static Matrix<T> FromGC<T>(int columnLength, int rowLength) where T : unmanaged, INumber<T>
+        public static Matrix<T> FromGC<T>(int columnLength, int rowLength)
+            where T : unmanaged, INumber<T>
         {
             return new GCMatrix<T>(columnLength, rowLength);
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static Matrix<T> FromHlsf<T>(int columnLength, int rowLength, hlsfAllocator allocator)
+        public static Matrix<T> FromHlsf<T>(
+                                int columnLength,
+                                int rowLength,
+                                hlsfAllocator allocator)
             where T : unmanaged, INumber<T>
         {
+            ArgumentNullException.ThrowIfNull(allocator, nameof(allocator));
             return new HlsfMatrix<T>(columnLength, rowLength, allocator);
         }
 

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<T>)));
+            long size = (long)KernelColumnLength * (long)KernelRowLength * (long)sizeof(MatrixKernel<T>);
+            this.handle = allocator.Alloc((uint)(size));
+            if (handle.AllocSize < size) {
+                throw new InsufficientMemoryException();
+            }
             k_buffer = (MatrixKernel<T>*)handle.Pointer;
         }
 
@@ -68,7 +72,9 @@ namespace StgSharp.Mathematics.Numeric
             }
         }
 
-        public override Matrix<T> this[ReadOnlySpan<MatrixIndex> column, ReadOnlySpan<MatrixIndex> row]
+        public override Matrix<T> this[
+                                  ReadOnlySpan<MatrixIndex> column,
+                                  ReadOnlySpan<MatrixIndex> row]
         {
             get => throw new NotImplementedException();
             set => throw new NotImplementedException();
@@ -106,14 +112,15 @@ namespace StgSharp.Mathematics.Numeric
         public override unsafe MatrixSegmentEnumeration<T> ColumnEnumeration
         {
             [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => new MatrixSegmentEnumeration<T>(k_buffer, KernelColumnLength, KernelColumnLength * KernelRowLength,
+            get => new MatrixSegmentEnumeration<T>(k_buffer, KernelColumnLength,
-                                                   1, KernelColumnLength);
+                                                   KernelColumnLength * KernelRowLength, 1,
+                                                   KernelColumnLength);
         }
 
         public override unsafe MatrixSegmentEnumeration<T> GetRowEnumeration()
         {
-            return new MatrixSegmentEnumeration<T>(k_buffer, 1, KernelColumnLength, KernelColumnLength,
+            return new MatrixSegmentEnumeration<T>(k_buffer, 1, KernelColumnLength,
-                                                   KernelRowLength);
+                                                   KernelColumnLength, KernelRowLength);
         }
 
         public override unsafe MatrixSegmentEnumeration<T> SequentialEnumeration

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<float> matrix,
+                           Matrix<float> p,
+                           Matrix<float> lower,
+                           Matrix<float> 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 { }
+        }
+
+    }
+}

StgSharp.Native/StgSharp.Native.vcxproj

@@ -152,6 +152,7 @@
     <ClCompile Include="sn-src\k_basicop.c" />
     <ClCompile Include="sn-src\dllmain.c" />
     <ClCompile Include="sn-src\k_bufop.c" />
+    <ClCompile Include="sn-src\k_pivot.c" />
     <ClCompile Include="sn-src\str_factorial.c" />
     <ClCompile Include="sn-src\k_fma.c" />
     <ClCompile Include="sn-src\str_hash.c" />

StgSharp.Native/StgSharp.Native.vcxproj.filters

@@ -54,6 +54,9 @@
     <ClCompile Include="sn-src\k_panel.c">
       <Filter>源文件</Filter>
     </ClCompile>
+    <ClCompile Include="sn-src\k_pivot.c">
+      <Filter>源文件</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="sn-head\gl.h">

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -7,6 +7,7 @@
 #include "sn_matkernel.h"
 #include <immintrin.h>
 #include <smmintrin.h>
+#include <stdint.h>
 
 #ifdef _MSC_VER
 #include <intrin.h>
@@ -203,6 +204,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);

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);
+}

StgSharpDebug/Program.cs

@@ -58,11 +58,11 @@ namespace StgSharp
                 unchecked {
                     hlsf = new((nuint)(64 * (ulong)10000 * 10000));
                 }
-                Matrix<float> m = Matrix.FromHlsf<float>(400, 400, hlsf);
+                Matrix<float> m = Matrix.FromHlsf<float>(4000, 4000, hlsf);
                 MatrixCompute.Fill(m, 1.0f);
-                Matrix<float> n = Matrix.FromHlsf<float>(400, 400, hlsf);
+                Matrix<float> n = Matrix.FromHlsf<float>(4000, 4000, hlsf);
                 MatrixCompute.Fill(n, 1.0f);
-                Matrix<float> p = Matrix.FromHlsf<float>(400, 400, hlsf);
+                Matrix<float> p = Matrix.FromHlsf<float>(4000, 4000, hlsf);
 
                 // MatrixCompute.Fill(p, 0.0f);
                 /*