add double compute kernel

StgSharp.Common/HighPerformance/MatrixPanel.cs

@@ -61,17 +61,9 @@ namespace StgSharp.HighPerformance
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static int ElementSideLength<T>() where T : unmanaged, INumber<T>
         {
-            return Math.Max(8, AlignmentBytes / sizeof(T));
+            return PanelSideLength<T>();
         }
 
-        // View an existing pointer; caller guarantees lifetime/alignment.
-        public static unsafe MatrixPanel<T>* FromPointer<T>(void* ptr)
-            where T : unmanaged, INumber<T>
-        {
-            return (MatrixPanel<T>*)Align64(ptr);
-        }
-
-        [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static int MaxElementCount<T>() where T : unmanaged, INumber<T>
         {
             int side = ElementSideLength<T>();
@@ -115,28 +107,43 @@ namespace StgSharp.HighPerformance
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         private static int KernelSideCountAvx2<T>() where T : unmanaged, INumber<T>
         {
-            const int minQ = 8;
+            return PanelSideFromBits<T>(256) / 4;
-            int lanes = 32 / sizeof(T);
-            int side = Math.Max(lanes, minQ);
-            return side / 4;
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         private static int KernelSideCountAvx512<T>() where T : unmanaged, INumber<T>
         {
-            const int minQ = 8;
+            return PanelSideFromBits<T>(512) / 4;
-            int lanes = 64 / sizeof(T);
-            int side = Math.Max(lanes, minQ);
-            return side / 4;
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         private static int KernelSideCountSse<T>() where T : unmanaged, INumber<T>
         {
-            const int minQ = 4;
+            return PanelSideFromBits<T>(128) / 4;
-            int lanes = 16 / sizeof(T);
+        }
-            int side = Math.Max(lanes, minQ);
+
-            return side / 4;
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int PanelSideFromBits<T>(int vectorBits) where T : unmanaged, INumber<T>
+        {
+            int elemBits = sizeof(T) * 8;
+
+            // Q = max(4, vectorBits/64, vectorBits/elemBits)
+            int by64 = vectorBits / 64;
+            int byElem = vectorBits / elemBits;
+            int side = Math.Max(4, Math.Max(by64, byElem));
+            return side;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int PanelSideLength<T>() where T : unmanaged, INumber<T>
+        {
+            return NativeIntrinsic.IntrinsicLevel switch
+            {
+                IntrinsicLevel.SSE => PanelSideFromBits<T>(128),
+                IntrinsicLevel.AVX2 => PanelSideFromBits<T>(256),
+                IntrinsicLevel.AVX512 => PanelSideFromBits<T>(512),
+                _ => 0
+            };
         }
 
     }

StgSharp.Common/Internal/Intrinsic/IntrinsicContext.cs

@@ -39,20 +39,29 @@ namespace StgSharp.Internal.Intrinsic
     {
 
         public delegate* unmanaged[Cdecl]<char*, int, int> city_hash_simplify;              //city_hash_simplify
-        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> f32_buffer_add;             //f32_buffer_add
-        public delegate* unmanaged[Cdecl]<void*, ScalarPacket*,long, void> f32_buffer_fill;                //f32_buffer_fill
-        public delegate* unmanaged[Cdecl]<void*, void*, ScalarPacket*, long, void> f32_buffer_scalar_mul;                //f32_buffer_s_m
-        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> f32_buffer_sub;             //f32_kernel_sub
-        public delegate* unmanaged[Cdecl]<void*, void*, void> f32_buffer_transpose;                //f32_kernel_transpose
-        public delegate* unmanaged[Cdecl]<MatrixPanel*, MatrixKernel*,  int, int, int, int, void> f32_build_panel;              //f32_build_panel
-        public delegate* unmanaged[Cdecl]<MatrixPanel*, void> f32_clear_panel;
-        public delegate* unmanaged[Cdecl]<Vector4*, Vector4*, void> f32_normalize_3;                //f32_normalize_3
-        public delegate* unmanaged[Cdecl]<void*, void*, void*, void> f32_panel_fma;               //f32_panel_fma
-        public delegate* unmanaged[Cdecl]<void*, ulong, ulong, ulong, void> f32_pivot;               //f32_panel_mul
-        public delegate* unmanaged[Cdecl]<void*, void*, int, int, int, int, void> f32_store_panel; //f32_store_panel
         public delegate* unmanaged[Cdecl]<int, ulong> factorial_simd;               //factorial_simd
         public delegate* unmanaged[Cdecl]<char*, int, int, int> index_pair;             //index_pair
-        public delegate* unmanaged[Cdecl]<void*, void*, ScalarPacket*, long, void> variant ;              //variant 
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> mat_f32_buffer_add;             //f32_buffer_add
+        public delegate* unmanaged[Cdecl]<void*, ScalarPacket*, long, void> mat_f32_buffer_fill;                //f32_buffer_fill
+        public delegate* unmanaged[Cdecl]<void*, void*, ScalarPacket*, long, void> mat_f32_buffer_scalar_mul;                //f32_buffer_s_m
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> mat_f32_buffer_sub;             //f32_kernel_sub
+        public delegate* unmanaged[Cdecl]<void*, void*, void> mat_f32_buffer_transpose;                //f32_kernel_transpose
+        public delegate* unmanaged[Cdecl]<MatrixPanel*, MatrixKernel*, int, int, int, int, void> mat_f32_build_panel;              //f32_build_panel
+        public delegate* unmanaged[Cdecl]<MatrixPanel*, void> mat_f32_clear_panel;                      //f32_clear_panel
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, void> mat_f32_panel_fma;               //f32_panel_fma
+        public delegate* unmanaged[Cdecl]<void*, ulong, ulong, ulong, void> mat_f32_pivot;               //f32_panel_mul
+        public delegate* unmanaged[Cdecl]<void*, void*, int, int, int, int, void> mat_f32_store_panel; //f32_store_panel
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> mat_f64_buffer_add;             //f32_buffer_add
+        public delegate* unmanaged[Cdecl]<void*, ScalarPacket*, long, void> mat_f64_buffer_fill;                //f32_buffer_fill
+        public delegate* unmanaged[Cdecl]<void*, void*, ScalarPacket*, long, void> mat_f64_buffer_scalar_mul;                //f32_buffer_s_m
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, long, void> mat_f64_buffer_sub;             //f32_kernel_sub
+        public delegate* unmanaged[Cdecl]<void*, void*, void> mat_f64_buffer_transpose;                //f32_kernel_transpose
+        public delegate* unmanaged[Cdecl]<MatrixPanel*, MatrixKernel*, int, int, int, int, void> mat_f64_build_panel;              //f32_build_panel
+        public delegate* unmanaged[Cdecl]<MatrixPanel*, void> mat_f64_clear_panel;                      //f32_clear_panel
+        public delegate* unmanaged[Cdecl]<void*, void*, void*, void> mat_f64_panel_fma;               //f32_panel_fma
+        public delegate* unmanaged[Cdecl]<void*, ulong, ulong, ulong, void> mat_f64_pivot;               //f32_panel_mul
+        public delegate* unmanaged[Cdecl]<void*, void*, int, int, int, int, void> mat_f64_store_panel; //f32_store_panel
+        public delegate* unmanaged[Cdecl]<Vector4*, Vector4*, void> vec_f32_normalize_3;                //f32_normalize_3
 
         public IntrinsicContext()
         {

StgSharp.Common/Mathematics/Graphic/Matrix44.cs

@@ -169,7 +169,7 @@ namespace StgSharp.Mathematics.Graphic
         internal unsafe void InternalTranspose(ref Matrix44 transpose)
         {
             fixed (MatrixKernel<float>* source = &this.kernel, target = &transpose.kernel) {
-                NativeIntrinsic.Context.f32_buffer_transpose(source, target);
+                NativeIntrinsic.Context.mat_f32_buffer_transpose(source, target);
             }
         }
 
@@ -177,7 +177,7 @@ namespace StgSharp.Mathematics.Graphic
         public static unsafe Matrix44 operator -(Matrix44 left, Matrix44 right)
         {
             Matrix44 ret = new();
-            NativeIntrinsic.Context.f32_buffer_sub(&left, &right, &ret, 1);
+            NativeIntrinsic.Context.mat_f32_buffer_sub(&left, &right, &ret, 1);
             return ret;
         }
 
@@ -219,7 +219,7 @@ namespace StgSharp.Mathematics.Graphic
         public static unsafe Matrix44 operator +(Matrix44 left, Matrix44 right)
         {
             Matrix44 ret = new();
-            NativeIntrinsic.Context.f32_buffer_add(&left, &right, &ret, 1);
+            NativeIntrinsic.Context.mat_f32_buffer_add(&left, &right, &ret, 1);
             return ret;
         }
 

StgSharp.Common/Mathematics/Linear/Matrix.cs

@@ -50,14 +50,17 @@ namespace StgSharp.Mathematics.Numeric
                           MatrixSize size,
                           MatrixStorage<T> buffer
 
-            )
+        )
         {
             _boarder = boarder;
             _size = size;
             _buffer = buffer;
         }
 
-        public T this[int x, int y]
+        public T this[
+                 int x,
+                 int y
+        ]
         {
             get
             {
@@ -108,7 +111,8 @@ namespace StgSharp.Mathematics.Numeric
         /// </exception>
         public unsafe Matrix<T> this[
                                 ReadOnlySpan<MatrixIndex> column,
-                                ReadOnlySpan<MatrixIndex> row]
+                                ReadOnlySpan<MatrixIndex> row
+        ]
         {
             get { throw new NotImplementedException(); }
             set { throw new NotImplementedException(); }
@@ -137,7 +141,10 @@ namespace StgSharp.Mathematics.Numeric
         protected internal int KernelRowLength => _size.KernelRowLength;
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public unsafe ref MatrixKernel<T> GetMatrixKernel(int x, int y)
+        public unsafe ref MatrixKernel<T> GetMatrixKernel(
+                                          int x,
+                                          int y
+        )
         {
             MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength,
                                                     x, y);
@@ -147,7 +154,10 @@ namespace StgSharp.Mathematics.Numeric
         }
 
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        internal unsafe ref T GetElementUnsafe(int x, int y)
+        internal unsafe ref T GetElementUnsafe(
+                              int x,
+                              int y
+        )
         {
             // Caller guarantees: indices in range, layout is dense square.
             int k_c = y >> 2;
@@ -161,5 +171,37 @@ namespace StgSharp.Mathematics.Numeric
             return ref kernel[c, r];
         }
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe MatrixKernel<T>* GetKernelAddressUnsafe(
+                                         int x,
+                                         int y
+        )
+        {
+            // Caller guarantees: indices in range, layout is dense square.
+            int k_c = y >> 2;
+            int k_r = x >> 2;
+
+            MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength,
+                                                    k_c, k_r);
+            return Buffer + _boarder.GetIndexOffset(k_c, k_r, in enumerator);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe ref MatrixKernel<T> GetKernelUnsafe(
+                                            int x,
+                                            int y
+        )
+        {
+            // 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;
+
+            MatrixKernelEnumerator enumerator = new(_size.KernelColumnLength, _size.KernelRowLength,
+                                                    k_c, k_r);
+            return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in enumerator)];
+        }
+
     }
 }

StgSharp.Mathematics/Numeric/MatrixCompute.Additional.cs

@@ -30,6 +30,7 @@ using StgSharp.Internal;
 using System;
 using System.Collections.Generic;
 using System.Linq;
+using System.Numerics;
 using System.Runtime.CompilerServices;
 using System.Text;
 using System.Threading.Tasks;
@@ -39,45 +40,84 @@ namespace StgSharp.Mathematics.Numeric
     public partial class MatrixCompute
     {
 
-        internal static unsafe void ComputeFmaPanelF32(MatrixParallelTaskPackage* p)
+        internal static unsafe void ComputeFmaPanelF32(
+                                    MatrixParallelTaskPackage* p
+        )
         {
             int size = p->ElementSize;
             MatrixKernel* leftBuffer = (MatrixKernel*)p->Left;
             MatrixKernel* rightBuffer = (MatrixKernel*)p->Right;
             MatrixKernel* ansBuffer = (MatrixKernel*)p->Answer;
-            MatrixPanel<float>* leftPanel = MatrixPanel.Create<float>();
-            MatrixPanel<float>* rightPanel = MatrixPanel.Create<float>();
-            MatrixPanel<float>* ansPanel = MatrixPanel.Create<float>();
             ref long offsetRef = ref Unsafe.As<int, long>(ref p->PrimTileOffset);
             ref long lengthRef = ref Unsafe.As<int, long>(ref p->PrimCount);
-            int ansWidth = p-> AnsPrimOffset;
+            int ansWidth = p->AnsPrimOffset;
             int ansHeight = p->AnsSecOffset;
             int commonK = p->LeftPrimOffset;
             int blocksPerRow = ansWidth / size;
-
+            if (NativeIntrinsic.Context.mat_f32_build_panel == null)
-            for (long c = offsetRef; c < offsetRef + lengthRef; c++)
             {
-                int i = (int)(c / blocksPerRow) * size;
+                // in this occasion size must be 1
-                int j = (int)(c % blocksPerRow) * size;
+                for (long c = offsetRef; c < offsetRef + lengthRef; c++)
-                NativeIntrinsic.Context.f32_clear_panel((MatrixPanel*)ansPanel);
-                for (int k = 0; k < commonK; k += size)
                 {
-                    NativeIntrinsic.Context
+                    int i = (int)(c / blocksPerRow);
-                                   .f32_build_panel((MatrixPanel*)leftPanel, leftBuffer, commonK,
+                    int j = (int)(c % blocksPerRow);
-                                                    ansHeight, k, i);
+                    MatrixKernel<float>* ansPanel = GetKernelAddressUnsafe((MatrixKernel<float>*)ansBuffer,
-                    NativeIntrinsic.Context
+                                                                           ansHeight, i, j);
-                                   .f32_build_panel((MatrixPanel*)rightPanel, rightBuffer, ansWidth,
+                    NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel);
-                                                    commonK, j, k);
+                    for (int k = 0; k < commonK; k += size) {
-                    NativeIntrinsic.Context
+                        NativeIntrinsic.Context
-                                   .f32_panel_fma((MatrixPanel*)leftPanel, (MatrixPanel*)rightPanel, (MatrixPanel*)ansPanel);
+                                       .mat_f32_panel_fma(
+                                        (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel<float>*)leftBuffer,
+                                                                             ansHeight, k, j), 
+                                        (MatrixPanel*)GetKernelAddressUnsafe((MatrixKernel<float>*)rightBuffer,
+                                                                             commonK, j, k), 
+                                        (MatrixPanel*)ansPanel);
+                    }
                 }
-                NativeIntrinsic.Context
+            } else
-                               .f32_store_panel((MatrixPanel*)ansPanel, ansBuffer, ansWidth,
+            {
-                                                ansHeight, j, i);
+                MatrixPanel<float>* leftPanel = MatrixPanel.Create<float>();
+                MatrixPanel<float>* rightPanel = MatrixPanel.Create<float>();
+                MatrixPanel<float>* ansPanel = MatrixPanel.Create<float>();
+                for (long c = offsetRef; c < offsetRef + lengthRef; c++)
+                {
+                    int i = (int)(c / blocksPerRow) * size;
+                    int j = (int)(c % blocksPerRow) * size;
+                    NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel);
+                    for (int k = 0; k < commonK; k += size)
+                    {
+                        NativeIntrinsic.Context
+                                       .mat_f32_build_panel((MatrixPanel*)leftPanel, leftBuffer,
+                                                            commonK, ansHeight, k, i);
+                        NativeIntrinsic.Context
+                                       .mat_f32_build_panel((MatrixPanel*)rightPanel, rightBuffer,
+                                                            ansWidth, commonK, j, k);
+                        NativeIntrinsic.Context
+                                       .mat_f32_panel_fma((MatrixPanel*)leftPanel, (MatrixPanel*)rightPanel, (MatrixPanel*)ansPanel);
+                    }
+                    NativeIntrinsic.Context
+                                   .mat_f32_store_panel((MatrixPanel*)ansPanel, ansBuffer, ansWidth,
+                                                        ansHeight, j, i);
+                }
+                MatrixPanel.Destroy(leftPanel);
+                MatrixPanel.Destroy(rightPanel);
+                MatrixPanel.Destroy(ansPanel);
             }
