redesign matrix design

StgSharp.Common/HighPerformance/MatrixKernel.cs

@@ -44,6 +44,8 @@ namespace StgSharp.HighPerformance
     public unsafe struct MatrixKernel<T> where T : unmanaged, INumber<T>
     {
 
+        public static readonly int Size = 16 * sizeof(T);
+
         private T c0r0;
 
         public MatrixKernel()
@@ -58,10 +60,10 @@ namespace StgSharp.HighPerformance
             }
         }
 
-        public ref T this[int column, int row]
+        public ref T this[int row, int column]
         {
             [MethodImpl(MethodImplOptions.AggressiveInlining)]
-            get => ref Unsafe.Add(ref c0r0, (column * 4) + row);
+            get => ref Unsafe.Add(ref c0r0, (row * 4) + column);
         }
 
     }

StgSharp.Common/Mathematics/Linear/CoordinationNaming.md

@@ -0,0 +1,88 @@
+# Matrix Coordinate & Naming Conventions (StgSharp)
+
+This document defines the coordinate terminology and naming rules used across `Matrix<T>`, `MatrixKernel<T>`, and packed/triangular storage layouts.  
+Goal: eliminate ambiguity around `column`/`row` names (e.g., “index **of** a column” vs “index **along** a column”).
+
+---
+
+## 1) The key rule: encode prepositions into identifiers (`Of` vs `In`)
+
+### `...OfX` — “which X” (object identity / ordinal)
+Use **`Of`** when the value identifies **which object** (which column, which row, which kernel, etc.).
+
+Examples:
+- `columnOfElement`: the element is located in **which column**
+- `rowOfKernel`: the kernel block is located in **which row**
+
+### `...InX` / `...WithinX` — “inside X” (local coordinate)
+Use **`In`** (or `Within`) when the value is a **local coordinate inside** another object.
+
+Examples:
+- `columnInKernel`: local column coordinate **inside** a kernel block
+- `rowInKernel`: local row coordinate **inside** a kernel block
+
+> Quick mental model: **`Of` = identity**, **`In` = inside**.
+
+---
+
+## 2) Coordinate levels: `Element`, `KernelGrid`, `InKernel`
+
+StgSharp commonly uses three coordinate levels:
+
+1. **Element-level**: final matrix element coordinates (`Matrix<T>[col,row]`)
+2. **Kernel-grid level**: coordinates of a kernel block in the kernel grid
+3. **In-kernel level**: local element coordinates within a kernel (e.g., 4×4)
+
+Recommended names:
+
+| Level | Column name | Row name | Meaning |
+|---|---|---|---|
+| Element | `columnOfElement` | `rowOfElement` | element located in which column/row (0-based) |
+| Kernel grid | `columnOfKernel` | `rowOfKernel` | kernel block located in which kernel-column/kernel-row (0-based) |
+| In-kernel | `columnInKernel` | `rowInKernel` | local coordinate inside the kernel (0..KernelSize-1) |
+
+---
+
+## 3) Index base: 0-based (mandatory)
+
+Unless explicitly documented otherwise, all coordinates are **0-based**:
+
+- `columnOfElement ∈ [0, ColumnLength-1]`
+- `rowOfElement ∈ [0, RowLength-1]`
+- `columnOfKernel ∈ [0, KernelColumnLength-1]`
+- `rowOfKernel ∈ [0, KernelRowLength-1]`
+- `columnInKernel ∈ [0, KernelSize-1]`
+- `rowInKernel ∈ [0, KernelSize-1]`
+
+---
+
+## 4) Mapping element coordinates to kernel coordinates (KernelSize = 4 example)
+
+For a fixed kernel size `KernelSize = 4`:
+
+```csharp
+int columnOfKernel = columnOfElement / 4;
+int rowOfKernel    = rowOfElement    / 4;
+
+int columnInKernel = columnOfElement % 4;
+int rowInKernel    = rowOfElement    % 4;
+```
+
+If `columnOfElement` and `rowOfElement` are guaranteed non-negative, bit operations can be used:
+
+```csharp
+int columnOfKernel = columnOfElement >> 2;
+int rowOfKernel    = rowOfElement    >> 2;
+
+int columnInKernel = columnOfElement & 3;
+int rowInKernel    = rowOfElement    & 3;
+```
+
+---
+
+## 5) Packed storage: distinguish `Index` from `Offset`
+
+Packed layouts do **not** store “empty blocks”; stored blocks are placed contiguously in a packed buffer.
+
+### Use the correct terms
+- **`packedKernelIndex`**: “the N-th kernel