-            MatrixPanel.Destroy(leftPanel);
+        }
-            MatrixPanel.Destroy(rightPanel);
+
-            MatrixPanel.Destroy(ansPanel);
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static unsafe MatrixKernel<T>* GetKernelAddressUnsafe<T>(
+                                                MatrixKernel<T>* buffer,
+                                                int colLength,
+                                                int x,
+                                                int y
+        ) where T : unmanaged, INumber<T>
+        {
+            // Caller guarantees: indices in range, layout is dense square.
+            int k_c = y >> 2;
+            int k_r = x >> 2;
+
+            return buffer + colLength * x + y;
         }
 
     }

StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs

@@ -77,9 +77,9 @@ namespace StgSharp.Mathematics.Numeric
                 for (int kBlock = 0; kBlock < n; kBlock += 4)
                 {
                     NativeIntrinsic.Context
-                                   .f32_pivot(matrix.Buffer, (ulong)(kBlock >> 2),
+                                   .mat_f32_pivot(matrix.Buffer, (ulong)(kBlock >> 2),
-                                              (ulong)(matrix.KernelRowLength - (kBlock >> 2)),
+                                                  (ulong)(matrix.KernelRowLength - (kBlock >> 2)),
-                                              (ulong)pivots);
+                                                  (ulong)pivots);
                     int colsInBlock = Math.Min(4, n - kBlock);
                     for (int offset = 0; offset < colsInBlock; offset++)
                     {

StgSharp.Mathematics/Numeric/MatrixCompute.cs

@@ -56,7 +56,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
+                NativeIntrinsic.Context.mat_f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
 
                 return;
             } else if (count <= 1024)
@@ -66,7 +66,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
+                NativeIntrinsic.Context.mat_f32_buffer_add(leftPtr, rightPtr, ansPtr, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
@@ -79,7 +79,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.f32_buffer_add;
+            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_add;
             MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
@@ -94,12 +94,12 @@ namespace StgSharp.Mathematics.Numeric
             scalar->Data<float>(0) = value;
             if (count <= 256)
             {
-                NativeIntrinsic.Context.f32_buffer_fill(ans.Buffer, scalar, count);
+                NativeIntrinsic.Context.mat_f32_buffer_fill(ans.Buffer, scalar, count);
                 return;
             } else if (count <= 1024)
             {
                 MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
-                NativeIntrinsic.Context.f32_buffer_fill(ans.Buffer, scalar, count);
+                NativeIntrinsic.Context.mat_f32_buffer_fill(ans.Buffer, scalar, count);
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
             }
@@ -109,7 +109,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = ans.KernelColumnLength;
             package->SecCount = ans.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.f32_buffer_fill;
+            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_fill;
             package->Scalar = scalar;
             MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
@@ -135,43 +135,65 @@ namespace StgSharp.Mathematics.Numeric
                 if (count > 256) {
                     pool = MatrixParallel.LeadParallel(1);
                 }
-                MatrixPanel<float>* leftPanel = MatrixPanel.Create<float>();
+                if (NativeIntrinsic.Context.mat_f32_build_panel == null)
-                MatrixPanel<float>* rightPanel = MatrixPanel.Create<float>();
-                MatrixPanel<float>* ansPanel = MatrixPanel.Create<float>();
-                for (int i = 0; i < ans.KernelRowLength; i += size)
                 {
-                    for (int j = 0; j < ans.KernelColumnLength; j += size)
+                    for (int i = 0; i < ans.KernelRowLength; i += size)
                     {
-                        NativeIntrinsic.Context.f32_clear_panel((MatrixPanel*)ansPanel);
+                        for (int j = 0; j < ans.KernelColumnLength; j += size)
-
-                        // Console.WriteLine($"{i}{j}");
-                        for (int k = 0; k < left.KernelColumnLength; k += size)
                         {
-                            NativeIntrinsic.Context
+                            MatrixKernel<float>* ansPanel = ans.GetKernelAddressUnsafe(i, j);
-                                           .f32_build_panel((MatrixPanel*)leftPanel,
+                            NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel);
-                                                            (MatrixKernel*)left.Buffer,
+
-                                                            left.KernelColumnLength,
+                            // Console.WriteLine($"{i}{j}");
-                                                            left.KernelRowLength, i, k);
+                            for (int k = 0; k < left.KernelColumnLength; k += size)
-                            NativeIntrinsic.Context
+                            {
-                                           .f32_build_panel((MatrixPanel*)rightPanel,
+                                NativeIntrinsic.Context
-                                                            (MatrixKernel*)right.Buffer,
+                                               .mat_f32_panel_fma((MatrixPanel*)left.GetKernelAddressUnsafe(i, k),
-                                                            right.KernelColumnLength,
+                                                                (MatrixPanel*)right.GetKernelAddressUnsafe(k, j),
-                                                            right.KernelRowLength, k, j);
+                                                                (MatrixPanel*)ansPanel);
-                            NativeIntrinsic.Context
+                            }
-                                           .f32_panel_fma((MatrixPanel*)leftPanel,
-                                                            (MatrixPanel*)rightPanel,
-                                                            (MatrixPanel*)ansPanel);
                         }
-                        NativeIntrinsic.Context
-                                       .f32_store_panel((MatrixPanel*)ansPanel,
-                                                        (MatrixKernel*)ans.Buffer,
-                                                        ans.KernelColumnLength, ans.KernelRowLength,
-                                                        i, j);
                     }
+                } else
+                {
+                    MatrixPanel<float>* leftPanel = MatrixPanel.Create<float>();
+                    MatrixPanel<float>* rightPanel = MatrixPanel.Create<float>();
+                    MatrixPanel<float>* ansPanel = MatrixPanel.Create<float>();
+                    for (int i = 0; i < ans.KernelRowLength; i += size)
+                    {
+                        for (int j = 0; j < ans.KernelColumnLength; j += size)
+                        {
+                            NativeIntrinsic.Context.mat_f32_clear_panel((MatrixPanel*)ansPanel);
+
+                            // Console.WriteLine($"{i}{j}");
+                            for (int k = 0; k < left.KernelColumnLength; k += size)
+                            {
+                                NativeIntrinsic.Context
+                                               .mat_f32_build_panel((MatrixPanel*)leftPanel,
+                                                                    (MatrixKernel*)left.Buffer,
+                                                                    left.KernelColumnLength,
+                                                                    left.KernelRowLength, i, k);
+                                NativeIntrinsic.Context
+                                               .mat_f32_build_panel((MatrixPanel*)rightPanel,
+                                                                    (MatrixKernel*)right.Buffer,
+                                                                    right.KernelColumnLength,
+                                                                    right.KernelRowLength, k, j);
+                                NativeIntrinsic.Context
+                                               .mat_f32_panel_fma((MatrixPanel*)leftPanel,
+                                                                (MatrixPanel*)rightPanel,
+                                                                (MatrixPanel*)ansPanel);
+                            }
+                            NativeIntrinsic.Context
+                                           .mat_f32_store_panel((MatrixPanel*)ansPanel,
+                                                                (MatrixKernel*)ans.Buffer,
+                                                                ans.KernelColumnLength,
+                                                                ans.KernelRowLength, i, j);
+                        }
+                    }
+                    MatrixPanel.Destroy(leftPanel);
+                    MatrixPanel.Destroy(rightPanel);
+                    MatrixPanel.Destroy(ansPanel);
                 }
-                MatrixPanel.Destroy(leftPanel);
-                MatrixPanel.Destroy(rightPanel);
-                MatrixPanel.Destroy(ansPanel);
                 if (count > 256) {
                     MatrixParallel.ReturnParallelResources(ref pool!);
                 }
@@ -211,7 +233,7 @@ namespace StgSharp.Mathematics.Numeric
             {
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
+                NativeIntrinsic.Context.mat_f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
 
                 return;
             } else if (count <= 1024)
@@ -220,7 +242,7 @@ namespace StgSharp.Mathematics.Numeric
 
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
+                NativeIntrinsic.Context.mat_f32_buffer_scalar_mul(leftPtr, ansPtr, scalar, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
@@ -233,7 +255,7 @@ namespace StgSharp.Mathematics.Numeric
             package->SecCount = left.KernelRowLength;
             package->Scalar = scalar;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
-            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.f32_buffer_scalar_mul;
+            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_scalar_mul;
             MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);
@@ -261,7 +283,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
+                NativeIntrinsic.Context.mat_f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
 
                 return;
             } else if (count <= 1024)