StgSharp.Common/Mathematics/Linear/Matrix.Create.cs

@@ -25,10 +25,12 @@
 //     
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
+using StgSharp.HighPerformance;
 using StgSharp.HighPerformance.Memory;
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Runtime.CompilerServices;
 using System.Text;
 using System.Threading.Tasks;
 
@@ -37,21 +39,72 @@ namespace StgSharp.Mathematics.Numeric
     public unsafe partial class Matrix<T>
     {
 
+        private const int KernelSideLength = 4;
+
         public static Matrix<T> Create(
                                 int columnLength,
                                 int rowLength,
+                                MatrixLayout layout,
                                 HybridLayerSegregatedFitAllocator allocator)
         {
-            Matrix<T> matrix = new Matrix<T>();
+            ArgumentNullException.ThrowIfNull(allocator);
-            int tSize = sizeof(T);
+
-            int totalSize = columnLength * rowLength * tSize;
+            ArgumentOutOfRangeException.ThrowIfNegativeOrZero(columnLength);
-            HybridLayerSegregatedFitAllocationHandle handle = allocator.Alloc((uint)totalSize);
+
-            MatrixStorageHLSF<T> storage = new MatrixStorageHLSF<T>(handle);
+            ArgumentOutOfRangeException.ThrowIfNegativeOrZero(rowLength);
-            storage._handle = handle;
+
+            int kernelColumns = DivideRoundUp(columnLength, KernelSideLength);
+            int kernelRows = DivideRoundUp(rowLength, KernelSideLength);
+
+            if (layout is MatrixLayout.UpperTriangle or MatrixLayout.LowerTriangle)
+            {
+                if (kernelColumns != kernelRows) {
+                    throw new ArgumentException("Triangular layouts require square matrices.");
+                }
+            }
+
+            long kernelCount = layout switch
+            {
+                MatrixLayout.DenseSquare => (long)kernelColumns * kernelRows,
+                MatrixLayout.UpperTriangle => TriangularKernelCount(kernelColumns),
+                MatrixLayout.LowerTriangle => TriangularKernelCount(kernelRows),
+                _ => throw new ArgumentOutOfRangeException(nameof(layout))
+            };
+
+            long requiredBytes = kernelCount * MatrixKernel<T>.Size;
+            if (requiredBytes <= 0 || requiredBytes > uint.MaxValue) {
+                throw new ArgumentOutOfRangeException(nameof(requiredBytes), "Requested matrix exceeds allocator capacity.");
+            }
+
+            HybridLayerSegregatedFitAllocationHandle handle = allocator.Alloc((uint)requiredBytes);
+            MatrixStorageHLSF<T> storage = new MatrixStorageHLSF<T>(allocator, handle);
+
+            Matrix<T> matrix = new Matrix<T>
+            {
+                Layout = layout,
+                KernelColumnLength = kernelColumns,
+                KernelRowLength = kernelRows,
+            };
+
             matrix._buffer = storage;
-            matrix._enumeration = new MatrixMetadata(columnLength, rowLength);
+            matrix._boarder = MatrixKernelBoarder.Create(layout);
+            matrix._size = new MatrixSize(columnLength, rowLength, kernelColumns, kernelRows);
+
             return matrix;
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int DivideRoundUp(int value, int divisor)
+        {
+            return (value + (divisor - 1)) / divisor;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static long TriangularKernelCount(int order)
+        {
+            long n = order;
+            return (n * (n + 1L)) / 2L;
+        }
+
     }
 }

StgSharp.Common/Mathematics/Linear/Matrix.cs

@@ -30,7 +30,6 @@ using StgSharp.Mathematics.Numeric;
 using System;
 using System.Numerics;
 using System.Runtime.CompilerServices;
-using System.Runtime.InteropServices;
 
 // using hlsfHandle = global::StgSharp.HighPerformance.Memory.HybridLayerSegregatedFitAllocationHandle;
 
@@ -42,17 +41,61 @@ namespace StgSharp.Mathematics.Numeric
         private protected static readonly int t_count = 16 / sizeof(T);
 
         private protected static readonly int t_size = sizeof(T);
-        private MatrixMetadata _enumeration;
+        private MatrixKernelBoarder _boarder;
-
+        private MatrixSize _size;
         private MatrixStorage<T> _buffer;
 
         private protected Matrix() { }
 
-        public ref T this[int column, int row]
+        public T this[int x, int y]
         {
-            get => throw new NotImplementedException();
+            get
+            {
+                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,
+                    k_r = x / 4,
+                    c = y % 4,
+                    r = x % 4;
+                ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_c, k_r);
+                return Unsafe.IsNullRef(ref kernel) ? T.Zero : kernel[c, r];
+            }
+            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,
+                    k_r = x / 4,
+                    c = y % 4,
+                    r = x % 4;
+                ref MatrixKernel<T> kernel = ref GetMatrixKernel(k_c, k_r);
+                if (!Unsafe.IsNullRef(ref kernel)) {
+                    kernel[c, r] = value;
+                }
+            }
         }
 