@@ -271,7 +293,7 @@ namespace StgSharp.Mathematics.Numeric
                 MatrixKernel<float>* leftPtr = left.Buffer;
                 MatrixKernel<float>* rightPtr = right.Buffer;
                 MatrixKernel<float>* ansPtr = ans.Buffer;
-                NativeIntrinsic.Context.f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
+                NativeIntrinsic.Context.mat_f32_buffer_sub(leftPtr, rightPtr, ansPtr, count);
 
                 MatrixParallel.ReturnParallelResources(ref pool);
                 return;
@@ -284,7 +306,7 @@ namespace StgSharp.Mathematics.Numeric
             package->PrimCount = left.KernelColumnLength;
             package->SecCount = left.KernelRowLength;
             package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
-            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.f32_buffer_sub;
+            package->ComputeHandle = (IntPtr)NativeIntrinsic.Context.mat_f32_buffer_sub;
             MatrixParallelHandle taskGroup = MatrixParallel.ScheduleTask(package);
 
             MatrixParallel.LaunchParallel(taskGroup, SleepMode.DeepSleep);

StgSharp.Native/sn-head/sn_intrinsic.h

@@ -215,10 +215,18 @@ DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, add);
 DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, add);
 DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, add);
 
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, add);
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, add);
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, add);
+
 DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, sse, sub);
 DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, avx, sub);
 DECLARE_BUF_PROC_LEFT_RIGHT_ANS(float, 512, sub);
 
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, sub);
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, sub);
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 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);
@@ -227,10 +235,18 @@ DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, sse, scalar_mul);
 DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, avx, scalar_mul);
 DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(float, 512, scalar_mul);
 
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, sse, scalar_mul);
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, avx, scalar_mul);
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, 512, scalar_mul);
+
 DECLARE_BUF_PROC_ANS_SCALAR(float, sse, fill);
 DECLARE_BUF_PROC_ANS_SCALAR(float, avx, fill);
 DECLARE_BUF_PROC_ANS_SCALAR(float, 512, fill);
 
+DECLARE_BUF_PROC_ANS_SCALAR(double, sse, fill);
+DECLARE_BUF_PROC_ANS_SCALAR(double, avx, fill);
+DECLARE_BUF_PROC_ANS_SCALAR(double, 512, fill);
+
 #pragma endregion
 
 #pragma region fma
@@ -239,6 +255,10 @@ DECLARE_PNL_PROC_ANS(float, sse, clear);
 DECLARE_PNL_PROC_ANS(float, avx, clear);
 DECLARE_PNL_PROC_ANS(float, 512, clear);
 
+DECLARE_PNL_PROC_ANS(double, sse, clear);
+DECLARE_PNL_PROC_ANS(double, avx, clear);
+DECLARE_PNL_PROC_ANS(double, 512, clear);
+
 DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, sse, fma);
 DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx, fma);
 DECLARE_PNL_PROC_LEFT_RIGHT_ANS(float, avx_fma, fma);
@@ -278,20 +298,33 @@ INTERNAL int SN_DECL index_pair_sse(short const *str, uint32_t target, int lengt
 // Mirrors C# IntrinsicContext field order so managed/unmanaged layouts stay in sync.
 typedef struct sn_intrinsic {
         HASH city_hash_simplify;
-        BUF_PROC_LRA f32_buffer_add;
-        BUF_PROC_AS f32_buffer_fill;
-        BUF_PROC_RAS f32_buffer_scalar_mul;
-        BUF_PROC_LRA f32_buffer_sub;
-        KER_PROC_RA f32_buffer_transpose;
-        PANEL_PROC f32_build_panel;
-        PANEL_PROC f32_clear_panel;
-        VECTORNORMALIZEPROC f32_normalize_3;
-        PNL_PROC_LRA f32_panel_fma;
-        PIVOTPROC f32_pivot;
-        PANEL_PROC f32_store_panel;
         FACTORIALROC factorial_simd;
         INDEX_PAIR index_pair;
-        BUF_PROC_RAS variant;
+
+        BUF_PROC_LRA mat_f32_buffer_add;
+        BUF_PROC_AS mat_f32_buffer_fill;
+        BUF_PROC_RAS mat_f32_buffer_scalar_mul;
+        BUF_PROC_LRA mat_f32_buffer_sub;
+        KER_PROC_RA mat_f32_buffer_transpose;
+        PANEL_PROC mat_f32_build_panel;
+        PANEL_PROC mat_f32_clear_panel;
+        PNL_PROC_LRA mat_f32_panel_fma;
+        PIVOTPROC mat_f32_pivot;
+        PANEL_PROC mat_f32_store_panel;
+
+        // f64 slots (currently unimplemented; filled with NULL in loader)
+        BUF_PROC_LRA mat_f64_buffer_add;
+        BUF_PROC_AS mat_f64_buffer_fill;
+        BUF_PROC_RAS mat_f64_buffer_scalar_mul;
+        BUF_PROC_LRA mat_f64_buffer_sub;
+        KER_PROC_RA mat_f64_buffer_transpose;
+        PANEL_PROC mat_f64_build_panel;
+        PANEL_PROC mat_f64_clear_panel;
+        PNL_PROC_LRA mat_f64_panel_fma;
+        PIVOTPROC mat_f64_pivot;
+        PANEL_PROC mat_f64_store_panel;
+
+        VECTORNORMALIZEPROC vec_f32_normalize_3;
 } sn_intrinsic;
 
 typedef enum most_advanced_instruction {

StgSharp.Native/sn-head/sn_matkernel.h

@@ -12,20 +12,28 @@
 #define YMMCOUNT(T) (sizeof(T) / 2)
 #define ZMMCOUNT(T) (sizeof(T) / 4)
 
-#define MAT_KERNEL(T)                        \
+#define MAT_KERNEL(T)                       \
-        union {                              \
+        union /* MAT_KERNEL_##T## /**/ {    \
-                __m128 xmm[XMMCOUNT(T) / 4]; \
+                __m128 f32_x[XMMCOUNT(T)];  \
-                __m256 ymm[YMMCOUNT(T) / 4]; \
+                __m256 f32_y[YMMCOUNT(T)];  \
-                __m512 zmm[ZMMCOUNT(T) / 4]; \
+                __m512 f32_z[ZMMCOUNT(T)];  \
-                T m[4][4];                   \
+                                            \
+                __m128d f64_x[XMMCOUNT(T)]; \
+                __m256d f64_y[YMMCOUNT(T)]; \
+                __m512d f64_z[ZMMCOUNT(T)]; \
+                T m[4][4];                  \
         }
 
-#define VEC(T)                               \
+#define VEC(T)                                  \
-        union {                              \
+        union {                                 \
-                __m128 xmm[XMMCOUNT(T) / 4]; \
+                __m128 f32_x[XMMCOUNT(T) / 4];  \
-                __m256 ymm[YMMCOUNT(T) / 4]; \
+                __m256 f32_y[YMMCOUNT(T) / 4];  \
-                __m512 zmm[ZMMCOUNT(T) / 4]; \
+                __m512 f32_z[ZMMCOUNT(T) / 4];  \
-                T v[4];                      \
+                                                \
+                __m128d f64_x[XMMCOUNT(T) / 4]; \
+                __m256d f64_y[YMMCOUNT(T) / 4]; \
+                __m512d f64_z[ZMMCOUNT(T) / 4]; \
+                T v[4];                         \
         }
 
 /* Provide a named alias for the zero kernel type to ensure consistent linkage in C++. */
@@ -35,6 +43,8 @@ typedef MAT_KERNEL(uint64_t) sn_zero_kernel_u64;
 #define SN_LANES(T_VEC, T) (sizeof(T_VEC) / sizeof(T))
 #define SN_MAX(a, b) (((a)) > ((b)) ? (a) : (b))
 #define SN_CEIL_DIV(a, b) (((a) + (b) - 1) / (b))
+#define SN_BITS(x) ((int)(sizeof(x) * 8))
+#define PANEL_MINQ(T_VEC, T) SN_MAX(4, SN_MAX(SN_BITS(T_VEC) / 64, SN_BITS(T_VEC) / SN_BITS(T)))
 
 /* Panel geometry helpers (all compile-time): */
 #define PANEL_Q(T, T_VEC, MINQ) SN_MAX(SN_LANES(T_VEC, T), (MINQ))
@@ -60,26 +70,25 @@ typedef MAT_KERNEL(uint64_t) sn_zero_kernel_u64;
 
 /* ISA presets with minimum square side:
         - SSE  : MINQ=4   (avoid 2x2 for double)
-        - AVX2 : MINQ=8
+        - AVX2 : MINQ computed as max(4, vector_bits/64, vector_bits/element_bits)
-        - AVX-512: MINQ=8 (16 for float/i32 naturally, so this is a no-op there)
+        - AVX-512: same rule (float/i32 naturally expand to 16)
 */
 
 /* SSE panel layout: column-major buffer holding 4 logical columns with MINQ=4.
         Columns: each column packs up to four scalars so ALIGN equals sizeof(__m128).
         Alignment: every column boundary lands on a 16-byte multiple for XMM loads. */
-#define MAT_PANEL_sse(T) SN_ALIGNAS(16) MAT_PANEL_SQ_MIN_NOALIGN(T, __m128, 4)
+#define MAT_PANEL_sse(T) SN_ALIGNAS(16) MAT_PANEL_SQ_MIN_NOALIGN(T, __m128, PANEL_MINQ(__m128, T))
-/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ=8.
+/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ per PANEL_MINQ.
-        Columns: each column stores eight elements via two 256-bit chunks when T is 32-bit.
+        Columns: each column stores up to lanes elements via contiguous chunks.
         Alignment: column starts are 32-byte aligned so AVX2 loads/stores remain contiguous. */
-#define MAT_PANEL_avx(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, 8)
+#define MAT_PANEL_avx(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, PANEL_MINQ(__m256, T))
-/* AVX2 panel layout: column-major buffer holding 8 logical columns with MINQ=8.
+/* AVX2/FMA panel layout shares the same geometry as AVX2. */
-        Columns: each column stores eight elements via two 256-bit chunks when T is 32-bit.
+#define MAT_PANEL_avx_fma(T) \
-        Alignment: column starts are 32-byte aligned so AVX2 loads/stores remain contiguous. */
+        SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, PANEL_MINQ(__m256, T))
-#define MAT_PANEL_avx_fma(T) SN_ALIGNAS(32) MAT_PANEL_SQ_MIN_NOALIGN(T, __m256, 8)
 /* AVX-512 panel layout: column-major buffer holding 8 or more columns depending on T.
-        Columns: each column expands to ceil(Q/lanes) ZMM chunks; default MINQ=8 keeps doublse-wide tiles.
+        Columns: each column expands to ceil(Q/lanes) ZMM chunks.
         Alignment: column starts sit on 64-byte boundaries to match AVX-512 streaming access. */
-#define MAT_PANEL_512(T) SN_ALIGNAS(64) MAT_PANEL_SQ_MIN_NOALIGN(T, __m512, 8)
+#define MAT_PANEL_512(T) SN_ALIGNAS(64) MAT_PANEL_SQ_MIN_NOALIGN(T, __m512, PANEL_MINQ(__m512, T))
 
 /* Global zero kernel declaration (defined in dllmain.c), aligned to 64 bytes */
 #ifdef __cplusplus
@@ -117,11 +126,14 @@ extern __declspec(align(64)) sn_zero_kernel_u64 SN_ZERO_KERNEL;
                                                          int col_index, int row_index)
 
 /* Declarations for AVX float variant */
-DECLARE_BUILD_PANEL(sse, float);
 DECLARE_BUILD_PANEL(avx, float);
 DECLARE_BUILD_PANEL(512, float);
-DECLARE_STORE_PANEL(sse, float);
+
 DECLARE_STORE_PANEL(avx, float);
 DECLARE_STORE_PANEL(512, float);
 
+DECLARE_BUILD_PANEL(512, double);
+
+DECLARE_STORE_PANEL(512, double);
+
 #endif /* SN_MATKERNEL */

StgSharp.Native/sn-src/k_basicop.c

@@ -25,28 +25,53 @@ DECLARE_PNL_PROC_ANS(float, 512, clear)
         }
 }
 
+DECLARE_PNL_PROC_ANS(double, sse, clear)
+{
+        register __m128 zero = _mm_setzero_ps();
+        PANEL_COL_VEC((*ans), 0, 0) = zero;
+        PANEL_COL_VEC((*ans), 1, 0) = zero;
+        PANEL_COL_VEC((*ans), 2, 0) = zero;
+        PANEL_COL_VEC((*ans), 3, 0) = zero;
+}
+
+DECLARE_PNL_PROC_ANS(double, avx, clear)
+{
+        register __m256 zero = _mm256_setzero_ps();
+        for (size_t i = 0; i < 8; i++) {
+                PANEL_COL_VEC((*ans), i, 0) = zero;
+        }
+}
+
+DECLARE_PNL_PROC_ANS(double, 512, clear)
+{
+        register __m512 zero = _mm512_setzero_ps();
+        for (size_t i = 0; i < 16; i++) {
+                PANEL_COL_VEC((*ans), i, 0) = zero;
+        }
+}
+
 // Matrix transpose implementations
 DECLARE_KER_PROC_RIGHT_ANS(float, sse, transpose)
 {
-        register __m128 col0 = right->xmm[0];
+        register __m128 col0 = right->f32_x[0];
-        register __m128 col1 = right->xmm[1];
+        register __m128 col1 = right->f32_x[1];
-        register __m128 col2 = right->xmm[2];
+        register __m128 col2 = right->f32_x[2];
-        register __m128 col3 = right->xmm[3];
+        register __m128 col3 = right->f32_x[3];
         register __m128 tmp0 = _mm_unpacklo_ps(col0, col1);
         register __m128 tmp1 = _mm_unpackhi_ps(col0, col1);
         register __m128 tmp2 = _mm_unpacklo_ps(col2, col3);
         register __m128 tmp3 = _mm_unpackhi_ps(col2, col3);
 
-        ans->xmm[0] = _mm_movelh_ps(tmp0, tmp2);
+        ans->f32_x[0] = _mm_movelh_ps(tmp0, tmp2);
-        ans->xmm[1] = _mm_movehl_ps(tmp2, tmp0);
+        ans->f32_x[1] = _mm_movehl_ps(tmp2, tmp0);
-        ans->xmm[2] = _mm_movelh_ps(tmp1, tmp3);
+        ans->f32_x[2] = _mm_movelh_ps(tmp1, tmp3);
-        ans->xmm[3] = _mm_movehl_ps(tmp3, tmp1);
+        ans->f32_x[3] = _mm_movehl_ps(tmp3, tmp1);
 }
 
 DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose)
 {
-        register __m256 col_low = right->ymm[0];
+        register __m256 col_low = right->f32_y[0];
-        register __m256 col_high = right->ymm[1];
+        register __m256 col_high = right->f32_y[1];
 
         register __m256 unpack_low = _mm256_unpacklo_ps(col_low, col_high);
         register __m256 unpack_high = _mm256_unpackhi_ps(col_low, col_high);
@@ -54,16 +79,16 @@ DECLARE_KER_PROC_RIGHT_ANS(float, avx, transpose)
         col_low = _mm256_permute2f128_ps(unpack_low, unpack_high, 0x20);
         col_high = _mm256_permute2f128_ps(unpack_low, unpack_high, 0x31);
 
-        ans->ymm[0] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(2, 0, 2, 0));
+        ans->f32_y[0] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(2, 0, 2, 0));
-        ans->ymm[1] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(3, 1, 3, 1));
+        ans->f32_y[1] = _mm256_shuffle_ps(col_high, col_low, _MM_SHUFFLE(3, 1, 3, 1));
 }
 
 DECLARE_KER_PROC_RIGHT_ANS(float, 512, transpose)
 {
-        register __m512 matrix = right->zmm[0];
+        register __m512 matrix = right->f32_z[0];
 
         const __m512i transpose_indices =
                 _mm512_setr_epi32(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15);
 
-        ans->zmm[0] = _mm512_permutexvar_ps(transpose_indices, matrix);
+        ans->f32_z[0] = _mm512_permutexvar_ps(transpose_indices, matrix);
 }