+        /// <summary>
+        ///   Matlab style indexer to generate a new matrix
+        /// </summary>
+        /// <param name="column">
+        ///
+        /// </param>
+        /// <param name="row">
+        ///
+        /// </param>
+        /// <returns>
+        ///
+        /// </returns>
+        /// <exception cref="NotImplementedException">
+        ///
+        /// </exception>
         public unsafe Matrix<T> this[
                                 ReadOnlySpan<MatrixIndex> column,
                                 ReadOnlySpan<MatrixIndex> row]
@@ -64,34 +107,33 @@ namespace StgSharp.Mathematics.Numeric
         /// <summary>
         ///   Length of the matrix in columns.
         /// </summary>
-        public int ColumnLength { get; init; }
+        public int ColumnLength => _size.ElementColumnLength;
 
         /// <summary>
         ///   Length of the matrix in rows.
         /// </summary>
-        public int RowLength { get; init; }
+        public int RowLength => _size.ElementRowLength;
 
         public MatrixLayout Layout { get; init; }
 
-        internal MatrixKernel<T>* Buffer { get; }
+        internal MatrixKernel<T>* Buffer
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => (MatrixKernel<T>*)_buffer.BufferPointer;
+        }
 
         protected internal int KernelColumnLength { get; init; }
 
         protected internal int KernelRowLength { get; init; }
 
-        public ref MatrixKernel<T> GetMatrixKernel(int column, int row)
+        public unsafe ref MatrixKernel<T> GetMatrixKernel(int x, int y)
         {
-            if (_enumeration.GetCurrentKernelIndex) { } else { }
+            fixed (MatrixSize *msPtr = &_size)
-        }
+            {
-
+                MatrixKernelEnumerator* e = (MatrixKernelEnumerator*)msPtr;
-        private protected void AssertIndex(int column, int row)
+                return ref _boarder.IsInBuffer(x, y, e) ?
-        {
+                           ref _buffer[_boarder.GetIndexOffset(x, y, e)] :
-            if (column < 0 || column >= ColumnLength) {
+                           ref Unsafe.NullRef<MatrixKernel<T>>();
-                throw new ArgumentOutOfRangeException(nameof(column));
-            }
-
-            if (row < 0 || row >= RowLength) {
-                throw new ArgumentOutOfRangeException(nameof(row));
             }
         }
 

StgSharp.Common/Mathematics/Linear/MatrixMetadata.cs

@@ -46,7 +46,7 @@ namespace StgSharp.Mathematics.Numeric
     }
 
     [StructLayout(LayoutKind.Sequential)]
-    public unsafe struct MatrixKernelEnumeration
+    public unsafe struct MatrixKernelBoarder
     {
 
         public readonly delegate* unmanaged[Cdecl]
@@ -55,9 +55,9 @@ namespace StgSharp.Mathematics.Numeric
         public readonly delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> IsInBuffer;
 
         [Obsolete("Do not use constructor, use static method instead", true)]
-        public MatrixKernelEnumeration() { }
+        public MatrixKernelBoarder() { }
 
-        public MatrixKernelEnumeration(
+        public MatrixKernelBoarder(
                delegate* unmanaged[Cdecl]
            <int, int, MatrixKernelEnumerator*, long> getOffset,
                delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> inBuffer)
@@ -66,24 +66,128 @@ namespace StgSharp.Mathematics.Numeric
             IsInBuffer = inBuffer;
         }
 
-        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> IsInBufferDenseSquarePtr { get; } = (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool>)&IsInBufferDenseSquare;
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long> GetIndexOffsetDenseSquarePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long>)&GetIndexOffsetDenseSquare;
 