StgSharp.Native/sn-src/k_bufop.c

@@ -4,10 +4,10 @@ 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].f32_x[0] = vec;
-                ans[i].xmm[1] = vec;
+                ans[i].f32_x[1] = vec;
-                ans[i].xmm[2] = vec;
+                ans[i].f32_x[2] = vec;
-                ans[i].xmm[3] = vec;
+                ans[i].f32_x[3] = vec;
         }
 }
 
@@ -15,8 +15,8 @@ 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].f32_y[0] = vec;
-                ans[i].ymm[1] = vec;
+                ans[i].f32_y[1] = vec;
         }
 }
 
@@ -24,22 +24,57 @@ 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;
+                ans[i].f32_z[0] = vec;
+        }
+}
+
+DECLARE_BUF_PROC_ANS_SCALAR(double, sse, fill)
+{
+        register __m128d vec = _mm_set1_pd(((double *)(scalar->data))[0]);
+        for (size_t i = 0; i < count; i++) {
+                ans[i].f64_x[0] = vec;
+                ans[i].f64_x[1] = vec;
+                ans[i].f64_x[2] = vec;
+                ans[i].f64_x[3] = vec;
+                ans[i].f64_x[4] = vec;
+                ans[i].f64_x[5] = vec;
+                ans[i].f64_x[6] = vec;
+                ans[i].f64_x[7] = vec;
+        }
+}
+
+DECLARE_BUF_PROC_ANS_SCALAR(double, avx, fill)
+{
+        register __m256d vec = _mm256_broadcast_sd(((double *)(scalar->data)));
+        for (size_t i = 0; i < count; i++) {
+                ans[i].f64_y[0] = vec;
+                ans[i].f64_y[1] = vec;
+                ans[i].f64_y[2] = vec;
+                ans[i].f64_y[3] = vec;
+        }
+}
+
+DECLARE_BUF_PROC_ANS_SCALAR(double, 512, fill)
+{
+        register __m512d vec = _mm512_broadcastsd_pd(_mm_loadu_pd(scalar));
+        for (size_t i = 0; i < count; i++) {
+                ans[i].f64_z[0] = vec;
+                ans[i].f64_z[1] = 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 c0 = _mm_add_ps(left[i].f32_x[0], right[i].f32_x[0]);
-                register __m128 c1 = _mm_add_ps(left[i].xmm[1], right[i].xmm[1]);
+                register __m128 c1 = _mm_add_ps(left[i].f32_x[1], right[i].f32_x[1]);
-                register __m128 c2 = _mm_add_ps(left[i].xmm[2], right[i].xmm[2]);
+                register __m128 c2 = _mm_add_ps(left[i].f32_x[2], right[i].f32_x[2]);
-                register __m128 c3 = _mm_add_ps(left[i].xmm[3], right[i].xmm[3]);
+                register __m128 c3 = _mm_add_ps(left[i].f32_x[3], right[i].f32_x[3]);
 
-                ans[i].xmm[0] = c0;
+                ans[i].f32_x[0] = c0;
-                ans[i].xmm[1] = c1;
+                ans[i].f32_x[1] = c1;
-                ans[i].xmm[2] = c2;
+                ans[i].f32_x[2] = c2;
-                ans[i].xmm[3] = c3;
+                ans[i].f32_x[3] = c3;
         }
 }
 
@@ -47,51 +82,147 @@ 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 c0 = _mm256_add_ps(left[i].f32_y[0], right[i].f32_y[0]);
-                register __m256 c1 = _mm256_add_ps(left[i].ymm[1], right[i].ymm[1]);
+                register __m256 c1 = _mm256_add_ps(left[i].f32_y[1], right[i].f32_y[1]);
-                ans[i].ymm[0] = c0;
+                ans[i].f32_y[0] = c0;
-                ans[i].ymm[1] = c1;
+                ans[i].f32_y[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]);
+                register __m512 c0 = _mm512_add_ps(left[i].f32_z[0], right[i].f32_z[0]);
-                ans[i].zmm[0] = c0;
+                ans[i].f32_z[0] = c0;
+        }
+}
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, add)
+{
+        for (size_t i = 0; i < count; i++) {
+                register __m128d c0 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[0]),
+                                                 _mm_castps_pd(right[i].f32_x[0]));
+                register __m128d c1 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[1]),
+                                                 _mm_castps_pd(right[i].f32_x[1]));
+                register __m128d c2 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[2]),
+                                                 _mm_castps_pd(right[i].f32_x[2]));
+                register __m128d c3 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[3]),
+                                                 _mm_castps_pd(right[i].f32_x[3]));
+                register __m128d c4 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[4]),
+                                                 _mm_castps_pd(right[i].f32_x[4]));
+                register __m128d c5 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[5]),
+                                                 _mm_castps_pd(right[i].f32_x[5]));
+                register __m128d c6 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[6]),
+                                                 _mm_castps_pd(right[i].f32_x[6]));
+                register __m128d c7 = _mm_add_pd(_mm_castps_pd(left[i].f32_x[7]),
+                                                 _mm_castps_pd(right[i].f32_x[7]));
+
+                ans[i].f32_x[0] = _mm_castpd_ps(c0);
+                ans[i].f32_x[1] = _mm_castpd_ps(c1);
+                ans[i].f32_x[2] = _mm_castpd_ps(c2);
+                ans[i].f32_x[3] = _mm_castpd_ps(c3);
+                ans[i].f32_x[4] = _mm_castpd_ps(c4);
+                ans[i].f32_x[5] = _mm_castpd_ps(c5);
+                ans[i].f32_x[6] = _mm_castpd_ps(c6);
+                ans[i].f32_x[7] = _mm_castpd_ps(c7);
+        }
+}
+
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, add)
+{
+        for (size_t i = 0; i < count; i++) {
+                size_t s = sizeof(MAT_KERNEL(float));
+                register __m256d c0 = _mm256_add_pd(left[i].f64_y[0], right[i].f64_y[0]);
+                register __m256d c1 = _mm256_add_pd(left[i].f64_y[1], right[i].f64_y[1]);
+                register __m256d c2 = _mm256_add_pd(left[i].f64_y[0], right[i].f64_y[0]);
+                register __m256d c3 = _mm256_add_pd(left[i].f64_y[1], right[i].f64_y[1]);
+                ans[i].f64_y[0] = c0;
+                ans[i].f64_y[1] = c1;
+                ans[i].f64_y[2] = c2;
+                ans[i].f64_y[3] = c3;
+        }
+}
+
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, add)
+{
+        for (size_t i = 0; i < count; i++) {
+                register __m512d c0 = _mm512_add_pd(left[i].f64_z[0], right[i].f64_z[0]);
+                register __m512d c1 = _mm512_add_pd(left[i].f64_z[1], right[i].f64_z[1]);
+                ans[i].f64_z[0] = c0;
+                ans[i].f64_z[1] = c1;
         }
 }
 
 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 c0 = _mm_sub_ps(left[i].f32_x[0], right[i].f32_x[0]);
-                register __m128 c1 = _mm_sub_ps(left[i].xmm[1], right[i].xmm[1]);
+                register __m128 c1 = _mm_sub_ps(left[i].f32_x[1], right[i].f32_x[1]);
-                register __m128 c2 = _mm_sub_ps(left[i].xmm[2], right[i].xmm[2]);
+                register __m128 c2 = _mm_sub_ps(left[i].f32_x[2], right[i].f32_x[2]);
-                register __m128 c3 = _mm_sub_ps(left[i].xmm[3], right[i].xmm[3]);
+                register __m128 c3 = _mm_sub_ps(left[i].f32_x[3], right[i].f32_x[3]);
 
-                ans[i].xmm[0] = c0;
+                ans[i].f32_x[0] = c0;
-                ans[i].xmm[1] = c1;
+                ans[i].f32_x[1] = c1;
-                ans[i].xmm[2] = c2;
+                ans[i].f32_x[2] = c2;
-                ans[i].xmm[3] = c3;
+                ans[i].f32_x[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 c0 = _mm256_sub_ps(left[i].f32_y[0], right[i].f32_y[0]);
-                register __m256 c1 = _mm256_sub_ps(left[i].ymm[1], right[i].ymm[1]);
+                register __m256 c1 = _mm256_sub_ps(left[i].f32_y[1], right[i].f32_y[1]);
-                ans[i].ymm[0] = c0;
+                ans[i].f32_y[0] = c0;
-                ans[i].ymm[1] = c1;
+                ans[i].f32_y[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]);
+                register __m512 c0 = _mm512_sub_ps(left[i].f32_z[0], right[i].f32_z[0]);
-                ans[i].zmm[0] = c0;
+                ans[i].f32_z[0] = c0;
+        }
+}
+
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, sse, sub)
+{
+        for (size_t i = 0; i < count; i++) {
+                register __m128d c0 = _mm_sub_pd(left[i].f64_x[0], right[i].f64_x[0]);
+                register __m128d c1 = _mm_sub_pd(left[i].f64_x[1], right[i].f64_x[1]);
+                register __m128d c2 = _mm_sub_pd(left[i].f64_x[2], right[i].f64_x[2]);
+                register __m128d c3 = _mm_sub_pd(left[i].f64_x[3], right[i].f64_x[3]);
+                register __m128d c4 = _mm_sub_pd(left[i].f64_x[4], right[i].f64_x[4]);
+                register __m128d c5 = _mm_sub_pd(left[i].f64_x[5], right[i].f64_x[5]);
+                register __m128d c6 = _mm_sub_pd(left[i].f64_x[6], right[i].f64_x[6]);
+                register __m128d c7 = _mm_sub_pd(left[i].f64_x[7], right[i].f64_x[7]);
+
+                ans[i].f64_x[0] = c0;
+                ans[i].f64_x[1] = c1;
+                ans[i].f64_x[2] = c2;
+                ans[i].f64_x[3] = c3;
+                ans[i].f64_x[4] = c4;
+                ans[i].f64_x[5] = c5;
+                ans[i].f64_x[6] = c6;
+                ans[i].f64_x[7] = c7;
+        }
+}
+
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, avx, sub)
+{
+        for (size_t i = 0; i < count; i++) {
+                register __m256 c0 = _mm256_sub_ps(left[i].f32_y[0], right[i].f32_y[0]);
+                register __m256 c1 = _mm256_sub_ps(left[i].f32_y[1], right[i].f32_y[1]);
+                ans[i].f32_y[0] = c0;
+                ans[i].f32_y[1] = c1;
+        }
+}
+
+DECLARE_BUF_PROC_LEFT_RIGHT_ANS(double, 512, sub)
+{
+        for (size_t i = 0; i < count; i++) {
+                register __m512 c0 = _mm512_sub_ps(left[i].f32_z[0], right[i].f32_z[0]);
+                ans[i].f32_z[0] = c0;
         }
 }
 