-        public static MatrixKernelEnumeration Create(MatrixLayout layout)
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long> GetIndexOffsetUpperTrianglePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long>)&GetIndexOffsetUpperTriangle;
+
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long> GetIndexOffsetLowerTrianglePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, long>)&GetIndexOffsetLowerTriangle;
+
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> IsInBufferDenseSquarePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool>)&IsInBufferDenseSquare;
+
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> IsInBufferUpperTrianglePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool>)&IsInBufferUpperTriangle;
+
+        private static delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool> IsInBufferLowerTrianglePtr { get; } =
+            (delegate* unmanaged[Cdecl]<int, int, MatrixKernelEnumerator*, bool>)&IsInBufferLowerTriangle;
+
+        public static MatrixKernelBoarder Create(MatrixLayout layout)
         {
             return layout switch
             {
-                MatrixLayout.DenseSquare => new MatrixKernelEnumeration(null, IsInBufferDenseSquarePtr),
+                MatrixLayout.DenseSquare => new MatrixKernelBoarder(GetIndexOffsetDenseSquarePtr, IsInBufferDenseSquarePtr),
 
-                MatrixLayout.UpperTriangle => throw new NotImplementedException(),
+                MatrixLayout.UpperTriangle => new MatrixKernelBoarder(GetIndexOffsetUpperTrianglePtr, IsInBufferUpperTrianglePtr),
+                MatrixLayout.LowerTriangle => new MatrixKernelBoarder(GetIndexOffsetLowerTrianglePtr, IsInBufferLowerTrianglePtr),
+                _ => throw new ArgumentOutOfRangeException(nameof(layout)),
             };
         }
 
-        private static bool IsInBufferDenseSquare(
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
-                            int col,
+        private static long GetIndexOffsetDenseSquare(
-                            int row,
+                            int x_kernel,
+                            int y_kernel,
                             MatrixKernelEnumerator* enumerator)
         {
-            return true;
+            long x = x_kernel;
+            long y = y_kernel;
+            long rows = enumerator->KernelRowLength;
+            return x * rows + y;
+        }
+
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
+        private static long GetIndexOffsetLowerTriangle(
+                            int x_kernel,
+                            int y_kernel,
+                            MatrixKernelEnumerator* enumerator)
+        {
+            long order = GetTriangularOrder(enumerator);
+            long x = x_kernel;
+            long y = y_kernel;
+
+            long prefixColumns = x * order - (x * (x - 1L)) / 2L;
+            return prefixColumns + (y - x);
+        }
+
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
+        private static long GetIndexOffsetUpperTriangle(
+                            int x_kernel,
+                            int y_kernel,
+                            MatrixKernelEnumerator* enumerator)
+        {
+            long x = x_kernel;
+            long y = y_kernel;
+
+            long prefixColumns = (x * (x + 1L)) / 2L;
+            return prefixColumns + y;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static long GetTriangularOrder(MatrixKernelEnumerator* enumerator)
+        {
+            return Math.Min(enumerator->KernelColumnLength, enumerator->KernelRowLength);
+        }
+
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
+        private static bool IsInBufferDenseSquare(
+                            int x_kernel,
+                            int y_kernel,
+                            MatrixKernelEnumerator* enumerator)
+        {
+            return IsWithinKernelDomain(enumerator, x_kernel, y_kernel);
+        }
+
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
+        private static bool IsInBufferLowerTriangle(
+                            int x_kernel,
+                            int y_kernel,
+                            MatrixKernelEnumerator* enumerator)
+        {
+            return IsWithinTriangularDomain(enumerator, x_kernel, y_kernel) && y_kernel >= x_kernel;
+        }
+
+        [UnmanagedCallersOnly(CallConvs =[ typeof(CallConvCdecl) ])]
+        private static bool IsInBufferUpperTriangle(
+                            int x_kernel,
+                            int y_kernel,
+                            MatrixKernelEnumerator* enumerator)
+        {
+            return IsWithinTriangularDomain(enumerator, x_kernel, y_kernel) && y_kernel <= x_kernel;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static bool IsWithinKernelDomain(
+                            MatrixKernelEnumerator* enumerator,
+                            int x_kernel,
+                            int y_kernel)
+        {
+            return x_kernel >= 0 &&
+                   y_kernel >= 0 &&
+                   x_kernel < enumerator->KernelColumnLength &&
+                   y_kernel < enumerator->KernelRowLength;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static bool IsWithinTriangularDomain(
+                            MatrixKernelEnumerator* enumerator,
+                            int x_kernel,
+                            int y_kernel)
+        {
+            long order = GetTriangularOrder(enumerator);
+            return x_kernel >= 0 && y_kernel >= 0 && x_kernel < order && y_kernel < order;
         }
 
     }
@@ -98,6 +202,18 @@ namespace StgSharp.Mathematics.Numeric
         [FieldOffset(8)] public readonly int KernelColumnLength;
         [FieldOffset(12)] public readonly int KernelRowLength;
 
+        public MatrixSize(
+               int elementColumnLength,
+               int elementRowLength,
+               int kernelColumnLength,
+               int kernelRowLength)
+        {
+            ElementColumnLength = elementColumnLength;
+            ElementRowLength = elementRowLength;
+            KernelColumnLength = kernelColumnLength;
+            KernelRowLength = kernelRowLength;
+        }
+
     }
 
     [StructLayout(LayoutKind.Explicit)]

StgSharp.Common/Mathematics/Linear/MatrixStorage.cs

@@ -25,11 +25,13 @@
 //     
 // -----------------------------------------------------------------------
 // -----------------------------------------------------------------------
+using StgSharp.HighPerformance;
 using StgSharp.HighPerformance.Memory;
 using System;
 using System.Collections.Generic;
 using System.Linq;
 using System.Numerics;
+using System.Runtime.CompilerServices;
 using System.Text;
 using System.Threading.Tasks;
 
@@ -41,6 +43,11 @@ namespace StgSharp.Mathematics.Numeric
     public abstract class MatrixStorage<T> : IDisposable where T : unmanaged, INumber<T>
     {
 
+        internal unsafe ref MatrixKernel<T> this[long index]
+        {
+            get { return ref Unsafe.AsRef<MatrixKernel<T>>(BufferPointer + MatrixKernel<T>.Size * index); }
+        }
+
         protected internal abstract unsafe T* BufferPointer { get; }
 
         public abstract void Dispose();

StgSharp.Native/StgSharp.Native.vcxproj

@@ -151,7 +151,7 @@
     <ClCompile Include="sn-src\k_panel.c" />
     <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_kerop.c" />
     <ClCompile Include="sn-src\k_pivot.c" />
     <ClCompile Include="sn-src\str_factorial.c" />
     <ClCompile Include="sn-src\k_fma.c" />

StgSharp.Native/StgSharp.Native.vcxproj.filters

@@ -48,7 +48,7 @@
     <ClCompile Include="sn-src\load_intrinsic.c">
       <Filter>源文件</Filter>
     </ClCompile>
-    <ClCompile Include="sn-src\k_bufop.c">
+    <ClCompile Include="sn-src\k_kerop.c">
       <Filter>源文件</Filter>
     </ClCompile>
     <ClCompile Include="sn-src\k_panel.c">

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -96,40 +96,6 @@ static_assert(sizeof(matrix_task) == 128, "matrix_parallel_task_package size mis
 _Static_assert(sizeof(MAT_TASK(float)) == 128, "matrix_parallel_task_package size mismatch");
 #endif
 
-#pragma region matrix buffer operation function type
-
-typedef void(SN_DECL *BUF_PROC_AS)(mat_kernel *ans, scalar_pack *scalar, int count);
-typedef void(SN_DECL *BUF_PROC_RAS)(mat_kernel *right, mat_kernel *ans, scalar_pack *scalar,
-                                    int count);
-typedef void(SN_DECL *BUF_PROC_LRAS)(mat_kernel *left, mat_kernel *right, mat_kernel *ans,
-                                     scalar_pack *scalar, int count);
-typedef void(SN_DECL *BUF_PROC_RA)(mat_kernel *right, mat_kernel *ans, int count);
-typedef void(SN_DECL *BUF_PROC_LRA)(mat_kernel *left, mat_kernel *right, mat_kernel *ans,
-                                    int count);
-
-#define BUF_PROC(T_type, T_arch, T_op_name) T_type##_buf_##T_op_name##_##T_arch
-
-#define DECLARE_BUF_PROC_ANS_SCALAR(T_type, T_arch, T_op_name)     \
-        INTERNAL void SN_DECL BUF_PROC(T_type, T_arch, T_op_name)( \
-                MAT_KERNEL(T_type) *restrict ans, scalar_pack const *scalar, size_t count)
-#define DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(T_type, T_arch, T_op_name)               \
-        INTERNAL void SN_DECL BUF_PROC(T_type, T_arch, T_op_name)(                 \
-                MAT_KERNEL(T_type) const *right, MAT_KERNEL(T_type) *restrict ans, \
-                scalar_pack const *scalar, int count)
-#define DECLARE_BUF_PROC_LEFT_RIGHT_ANS_SCALAR(T_type, T_arch, T_op_name)        \
-        INTERNAL void SN_DECL BUF_PROC(T_type, T_arch, T_op_name)(               \
-                MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \
-                MAT_KERNEL(T_type) *restrict ans, scalar_pack const *scalar, size_t count)
-#define DECLARE_BUF_PROC_RIGHT_ANS(T_type, T_arch, T_op_name)      \
-        INTERNAL void SN_DECL BUF_PROC(T_type, T_arch, T_op_name)( \
-                MAT_KERNEL(T_type) const *right, MAT_KERNEL(T_type) *restrict ans, size_t count)
-#define DECLARE_BUF_PROC_LEFT_RIGHT_ANS(T_type, T_arch, T_op_name)               \
-        INTERNAL void SN_DECL BUF_PROC(T_type, T_arch, T_op_name)(               \
-                MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \
-                MAT_KERNEL(T_type) *restrict ans, size_t count)
-
-#pragma endregion
-
 #pragma region matrix panel function type
 
 typedef void(SN_DECL *PANEL_PROC)(mat_panel *panel, mat_kernel *kernel, int col_length,
@@ -176,6 +142,9 @@ typedef void(SN_DECL *KER_PROC_LAS)(mat_kernel *matrix, mat_kernel *ans, scalar_
 
 #define KER_PROC(T_type, T_arch, T_op_name) T_type##_kernel_##T_op_name##_##T_arch
 
+#define DECLARE_KER_PROC_ANS_SCALAR(T_type, T_arch, T_op_name)     \
+        INTERNAL void SN_DECL KER_PROC(T_type, T_arch, T_op_name)( \
+                MAT_KERNEL(T_type) *restrict ans, scalar_pack const *scalar)
 #define DECLARE_KER_PROC_LEFT_RIGHT_ANS_SCALAR(T_type, T_arch, T_op_name)        \
         INTERNAL void SN_DECL KER_PROC(T_type, T_arch, T_op_name)(               \
                 MAT_KERNEL(T_type) const *left, MAT_KERNEL(T_type) const *right, \
@@ -216,25 +185,25 @@ typedef uint64_t(SN_DECL *FACTORIALROC)(int n);
 
 INTERNAL void f32_normalize(VEC(float) * source, VEC(float) * target);
 
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, add);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, sse, add);
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, add);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, avx, add);
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, add);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, 512, add);
 
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, sub);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, sse, sub);
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, sub);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, avx, sub);
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, sub);
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, 512, sub);
 
 DECLARE_KER_PROC_RIGHT_ANS(float, sse, transpose);
 DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose);
 DECLARE_KER_PROC_RIGHT_ANS(float, 512, transpose);
 
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul);
+DECLARE_KER_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul);
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul);
+DECLARE_KER_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul);
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul);
+DECLARE_KER_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul);
 