@@ -99,15 +230,15 @@ 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 c0 = _mm_mul_ps(right[i].f32_x[0], vec);
-                register __m128 c1 = _mm_mul_ps(right[i].xmm[1], vec);
+                register __m128 c1 = _mm_mul_ps(right[i].f32_x[1], vec);
-                register __m128 c2 = _mm_mul_ps(right[i].xmm[2], vec);
+                register __m128 c2 = _mm_mul_ps(right[i].f32_x[2], vec);
-                register __m128 c3 = _mm_mul_ps(right[i].xmm[3], vec);
+                register __m128 c3 = _mm_mul_ps(right[i].f32_x[3], vec);
 
-                ans[i].xmm[0] = c0;
+                ans[i].f32_x[0] = c0;
-                ans[i].xmm[1] = c1;
+                ans[i].f32_x[1] = c1;
-                ans[i].xmm[2] = c2;
+                ans[i].f32_x[2] = c2;
-                ans[i].xmm[3] = c3;
+                ans[i].f32_x[3] = c3;
         }
 }
 
@@ -115,11 +246,11 @@ 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 c0 = _mm256_mul_ps(right[i].f32_y[0], vec);
-                register __m256 c1 = _mm256_mul_ps(right[i].ymm[1], vec);
+                register __m256 c1 = _mm256_mul_ps(right[i].f32_y[1], vec);
 
-                ans[i].ymm[0] = c0;
+                ans[i].f32_y[0] = c0;
-                ans[i].ymm[1] = c1;
+                ans[i].f32_y[1] = c1;
         }
 }
 
@@ -127,8 +258,60 @@ 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);
+                register __m512 c0 = _mm512_mul_ps(right[i].f32_z[0], vec);
 
-                ans[i].zmm[0] = c0;
+                ans[i].f32_z[0] = c0;
+        }
+}
+
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, sse, scalar_mul)
+{
+        register __m128d vec = _mm_set1_pd(((double *)(scalar->data))[0]);
+        for (size_t i = 0; i < count; i++) {
+                register __m128d c0 = _mm_mul_pd(right[i].f64_x[0], vec);
+                register __m128d c1 = _mm_mul_pd(right[i].f64_x[1], vec);
+                register __m128d c2 = _mm_mul_pd(right[i].f64_x[2], vec);
+                register __m128d c3 = _mm_mul_pd(right[i].f64_x[3], vec);
+                register __m128d c4 = _mm_mul_pd(right[i].f64_x[4], vec);
+                register __m128d c5 = _mm_mul_pd(right[i].f64_x[5], vec);
+                register __m128d c6 = _mm_mul_pd(right[i].f64_x[6], vec);
+                register __m128d c7 = _mm_mul_pd(right[i].f64_x[7], vec);
+
+                ans[i].f64_x[0] = c0;
+                ans[i].f64_x[1] = c1;
+                ans[i].f64_x[2] = c2;
+                ans[i].f64_x[3] = c3;
+                ans[i].f64_x[4] = c4;
+                ans[i].f64_x[5] = c5;
+                ans[i].f64_x[6] = c6;
+                ans[i].f64_x[7] = c7;
+        }
+}
+
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, avx, scalar_mul)
+{
+        register __m256d vec = _mm256_broadcast_sd(((double *)(scalar->data)));
+        for (size_t i = 0; i < count; i++) {
+                register __m256d c0 = _mm256_mul_pd(right[i].f64_y[0], vec);
+                register __m256d c1 = _mm256_mul_pd(right[i].f64_y[1], vec);
+                register __m256d c2 = _mm256_mul_pd(right[i].f64_y[2], vec);
+                register __m256d c3 = _mm256_mul_pd(right[i].f64_y[3], vec);
+
+                ans[i].f64_y[0] = c0;
+                ans[i].f64_y[1] = c1;
+                ans[i].f64_y[2] = c2;
+                ans[i].f64_y[3] = c3;
+        }
+}
+
+DECLARE_BUF_PROC_RIGHT_ANS_SCALAR(double, 512, scalar_mul)
+{
+        register __m512d vec = _mm512_broadcastsd_pd(_mm_load_pd(scalar));
+        for (size_t i = 0; i < count; i++) {
+                register __m512d c0 = _mm512_mul_pd(right[i].f64_z[0], vec);
+                register __m512d c1 = _mm512_mul_pd(right[i].f64_z[1], vec);
+
+                ans[i].f64_z[0] = c0;
+                ans[i].f64_z[1] = c1;
         }
 }

StgSharp.Native/sn-src/k_panel.c

@@ -31,8 +31,8 @@ DECLARE_BUILD_PANEL(avx, float)
                 const MAT_KERNEL(float) *src_bot = is_right_half ? k_rb : k_lb;
 
                 /* Vector loads: each kernel column is a compact __m128 */
-                top128 = src_top->xmm[kc_col];
+                top128 = src_top->f32_x[kc_col];
-                bot128 = src_bot->xmm[kc_col];
+                bot128 = src_bot->f32_x[kc_col];
 
                 /* Combine into an AVX column: low 128 = top, high 128 = bottom */
                 col256 = _mm256_set_m128(bot128, top128);
@@ -41,15 +41,6 @@ DECLARE_BUILD_PANEL(avx, float)
         }
 }
 
-DECLARE_BUILD_PANEL(sse, float)
-{
-        /* Single kernel -> 4 columns, each is one __m128 column vector. */
-        const MAT_KERNEL(float) *src = matrix + (row_index * col_length + col_index);
-
-        for (int pc = 0; pc < 4; ++pc) {
-                PANEL_COL_VEC((*panel), pc, 0) = src->xmm[pc];
-        }
-}
 
 DECLARE_BUILD_PANEL(512, float)
 {
@@ -71,10 +62,10 @@ DECLARE_BUILD_PANEL(512, float)
                 const int blk_col = pc >> 2; /* which kernel column block (0..3) */
                 const int kc_col = pc & 3; /* column inside a 4x4 kernel */
 
-                const __m128 r0 = k[blk_col][0]->xmm[kc_col];
+                const __m128 r0 = k[blk_col][0]->f32_x[kc_col];
-                const __m128 r1 = k[blk_col][1]->xmm[kc_col];
+                const __m128 r1 = k[blk_col][1]->f32_x[kc_col];
-                const __m128 r2 = k[blk_col][2]->xmm[kc_col];
+                const __m128 r2 = k[blk_col][2]->f32_x[kc_col];
-                const __m128 r3 = k[blk_col][3]->xmm[kc_col];
+                const __m128 r3 = k[blk_col][3]->f32_x[kc_col];
 
                 const __m256 top = _mm256_set_m128(r1, r0);
                 const __m256 bot = _mm256_set_m128(r3, r2);
@@ -86,6 +77,43 @@ DECLARE_BUILD_PANEL(512, float)
         }
 }
 