-DECLARE_BUF_PROC_ANS_SCALAR(float, sse, fill);
+DECLARE_KER_PROC_ANS_SCALAR(float, sse, fill);
-DECLARE_BUF_PROC_ANS_SCALAR(float, avx, fill);
+DECLARE_KER_PROC_ANS_SCALAR(float, avx, fill);
-DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill);
+DECLARE_KER_PROC_ANS_SCALAR(float, 512, fill);
 
 #pragma endregion
 
@@ -280,6 +249,18 @@ INTERNAL int SN_DECL index_pair_sse(short const *str, uint32_t target, int lengt
 
 #pragma endregion
 
+typedef struct sn_mat_enumerator {
+        int CurrentColumn;
+        int CurrentRow;
+        int KernelColumnLength;
+        int KernelRowLength;
+} sn_mat_enumerator;
+
+typedef struct sn_mat_boarder {
+        int (*SN_DECL is_in_buffer)(int col, int row, sn_mat_enumerator *enumerator);
+        long (*SN_DECL get_index_offset)(int col, int row, sn_mat_enumerator *wnumerator);
+} sn_mat_boarder;
+
 // Mirrors C# IntrinsicContext field order so managed/unmanaged layouts stay in sync.
 typedef struct sn_intrinsic {
         HASH city_hash_simplify;

StgSharp.Native/sn-src/k_bufop.c

@@ -1,134 +0,0 @@
-#include "sn_intrinsic.h"
-
-DECLARE_BUF_PROC_ANS_SCALAR(float, sse, fill)
-{
-        register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]);
-        for (size_t i = 0; i < count; i++) {
-                ans[i].xmm[0] = vec;
-                ans[i].xmm[1] = vec;
-                ans[i].xmm[2] = vec;
-                ans[i].xmm[3] = vec;
-        }
-}
-
-DECLARE_BUF_PROC_ANS_SCALAR(float, avx, fill)
-{
-        register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data)));
-        for (size_t i = 0; i < count; i++) {
-                ans[i].ymm[0] = vec;
-                ans[i].ymm[1] = vec;
-        }
-}
-
-DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill)
-{
-        register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]);
-        for (size_t i = 0; i < count; i++) {
-                ans[i].zmm[0] = vec;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, add)
-{
-        for (size_t i = 0; i < count; i++) {
-                register __m128 c0 = _mm_add_ps(left[i].xmm[0], right[i].xmm[0]);
-                register __m128 c1 = _mm_add_ps(left[i].xmm[1], right[i].xmm[1]);
-                register __m128 c2 = _mm_add_ps(left[i].xmm[2], right[i].xmm[2]);
-                register __m128 c3 = _mm_add_ps(left[i].xmm[3], right[i].xmm[3]);
-
-                ans[i].xmm[0] = c0;
-                ans[i].xmm[1] = c1;
-                ans[i].xmm[2] = c2;
-                ans[i].xmm[3] = c3;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, add)
-{
-        for (size_t i = 0; i < count; i++) {
-                size_t s = sizeof(MAT_KERNEL(float));
-                register __m256 c0 = _mm256_add_ps(left[i].ymm[0], right[i].ymm[0]);
-                register __m256 c1 = _mm256_add_ps(left[i].ymm[1], right[i].ymm[1]);
-                ans[i].ymm[0] = c0;
-                ans[i].ymm[1] = c1;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, add)
-{
-        for (size_t i = 0; i < count; i++) {
-                register __m512 c0 = _mm512_add_ps(left[i].zmm[0], right[i].zmm[0]);
-                ans[i].zmm[0] = c0;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, sub)
-{
-        for (size_t i = 0; i < count; i++) {
-                register __m128 c0 = _mm_sub_ps(left[i].xmm[0], right[i].xmm[0]);
-                register __m128 c1 = _mm_sub_ps(left[i].xmm[1], right[i].xmm[1]);
-                register __m128 c2 = _mm_sub_ps(left[i].xmm[2], right[i].xmm[2]);
-                register __m128 c3 = _mm_sub_ps(left[i].xmm[3], right[i].xmm[3]);
-
-                ans[i].xmm[0] = c0;
-                ans[i].xmm[1] = c1;
-                ans[i].xmm[2] = c2;
-                ans[i].xmm[3] = c3;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, sub)
-{
-        for (size_t i = 0; i < count; i++) {
-                register __m256 c0 = _mm256_sub_ps(left[i].ymm[0], right[i].ymm[0]);
-                register __m256 c1 = _mm256_sub_ps(left[i].ymm[1], right[i].ymm[1]);
-                ans[i].ymm[0] = c0;
-                ans[i].ymm[1] = c1;
-        }
-}
-
-DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, sub)
-{
-        for (size_t i = 0; i < count; i++) {
-                register __m512 c0 = _mm512_sub_ps(left[i].zmm[0], right[i].zmm[0]);
-                ans[i].zmm[0] = c0;
-        }
-}
-
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul)
-{
-        register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]);
-        for (size_t i = 0; i < count; i++) {
-                register __m128 c0 = _mm_mul_ps(right[i].xmm[0], vec);
-                register __m128 c1 = _mm_mul_ps(right[i].xmm[1], vec);
-                register __m128 c2 = _mm_mul_ps(right[i].xmm[2], vec);
-                register __m128 c3 = _mm_mul_ps(right[i].xmm[3], vec);
-
-                ans[i].xmm[0] = c0;
-                ans[i].xmm[1] = c1;
-                ans[i].xmm[2] = c2;
-                ans[i].xmm[3] = c3;
-        }
-}
-
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul)
-{
-        register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data)));
-        for (size_t i = 0; i < count; i++) {
-                register __m256 c0 = _mm256_mul_ps(right[i].ymm[0], vec);
-                register __m256 c1 = _mm256_mul_ps(right[i].ymm[1], vec);
-
-                ans[i].ymm[0] = c0;
-                ans[i].ymm[1] = c1;
-        }
-}
-
-DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul)
-{
-        register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]);
-        for (size_t i = 0; i < count; i++) {
-                register __m512 c0 = _mm512_mul_ps(right[i].zmm[0], vec);
-
-                ans[i].zmm[0] = c0;
-        }
-}

StgSharp.Native/sn-src/k_kerop.c

@@ -0,0 +1,113 @@
+#include "sn_intrinsic.h"
+
+DECLARE_KER_PROC_ANS_SCALAR(float, sse, fill)
+{
+        register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]);
+        ans->xmm[0] = vec;
+        ans->xmm[1] = vec;
+        ans->xmm[2] = vec;
+        ans->xmm[3] = vec;
+}
+
+DECLARE_KER_PROC_ANS_SCALAR(float, avx, fill)
+{
+        register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data)));
+        ans->ymm[0] = vec;
+        ans->ymm[1] = vec;
+}
+
+DECLARE_KER_PROC_ANS_SCALAR(float, 512, fill)
+{
+        register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]);
+        ans[i].zmm[0] = vec;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, sse, add)
+{
+        register __m128 c0 = _mm_add_ps(left->xmm[0], right->xmm[0]);
+        register __m128 c1 = _mm_add_ps(left->xmm[1], right->xmm[1]);
+        register __m128 c2 = _mm_add_ps(left->xmm[2], right->xmm[2]);
+        register __m128 c3 = _mm_add_ps(left->xmm[3], right->xmm[3]);
+
+        ans->xmm[0] = c0;
+        ans->xmm[1] = c1;
+        ans->xmm[2] = c2;
+        ans->xmm[3] = c3;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, avx, add)
+{
+        register __m256 c0 = _mm256_add_ps(left->ymm[0], right->ymm[0]);
+        register __m256 c1 = _mm256_add_ps(left->ymm[1], right->ymm[1]);
+        ans->ymm[0] = c0;
+        ans->ymm[1] = c1;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, 512, add)
+{
+        register __m512 c0 = _mm512_add_ps(left->zmm[0], right->zmm[0]);
+        ans->zmm[0] = c0;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, sse, sub)
+{
+        register __m128 c0 = _mm_sub_ps(left->xmm[0], right->xmm[0]);
+        register __m128 c1 = _mm_sub_ps(left->xmm[1], right->xmm[1]);
+        register __m128 c2 = _mm_sub_ps(left->xmm[2], right->xmm[2]);
+        register __m128 c3 = _mm_sub_ps(left->xmm[3], right->xmm[3]);
+
+        ans->xmm[0] = c0;
+        ans->xmm[1] = c1;
+        ans->xmm[2] = c2;
+        ans->xmm[3] = c3;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, avx, sub)
+{
+        register __m256 c0 = _mm256_sub_ps(left->ymm[0], right->ymm[0]);
+        register __m256 c1 = _mm256_sub_ps(left->ymm[1], right->ymm[1]);
+        ans->ymm[0] = c0;
+        ans->ymm[1] = c1;
+}
+
+DECLARE_KER_PROC_LEFT_RIGHT_ANS(float, 512, sub)
+{
+        register __m512 c0 = _mm512_sub_ps(left->zmm[0], right->zmm[0]);
+        ans->zmm[0] = c0;
+}
+
+DECLARE_KER_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul)
+{
+        register __m128 vec = _mm_set_ps1(((float *)(scalar->data))[0]);
+        register __m128 c0 = _mm_mul_ps(right->xmm[0], vec);
+        register __m128 c1 = _mm_mul_ps(right->xmm[1], vec);
+        register __m128 c2 = _mm_mul_ps(right->xmm[2], vec);
+        register __m128 c3 = _mm_mul_ps(right->xmm[3], vec);
+
+        ans->xmm[0] = c0;
+        ans->xmm[1] = c1;
+        ans->xmm[2] = c2;
+        ans->xmm[3] = c3;
+}
+
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul)
+{
+        register __m256 vec = _mm256_broadcast_ss(((float *)(scalar->data)));
+        for (size_t i = 0; i < count; i++) {
+                register __m256 c0 = _mm256_mul_ps(right[i].ymm[0], vec);
+                register __m256 c1 = _mm256_mul_ps(right[i].ymm[1], vec);
+
+                ans[i].ymm[0] = c0;
+                ans[i].ymm[1] = c1;
+        }
+}
+
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul)
+{
+        register __m512 vec = _mm512_broadcastss_ps(((__m128 *)(scalar->data))[0]);
+        for (size_t i = 0; i < count; i++) {
+                register __m512 c0 = _mm512_mul_ps(right[i].zmm[0], vec);
+
+                ans[i].zmm[0] = c0;
+        }
+}

StgSharpDebug/Program.cs

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