+
+DECLARE_BUILD_PANEL(512, double)
+{
+        const int Q = 8; /* max(4, lanes) for double with ZMM -> 8 columns */
+
+        const int has_right = (col_index + 1) < col_length;
+        const int has_bottom = (row_index + 1) < row_length;
+
+        const MAT_KERNEL(double) *const k_lt = matrix + (row_index * col_length + col_index);
+        const MAT_KERNEL(double) *const k_rt =
+                has_right ? (matrix + (row_index * col_length + (col_index + 1))) :
+                            (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL;
+        const MAT_KERNEL(double) *const k_lb =
+                has_bottom ? (matrix + ((row_index + 1) * col_length + col_index)) :
+                             (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL;
+        const MAT_KERNEL(double) *const k_rb =
+                (has_right && has_bottom) ?
+                        (matrix + ((row_index + 1) * col_length + (col_index + 1))) :
+                        (const MAT_KERNEL(double) *)&SN_ZERO_KERNEL;
+
+        for (int pc = 0; pc < Q; ++pc) {
+                const int is_right_half = (pc >= 4);
+                const int kc_col = pc & 3; /* column index inside 4x4 kernel */
+
+                const MAT_KERNEL(double) *src_top = is_right_half ? k_rt : k_lt;
+                const MAT_KERNEL(double) *src_bot = is_right_half ? k_rb : k_lb;
+
+                const __m256d top256 = src_top->f64_y[kc_col];
+                const __m256d bot256 = src_bot->f64_y[kc_col];
+
+                __m512d col = _mm512_castpd256_pd512(top256);
+                col = _mm512_insertf64x4(col, bot256, 1);
+
+                PANEL_COL_VEC((*panel), pc, 0) = _mm512_castpd_ps(col);
+        }
+}
+
 DECLARE_STORE_PANEL(avx, float)
 {
         const int Q = 8; /* fixed for AVX float panel */
@@ -124,8 +152,8 @@ DECLARE_STORE_PANEL(avx, float)
                 bot128 = _mm256_extractf128_ps(col256, 1); /* high 128 */
 
                 /* Store back into kernels directly (out-of-bounds writes go to dummy) */
-                k_rt->xmm[kc_col] = top128;
+                k_rt->f32_x[kc_col] = top128;
-                k_rb->xmm[kc_col] = bot128;
+                k_rb->f32_x[kc_col] = bot128;
         }
         for (; i < Q; i++) {
                 const int kc_col = i - 4; /* column index inside 4x4 kernel */
@@ -137,20 +165,11 @@ DECLARE_STORE_PANEL(avx, float)
                 bot128 = _mm256_extractf128_ps(col256, 1); /* high 128 */
 
                 /* Store back into kernels directly (out-of-bounds writes go to dummy) */
-                k_lt->xmm[kc_col] = top128;
+                k_lt->f32_x[kc_col] = top128;
-                k_lb->xmm[kc_col] = bot128;
+                k_lb->f32_x[kc_col] = bot128;
         }
 }
 
-DECLARE_STORE_PANEL(sse, float)
-{
-        MAT_KERNEL(float) *dst = matrix + (row_index * col_length + col_index);
-
-        for (int pc = 0; pc < 4; ++pc) {
-                const __m128 col = PANEL_COL_VEC((*panel), pc, 0);
-                dst->xmm[pc] = col;
-        }
-}
 
 DECLARE_STORE_PANEL(512, float)
 {
@@ -181,9 +200,99 @@ DECLARE_STORE_PANEL(512, float)
                 const __m128 r2 = _mm256_castps256_ps128(bot);
                 const __m128 r3 = _mm256_extractf128_ps(bot, 1);
 
-                k[blk_col][0]->xmm[kc_col] = r0;
+                k[blk_col][0]->f32_x[kc_col] = r0;
-                k[blk_col][1]->xmm[kc_col] = r1;
+                k[blk_col][1]->f32_x[kc_col] = r1;
-                k[blk_col][2]->xmm[kc_col] = r2;
+                k[blk_col][2]->f32_x[kc_col] = r2;
-                k[blk_col][3]->xmm[kc_col] = r3;
+                k[blk_col][3]->f32_x[kc_col] = r3;
+        }
+}
+
+DECLARE_STORE_PANEL(sse, double)
+{
+        MAT_KERNEL(double) *dst = matrix + (row_index * col_length + col_index);
+
+        for (int pc = 0; pc < 4; ++pc) {
+                const int base = pc * 2;
+                dst->f64_x[base] = _mm_castps_pd(PANEL_COL_VEC((*panel), pc, 0));
+                dst->f64_x[base + 1] = _mm_castps_pd(PANEL_COL_VEC((*panel), pc, 1));
+        }
+}
+
+DECLARE_STORE_PANEL(avx, double)
+{
+        const int Q = 8; /* fixed for AVX double panel */
+
+        MAT_KERNEL(double) zero;
+
+        const int has_right = (col_index + 1) < col_length;
+        const int has_bottom = (row_index + 1) < row_length;
+
+        MAT_KERNEL(double) *const k_lt = (col_index < col_length && row_index < row_length) ?
+                                                 (matrix + (row_index * col_length + col_index)) :
+                                                 &zero;
+        MAT_KERNEL(double) *const k_rt =
+                has_right && (row_index < row_length) ?
+                        (matrix + (row_index * col_length + (col_index + 1))) :
+                        &zero;
+        MAT_KERNEL(double) *const k_lb =
+                has_bottom && (col_index < col_length) ?
+                        (matrix + ((row_index + 1) * col_length + col_index)) :
+                        &zero;
+        MAT_KERNEL(double) *const k_rb =
+                (has_right && has_bottom) ?
+                        (matrix + ((row_index + 1) * col_length + (col_index + 1))) :
+                        &zero;
+
+        for (int i = 0; i < Q; ++i) {
+                const int kc_col = (i < 4) ? i : (i - 4);
+                const __m256 top_chunk = PANEL_COL_VEC((*panel), i, 0);
+                const __m256 bot_chunk = PANEL_COL_VEC((*panel), i, 1);
+
+                MAT_KERNEL(double) *const top_dst = (i < 4) ? k_rt : k_lt;
+                MAT_KERNEL(double) *const bot_dst = (i < 4) ? k_rb : k_lb;
+
+                top_dst->f64_y[kc_col] = _mm256_castps_pd(top_chunk);
+                bot_dst->f64_y[kc_col] = _mm256_castps_pd(bot_chunk);
+        }
+}
+
+DECLARE_STORE_PANEL(512, double)
+{
+        const int Q = 8; /* max(4, lanes) -> 8 columns, one ZMM per column */
+
+        MAT_KERNEL(double) zero;
+
+        const int has_right = (col_index + 1) < col_length;
+        const int has_bottom = (row_index + 1) < row_length;
+
+        MAT_KERNEL(double) *const k_lt = (col_index < col_length && row_index < row_length) ?
+                                                 (matrix + (row_index * col_length + col_index)) :
+                                                 &zero;
+        MAT_KERNEL(double) *const k_rt =
+                has_right && (row_index < row_length) ?
+                        (matrix + (row_index * col_length + (col_index + 1))) :
+                        &zero;
+        MAT_KERNEL(double) *const k_lb =
+                has_bottom && (col_index < col_length) ?
+                        (matrix + ((row_index + 1) * col_length + col_index)) :
+                        &zero;
+        MAT_KERNEL(double) *const k_rb =
+                (has_right && has_bottom) ?
+                        (matrix + ((row_index + 1) * col_length + (col_index + 1))) :
+                        &zero;
+
+        for (int i = 0; i < Q; ++i) {
+                const int kc_col = (i < 4) ? i : (i - 4);
+                const __m512 col = PANEL_COL_VEC((*panel), i, 0);
+                const __m512d col_pd = _mm512_castps_pd(col);
+
+                const __m256d top256 = _mm512_castpd512_pd256(col_pd);
+                const __m256d bot256 = _mm512_extractf64x4_pd(col_pd, 1);
+
+                MAT_KERNEL(double) *const top_dst = (i < 4) ? k_rt : k_lt;
+                MAT_KERNEL(double) *const bot_dst = (i < 4) ? k_rb : k_lb;
+
+                top_dst->f64_y[kc_col] = top256;
+                bot_dst->f64_y[kc_col] = bot256;
         }
 }

StgSharp.Native/sn-src/k_pivot.c

@@ -2,8 +2,8 @@
 
 // Hybrid implementation: SIMD for Abs, Scalar for Max
 // Assumes Column-Major layout for both MAT_KERNEL content and Kernel sequence
-static void pivot_float(MAT_KERNEL(float) const *source, size_t col_begin, size_t col_length,
+INTERNAL void SN_DECL pivot_float(MAT_KERNEL(float) const *source, size_t col_begin,
-                        size_t *restrict ans)
+                                  size_t col_length, size_t *restrict ans)
 {
         float max_v[4] = { -1.0f, -1.0f, -1.0f, -1.0f };
 

StgSharp.Native/sn-src/load_intrinsic.c

@@ -65,17 +65,30 @@ void load_intrinsic_sse(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
-        context->f32_store_panel = STORE_PANEL(sse, float);
+
-        context->f32_clear_panel = PNL_PROC(float, sse, clear);
+        context->mat_f32_buffer_add = BUF_PROC(float, sse, add);
-        context->f32_build_panel = BUILD_PANEL(sse, float);
+        context->mat_f32_buffer_fill = BUF_PROC(float, sse, fill);
-        context->f32_buffer_add = BUF_PROC(float, sse, add);
+        context->mat_f32_buffer_scalar_mul = BUF_PROC(float, sse, scalar_mul);
-        context->f32_panel_fma = PNL_PROC(float, sse, fma);
+        context->mat_f32_buffer_sub = BUF_PROC(float, sse, sub);
-        context->f32_buffer_fill = BUF_PROC(float, sse, fill);
+        context->mat_f32_buffer_transpose = KER_PROC(float, sse, transpose);
-        context->variant = BUF_PROC(float, sse, scalar_mul);
+        context->mat_f32_build_panel = NULL;
-        context->f32_buffer_sub = BUF_PROC(float, sse, sub);
+        context->mat_f32_clear_panel = PNL_PROC(float, sse, clear);
-        context->f32_buffer_transpose = KER_PROC(float, sse, transpose);
+        context->mat_f32_panel_fma = PNL_PROC(float, sse, fma);
-        context->f32_buffer_scalar_mul = BUF_PROC(float, sse, scalar_mul);
+        context->mat_f32_pivot = pivot_float;
-        context->f32_normalize_3 = f32_normalize;
+        context->mat_f32_store_panel = NULL;
+
+        context->mat_f64_buffer_add = NULL;
+        context->mat_f64_buffer_fill = NULL;
+        context->mat_f64_buffer_scalar_mul = NULL;
+        context->mat_f64_buffer_sub = NULL;
+        context->mat_f64_buffer_transpose = NULL;
+        context->mat_f64_build_panel = NULL;
+        context->mat_f64_clear_panel = NULL;
+        context->mat_f64_panel_fma = NULL;
+        context->mat_f64_pivot = NULL;
+        context->mat_f64_store_panel = NULL;
+
+        context->vec_f32_normalize_3 = f32_normalize;
 }
 
 void load_intrinsic_avx(sn_intrinsic *context)
@@ -83,17 +96,30 @@ void load_intrinsic_avx(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
-        context->f32_store_panel = STORE_PANEL(avx, float);
+
-        context->f32_clear_panel = PNL_PROC(float, avx, clear);
+        context->mat_f32_buffer_add = BUF_PROC(float, avx, add);
-        context->f32_build_panel = BUILD_PANEL(avx, float);
+        context->mat_f32_buffer_fill = BUF_PROC(float, avx, fill);
-        context->f32_buffer_add = BUF_PROC(float, avx, add);
+        context->mat_f32_buffer_scalar_mul = BUF_PROC(float, avx, scalar_mul);
-        context->f32_panel_fma = PNL_PROC(float, avx, fma);
+        context->mat_f32_buffer_sub = BUF_PROC(float, avx, sub);
-        context->f32_buffer_fill = BUF_PROC(float, avx, fill);
+        context->mat_f32_buffer_transpose = KER_PROC(float, avx, transpose);
-        context->variant = BUF_PROC(float, avx, scalar_mul);
+        context->mat_f32_build_panel = BUILD_PANEL(avx, float);
-        context->f32_buffer_sub = BUF_PROC(float, avx, sub);
+        context->mat_f32_clear_panel = PNL_PROC(float, avx, clear);
-        context->f32_buffer_transpose = KER_PROC(float, avx, transpose);
+        context->mat_f32_panel_fma = PNL_PROC(float, avx, fma);
-        context->f32_buffer_scalar_mul = BUF_PROC(float, avx, scalar_mul);
+        context->mat_f32_pivot = pivot_float;
-        context->f32_normalize_3 = f32_normalize;
+        context->mat_f32_store_panel = STORE_PANEL(avx, float);
+
+        context->mat_f64_buffer_add = BUF_PROC(double, avx, add);
+        context->mat_f64_buffer_fill = BUF_PROC(double, avx, fill);
+        context->mat_f64_buffer_scalar_mul = BUF_PROC(double, avx, scalar_mul);
+        context->mat_f64_buffer_sub = BUF_PROC(double, avx, sub);
+        context->mat_f64_buffer_transpose = NULL;
+        context->mat_f64_build_panel = NULL;
+        context->mat_f64_clear_panel = PNL_PROC(double, avx, clear);
+        context->mat_f64_panel_fma = NULL;
+        context->mat_f64_pivot = NULL;
+        context->mat_f64_store_panel = NULL;
+
+        context->vec_f32_normalize_3 = f32_normalize;
 }
 
 void load_intrinsic_512(sn_intrinsic *context)
@@ -101,15 +127,28 @@ void load_intrinsic_512(sn_intrinsic *context)
         context->city_hash_simplify = city_hash_simplify_sse;
         context->factorial_simd = factorial_simd_sse;
         context->index_pair = index_pair_sse;
-        context->f32_store_panel = STORE_PANEL(512, float);
+
-        context->f32_clear_panel = PNL_PROC(float, 512, clear);
+        context->mat_f32_buffer_add = BUF_PROC(float, 512, add);
-        context->f32_build_panel = BUILD_PANEL(512, float);
+        context->mat_f32_buffer_fill = BUF_PROC(float, 512, fill);
-        context->f32_buffer_add = BUF_PROC(float, 512, add);
+        context->mat_f32_buffer_scalar_mul = BUF_PROC(float, 512, scalar_mul);
-        context->f32_panel_fma = PNL_PROC(float, 512, fma);
+        context->mat_f32_buffer_sub = BUF_PROC(float, 512, sub);
-        context->f32_buffer_fill = BUF_PROC(float, 512, fill);
+        context->mat_f32_buffer_transpose = KER_PROC(float, 512, transpose);
-        context->variant = BUF_PROC(float, 512, scalar_mul);
+        context->mat_f32_build_panel = BUILD_PANEL(512, float);
-        context->f32_buffer_sub = BUF_PROC(float, 512, sub);
+        context->mat_f32_clear_panel = PNL_PROC(float, 512, clear);
-        context->f32_buffer_transpose = KER_PROC(float, 512, transpose);
+        context->mat_f32_panel_fma = PNL_PROC(float, 512, fma);
-        context->f32_buffer_scalar_mul = BUF_PROC(float, 512, scalar_mul);
+        context->mat_f32_pivot = pivot_float;
-        context->f32_normalize_3 = f32_normalize;
+        context->mat_f32_store_panel = STORE_PANEL(512, float);
+
+        context->mat_f64_buffer_add = NULL;
+        context->mat_f64_buffer_fill = NULL;
+        context->mat_f64_buffer_scalar_mul = NULL;
+        context->mat_f64_buffer_sub = NULL;
+        context->mat_f64_buffer_transpose = NULL;
+        context->mat_f64_build_panel = NULL;
+        context->mat_f64_clear_panel = NULL;
+        context->mat_f64_panel_fma = NULL;
+        context->mat_f64_pivot = NULL;
+        context->mat_f64_store_panel = NULL;
+
+        context->vec_f32_normalize_3 = f32_normalize;
 }

StgSharp.Native/sn-src/vector.c

@@ -3,7 +3,7 @@
 
 INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *restrict target)
 {
-        register __m128 s = _mm_load_ps(&source->xmm);
+        register __m128 s = _mm_load_ps(&source->f32_x);
         register __m128 a_length_sq = _mm_mul_ps(s, s);
         a_length_sq = _mm_add_ps(a_length_sq,
                                  _mm_shuffle_ps(a_length_sq, a_length_sq, _MM_SHUFFLE(2, 3, 0, 1)));
@@ -13,13 +13,13 @@ INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *res
         __m128 a_inv_length = _mm_rsqrt_ps(a_length_sq);
         s = _mm_mul_ps(s, a_inv_length);
 
-        __m128 dot_product = _mm_mul_ps(target->xmm[0], s);
+        __m128 dot_product = _mm_mul_ps(target->f32_x[0], s);
         dot_product = _mm_add_ps(dot_product,
                                  _mm_shuffle_ps(dot_product, dot_product, _MM_SHUFFLE(2, 3, 0, 1)));
         dot_product = _mm_add_ss(dot_product,
                                  _mm_shuffle_ps(dot_product, dot_product, _MM_SHUFFLE(0, 0, 3, 2)));
 
-        __m128 b_orthogonal = _mm_sub_ps(target->xmm[0], _mm_mul_ps(s, dot_product));
+        __m128 b_orthogonal = _mm_sub_ps(target->f32_x[0], _mm_mul_ps(s, dot_product));
 
         __m128 b_length_sq = _mm_mul_ps(b_orthogonal, b_orthogonal);
         b_length_sq = _mm_add_ps(b_length_sq,
@@ -28,5 +28,5 @@ INTERNAL void SN_DECL f32_normalize(VEC(float) *restrict source, VEC(float) *res
                                  _mm_shuffle_ps(b_length_sq, b_length_sq, _MM_SHUFFLE(0, 0, 3, 2)));
 
         __m128 b_inv_length = _mm_rsqrt_ps(b_length_sq);
-        target->xmm[0] = _mm_mul_ps(b_orthogonal, b_inv_length);
+        target->f32_x[0] = _mm_mul_ps(b_orthogonal, b_inv_length);
 }

StgSharp.Native/sn-src/vector_dot.c

@@ -6,12 +6,12 @@ static const int avxdot_indices[] = { 0, 4, 1, 5, 2, 3, 6, 7 };
 INTERNAL void FORCEINLINE SN_DECL dot_41_sse(MAT_KERNEL(float) const *transpose,
                                              __m128 const *vector, __m128 *ans)
 {
-        *ans = _mm_dp_ps(transpose->xmm[0], *vector, 0B11110001);
+        *ans = _mm_dp_ps(transpose->f32_x[0], *vector, 0B11110001);
-        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], *vector, 0B11110001),
+        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], *vector, 0B11110001),
                              SSE_INSERT_SIMPLE(1, 1));
-        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], *vector, 0B11110001),
+        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], *vector, 0B11110001),
                              SSE_INSERT_SIMPLE(1, 3));
-        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], *vector, 0B11110001),
+        *ans = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], *vector, 0B11110001),
                              SSE_INSERT_SIMPLE(1, 4));
 }
 
@@ -19,8 +19,8 @@ INTERNAL void SN_DECL dot_41_avx(MAT_KERNEL(float) const *transpose, __m128 cons
                                  __m128 *ans)
 {
         __m256 _double = _mm256_set_m128(vector[0], vector[0]);
-        __m256 tmp = _mm256_dp_ps(transpose->ymm[0], _double, 0B11110001);
+        __m256 tmp = _mm256_dp_ps(transpose->f32_y[0], _double, 0B11110001);
-        tmp = _mm256_blend_ps(tmp, _mm256_dp_ps(transpose->ymm[0], _double, 0B11110010),
+        tmp = _mm256_blend_ps(tmp, _mm256_dp_ps(transpose->f32_y[0], _double, 0B11110010),
                               0B00100010);
         *ans = _mm256_castps256_ps128(
                 _mm256_permutevar8x32_ps(tmp, _mm256_loadu_epi32(avxdot_indices)));
@@ -29,34 +29,34 @@ INTERNAL void SN_DECL dot_41_avx(MAT_KERNEL(float) const *transpose, __m128 cons
 INTERNAL void SN_DECL dot_42_sse(MAT_KERNEL(float) const *transpose, __m128 const *vector,
                                  __m128 *ans)
 {
-        ans[0] = _mm_dp_ps(transpose->xmm[0], vector[0], 0B11110001);
+        ans[0] = _mm_dp_ps(transpose->f32_x[0], vector[0], 0B11110001);
-        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], vector[0], 0B11110001),
+        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], vector[0], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 2));
-        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], vector[0], 0B11110001),
+        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], vector[0], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 3));
-        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], vector[0], 0B11110001),
+        ans[0] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], vector[0], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 4));
 
-        ans[1] = _mm_dp_ps(transpose->xmm[0], vector[1], 0B11110001);
+        ans[1] = _mm_dp_ps(transpose->f32_x[0], vector[1], 0B11110001);
-        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[1], vector[1], 0B11110001),
+        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[1], vector[1], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 2));
-        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[2], vector[1], 0B11110001),
+        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[2], vector[1], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 3));
-        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->xmm[3], vector[1], 0B11110001),
+        ans[1] = _mm_insert_ps(*ans, _mm_dp_ps(transpose->f32_x[3], vector[1], 0B11110001),
                                SSE_INSERT_SIMPLE(1, 4));
 }
 
 INTERNAL void SN_DECL dot_42_avx(MAT_KERNEL(float) const *transpose, __m128 const *vector,
                                  __m128 *ans)
 {
-        __m256 temp = transpose->ymm[0];
+        __m256 temp = transpose->f32_y[0];
         __m256 vector_stream = _mm256_set_m128(vector[1], vector[0]);
 
         __m256 result = _mm256_dp_ps(temp, vector_stream, 0B11110001);
         temp = _mm256_permute2f128_ps(temp, temp, 0B00000001);
         result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11110010), 0B00100010);
 
-        temp = transpose->ymm[1];
+        temp = transpose->f32_y[1];
         result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11110100), 0B01000100);
         temp = _mm256_permute2f128_ps(temp, temp, 0B00000001);
         result = _mm256_blend_ps(result, _mm256_dp_ps(temp, vector_stream, 0B11111000), 0B10001000);