diff --git a/src/AiDotNet.Tensors/Engines/AiDotNetEngine.cs b/src/AiDotNet.Tensors/Engines/AiDotNetEngine.cs
index 47e4342706..30e7b05871 100644
--- a/src/AiDotNet.Tensors/Engines/AiDotNetEngine.cs
+++ b/src/AiDotNet.Tensors/Engines/AiDotNetEngine.cs
@@ -1,3 +1,7 @@
+using System;
+using System.Diagnostics;
+using AiDotNet.Tensors.Interfaces;
+
 namespace AiDotNet.Tensors.Engines;
 
 /// <summary>
@@ -116,18 +120,18 @@ public static bool AutoDetectAndConfigureGpu()
             if (gpuEngine.SupportsGpu)
             {
                 Current = gpuEngine;
-                Console.WriteLine($"[AiDotNet] GPU acceleration enabled: {gpuEngine.Name}");
+                Trace.WriteLine($"[AiDotNet] GPU acceleration enabled: {gpuEngine.Name}");
                 return true;
             }
 
             gpuEngine.Dispose();
-            Console.WriteLine("[AiDotNet] GPU not available, using CPU");
+            Trace.WriteLine("[AiDotNet] GPU not available, using CPU");
             return false;
         }
         catch (Exception ex)
         {
-            Console.WriteLine($"[AiDotNet] Failed to initialize DirectGpu: {ex.Message}");
-            Console.WriteLine("[AiDotNet] Falling back to CPU");
+            Trace.WriteLine($"[AiDotNet] Failed to initialize DirectGpu: {ex.Message}");
+            Trace.WriteLine("[AiDotNet] Falling back to CPU");
             return false;
         }
     }
@@ -143,7 +147,7 @@ public static bool AutoDetectAndConfigureGpu()
     public static void ResetToCpu()
     {
         Current = new CpuEngine();
-        Console.WriteLine("[AiDotNet] Reset to CPU engine");
+        Trace.WriteLine("[AiDotNet] Reset to CPU engine");
     }
 
     /// <summary>
diff --git a/src/AiDotNet.Tensors/Engines/CpuEngine.cs b/src/AiDotNet.Tensors/Engines/CpuEngine.cs
index 483fa82ee4..e9f919d869 100644
--- a/src/AiDotNet.Tensors/Engines/CpuEngine.cs
+++ b/src/AiDotNet.Tensors/Engines/CpuEngine.cs
@@ -1,4 +1,5 @@
 using System;
+using System.Diagnostics;
 using System.Runtime.CompilerServices;
 using AiDotNet.Tensors.Helpers;
 using AiDotNet.Tensors.Interfaces;
@@ -32,15 +33,39 @@ namespace AiDotNet.Tensors.Engines;
 /// - You're using custom numeric types
 /// </para>
 /// </remarks>
+/// <summary>
+/// CPU-based execution engine using INumericOperations for type-generic operations.
+/// </summary>
+/// <remarks>
+/// <para>
+/// CpuEngine provides the default execution backend for AiDotNet. It works with
+/// any numeric type that implements INumericOperations{T}, including decimal,
+/// BigInteger, and custom numeric types.
+/// </para>
+/// <para><b>For Beginners:</b> This is the standard, "always works" mode. 
+/// It uses your computer's main processor (CPU) to do the math. While not as 
+/// fast as a graphics card (GPU) for huge problems, it is very reliable 
+/// and works on every computer without any extra setup.</para>
+/// </remarks>
 public class CpuEngine : IEngine
 {
+    private const int CpuMatMulTileSizeFloat = 64;
+    private const int CpuMatMulTileSizeDouble = 32;
+    private const long CpuMatMulParallelThresholdOps = 8_000_000;
+    private static readonly bool CpuMatMulTraceEnabled =
+        Environment.GetEnvironmentVariable("AIDOTNET_CPU_MATMUL_TRACE") == "1";
+    private static readonly bool CpuMatMulSingleThread =
+        Environment.GetEnvironmentVariable("AIDOTNET_CPU_MATMUL_SINGLE_THREAD") == "1";
+
     /// <inheritdoc/>
     public string Name => "CPU Engine";
 
     /// <inheritdoc/>
     public bool SupportsGpu => false;
 
-    /// <inheritdoc/>
+    /// <summary>
+    /// Gets the direct GPU engine if available, for potential offloading.
+    /// </summary>
     public DirectGpu.DirectGpuEngine? DirectGpu => Engine.DirectGpu;
 
     /// <inheritdoc/>
@@ -1412,8 +1437,22 @@ public Matrix<T> MatrixMultiply<T>(Matrix<T> a, Matrix<T> b)
                 $"First matrix columns ({a.Columns}) must equal second matrix rows ({b.Rows}).");
         }
 
+#if NET6_0_OR_GREATER
+        if (typeof(T) == typeof(float) && a is Matrix<float> floatA && b is Matrix<float> floatB)
+        {
+            var result = MatrixMultiplyFloat(floatA, floatB);
+            return Unsafe.As<Matrix<float>, Matrix<T>>(ref result);
+        }
+
+        if (typeof(T) == typeof(double) && a is Matrix<double> doubleA && b is Matrix<double> doubleB)
+        {
+            var result = MatrixMultiplyDouble(doubleA, doubleB);
+            return Unsafe.As<Matrix<double>, Matrix<T>>(ref result);
+        }
+#endif
+
         var numOps = MathHelper.GetNumericOperations<T>();
-        var result = new Matrix<T>(a.Rows, b.Columns);
+        var genericResult = new Matrix<T>(a.Rows, b.Columns);
 
         // Standard O(n³) matrix multiplication
         for (int i = 0; i < a.Rows; i++)
@@ -1425,13 +1464,211 @@ public Matrix<T> MatrixMultiply<T>(Matrix<T> a, Matrix<T> b)
                 {
                     sum = numOps.Add(sum, numOps.Multiply(a[i, k], b[k, j]));
                 }
-                result[i, j] = sum;
+                genericResult[i, j] = sum;
             }
         }
 
+        return genericResult;
+    }
+
+#if NET6_0_OR_GREATER
+    private static Matrix<float> MatrixMultiplyFloat(Matrix<float> a, Matrix<float> b)
+    {
+        int m = a.Rows;
+        int k = a.Columns;
+        int n = b.Columns;
+        var result = new Matrix<float>(m, n);
+
+        if (m == 0 || n == 0 || k == 0)
+            return result;
+
+        int tileSize = GetCpuMatMulTileSize(m, n, k, CpuMatMulTileSizeFloat);
+        bool useParallel = ShouldParallelizeMatMul(m, n, k);
+
+        if (CpuMatMulSingleThread)
+            useParallel = false;
+
+        if (CpuMatMulTraceEnabled)
+        {
+            Trace.WriteLine($"[CpuMatMul] float {m}x{k}x{n} tile={tileSize} parallel={useParallel}");
+        }
+
+        var aSpan = a.AsSpan();
+        var bSpan = b.AsSpan();
+        var cSpan = result.AsWritableSpan();
+        cSpan.Clear();
+
+        if (useParallel)
+        {
+            int blockCount = (m + tileSize - 1) / tileSize;
+            Parallel.For(0, blockCount, block =>
+            {
+                int iStart = block * tileSize;
+                int iEnd = Math.Min(iStart + tileSize, m);
+                var aLocal = a.AsSpan();
+                var bLocal = b.AsSpan();
+                var cLocal = result.AsWritableSpan();
+                MultiplyFloatBlock(aLocal, bLocal, cLocal, m, k, n, iStart, iEnd, tileSize);
+            });
+        }
+        else
+        {
+            MultiplyFloatBlock(aSpan, bSpan, cSpan, m, k, n, 0, m, tileSize);
+        }
+
+        return result;
+    }
+
+    private static Matrix<double> MatrixMultiplyDouble(Matrix<double> a, Matrix<double> b)
+    {
+        int m = a.Rows;
+        int k = a.Columns;
+        int n = b.Columns;
+        var result = new Matrix<double>(m, n);
+
+        if (m == 0 || n == 0 || k == 0)
+            return result;
+
+        int tileSize = GetCpuMatMulTileSize(m, n, k, CpuMatMulTileSizeDouble);
+        bool useParallel = ShouldParallelizeMatMul(m, n, k);
+
+        if (CpuMatMulSingleThread)
+            useParallel = false;
+
+        if (CpuMatMulTraceEnabled)
+        {
+            Trace.WriteLine($"[CpuMatMul] double {m}x{k}x{n} tile={tileSize} parallel={useParallel}");
+        }
+
+        var aSpan = a.AsSpan();
+        var bSpan = b.AsSpan();
+        var cSpan = result.AsWritableSpan();
+        cSpan.Clear();
+
+        if (useParallel)
+        {
+            int blockCount = (m + tileSize - 1) / tileSize;
+            Parallel.For(0, blockCount, block =>
+            {
+                int iStart = block * tileSize;
+                int iEnd = Math.Min(iStart + tileSize, m);
+                var aLocal = a.AsSpan();
+                var bLocal = b.AsSpan();
+                var cLocal = result.AsWritableSpan();
+                MultiplyDoubleBlock(aLocal, bLocal, cLocal, m, k, n, iStart, iEnd, tileSize);
+            });
+        }
+        else
+        {
+            MultiplyDoubleBlock(aSpan, bSpan, cSpan, m, k, n, 0, m, tileSize);
+        }
+
         return result;
     }
 
+    private static int GetCpuMatMulTileSize(int m, int n, int k, int defaultTile)
+    {
+        int tile = defaultTile;
+        if (m < tile || n < tile || k < tile)
+        {
+            tile = Math.Max(16, defaultTile / 2);
+        }
+
+        return tile;
+    }
+
+    private static bool ShouldParallelizeMatMul(int m, int n, int k)
+    {
+        if (Environment.ProcessorCount <= 1)
+            return false;
+
+        if (m <= 0 || n <= 0 || k <= 0)
+            return false;
+
+        // Use double to avoid 64-bit overflow during operation count estimation
+        double ops = (double)m * n * k;
+        return ops >= CpuMatMulParallelThresholdOps;
+    }
+
+    private static void MultiplyFloatBlock(
+        ReadOnlySpan<float> a,
+        ReadOnlySpan<float> b,
+        Span<float> c,
+        int m,
+        int k,
+        int n,
+        int iStart,
+        int iEnd,
+        int tileSize)
+    {
+        for (int k0 = 0; k0 < k; k0 += tileSize)
+        {
+            int kEnd = Math.Min(k0 + tileSize, k);
+
+            for (int j0 = 0; j0 < n; j0 += tileSize)
+            {
+                int jEnd = Math.Min(j0 + tileSize, n);
+
+                for (int i = iStart; i < iEnd; i++)
+                {
+                    int rowOffset = i * n;
+                    int aRowOffset = i * k;
+                    for (int kk = k0; kk < kEnd; kk++)
+                    {
+                        float aik = a[aRowOffset + kk];
+                        int bRowOffset = kk * n;
+                        SimdVector.MatMulInnerLoopFloat(
+                            aik,
+                            b.Slice(bRowOffset + j0, jEnd - j0),
+                            c.Slice(rowOffset + j0, jEnd - j0),
+                            0,
+                            jEnd - j0);
+                    }
+                }
+            }
+        }
+    }
+
+    private static void MultiplyDoubleBlock(
+        ReadOnlySpan<double> a,
+        ReadOnlySpan<double> b,
+        Span<double> c,
+        int m,
+        int k,
+        int n,
+        int iStart,
+        int iEnd,
+        int tileSize)
+    {
+        for (int k0 = 0; k0 < k; k0 += tileSize)
+        {
+            int kEnd = Math.Min(k0 + tileSize, k);
+
+            for (int j0 = 0; j0 < n; j0 += tileSize)
+            {
+                int jEnd = Math.Min(j0 + tileSize, n);
+
+                for (int i = iStart; i < iEnd; i++)
+                {
+                    int rowOffset = i * n;
+                    int aRowOffset = i * k;
+                    for (int kk = k0; kk < kEnd; kk++)
+                    {
+                        double aik = a[aRowOffset + kk];
+                        int bRowOffset = kk * n;
+                        SimdVector.MatMulInnerLoopDouble(
+                            aik,
+                            b.Slice(bRowOffset + j0, jEnd - j0),
+                            c.Slice(rowOffset + j0, jEnd - j0),
+                            0,
+                            jEnd - j0);
+                    }
+                }
+            }
+        }
+    }
+#endif
+
     /// <inheritdoc/>
     public Vector<T> MatrixVectorMultiply<T>(Matrix<T> matrix, Vector<T> vector)
     {
@@ -1448,14 +1685,25 @@ public Vector<T> MatrixVectorMultiply<T>(Matrix<T> matrix, Vector<T> vector)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Vector<T>(matrix.Rows);
 
-        for (int i = 0; i < matrix.Rows; i++)
+        if ((long)matrix.Rows * matrix.Columns > 10000)
         {
-            T sum = numOps.Zero;
-            for (int j = 0; j < matrix.Columns; j++)
+            var vectorData = vector.Data;
+            var resultData = result.Data;
+            Parallel.For(0, matrix.Rows, i =>
+            {
+                var rowSpan = matrix.GetRowReadOnlySpan(i);
+                var vectorSpan = new ReadOnlySpan<T>(vectorData);
+                resultData[i] = numOps.Dot(rowSpan, vectorSpan);
+            });
+        }
+        else
+        {
+            var vectorSpan = vector.AsSpan();
+            var resultSpan = result.AsWritableSpan();
+            for (int i = 0; i < matrix.Rows; i++)
             {
-                sum = numOps.Add(sum, numOps.Multiply(matrix[i, j], vector[j]));
+                resultSpan[i] = numOps.Dot(matrix.GetRowReadOnlySpan(i), vectorSpan);
             }
-            result[i] = sum;
         }
 
         return result;
@@ -1494,13 +1742,7 @@ public Matrix<T> MatrixAdd<T>(Matrix<T> a, Matrix<T> b)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Matrix<T>(a.Rows, a.Columns);
 
-        for (int i = 0; i < a.Rows; i++)
-        {
-            for (int j = 0; j < a.Columns; j++)
-            {
-                result[i, j] = numOps.Add(a[i, j], b[i, j]);
-            }
-        }
+        numOps.Add(a.AsSpan(), b.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -1513,13 +1755,7 @@ public Matrix<T> MatrixMultiplyScalar<T>(Matrix<T> matrix, T scalar)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Matrix<T>(matrix.Rows, matrix.Columns);
 
-        for (int i = 0; i < matrix.Rows; i++)
-        {
-            for (int j = 0; j < matrix.Columns; j++)
-            {
-                result[i, j] = numOps.Multiply(matrix[i, j], scalar);
-            }
-        }
+        numOps.MultiplyScalar(matrix.AsSpan(), scalar, result.AsWritableSpan());
 
         return result;
     }
@@ -1531,16 +1767,9 @@ public Matrix<T> MatrixSubtract<T>(Matrix<T> a, Matrix<T> b)
         if (a.Rows != b.Rows || a.Columns != b.Columns)
             throw new ArgumentException("Matrix dimensions must match for subtraction");
 
+        var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Matrix<T>(a.Rows, a.Columns);
-
-        // VECTORIZED: Use existing Vector Subtract operation on each row
-        for (int i = 0; i < a.Rows; i++)
-        {
-            var rowA = a.GetRow(i);
-            var rowB = b.GetRow(i);
-            var diffRow = Subtract(rowA, rowB); // Reuse vectorized Vector Subtract
-            result.SetRow(i, diffRow);
-        }
+        numOps.Subtract(a.AsSpan(), b.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -1550,17 +1779,8 @@ public T MatrixSumOfSquares<T>(Matrix<T> matrix)
         if (matrix == null) throw new ArgumentNullException(nameof(matrix));
 
         var numOps = MathHelper.GetNumericOperations<T>();
-        T sum = numOps.Zero;
-
-        // VECTORIZED: Use existing DotProduct operation on each row
-        for (int i = 0; i < matrix.Rows; i++)
-        {
-            var row = matrix.GetRow(i);
-            T rowSumSquares = DotProduct(row, row); // row · row = sum of squares for row
-            sum = numOps.Add(sum, rowSumSquares);
-        }
-
-        return sum;
+        var dataSpan = matrix.AsSpan();
+        return numOps.Dot(dataSpan, dataSpan);
     }
 
     public void SwapColumns<T>(Matrix<T> matrix, int col1, int col2)
@@ -1593,39 +1813,27 @@ public Matrix<T> OuterProduct<T>(Vector<T> a, Vector<T> b)
         if (a == null) throw new ArgumentNullException(nameof(a));
         if (b == null) throw new ArgumentNullException(nameof(b));
 
+        var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Matrix<T>(a.Length, b.Length);
-        var aArray = a.ToArray();
-        var bArray = b.ToArray();
 
-        // Use SIMD-optimized TensorPrimitives for float type
-        if (typeof(T) == typeof(float) && bArray.Length >= 16)
+        if ((long)a.Length * b.Length > 10000)
         {
-            var bFloat = (float[])(object)bArray;
-            var aFloat = (float[])(object)aArray;
-
-            for (int i = 0; i < aFloat.Length; i++)
+            var aData = a.Data;
+            var bData = b.Data;
+            Parallel.For(0, aData.Length, i =>
             {
-                var rowData = new float[bFloat.Length];
-                // SIMD vectorized: multiply vector b by scalar a[i]
-                TensorPrimitives.Multiply(bFloat, aFloat[i], rowData);
-
-                // Copy result to matrix
-                for (int j = 0; j < bFloat.Length; j++)
-                {
-                    result[i, j] = (T)(object)rowData[j];
-                }
-            }
+                var rowSpan = result.GetRowSpan(i);
+                var bSpan = new ReadOnlySpan<T>(bData);
+                numOps.MultiplyScalar(bSpan, aData[i], rowSpan);
+            });
         }
         else
         {
-            // Fallback using NumOps
-            var numOps = MathHelper.GetNumericOperations<T>();
-            for (int i = 0; i < aArray.Length; i++)
+            var aData = a.Data;
+            var bSpan = b.AsSpan();
+            for (int i = 0; i < aData.Length; i++)
             {
-                for (int j = 0; j < bArray.Length; j++)
-                {
-                    result[i, j] = numOps.Multiply(aArray[i], bArray[j]);
-                }
+                numOps.MultiplyScalar(bSpan, aData[i], result.GetRowSpan(i));
             }
         }
 
@@ -1655,13 +1863,10 @@ public Vector<T> GetRow<T>(Matrix<T> matrix, int rowIndex)
             throw new ArgumentOutOfRangeException(nameof(rowIndex),
                 $"Row index {rowIndex} is out of range. Valid range is 0 to {matrix.Rows - 1}.");
 
-        // Row access is contiguous - can use direct array copy
-        var result = new T[matrix.Columns];
-        for (int j = 0; j < matrix.Columns; j++)
-        {
-            result[j] = matrix[rowIndex, j];
-        }
-        return new Vector<T>(result);
+        var numOps = MathHelper.GetNumericOperations<T>();
+        var result = new Vector<T>(matrix.Columns);
+        numOps.Copy(matrix.GetRowReadOnlySpan(rowIndex), result.AsWritableSpan());
+        return result;
     }
 
     public void SetColumn<T>(Matrix<T> matrix, int columnIndex, Vector<T> values)
@@ -1696,12 +1901,8 @@ public void SetRow<T>(Matrix<T> matrix, int rowIndex, Vector<T> values)
                 $"Values vector length ({values.Length}) must match matrix columns ({matrix.Columns}).",
                 nameof(values));
 
-        // Row access is contiguous - direct assignment
-        var valuesArray = values.ToArray();
-        for (int j = 0; j < matrix.Columns; j++)
-        {
-            matrix[rowIndex, j] = valuesArray[j];
-        }
+        var numOps = MathHelper.GetNumericOperations<T>();
+        numOps.Copy(values.AsSpan(), matrix.GetRowSpan(rowIndex));
     }
 
     #endregion
@@ -1850,10 +2051,7 @@ public Tensor<T> TensorAdd<T>(Tensor<T> a, Tensor<T> b)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(a.Shape);
 
-        for (int i = 0; i < a.Length; i++)
-        {
-            result.SetFlat(i, numOps.Add(a.GetFlat(i), b.GetFlat(i)));
-        }
+        numOps.Add(a.AsSpan(), b.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -1919,35 +2117,12 @@ public Tensor<T> TensorAddMany<T>(params Tensor<T>[] tensors)
 
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(referenceShape);
-        int length = tensors[0].Length;
+        var resultSpan = result.AsWritableSpan();
 
-        // Single-pass addition: accumulate all tensors element by element
-        // This avoids n-1 intermediate allocations from chained binary additions
-        if (length > 10000)
+        numOps.Copy(tensors[0].AsSpan(), resultSpan);
+        for (int t = 1; t < tensors.Length; t++)
         {
-            // Parallel execution for large tensors
-            Parallel.For(0, length, i =>
-            {
-                T sum = numOps.Zero;
-                for (int t = 0; t < tensors.Length; t++)
-                {
-                    sum = numOps.Add(sum, tensors[t].GetFlat(i));
-                }
-                result.SetFlat(i, sum);
-            });
-        }
-        else
-        {
-            // Sequential execution for smaller tensors (avoids parallel overhead)
-            for (int i = 0; i < length; i++)
-            {
-                T sum = numOps.Zero;
-                for (int t = 0; t < tensors.Length; t++)
-                {
-                    sum = numOps.Add(sum, tensors[t].GetFlat(i));
-                }
-                result.SetFlat(i, sum);
-            }
+            numOps.Add(resultSpan, tensors[t].AsSpan(), resultSpan);
         }
 
         return result;
@@ -1967,10 +2142,7 @@ public Tensor<T> TensorSubtract<T>(Tensor<T> a, Tensor<T> b)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(a.Shape);
 
-        for (int i = 0; i < a.Length; i++)
-        {
-            result.SetFlat(i, numOps.Subtract(a.GetFlat(i), b.GetFlat(i)));
-        }
+        numOps.Subtract(a.AsSpan(), b.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -1989,10 +2161,7 @@ public Tensor<T> TensorMultiply<T>(Tensor<T> a, Tensor<T> b)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(a.Shape);
 
-        for (int i = 0; i < a.Length; i++)
-        {
-            result.SetFlat(i, numOps.Multiply(a.GetFlat(i), b.GetFlat(i)));
-        }
+        numOps.Multiply(a.AsSpan(), b.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2018,35 +2187,12 @@ public Tensor<T> TensorMultiplyMany<T>(params Tensor<T>[] tensors)
 
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(referenceShape);
-        int length = tensors[0].Length;
+        var resultSpan = result.AsWritableSpan();
 
-        // Single-pass multiplication: accumulate all tensors element by element
-        // This avoids n-1 intermediate allocations from chained binary multiplications
-        if (length > 10000)
-        {
-            // Parallel execution for large tensors
-            Parallel.For(0, length, i =>
-            {
-                T product = numOps.One;
-                for (int t = 0; t < tensors.Length; t++)
-                {
-                    product = numOps.Multiply(product, tensors[t].GetFlat(i));
-                }
-                result.SetFlat(i, product);
-            });
-        }
-        else
+        numOps.Copy(tensors[0].AsSpan(), resultSpan);
+        for (int t = 1; t < tensors.Length; t++)
         {
-            // Sequential execution for smaller tensors (avoids parallel overhead)
-            for (int i = 0; i < length; i++)
-            {
-                T product = numOps.One;
-                for (int t = 0; t < tensors.Length; t++)
-                {
-                    product = numOps.Multiply(product, tensors[t].GetFlat(i));
-                }
-                result.SetFlat(i, product);
-            }
+            numOps.Multiply(resultSpan, tensors[t].AsSpan(), resultSpan);
         }
 
         return result;
@@ -2060,10 +2206,7 @@ public Tensor<T> TensorMultiplyScalar<T>(Tensor<T> tensor, T scalar)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        for (int i = 0; i < tensor.Length; i++)
-        {
-            result.SetFlat(i, numOps.Multiply(tensor.GetFlat(i), scalar));
-        }
+        numOps.MultiplyScalar(tensor.AsSpan(), scalar, result.AsWritableSpan());
 
         return result;
     }
@@ -2082,15 +2225,18 @@ public Tensor<T> TensorDivide<T>(Tensor<T> a, Tensor<T> b)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(a.Shape);
 
-        for (int i = 0; i < a.Length; i++)
+        var aSpan = a.AsSpan();
+        var bSpan = b.AsSpan();
+        var resultSpan = result.AsWritableSpan();
+        bool checkZero = !MathHelper.IsFloatingPoint<T>();
+        for (int i = 0; i < bSpan.Length; i++)
         {
-            // Check for division by zero
-            if (numOps.Equals(b.GetFlat(i), numOps.Zero))
+            var divisor = bSpan[i];
+            if (checkZero && numOps.Equals(divisor, numOps.Zero))
             {
                 throw new DivideByZeroException($"Division by zero at index {i}");
             }
-
-            result.SetFlat(i, numOps.Divide(a.GetFlat(i), b.GetFlat(i)));
+            resultSpan[i] = numOps.Divide(aSpan[i], divisor);
         }
 
         return result;
@@ -2343,20 +2489,7 @@ public Tensor<T> TensorLog<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Log(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Log(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Log(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2369,20 +2502,7 @@ public Tensor<T> TensorExp<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Exp(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Exp(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Exp(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2395,20 +2515,7 @@ public Tensor<T> TensorSqrt<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Sqrt(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Sqrt(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Sqrt(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2421,20 +2528,7 @@ public Tensor<T> TensorAbs<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Abs(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Abs(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Abs(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2447,20 +2541,7 @@ public Tensor<T> TensorNegate<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Negate(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Negate(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Negate(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2473,20 +2554,7 @@ public Tensor<T> TensorPower<T>(Tensor<T> tensor, T exponent)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Power(tensor.GetFlat(i), exponent));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Power(tensor.GetFlat(i), exponent));
-            }
-        }
+        numOps.Pow(tensor.AsSpan(), exponent, result.AsWritableSpan());
 
         return result;
     }
@@ -2502,18 +2570,22 @@ public Tensor<T> TensorPower<T>(Tensor<T> bases, Tensor<T> exponents)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(bases.Shape);
 
+        var baseData = bases.Data;
+        var expData = exponents.Data;
+        var resultData = result.Data;
+
         if (bases.Length > 10000)
         {
             Parallel.For(0, bases.Length, i =>
             {
-                result.SetFlat(i, numOps.Power(bases.GetFlat(i), exponents.GetFlat(i)));
+                resultData[i] = numOps.Power(baseData[i], expData[i]);
             });
         }
         else
         {
             for (int i = 0; i < bases.Length; i++)
             {
-                result.SetFlat(i, numOps.Power(bases.GetFlat(i), exponents.GetFlat(i)));
+                resultData[i] = numOps.Power(baseData[i], expData[i]);
             }
         }
 
@@ -2528,20 +2600,7 @@ public Tensor<T> TensorFloor<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Floor(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Floor(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Floor(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2554,20 +2613,7 @@ public Tensor<T> TensorCeiling<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Ceiling(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Ceiling(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Ceiling(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2580,20 +2626,7 @@ public Tensor<T> TensorFrac<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Frac(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Frac(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Frac(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2606,20 +2639,7 @@ public Tensor<T> TensorSin<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Sin(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Sin(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Sin(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -2632,20 +2652,7 @@ public Tensor<T> TensorCos<T>(Tensor<T> tensor)
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
 
-        if (tensor.Length > 10000)
-        {
-            Parallel.For(0, tensor.Length, i =>
-            {
-                result.SetFlat(i, numOps.Cos(tensor.GetFlat(i)));
-            });
-        }
-        else
-        {
-            for (int i = 0; i < tensor.Length; i++)
-            {
-                result.SetFlat(i, numOps.Cos(tensor.GetFlat(i)));
-            }
-        }
+        numOps.Cos(tensor.AsSpan(), result.AsWritableSpan());
 
         return result;
     }
@@ -3018,14 +3025,7 @@ public T TensorSum<T>(Tensor<T> tensor)
         if (tensor == null) throw new ArgumentNullException(nameof(tensor));
 
         var numOps = MathHelper.GetNumericOperations<T>();
-        T sum = numOps.Zero;
-
-        for (int i = 0; i < tensor.Length; i++)
-        {
-            sum = numOps.Add(sum, tensor.GetFlat(i));
-        }
-
-        return sum;
+        return numOps.Sum(tensor.AsSpan());
     }
 
     /// <inheritdoc/>
@@ -3091,35 +3091,26 @@ public T TensorMaxValue<T>(Tensor<T> tensor)
         if (tensor.Length == 0) throw new ArgumentException("Cannot compute max of empty tensor.", nameof(tensor));
 
         var numOps = MathHelper.GetNumericOperations<T>();
-        T maxVal = tensor.GetFlat(0);
-
-        // Parallel reduction for large tensors
         if (tensor.Length > 10000)
         {
             int workerCount = Environment.ProcessorCount;
             var localMaxes = new T[workerCount];
             var hasValue = new bool[workerCount];
             int chunkSize = (tensor.Length + workerCount - 1) / workerCount;
+            var data = tensor.Data;
 
             Parallel.For(0, workerCount, threadIdx =>
             {
                 int start = threadIdx * chunkSize;
-                int end = Math.Min(start + chunkSize, tensor.Length);
-                if (start >= tensor.Length) return;
+                int end = Math.Min(start + chunkSize, data.Length);
+                if (start >= data.Length) return;
 
-                T localMax = tensor.GetFlat(start);
-                for (int i = start + 1; i < end; i++)
-                {
-                    var val = tensor.GetFlat(i);
-                    if (numOps.GreaterThan(val, localMax))
-                        localMax = val;
-                }
-                localMaxes[threadIdx] = localMax;
+                localMaxes[threadIdx] = numOps.Max(new ReadOnlySpan<T>(data, start, end - start));
                 hasValue[threadIdx] = true;
             });
 
-            // Combine only populated slots
             bool first = true;
+            T maxVal = numOps.Zero;
             for (int i = 0; i < workerCount; i++)
             {
                 if (!hasValue[i]) continue;
@@ -3133,18 +3124,11 @@ public T TensorMaxValue<T>(Tensor<T> tensor)
                     maxVal = localMaxes[i];
                 }
             }
-        }
-        else
-        {
-            for (int i = 1; i < tensor.Length; i++)
-            {
-                var val = tensor.GetFlat(i);
-                if (numOps.GreaterThan(val, maxVal))
-                    maxVal = val;
-            }
+
+            return maxVal;
         }
 
-        return maxVal;
+        return numOps.Max(tensor.AsSpan());
     }
 
     /// <inheritdoc/>
@@ -3154,35 +3138,26 @@ public T TensorMinValue<T>(Tensor<T> tensor)
         if (tensor.Length == 0) throw new ArgumentException("Cannot compute min of empty tensor.", nameof(tensor));
 
         var numOps = MathHelper.GetNumericOperations<T>();
-        T minVal = tensor.GetFlat(0);
-
-        // Parallel reduction for large tensors
         if (tensor.Length > 10000)
         {
             int workerCount = Environment.ProcessorCount;
             var localMins = new T[workerCount];
             var hasValue = new bool[workerCount];
             int chunkSize = (tensor.Length + workerCount - 1) / workerCount;
+            var data = tensor.Data;
 
             Parallel.For(0, workerCount, threadIdx =>
             {
                 int start = threadIdx * chunkSize;
-                int end = Math.Min(start + chunkSize, tensor.Length);
-                if (start >= tensor.Length) return;
+                int end = Math.Min(start + chunkSize, data.Length);
+                if (start >= data.Length) return;
 
-                T localMin = tensor.GetFlat(start);
-                for (int i = start + 1; i < end; i++)
-                {
-                    var val = tensor.GetFlat(i);
-                    if (numOps.LessThan(val, localMin))
-                        localMin = val;
-                }
-                localMins[threadIdx] = localMin;
+                localMins[threadIdx] = numOps.Min(new ReadOnlySpan<T>(data, start, end - start));
                 hasValue[threadIdx] = true;
             });
 
-            // Combine only populated slots
             bool first = true;
+            T minVal = numOps.Zero;
             for (int i = 0; i < workerCount; i++)
             {
                 if (!hasValue[i]) continue;
@@ -3196,18 +3171,11 @@ public T TensorMinValue<T>(Tensor<T> tensor)
                     minVal = localMins[i];
                 }
             }
-        }
-        else
-        {
-            for (int i = 1; i < tensor.Length; i++)
-            {
-                var val = tensor.GetFlat(i);
-                if (numOps.LessThan(val, minVal))
-                    minVal = val;
-            }
+
+            return minVal;
         }
 
-        return minVal;
+        return numOps.Min(tensor.AsSpan());
     }
 
     /// <inheritdoc/>
@@ -12285,18 +12253,8 @@ public T TensorSumOfSquares<T>(Tensor<T> tensor)
         if (tensor == null) throw new ArgumentNullException(nameof(tensor));
 
         var numOps = MathHelper.GetNumericOperations<T>();
-        T sum = numOps.Zero;
-        int length = tensor.Length;
-
-        // Use SIMD-friendly sequential access pattern
-        var data = tensor.ToArray();
-        for (int i = 0; i < length; i++)
-        {
-            T val = data[i];
-            sum = numOps.Add(sum, numOps.Multiply(val, val));
-        }
-
-        return sum;
+        var dataSpan = tensor.AsSpan();
+        return numOps.Dot(dataSpan, dataSpan);
     }
 
     /// <inheritdoc/>
@@ -13371,12 +13329,7 @@ public Tensor<T> TensorAddScalar<T>(Tensor<T> tensor, T scalar)
 
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
-        int totalElements = tensor.Length;
-
-        Parallel.For(0, totalElements, i =>
-        {
-            result.SetFlat(i, numOps.Add(tensor.GetFlat(i), scalar));
-        });
+        numOps.AddScalar(tensor.AsSpan(), scalar, result.AsWritableSpan());
 
         return result;
     }
@@ -13388,12 +13341,7 @@ public Tensor<T> TensorSubtractScalar<T>(Tensor<T> tensor, T scalar)
 
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
-        int totalElements = tensor.Length;
-
-        Parallel.For(0, totalElements, i =>
-        {
-            result.SetFlat(i, numOps.Subtract(tensor.GetFlat(i), scalar));
-        });
+        numOps.SubtractScalar(tensor.AsSpan(), scalar, result.AsWritableSpan());
 
         return result;
     }
@@ -13405,12 +13353,7 @@ public Tensor<T> TensorDivideScalar<T>(Tensor<T> tensor, T scalar)
 
         var numOps = MathHelper.GetNumericOperations<T>();
         var result = new Tensor<T>(tensor.Shape);
-        int totalElements = tensor.Length;
-
-        Parallel.For(0, totalElements, i =>
-        {
-            result.SetFlat(i, numOps.Divide(tensor.GetFlat(i), scalar));
-        });
+        numOps.DivideScalar(tensor.AsSpan(), scalar, result.AsWritableSpan());
 
         return result;
     }
@@ -15014,11 +14957,18 @@ public Tensor<T> FusedLinear<T>(Tensor<T> input, Tensor<T> weights, Tensor<T>? b
             // Use optimized fused operations for float type
             if (typeof(T) == typeof(float))
             {
-                // Cast arrays directly (boxing avoids generic constraint issues)
-                var inputArray = (float[])(object)input.Data;
-                var weightsArray = (float[])(object)weights.Data;
-                var biasArray = bias != null ? (float[])(object)bias.Data : null;
-                var outputArray = (float[])(object)result.Data;
+                var inputData = input.Data;
+                var weightsData = weights.Data;
+                var outputData = result.Data;
+                var inputArray = Unsafe.As<T[], float[]>(ref inputData);
+                var weightsArray = Unsafe.As<T[], float[]>(ref weightsData);
+                var outputArray = Unsafe.As<T[], float[]>(ref outputData);
+                float[]? biasArray = null;
+                if (bias is not null)
+                {
+                    var biasData = bias.Data;
+                    biasArray = Unsafe.As<T[], float[]>(ref biasData);
+                }
 
                 CpuFusedOperations.FusedGemmBiasActivation(
                     inputArray, weightsArray, biasArray, outputArray,
@@ -15030,10 +14980,18 @@ public Tensor<T> FusedLinear<T>(Tensor<T> input, Tensor<T> weights, Tensor<T>? b
             // Use optimized fused operations for double type
             if (typeof(T) == typeof(double))
             {
-                var inputArray = (double[])(object)input.Data;
-                var weightsArray = (double[])(object)weights.Data;
-                var biasArray = bias != null ? (double[])(object)bias.Data : null;
-                var outputArray = (double[])(object)result.Data;
+                var inputData = input.Data;
+                var weightsData = weights.Data;
+                var outputData = result.Data;
+                var inputArray = Unsafe.As<T[], double[]>(ref inputData);
+                var weightsArray = Unsafe.As<T[], double[]>(ref weightsData);
+                var outputArray = Unsafe.As<T[], double[]>(ref outputData);
+                double[]? biasArray = null;
+                if (bias is not null)
+                {
+                    var biasData = bias.Data;
+                    biasArray = Unsafe.As<T[], double[]>(ref biasData);
+                }
 
                 CpuFusedOperations.FusedGemmBiasActivation(
                     inputArray, weightsArray, biasArray, outputArray,
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/DirectGpuEngine.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/DirectGpuEngine.cs
index f733dc561e..997d26a199 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/DirectGpuEngine.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/DirectGpuEngine.cs
@@ -1,3 +1,4 @@
+using System.Diagnostics;
 using AiDotNet.Tensors.Engines.DirectGpu.CUDA;
 using AiDotNet.Tensors.Engines.DirectGpu.HIP;
 using AiDotNet.Tensors.Engines.DirectGpu.OpenCL;
@@ -13,19 +14,12 @@ namespace AiDotNet.Tensors.Engines.DirectGpu;
 /// <remarks>
 /// <para><b>Design Philosophy:</b></para>
 /// <para>
-/// This engine provides the fastest GPU path by:
-/// 1. Using float32-only kernels (optimal GPU performance)
-/// 2. Converting generic types at the boundary (preserves clean API)
-/// 3. Implementing optimizations CLBlast misses (tensor cores, fusion, double-buffering)
+/// This engine provides the fastest GPU path by using specialized math kernels 
+/// that run directly on your graphics card.
 /// </para>
-/// <para><b>Fallback Tiers:</b></para>
-/// <code>
-/// DirectGpuEngine (this) - custom optimized kernels
-///     ↓ fallback
-/// CLBlast - tuned but missing optimizations
-///     ↓ fallback
-/// CPU - always available
-/// </code>
+/// <para><b>For Beginners:</b> Think of this as a high-speed lane for math. 
+/// It bypasses the normal slow ways computers do math and uses the massive 
+/// parallel power of your graphics card to speed up AI calculations.</para>
 /// </remarks>
 public sealed class DirectGpuEngine : IDisposable
 {
@@ -35,6 +29,8 @@ public sealed class DirectGpuEngine : IDisposable
     private bool _disposed;
     private const string BackendOrderEnvVar = "AIDOTNET_DIRECTGPU_BACKENDS";
     private static readonly string[] DefaultBackendOrder = new[] { "cuda", "opencl", "hip" };
+    private static readonly bool GemmValidateEnabled =
+        Environment.GetEnvironmentVariable("AIDOTNET_GEMM_VALIDATE") == "1";
 
     /// <summary>
     /// Gets whether the direct GPU engine is available.
@@ -88,30 +84,28 @@ public sealed class DirectGpuEngine : IDisposable
     public long LocalMemoryBytes => _backend?.LocalMemoryBytes ?? 0;
 
     /// <summary>
-    /// Initializes the DirectGpuEngine, automatically selecting the best available backend.
+    /// Initializes a new instance of the <see cref="DirectGpuEngine"/> class.
     /// </summary>
     /// <remarks>
-    /// Backend selection order:
-    /// 1. CUDA (NVIDIA GPUs with NVRTC available)
-    /// 2. OpenCL (works on AMD, Intel, and NVIDIA)
-    /// 3. HIP (AMD GPUs with MFMA support - MI100/200/300, RDNA3)
-    /// </remarks>
+    /// <para><b>For Beginners:</b> This method searches your computer for a 
+/// graphics card it can use. It tries CUDA first (NVIDIA), then OpenCL (standard), 
+/// then HIP (AMD) until it finds a way to run fast math.</para>
+/// </remarks>
     public DirectGpuEngine()
     {
-        // Initialize fusion manager
+        // Initialize local state first
         _fusionManager = new KernelFusionManager();
-
-        Console.WriteLine("[DirectGpuEngine] Initializing GPU backends...");
+        Trace.WriteLine("[DirectGpuEngine] Initializing GPU backends...");
 
         var backendOrder = GetBackendOrderFromEnv();
         if (backendOrder.Count == 0)
         {
-            Console.WriteLine($"[DirectGpuEngine] Direct GPU backends disabled via {BackendOrderEnvVar}.");
+            Trace.WriteLine($"[DirectGpuEngine] Direct GPU backends disabled via {BackendOrderEnvVar}.");
             _isAvailable = false;
             return;
         }
 
-        Console.WriteLine($"[DirectGpuEngine] Backend order: {string.Join(", ", backendOrder)}");
+        Trace.WriteLine($"[DirectGpuEngine] Backend order: {string.Join(", ", backendOrder)}");
         foreach (var backendName in backendOrder)
         {
             var backend = TryCreateBackend(backendName);
@@ -123,7 +117,7 @@ public DirectGpuEngine()
             }
         }
 
-        Console.WriteLine("[DirectGpuEngine] No GPU backends available. Falling back to CPU.");
+        Trace.WriteLine("[DirectGpuEngine] No GPU backends available. Falling back to CPU.");
         _isAvailable = false;
     }
 
@@ -170,7 +164,7 @@ private static IReadOnlyList<string> GetBackendOrderFromEnv()
                 case "":
                     break;
                 default:
-                    Console.WriteLine($"[DirectGpuEngine] Unknown backend token '{token}' in {BackendOrderEnvVar}. Expected: cuda, opencl, hip, auto, none.");
+                    Trace.WriteLine($"[DirectGpuEngine] Unknown backend token '{token}' in {BackendOrderEnvVar}. Expected: cuda, opencl, hip, auto, none.");
                     break;
             }
         }
@@ -199,45 +193,45 @@ private static void AddBackend(List<string> backends, HashSet<string> seen, stri
     {
         try
         {
-            Console.WriteLine("[DirectGpuEngine] Checking CUDA availability...");
+            Trace.WriteLine("[DirectGpuEngine] Checking CUDA availability...");
             if (CudaBackend.IsCudaAvailable)
             {
-                Console.WriteLine("[DirectGpuEngine] Creating CUDA backend...");
+                Trace.WriteLine("[DirectGpuEngine] Creating CUDA backend...");
                 var cudaBackend = new CudaBackend();
-                Console.WriteLine($"[DirectGpuEngine] CUDA backend created, IsAvailable = {cudaBackend.IsAvailable}");
+                Trace.WriteLine($"[DirectGpuEngine] CUDA backend created, IsAvailable = {cudaBackend.IsAvailable}");
 
                 if (cudaBackend.IsAvailable)
                 {
-                    Console.WriteLine($"[DirectGpuEngine] SUCCESS: Using CUDA backend on {cudaBackend.DeviceName}");
+                    Trace.WriteLine($"[DirectGpuEngine] SUCCESS: Using CUDA backend on {cudaBackend.DeviceName}");
                     System.Diagnostics.Debug.WriteLine($"DirectGpuEngine: Using CUDA backend on {cudaBackend.DeviceName}");
                     return cudaBackend;
                 }
-                Console.WriteLine("[DirectGpuEngine] CUDA backend created but not available, disposing...");
+                Trace.WriteLine("[DirectGpuEngine] CUDA backend created but not available, disposing...");
                 cudaBackend.Dispose();
             }
             else
             {
-                Console.WriteLine("[DirectGpuEngine] CUDA is not available on this system");
+                Trace.WriteLine("[DirectGpuEngine] CUDA is not available on this system");
             }
         }
         catch (DllNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] CUDA library not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] CUDA library not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"CUDA library not found: {ex.Message}");
         }
         catch (TypeInitializationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] CUDA type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] CUDA type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
             System.Diagnostics.Debug.WriteLine($"CUDA type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
         }
         catch (InvalidOperationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] CUDA initialization failed: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] CUDA initialization failed: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"CUDA initialization failed: {ex.Message}");
         }
         catch (EntryPointNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] CUDA function not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] CUDA function not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"CUDA function not found: {ex.Message}");
         }
 
@@ -248,40 +242,40 @@ private static void AddBackend(List<string> backends, HashSet<string> seen, stri
     {
         try
         {
-            Console.WriteLine("[DirectGpuEngine] Checking OpenCL availability...");
-            Console.WriteLine($"[DirectGpuEngine] OpenClBackend.IsOpenClAvailable = {OpenClBackend.IsOpenClAvailable}");
+            Trace.WriteLine("[DirectGpuEngine] Checking OpenCL availability...");
+            Trace.WriteLine($"[DirectGpuEngine] OpenClBackend.IsOpenClAvailable = {OpenClBackend.IsOpenClAvailable}");
 
-            Console.WriteLine("[DirectGpuEngine] Creating OpenCL backend...");
+            Trace.WriteLine("[DirectGpuEngine] Creating OpenCL backend...");
             var openClBackend = new OpenClBackend();
-            Console.WriteLine($"[DirectGpuEngine] OpenCL backend created, IsAvailable = {openClBackend.IsAvailable}");
+            Trace.WriteLine($"[DirectGpuEngine] OpenCL backend created, IsAvailable = {openClBackend.IsAvailable}");
 
             if (openClBackend.IsAvailable)
             {
-                Console.WriteLine($"[DirectGpuEngine] SUCCESS: Using OpenCL backend on {openClBackend.DeviceName}");
+                Trace.WriteLine($"[DirectGpuEngine] SUCCESS: Using OpenCL backend on {openClBackend.DeviceName}");
                 System.Diagnostics.Debug.WriteLine($"DirectGpuEngine: Using OpenCL backend on {openClBackend.DeviceName}");
                 return openClBackend;
             }
-            Console.WriteLine("[DirectGpuEngine] OpenCL backend created but not available, disposing...");
+            Trace.WriteLine("[DirectGpuEngine] OpenCL backend created but not available, disposing...");
             openClBackend.Dispose();
         }
         catch (DllNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] OpenCL library not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] OpenCL library not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"OpenCL library not found: {ex.Message}");
         }
         catch (TypeInitializationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] OpenCL type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] OpenCL type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
             System.Diagnostics.Debug.WriteLine($"OpenCL type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
         }
         catch (InvalidOperationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] OpenCL initialization failed: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] OpenCL initialization failed: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"OpenCL initialization failed: {ex.Message}");
         }
         catch (EntryPointNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] OpenCL function not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] OpenCL function not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"OpenCL function not found: {ex.Message}");
         }
 
@@ -292,43 +286,43 @@ private static void AddBackend(List<string> backends, HashSet<string> seen, stri
     {
         try
         {
-            Console.WriteLine("[DirectGpuEngine] Checking HIP availability...");
+            Trace.WriteLine("[DirectGpuEngine] Checking HIP availability...");
             if (HipBackend.IsHipAvailable)
             {
-                Console.WriteLine("[DirectGpuEngine] HIP is available, creating HIP backend...");
+                Trace.WriteLine("[DirectGpuEngine] HIP is available, creating HIP backend...");
                 var hipBackend = new HipBackend();
                 if (hipBackend.IsAvailable)
                 {
-                    Console.WriteLine($"[DirectGpuEngine] SUCCESS: Using HIP backend with {hipBackend.Architecture} architecture");
+                    Trace.WriteLine($"[DirectGpuEngine] SUCCESS: Using HIP backend with {hipBackend.Architecture} architecture");
                     System.Diagnostics.Debug.WriteLine($"DirectGpuEngine: Using HIP backend with {hipBackend.Architecture} architecture");
                     return hipBackend;
                 }
-                Console.WriteLine("[DirectGpuEngine] HIP backend created but not available, disposing...");
+                Trace.WriteLine("[DirectGpuEngine] HIP backend created but not available, disposing...");
                 hipBackend.Dispose();
             }
             else
             {
-                Console.WriteLine("[DirectGpuEngine] HIP is not available on this system");
+                Trace.WriteLine("[DirectGpuEngine] HIP is not available on this system");
             }
         }
         catch (DllNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] HIP library not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] HIP library not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"HIP library not found: {ex.Message}");
         }
         catch (TypeInitializationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] HIP type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] HIP type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
             System.Diagnostics.Debug.WriteLine($"HIP type initialization failed: {ex.InnerException?.Message ?? ex.Message}");
         }
         catch (InvalidOperationException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] HIP initialization failed: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] HIP initialization failed: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"HIP initialization failed: {ex.Message}");
         }
         catch (EntryPointNotFoundException ex)
         {
-            Console.WriteLine($"[DirectGpuEngine] HIP function not found: {ex.Message}");
+            Trace.WriteLine($"[DirectGpuEngine] HIP function not found: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"HIP function not found: {ex.Message}");
         }
 
@@ -400,10 +394,22 @@ public static T[] FromFloatArray<T>(float[] data)
         // Download result
         float[] resultFloat = _backend.DownloadBuffer(bufferC);
 
+        if (GemmValidateEnabled && IsAnyNonFinite(resultFloat, out int badIndex))
+        {
+            Trace.WriteLine($"[DirectGpuEngine] GEMM produced non-finite values (first index {badIndex}). Falling back to CPU.");
+            return null;
+        }
+
         // Convert back to T
         return FromFloatArray<T>(resultFloat);
     }
 
+    private static bool IsAnyNonFinite(float[] data, out int badIndex)
+    {
+        var numOps = MathHelper.GetNumericOperations<float>();
+        return numOps.IsAnyNonFinite(data, out badIndex);
+    }
+
     /// <summary>
     /// Matrix multiplication with pre-allocated GPU weight buffer (for cached weights).
     /// </summary>
@@ -417,6 +423,11 @@ public static T[] FromFloatArray<T>(float[] data)
         using var bufferC = _backend.MatMul(bufferInput, cachedWeights, M, N, K);
 
         float[] resultFloat = _backend.DownloadBuffer(bufferC);
+        if (GemmValidateEnabled && IsAnyNonFinite(resultFloat, out int badIndex))
+        {
+            Trace.WriteLine($"[DirectGpuEngine] GEMM produced non-finite values (first index {badIndex}). Falling back to CPU.");
+            return null;
+        }
         return FromFloatArray<T>(resultFloat);
     }
 
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/GemmBenchmark.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/GemmBenchmark.cs
index 3f334f2736..eca3fd11f7 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/GemmBenchmark.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/GemmBenchmark.cs
@@ -132,7 +132,7 @@ public GemmBenchmarkResult[] RunSuite()
         {
             var result = BenchmarkSingle(size, size, size);
             results.Add(result);
-            Console.WriteLine(result);
+            Trace.WriteLine(result);
         }
 
         // Neural network layer sizes (typical Dense layer dimensions)
@@ -150,7 +150,7 @@ public GemmBenchmarkResult[] RunSuite()
         {
             var result = BenchmarkSingle(M, N, K);
             results.Add(result);
-            Console.WriteLine(result);
+            Trace.WriteLine(result);
         }
 
         return results.ToArray();
@@ -201,39 +201,39 @@ public string GenerateReport(GemmBenchmarkResult[] results)
     /// </summary>
     public static void QuickTest()
     {
-        Console.WriteLine("DirectGpu GEMM Quick Test");
-        Console.WriteLine("=========================");
+        Trace.WriteLine("DirectGpu GEMM Quick Test");
+        Trace.WriteLine("=========================");
 
         using var engine = new DirectGpuEngine();
 
         if (!engine.IsAvailable)
         {
-            Console.WriteLine("ERROR: DirectGpu not available.");
-            Console.WriteLine("Make sure OpenCL is installed and a compatible GPU is present.");
+            Trace.WriteLine("ERROR: DirectGpu not available.");
+            Trace.WriteLine("Make sure OpenCL is installed and a compatible GPU is present.");
             return;
         }
 
-        Console.WriteLine($"GPU: {engine.DeviceName}");
-        Console.WriteLine($"Vendor: {engine.DeviceVendor}");
-        Console.WriteLine($"Compute Units: {engine.ComputeUnits}");
-        Console.WriteLine($"Global Memory: {engine.GlobalMemoryGB:F1} GB");
-        Console.WriteLine();
+        Trace.WriteLine($"GPU: {engine.DeviceName}");
+        Trace.WriteLine($"Vendor: {engine.DeviceVendor}");
+        Trace.WriteLine($"Compute Units: {engine.ComputeUnits}");
+        Trace.WriteLine($"Global Memory: {engine.GlobalMemoryGB:F1} GB");
+        Trace.WriteLine("");
 
         var benchmark = new GemmBenchmark(engine);
 
         // Quick test with 2048x2048
-        Console.WriteLine("Running quick benchmark (2048x2048x2048)...");
+        Trace.WriteLine("Running quick benchmark (2048x2048x2048)...");
         var result = benchmark.BenchmarkSingle(2048, 2048, 2048);
 
-        Console.WriteLine();
-        Console.WriteLine($"Result: {result.GFlops:N0} GFLOPS ({result.TFlops:F2} TFLOPS)");
-        Console.WriteLine($"Target: 25,000 GFLOPS (25.0 TFLOPS)");
-        Console.WriteLine($"Status: {(result.MetTarget ? "TARGET MET!" : "Below target")}");
+        Trace.WriteLine("");
+        Trace.WriteLine($"Result: {result.GFlops:N0} GFLOPS ({result.TFlops:F2} TFLOPS)");
+        Trace.WriteLine($"Target: 25,000 GFLOPS (25.0 TFLOPS)");
+        Trace.WriteLine($"Status: {(result.MetTarget ? "TARGET MET!" : "Below target")}");
 
         if (result.GFlops > 0)
         {
-            Console.WriteLine();
-            Console.WriteLine($"vs CLBlast (2,500 GFLOPS): {result.GFlops / 2500:F1}x faster");
+            Trace.WriteLine("");
+            Trace.WriteLine($"vs CLBlast (2,500 GFLOPS): {result.GFlops / 2500:F1}x faster");
         }
     }
 }
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipBackend.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipBackend.cs
index 51ab9e8468..28295c5f19 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipBackend.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipBackend.cs
@@ -4,6 +4,7 @@
 
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Runtime.InteropServices;
 using AiDotNet.Tensors.Engines.DirectGpu.HIP.Kernels;
 using AiDotNet.Tensors.Engines.DirectGpu.Sparsity;
@@ -207,7 +208,7 @@ public HipBackend(int deviceIndex)
 
             if (!hasGemmKernel)
             {
-                Console.WriteLine("[HipBackend] No GEMM kernels compiled - backend not available");
+                Trace.WriteLine("[HipBackend] No GEMM kernels compiled - backend not available");
                 IsAvailable = false;
                 return;
             }
@@ -261,7 +262,7 @@ private void CompileKernels()
         string source = HipMfmaKernel.GetSource();
         string compileFlags = HipMfmaKernel.GetCompileFlags(_architecture);
 
-        Console.WriteLine($"[HipBackend] Compiling kernels for {_architecture} with flags: {compileFlags}");
+        Trace.WriteLine($"[HipBackend] Compiling kernels for {_architecture} with flags: {compileFlags}");
 
         try
         {
@@ -346,12 +347,12 @@ private void CompileKernels()
             CompileKernelModule(Kernels.HipSpatialTransformerKernels.GetSource(), "spatial_transformer", ref _spatialTransformerModule,
                 Kernels.HipSpatialTransformerKernels.GetKernelNames());
 
-            Console.WriteLine($"[HipBackend] Kernel compilation complete. Available kernels: {_kernelCache.Count}");
+            Trace.WriteLine($"[HipBackend] Kernel compilation complete. Available kernels: {_kernelCache.Count}");
             System.Diagnostics.Debug.WriteLine($"HIP kernels compiled successfully for {_architecture}. Total: {_kernelCache.Count}");
         }
         catch (Exception ex)
         {
-            Console.WriteLine($"[HipBackend] Kernel compilation EXCEPTION: {ex.GetType().Name}: {ex.Message}");
+            Trace.WriteLine($"[HipBackend] Kernel compilation EXCEPTION: {ex.GetType().Name}: {ex.Message}");
             System.Diagnostics.Debug.WriteLine($"HIP kernel compilation failed: {ex.Message}");
         }
     }
@@ -371,7 +372,7 @@ private void CompileKernelModule(string source, string moduleName, ref IntPtr mo
 
         if (rtcResult != HipRtcResult.Success)
         {
-            Console.WriteLine($"[HipBackend] Failed to create program for {moduleName}: {rtcResult}");
+            Trace.WriteLine($"[HipBackend] Failed to create program for {moduleName}: {rtcResult}");
             return;
         }
 
@@ -395,7 +396,7 @@ private void CompileKernelModule(string source, string moduleName, ref IntPtr mo
                 HipNativeBindings.hiprtcGetProgramLog(prog, logPtr);
                 string log = Marshal.PtrToStringAnsi(logPtr) ?? "";
                 Marshal.FreeHGlobal(logPtr);
-                Console.WriteLine($"[HipBackend] Compile failed for {moduleName}: {log}");
+                Trace.WriteLine($"[HipBackend] Compile failed for {moduleName}: {log}");
             }
             HipNativeBindings.hiprtcDestroyProgram(ref prog);
             return;
@@ -424,7 +425,7 @@ private void CompileKernelModule(string source, string moduleName, ref IntPtr mo
 
         if (hipResult != HipError.Success)
         {
-            Console.WriteLine($"[HipBackend] Failed to load module {moduleName}: {hipResult}");
+            Trace.WriteLine($"[HipBackend] Failed to load module {moduleName}: {hipResult}");
             return;
         }
 
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipNativeBindings.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipNativeBindings.cs
index 131020e67b..daeaf0090e 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipNativeBindings.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/HIP/HipNativeBindings.cs
@@ -3,6 +3,7 @@
 // Provides direct access to HIP for MFMA kernel execution.
 
 using System;
+using System.Diagnostics;
 using System.Runtime.InteropServices;
 
 namespace AiDotNet.Tensors.Engines.DirectGpu.HIP;
@@ -290,7 +291,7 @@ private static void InitializeDllSearchPath()
     private static void LogDiagnostic(string message)
     {
         if (EnableDiagnostics)
-            Console.WriteLine(message);
+            Trace.WriteLine(message);
     }
 
     /// <summary>
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/DynamicGemmKernel.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/DynamicGemmKernel.cs
index 35451cce84..c16894a67a 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/DynamicGemmKernel.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/DynamicGemmKernel.cs
@@ -158,7 +158,7 @@ public static int GetEffectiveVariant(int M, int N, int K)
             {
                 if (EnableDiagnostics)
                 {
-                    Console.WriteLine($"[DynamicGemm] Size {M}x{N}x{K} < {MinSwizzleSize}, falling back from variant {KernelVariant} to CLBlast baseline");
+                    Trace.WriteLine($"[DynamicGemm] Size {M}x{N}x{K} < {MinSwizzleSize}, falling back from variant {KernelVariant} to CLBlast baseline");
                 }
                 return 0; // Fall back to CLBlast baseline
             }
@@ -255,12 +255,12 @@ private static void LogDiag(string message)
                 catch
                 {
                     // Fall back to console on error
-                    Console.WriteLine(logLine);
+                    Trace.WriteLine(logLine);
                 }
             }
             else
             {
-                Console.WriteLine(logLine);
+                Trace.WriteLine(logLine);
             }
         }
     }
@@ -374,10 +374,10 @@ private DirectOpenClKernel GetKernelWithVariant(GemmConfig config, int variant)
             if (EnableDiagnostics)
             {
                 var lines = source.Split('\n');
-                Console.WriteLine($"[DynamicGemm] Kernel source (first 30 lines):");
+                Trace.WriteLine($"[DynamicGemm] Kernel source (first 30 lines):");
                 for (int i = 0; i < Math.Min(30, lines.Length); i++)
                 {
-                    Console.WriteLine($"  {i + 1:D3}: {lines[i].TrimEnd()}");
+                    Trace.WriteLine($"  {i + 1:D3}: {lines[i].TrimEnd()}");
                 }
             }
             throw;
@@ -509,9 +509,9 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
                     var source = ClBlastXgemmKernel.BuildSourceWithSwizzle(config, gemmK, 0x0F);
                     if (EnableDiagnostics)
                     {
-                        Console.WriteLine($"[DynamicGemm] SELECTED XOR SWIZZLE kernel: {config.KernelName} GEMMK={gemmK}");
-                        Console.WriteLine($"[DynamicGemm] Swizzle defines present: LDS_SWIZZLE_A={source.Contains("LDS_SWIZZLE_A(kg, mg)")}, LDS_STRIDE_A={source.Contains("LDS_STRIDE_A")}");
-                        Console.WriteLine($"[DynamicGemm] Original pattern present: alm[kg*(MWG/VWM)={source.Contains("alm[kg*(MWG/VWM)")}");
+                        Trace.WriteLine($"[DynamicGemm] SELECTED XOR SWIZZLE kernel: {config.KernelName} GEMMK={gemmK}");
+                        Trace.WriteLine($"[DynamicGemm] Swizzle defines present: LDS_SWIZZLE_A={source.Contains("LDS_SWIZZLE_A(kg, mg)")}, LDS_STRIDE_A={source.Contains("LDS_STRIDE_A")}");
+                        Trace.WriteLine($"[DynamicGemm] Original pattern present: alm[kg*(MWG/VWM)={source.Contains("alm[kg*(MWG/VWM)")}");
                     }
                     return source;
                 }
@@ -519,14 +519,14 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
                 case 2:
                     // RDNA1 optimized - XOR swizzle + Wave32 hints
                     if (EnableDiagnostics)
-                        Console.WriteLine($"[DynamicGemm] SELECTED RDNA1 OPTIMIZED kernel: {config.KernelName} GEMMK={gemmK}");
+                        Trace.WriteLine($"[DynamicGemm] SELECTED RDNA1 OPTIMIZED kernel: {config.KernelName} GEMMK={gemmK}");
                     return ClBlastXgemmKernel.BuildSourceOptimizedRdna1(config, gemmK, 0x0F, true);
 
                 case 0:
                 default:
                     // Original CLBlast baseline (default)
                     if (EnableDiagnostics)
-                        Console.WriteLine($"[DynamicGemm] SELECTED CLBlast BASELINE kernel: {config.KernelName} GEMMK={gemmK}");
+                        Trace.WriteLine($"[DynamicGemm] SELECTED CLBlast BASELINE kernel: {config.KernelName} GEMMK={gemmK}");
                     return ClBlastXgemmKernel.BuildSource(config, gemmK);
             }
         }
@@ -689,7 +689,7 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
 
         if (EnableDiagnostics && config.UseTrueVectorLDS)
         {
-            Console.WriteLine($"[DynamicGemm] UseTrueVectorLDS check: config.UseTrueVectorLDS={config.UseTrueVectorLDS}, VWM={VWM}>1={VWM > 1}, MWI%VWM={MWI % VWM}==0, NWI%VWN={NWI % VWN}==0, result={useTrueVectorLDS}");
+            Trace.WriteLine($"[DynamicGemm] UseTrueVectorLDS check: config.UseTrueVectorLDS={config.UseTrueVectorLDS}, VWM={VWM}>1={VWM > 1}, MWI%VWM={MWI % VWM}==0, NWI%VWN={NWI % VWN}==0, result={useTrueVectorLDS}");
         }
 
         // COOPERATIVE LOADING KERNEL (CLBlast-style MDIMA/NDIMB)
@@ -707,7 +707,7 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
             // Use TRUE CLBlast-style vectorized LDS kernel
             // This achieves maximum performance by using vector types throughout
             if (EnableDiagnostics)
-                Console.WriteLine($"[DynamicGemm] SELECTED TRUE VECTORIZED kernel: {config.KernelName} VWM={VWM} VWN={VWN} MWI={MWI} NWI={NWI}");
+                Trace.WriteLine($"[DynamicGemm] SELECTED TRUE VECTORIZED kernel: {config.KernelName} VWM={VWM} VWN={VWN} MWI={MWI} NWI={NWI}");
             GenerateCLBlastTrueVectorizedKernel(sb, MWI, NWI, VWM, VWN, KWI, KREG, config.UseColumnMajorA);
         }
         else if (useCooperativeLoading)
@@ -715,7 +715,7 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
             // Use cooperative loading kernel - MDIMA/NDIMB differ from MDIMC/NDIMC
             // This is how CLBlast achieves maximum memory bandwidth
             if (EnableDiagnostics)
-                Console.WriteLine($"[DynamicGemm] SELECTED COOPERATIVE kernel: {config.KernelName} MDIMA={MDIMA} NDIMB={NDIMB}");
+                Trace.WriteLine($"[DynamicGemm] SELECTED COOPERATIVE kernel: {config.KernelName} MDIMA={MDIMA} NDIMB={NDIMB}");
             GenerateCooperativeLoadingKernel(sb, MWI, NWI, VWM, VWN, KWI, KREG, MDIMA, NDIMB, MDIMC, NDIMC);
         }
         else if (isHighOccupancy)
@@ -723,21 +723,21 @@ private static string GenerateKernelSource(GemmConfig config, int variant)
             // Use high-occupancy kernel with TRUE double-buffering (ping-pong)
             // This hides 100% of memory latency by overlapping load and compute
             if (EnableDiagnostics)
-                Console.WriteLine($"[DynamicGemm] SELECTED HIGH-OCCUPANCY kernel: {config.KernelName}");
+                Trace.WriteLine($"[DynamicGemm] SELECTED HIGH-OCCUPANCY kernel: {config.KernelName}");
             GenerateHighOccupancyDoubleBufferedKernel(sb, MWI, NWI, VWN, KWG);
         }
         else if (KREG > 1 && (VWN > 1 || VWM > 1))
         {
             // Use vectorized kernel WITH KREG for CLBlast-style performance
             if (EnableDiagnostics)
-                Console.WriteLine($"[DynamicGemm] SELECTED KREG kernel: {config.KernelName} KREG={KREG}");
+                Trace.WriteLine($"[DynamicGemm] SELECTED KREG kernel: {config.KernelName} KREG={KREG}");
             GenerateVectorizedKernelWithKreg(sb, MWI, NWI, VWM, VWN, KWI, KREG, useSubgroups);
         }
         else if (VWN > 1 || VWM > 1)
         {
             // Use vectorized kernel WITHOUT KREG (simpler, often faster!)
             if (EnableDiagnostics)
-                Console.WriteLine($"[DynamicGemm] SELECTED VECTORIZED kernel: {config.KernelName} VWM={VWM} VWN={VWN}");
+                Trace.WriteLine($"[DynamicGemm] SELECTED VECTORIZED kernel: {config.KernelName} VWM={VWM} VWN={VWN}");
             GenerateVectorizedKernel(sb, MWI, NWI, VWM, VWN, KWI);
         }
         else
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/GemmAutoTuner.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/GemmAutoTuner.cs
index d9794483e4..da50092165 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/GemmAutoTuner.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/GemmAutoTuner.cs
@@ -50,6 +50,34 @@ public readonly struct GemmConfig
     public bool UseTrueVectorLDS { get; init; }  // Use vectorized LDS arrays instead of scalar
     public bool UseColumnMajorA { get; init; }   // Treat A as column-major (packed/transpose path)
 
+    public GemmConfig WithColumnMajorA(bool useColumnMajorA)
+    {
+        return new GemmConfig
+        {
+            TileM = TileM,
+            TileN = TileN,
+            TileK = TileK,
+            ThreadTileM = ThreadTileM,
+            ThreadTileN = ThreadTileN,
+            VectorWidthM = VectorWidthM,
+            VectorWidthN = VectorWidthN,
+            UseDoubleBuffering = UseDoubleBuffering,
+            UseVectorizedLoads = UseVectorizedLoads,
+            KernelName = KernelName,
+            KReg = KReg,
+            KUnroll = KUnroll,
+            UseSubgroupOps = UseSubgroupOps,
+            StrideM = StrideM,
+            StrideN = StrideN,
+            CacheA = CacheA,
+            CacheB = CacheB,
+            MdimaSize = MdimaSize,
+            NdimbSize = NdimbSize,
+            UseTrueVectorLDS = UseTrueVectorLDS,
+            UseColumnMajorA = useColumnMajorA
+        };
+    }
+
     /// <summary>
     /// Generates a unique cache key for this configuration.
     /// Used by DynamicGemmKernel to cache compiled kernels.
@@ -153,7 +181,7 @@ private static void LogDiag(string message)
             catch
             {
                 if (EnableDiagnostics && Logger == null)
-                    Console.WriteLine(logLine);
+                    Trace.WriteLine(logLine);
             }
         }
         else if (EnableDiagnostics && Logger == null)
@@ -180,7 +208,7 @@ private static void WriteConsoleWithColor(string message)
     {
         if (Console.IsOutputRedirected)
         {
-            Console.WriteLine(message);
+            Trace.WriteLine(message);
             return;
         }
 
@@ -189,12 +217,12 @@ private static void WriteConsoleWithColor(string message)
         {
             var previous = Console.ForegroundColor;
             Console.ForegroundColor = color.Value;
-            Console.WriteLine(message);
+            Trace.WriteLine(message);
             Console.ForegroundColor = previous;
         }
         else
         {
-            Console.WriteLine(message);
+            Trace.WriteLine(message);
         }
     }
 
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClBackend.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClBackend.cs
index d5a9dc6300..d55e787381 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClBackend.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClBackend.cs
@@ -27,6 +27,23 @@ namespace AiDotNet.Tensors.Engines.DirectGpu.OpenCL
     /// <item>Bank-conflict-free shared memory</item>
     /// </list>
     /// </remarks>
+    /// <summary>
+    /// OpenCL backend for direct GPU access on AMD, Intel, and NVIDIA GPUs.
+    /// Uses pure P/Invoke with no managed GPU runtime dependency.
+    /// </summary>
+    /// <remarks>
+    /// <para><b>Key Features:</b></para>
+    /// <list type="bullet">
+    /// <item>Works on ALL .NET versions (4.6.2, 4.7.1, net8.0, etc.)</item>
+    /// <item>No managed GPU runtime dependency - pure P/Invoke</item>
+    /// <item>Double-buffered GEMM for compute/memory overlap</item>
+    /// <item>Fused operations (GEMM+Bias+Activation)</item>
+    /// <item>Bank-conflict-free shared memory</item>
+    /// </list>
+    /// <para><b>For Beginners:</b> This is the "driver" that talks directly to your graphics card (GPU). 
+    /// It translates math problems (like multiplying giant tables of numbers) into a language 
+    /// the GPU understands. This is much faster than using just your computer's main processor (CPU).</para>
+    /// </remarks>
     public sealed class OpenClBackend : IAsyncGpuBackend
     {
         private DirectOpenClContext? _context;
@@ -60,13 +77,36 @@ public sealed class OpenClBackend : IAsyncGpuBackend
         private ClBlastXgemmDirectParameters _clblastDirectParams;
         private int _clblastMinIndirectSize;
 
-        public bool IsAvailable { get; }
+        /// <inheritdoc/>
+        public bool IsAvailable { get; } = false;
+
+        /// <inheritdoc/>
         public string BackendName => "OpenCL";
-        public string DeviceName { get; }
-        public string DeviceVendor { get; }
-        public int ComputeUnits { get; }
-        public long GlobalMemoryBytes { get; }
-        public long LocalMemoryBytes { get; }
+
+        /// <summary>
+        /// Gets the name of the GPU device.
+        /// </summary>
+        public string DeviceName { get; } = "None";
+
+        /// <summary>
+        /// Gets the vendor of the GPU device (e.g., AMD, NVIDIA).
+        /// </summary>
+        public string DeviceVendor { get; } = "None";
+
+        /// <summary>
+        /// Gets the number of compute units available on the device.
+        /// </summary>
+        public int ComputeUnits { get; } = 0;
+
+        /// <summary>
+        /// Gets the total amount of global memory in bytes.
+        /// </summary>
+        public long GlobalMemoryBytes { get; } = 0;
+
+        /// <summary>
+        /// Gets the amount of local (shared) memory in bytes.
+        /// </summary>
+        public long LocalMemoryBytes { get; } = 0;
 
         // IAsyncGpuBackend properties
         public bool SupportsMultiStream => true;
@@ -88,10 +128,24 @@ public sealed class OpenClBackend : IAsyncGpuBackend
         /// </summary>
         public static bool IsOpenClAvailable => DirectOpenClContext.IsAvailable;
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="OpenClBackend"/> class.
+        /// </summary>
+        /// <param name="logger">Optional logger for diagnostic output.</param>
         public OpenClBackend(ILogger? logger = null) : this(0, logger)
         {
         }
 
+        /// <summary>
+        /// Initializes a new instance of the <see cref="OpenClBackend"/> class for a specific device.
+        /// </summary>
+        /// <param name="deviceIndex">The index of the OpenCL device to use.</param>
+        /// <param name="logger">Optional logger for diagnostic output.</param>
+        /// <remarks>
+        /// <para><b>For Beginners:</b> This sets up the connection between your code and the GPU. 
+        /// It checks if the GPU is working, finds out how much memory it has, and prepares 
+        /// the "toolset" (kernels) the model will use for math.</para>
+        /// </remarks>
         public OpenClBackend(int deviceIndex, ILogger? logger = null)
         {
             _logger = logger;
@@ -109,9 +163,9 @@ public OpenClBackend(int deviceIndex, ILogger? logger = null)
 
             try
             {
-                Console.WriteLine($"[OpenClBackend] Creating DirectOpenClContext for device {deviceIndex}...");
+                Trace.WriteLine($"[OpenClBackend] Creating DirectOpenClContext for device {deviceIndex}...");
                 _context = new DirectOpenClContext(deviceIndex);
-                Console.WriteLine($"[OpenClBackend] Context created: Device={_context.DeviceName}, Vendor={_context.DeviceVendor}");
+                Trace.WriteLine($"[OpenClBackend] Context created: Device={_context.DeviceName}, Vendor={_context.DeviceVendor}");
 
                 IsAvailable = true;
                 DeviceName = _context.DeviceName;
@@ -122,7 +176,7 @@ public OpenClBackend(int deviceIndex, ILogger? logger = null)
                 string? envCUs = Environment.GetEnvironmentVariable("AIDOTNET_GPU_COMPUTE_UNITS");
                 if (int.TryParse(envCUs, out int overrideCUs) && overrideCUs > 0 && overrideCUs <= 256)
                 {
-                    Console.WriteLine($"[OpenClBackend] CU override: {detectedCUs} -> {overrideCUs} (via AIDOTNET_GPU_COMPUTE_UNITS)");
+                    Trace.WriteLine($"[OpenClBackend] CU override: {detectedCUs} -> {overrideCUs} (via AIDOTNET_GPU_COMPUTE_UNITS)");
                     ComputeUnits = overrideCUs;
                 }
                 else
@@ -140,35 +194,35 @@ public OpenClBackend(int deviceIndex, ILogger? logger = null)
                 _supportsSubgroups = _context.SupportsSubgroups;
 
                 // Print GPU capabilities for diagnostics
-                Console.WriteLine($"[OpenClBackend] GPU Capabilities:");
-                Console.WriteLine($"[OpenClBackend]   Compute Units: {ComputeUnits}");
-                Console.WriteLine($"[OpenClBackend]   Max Work Group Size: {_maxWorkGroupSize}");
+                Trace.WriteLine($"[OpenClBackend] GPU Capabilities:");
+                Trace.WriteLine($"[OpenClBackend]   Compute Units: {ComputeUnits}");
+                Trace.WriteLine($"[OpenClBackend]   Max Work Group Size: {_maxWorkGroupSize}");
                 if (_maxWorkItemSizes.Length >= 2)
                 {
-                    Console.WriteLine($"[OpenClBackend]   Max Work Item Sizes: [{string.Join(", ", _maxWorkItemSizes)}]");
+                    Trace.WriteLine($"[OpenClBackend]   Max Work Item Sizes: [{string.Join(", ", _maxWorkItemSizes)}]");
                 }
-                Console.WriteLine($"[OpenClBackend]   Local Memory: {LocalMemoryBytes / 1024} KB");
-                Console.WriteLine($"[OpenClBackend]   Supports FP16: {_supportsFp16}");
-                Console.WriteLine($"[OpenClBackend]   Supports Subgroups: {_supportsSubgroups}");
+                Trace.WriteLine($"[OpenClBackend]   Local Memory: {LocalMemoryBytes / 1024} KB");
+                Trace.WriteLine($"[OpenClBackend]   Supports FP16: {_supportsFp16}");
+                Trace.WriteLine($"[OpenClBackend]   Supports Subgroups: {_supportsSubgroups}");
 
                 // Initialize default stream wrapper
                 _defaultStream = new OpenClCommandQueue(this, _context.CommandQueue, _context.Context, _context.Device,
                     GpuStreamType.Default, _context.IsProfilingEnabled, ownsHandle: false);
-                Console.WriteLine("[OpenClBackend] Default command queue wrapper initialized.");
+                Trace.WriteLine("[OpenClBackend] Default command queue wrapper initialized.");
 
-                Console.WriteLine("[OpenClBackend] Compiling kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling kernels...");
                 CompileKernels();
-                Console.WriteLine($"[OpenClBackend] Kernels compiled successfully. Total: {_kernelCache.Count}");
+                Trace.WriteLine($"[OpenClBackend] Kernels compiled successfully. Total: {_kernelCache.Count}");
 
                 // Initialize dynamic kernel generator for Bayesian-optimized GEMM
                 _dynamicGemm = new DynamicGemmKernel(_context);
-                Console.WriteLine("[OpenClBackend] Dynamic GEMM kernel generator initialized.");
+                Trace.WriteLine("[OpenClBackend] Dynamic GEMM kernel generator initialized.");
             }
             catch (Exception ex)
             {
-                Console.WriteLine($"[OpenClBackend] Initialization FAILED: {ex.GetType().Name}: {ex.Message}");
+                Trace.WriteLine($"[OpenClBackend] Initialization FAILED: {ex.GetType().Name}: {ex.Message}");
                 if (ex.InnerException != null)
-                    Console.WriteLine($"[OpenClBackend] Inner: {ex.InnerException.Message}");
+                    Trace.WriteLine($"[OpenClBackend] Inner: {ex.InnerException.Message}");
                 System.Diagnostics.Debug.WriteLine($"OpenClBackend initialization failed: {ex.Message}");
                 IsAvailable = false;
                 DeviceName = "None";
@@ -193,7 +247,7 @@ private void CompileKernels()
             try
             {
                 // Compile GEMM kernels with aggressive optimizations
-                Console.WriteLine("[OpenClBackend] Compiling GEMM kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling GEMM kernels...");
                 var gemmProgram = new DirectOpenClProgram(_context, GemmKernel.GetSource());
                 gemmProgram.Build(optimizationFlags);
                 _programs.Add(gemmProgram);
@@ -201,10 +255,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, gemmProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] GEMM kernels: {string.Join(", ", GemmKernel.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] GEMM kernels: {string.Join(", ", GemmKernel.GetKernelNames())}");
 
                 // Compile activation kernels
-                Console.WriteLine("[OpenClBackend] Compiling activation kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling activation kernels...");
                 var activationProgram = new DirectOpenClProgram(_context, ActivationKernels.GetSource());
                 activationProgram.Build(optimizationFlags);
                 _programs.Add(activationProgram);
@@ -212,10 +266,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, activationProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Activation kernels: {string.Join(", ", ActivationKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Activation kernels: {string.Join(", ", ActivationKernels.GetKernelNames())}");
 
                 // Compile fused kernels
-                Console.WriteLine("[OpenClBackend] Compiling fused kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling fused kernels...");
                 var fusedProgram = new DirectOpenClProgram(_context, FusedKernels.GetSource());
                 fusedProgram.Build(optimizationFlags);
                 _programs.Add(fusedProgram);
@@ -223,10 +277,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, fusedProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Fused kernels: {string.Join(", ", FusedKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Fused kernels: {string.Join(", ", FusedKernels.GetKernelNames())}");
 
                 // Compile reduction kernels
-                Console.WriteLine("[OpenClBackend] Compiling reduction kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling reduction kernels...");
                 var reductionProgram = new DirectOpenClProgram(_context, ReductionKernels.GetSource());
                 reductionProgram.Build(optimizationFlags);
                 _programs.Add(reductionProgram);
@@ -234,10 +288,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, reductionProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Reduction kernels: {string.Join(", ", ReductionKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Reduction kernels: {string.Join(", ", ReductionKernels.GetKernelNames())}");
 
                 // Compile packing/pad-copy kernels
-                Console.WriteLine("[OpenClBackend] Compiling packing kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling packing kernels...");
                 var packingProgram = new DirectOpenClProgram(_context, PackingKernels.GetSource());
                 packingProgram.Build(optimizationFlags);
                 _programs.Add(packingProgram);
@@ -245,10 +299,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, packingProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Packing kernels: {string.Join(", ", PackingKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Packing kernels: {string.Join(", ", PackingKernels.GetKernelNames())}");
 
                 // Compile sparse GEMM kernels (2:4 structured sparsity)
-                Console.WriteLine("[OpenClBackend] Compiling sparse GEMM kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling sparse GEMM kernels...");
                 var sparseProgram = new DirectOpenClProgram(_context, SparseGemmKernels.GetSource());
                 sparseProgram.Build(optimizationFlags);
                 _programs.Add(sparseProgram);
@@ -256,10 +310,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, sparseProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Sparse GEMM kernels: {string.Join(", ", SparseGemmKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Sparse GEMM kernels: {string.Join(", ", SparseGemmKernels.GetKernelNames())}");
 
                 // Compile CSR sparse kernels (general sparsity for GNN)
-                Console.WriteLine("[OpenClBackend] Compiling CSR sparse kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling CSR sparse kernels...");
                 var csrSparseProgram = new DirectOpenClProgram(_context, CsrSparseKernels.GetSource());
                 csrSparseProgram.Build(optimizationFlags);
                 _programs.Add(csrSparseProgram);
@@ -267,10 +321,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, csrSparseProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] CSR sparse kernels: {string.Join(", ", CsrSparseKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] CSR sparse kernels: {string.Join(", ", CsrSparseKernels.GetKernelNames())}");
 
                 // Compile convolution kernels
-                Console.WriteLine("[OpenClBackend] Compiling convolution kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling convolution kernels...");
                 var convProgram = new DirectOpenClProgram(_context, ConvolutionKernels.GetSource());
                 convProgram.Build(optimizationFlags);
                 _programs.Add(convProgram);
@@ -278,10 +332,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, convProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Convolution kernels: {string.Join(", ", ConvolutionKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Convolution kernels: {string.Join(", ", ConvolutionKernels.GetKernelNames())}");
 
                 // Compile fused convolution kernels
-                Console.WriteLine("[OpenClBackend] Compiling fused convolution kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling fused convolution kernels...");
                 var fusedConvProgram = new DirectOpenClProgram(_context, FusedConvolutionKernels.GetSource());
                 fusedConvProgram.Build(optimizationFlags);
                 _programs.Add(fusedConvProgram);
@@ -289,10 +343,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, fusedConvProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Fused convolution kernels: {string.Join(", ", FusedConvolutionKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Fused convolution kernels: {string.Join(", ", FusedConvolutionKernels.GetKernelNames())}");
 
                 // Compile pooling kernels
-                Console.WriteLine("[OpenClBackend] Compiling pooling kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling pooling kernels...");
                 var poolProgram = new DirectOpenClProgram(_context, PoolingKernels.GetSource());
                 poolProgram.Build(optimizationFlags);
                 _programs.Add(poolProgram);
@@ -300,10 +354,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, poolProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Pooling kernels: {string.Join(", ", PoolingKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Pooling kernels: {string.Join(", ", PoolingKernels.GetKernelNames())}");
 
                 // Compile normalization kernels
-                Console.WriteLine("[OpenClBackend] Compiling normalization kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling normalization kernels...");
                 var normProgram = new DirectOpenClProgram(_context, NormalizationKernels.GetSource());
                 normProgram.Build(optimizationFlags);
                 _programs.Add(normProgram);
@@ -311,10 +365,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, normProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Normalization kernels: {string.Join(", ", NormalizationKernels.GetKernelNames())}");
+                Trace.WriteLine($"[OpenClBackend] Normalization kernels: {string.Join(", ", NormalizationKernels.GetKernelNames())}");
 
                 // Compile neural network kernels (activation gradients, loss, optimizers)
-                Console.WriteLine("[OpenClBackend] Compiling neural network kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling neural network kernels...");
                 var nnProgram = new DirectOpenClProgram(_context, NeuralNetKernels.GetSource());
                 nnProgram.Build(optimizationFlags);
                 _programs.Add(nnProgram);
@@ -322,7 +376,7 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, nnProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Neural network kernels compiled: {NeuralNetKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Neural network kernels compiled: {NeuralNetKernels.GetKernelNames().Length} kernels");
 
                 // Compile mixed precision kernels only if device actually supports FP16
                 // This is non-fatal - if compilation fails, we fall back to no FP16 support
@@ -331,7 +385,7 @@ private void CompileKernels()
                 {
                     try
                     {
-                        Console.WriteLine("[OpenClBackend] Compiling mixed precision kernels...");
+                        Trace.WriteLine("[OpenClBackend] Compiling mixed precision kernels...");
                         var mixedPrecisionSource = string.Join("\n\n",
                             MixedPrecisionKernels.ConvertFp32ToFp16,
                             MixedPrecisionKernels.ConvertFp16ToFp32,
@@ -347,23 +401,23 @@ private void CompileKernels()
                         _kernelCache["mixed_precision_backward"] = new DirectOpenClKernel(_context, mpProgram, "mixed_precision_backward");
                         _kernelCache["accumulate_gradient_fp32"] = new DirectOpenClKernel(_context, mpProgram, "accumulate_gradient_fp32");
                         _mixedPrecisionKernelsAvailable = true;
-                        Console.WriteLine("[OpenClBackend] Mixed precision kernels compiled: 5 kernels");
+                        Trace.WriteLine("[OpenClBackend] Mixed precision kernels compiled: 5 kernels");
                     }
                     catch (Exception ex)
                     {
                         // Mixed precision compilation failed - this is non-fatal
                         // Device may report FP16 support but have driver issues with these kernels
-                        Console.WriteLine($"[OpenClBackend] Warning: Mixed precision kernel compilation failed (non-fatal): {ex.Message}");
-                        Console.WriteLine("[OpenClBackend] Continuing without mixed precision support.");
+                        Trace.WriteLine($"[OpenClBackend] Warning: Mixed precision kernel compilation failed (non-fatal): {ex.Message}");
+                        Trace.WriteLine("[OpenClBackend] Continuing without mixed precision support.");
                     }
                 }
                 else
                 {
-                    Console.WriteLine("[OpenClBackend] Skipping mixed precision kernels (FP16 not supported on this device).");
+                    Trace.WriteLine("[OpenClBackend] Skipping mixed precision kernels (FP16 not supported on this device).");
                 }
 
                 // Compile attention kernels (FlashAttention, GQA, ScaledDotProduct)
-                Console.WriteLine("[OpenClBackend] Compiling attention kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling attention kernels...");
                 var attnProgram = new DirectOpenClProgram(_context, AttentionKernels.GetSource());
                 attnProgram.Build(optimizationFlags);
                 _programs.Add(attnProgram);
@@ -371,10 +425,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, attnProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Attention kernels compiled: {AttentionKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Attention kernels compiled: {AttentionKernels.GetKernelNames().Length} kernels");
 
                 // Compile FFT kernels
-                Console.WriteLine("[OpenClBackend] Compiling FFT kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling FFT kernels...");
                 var fftProgram = new DirectOpenClProgram(_context, FFTKernels.GetSource());
                 fftProgram.Build(optimizationFlags);
                 _programs.Add(fftProgram);
@@ -382,10 +436,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, fftProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] FFT kernels compiled: {FFTKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] FFT kernels compiled: {FFTKernels.GetKernelNames().Length} kernels");
 
                 // Compile spatial transformer kernels (TopK, AffineGrid, GridSample)
-                Console.WriteLine("[OpenClBackend] Compiling spatial transformer kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling spatial transformer kernels...");
                 var stProgram = new DirectOpenClProgram(_context, SpatialTransformerKernels.GetSource());
                 stProgram.Build(optimizationFlags);
                 _programs.Add(stProgram);
@@ -393,10 +447,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, stProgram, name);
                 }
-                Console.WriteLine("[OpenClBackend] Spatial transformer kernels compiled: 4 kernels");
+                Trace.WriteLine("[OpenClBackend] Spatial transformer kernels compiled: 4 kernels");
 
                 // Compile locally connected convolution kernels
-                Console.WriteLine("[OpenClBackend] Compiling locally connected kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling locally connected kernels...");
                 var locallyConnectedProgram = new DirectOpenClProgram(_context, LocallyConnectedKernels.GetSource());
                 locallyConnectedProgram.Build(optimizationFlags);
                 _programs.Add(locallyConnectedProgram);
@@ -404,10 +458,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, locallyConnectedProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Locally connected kernels compiled: {LocallyConnectedKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Locally connected kernels compiled: {LocallyConnectedKernels.GetKernelNames().Length} kernels");
 
                 // Compile deformable convolution kernels
-                Console.WriteLine("[OpenClBackend] Compiling deformable convolution kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling deformable convolution kernels...");
                 var deformableProgram = new DirectOpenClProgram(_context, DeformableConvolutionKernels.GetSource());
                 deformableProgram.Build(optimizationFlags);
                 _programs.Add(deformableProgram);
@@ -415,10 +469,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, deformableProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Deformable convolution kernels compiled: {DeformableConvolutionKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Deformable convolution kernels compiled: {DeformableConvolutionKernels.GetKernelNames().Length} kernels");
 
                 // Compile random number generation kernels
-                Console.WriteLine("[OpenClBackend] Compiling random number kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling random number kernels...");
                 var randomProgram = new DirectOpenClProgram(_context, RandomKernels.GetKernels());
                 randomProgram.Build(optimizationFlags);
                 _programs.Add(randomProgram);
@@ -426,10 +480,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, randomProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Random number kernels compiled: 2 kernels");
+                Trace.WriteLine($"[OpenClBackend] Random number kernels compiled: 2 kernels");
 
                 // Compile specialized kernels (hyperbolic geometry, octonion algebra, quantum computing)
-                Console.WriteLine("[OpenClBackend] Compiling specialized kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling specialized kernels...");
                 var specializedProgram = new DirectOpenClProgram(_context, SpecializedKernels.GetSource());
                 specializedProgram.Build(optimizationFlags);
                 _programs.Add(specializedProgram);
@@ -437,10 +491,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, specializedProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Specialized kernels compiled: {SpecializedKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Specialized kernels compiled: {SpecializedKernels.GetKernelNames().Length} kernels");
 
                 // Compile FP16 conversion kernels (half-precision float conversion)
-                Console.WriteLine("[OpenClBackend] Compiling FP16 conversion kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling FP16 conversion kernels...");
                 var fp16Program = new DirectOpenClProgram(_context, Fp16Kernels.GetSource());
                 fp16Program.Build(optimizationFlags);
                 _programs.Add(fp16Program);
@@ -448,10 +502,10 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, fp16Program, name);
                 }
-                Console.WriteLine($"[OpenClBackend] FP16 conversion kernels compiled: {Fp16Kernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] FP16 conversion kernels compiled: {Fp16Kernels.GetKernelNames().Length} kernels");
 
                 // Compile loss function kernels (MSE, BCE, CE, Huber, Focal, Triplet, etc.)
-                Console.WriteLine("[OpenClBackend] Compiling loss function kernels...");
+                Trace.WriteLine("[OpenClBackend] Compiling loss function kernels...");
                 var lossProgram = new DirectOpenClProgram(_context, LossKernels.GetSource());
                 lossProgram.Build(optimizationFlags);
                 _programs.Add(lossProgram);
@@ -459,7 +513,7 @@ private void CompileKernels()
                 {
                     _kernelCache[name] = new DirectOpenClKernel(_context, lossProgram, name);
                 }
-                Console.WriteLine($"[OpenClBackend] Loss function kernels compiled: {LossKernels.GetKernelNames().Length} kernels");
+                Trace.WriteLine($"[OpenClBackend] Loss function kernels compiled: {LossKernels.GetKernelNames().Length} kernels");
             }
             catch (Exception ex)
             {
@@ -654,12 +708,12 @@ public void Log<TState>(LogLevel logLevel, EventId eventId, TState state, Except
                 {
                     var previous = Console.ForegroundColor;
                     Console.ForegroundColor = color.Value;
-                    Console.WriteLine(message);
+                    Trace.WriteLine(message);
                     Console.ForegroundColor = previous;
                 }
                 else
                 {
-                    Console.WriteLine(message);
+                    Trace.WriteLine(message);
                 }
             }
 
@@ -797,7 +851,7 @@ private bool TryGetTunedConfig(int M, int N, int K, out GemmConfig config)
 
                     if (EnableTuningDiagnostics)
                     {
-                        Console.WriteLine($"[GEMM] Cached config invalid: {validationError}");
+                        Trace.WriteLine($"[GEMM] Cached config invalid: {validationError}");
                     }
                 }
 
@@ -862,7 +916,7 @@ private bool TryGetTunedConfig(int M, int N, int K, out GemmConfig config)
             {
                 if (EnableTuningDiagnostics)
                 {
-                    Console.WriteLine($"[GEMM] Tuning lookup failed: {ex.Message}");
+                    Trace.WriteLine($"[GEMM] Tuning lookup failed: {ex.Message}");
                 }
 
                 return false;
@@ -896,10 +950,11 @@ private bool TryGetClBlastBaselineConfig(out GemmConfig config)
                 _clblastTransposeParams = ClBlastTransposeDatabase.GetParameters(deviceInfo);
                 _clblastDirectParams = ClBlastXgemmDirectDatabase.GetParameters(deviceInfo);
                 _clblastMinIndirectSize = ClBlastGemmRoutineDatabase.GetXgemmMinIndirectSize(deviceInfo);
-                Console.WriteLine($"[OpenClBackend] CLBlast MinIndirectSize threshold: {_clblastMinIndirectSize} (use INDIRECT for M/N >= {_clblastMinIndirectSize})");
+                Trace.WriteLine($"[OpenClBackend] CLBlast MinIndirectSize threshold: {_clblastMinIndirectSize} (use INDIRECT for M/N >= {_clblastMinIndirectSize})");
 
                 if (ClBlastXgemmDatabase.TryGetConfig(deviceInfo, out var baseline))
                 {
+                    baseline = NormalizeRowMajorConfig(baseline);
                     _clblastBaselineConfig = baseline;
                     config = baseline;
                     return true;
@@ -929,15 +984,24 @@ private bool TryGetClBlastBaselineConfig(out GemmConfig config)
                     MdimaSize = 16,
                     NdimbSize = 8,
                     UseTrueVectorLDS = true,
-                    UseColumnMajorA = true,
+                    UseColumnMajorA = false,
                     KernelName = "clblast_baseline_k0"
                 };
+                defaultBaseline = NormalizeRowMajorConfig(defaultBaseline);
                 _clblastBaselineConfig = defaultBaseline;
                 config = defaultBaseline;
                 return true;
             }
         }
 
+        private static GemmConfig NormalizeRowMajorConfig(GemmConfig config)
+        {
+            if (!config.UseColumnMajorA)
+                return config;
+
+            return config.WithColumnMajorA(false);
+        }
+
         private bool EnsureClBlastPackingKernels()
         {
             if (_context == null)
@@ -1258,19 +1322,19 @@ private bool TryExecuteClBlastBaselineGemm(
             if (!useIndirectPath || forceDirect)
             {
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying DIRECT path (M/N < {_clblastMinIndirectSize} or forceDirect={forceDirect})");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying DIRECT path (M/N < {_clblastMinIndirectSize} or forceDirect={forceDirect})");
                 if (TryExecuteClBlastDirectGemm(A, B, C, M, N, K, alpha, beta))
                 {
                     if (traceEnabled)
-                        Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: DIRECT path executed");
+                        Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: DIRECT path executed");
                     return true;
                 }
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] DIRECT path failed, trying INDIRECT");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] DIRECT path failed, trying INDIRECT");
             }
             else if (traceEnabled)
             {
-                Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Skipping DIRECT path (M/N >= {_clblastMinIndirectSize}), using INDIRECT");
+                Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Skipping DIRECT path (M/N >= {_clblastMinIndirectSize}), using INDIRECT");
             }
 
             if (_dynamicGemm == null)
@@ -1288,7 +1352,7 @@ private bool TryExecuteClBlastBaselineGemm(
             if (config.TileM <= 0 || config.TileN <= 0 || kUnit <= 0)
                 return false;
 
-            bool kernel1 = IsClBlastBaselineKernel1(config);
+            bool kernel1 = IsClBlastBaselineKernel(config);
 
             // ROW-MAJOR SWAP TRICK (Zero-Copy Optimization):
             // For row-major GEMM C = A × B where A(M×K), B(K×N), C(M×N):
@@ -1374,13 +1438,13 @@ private bool TryExecuteClBlastBaselineGemm(
                         if (timingEnabled) { Synchronize(); packBTime = sw!.ElapsedTicks; }
                     }
 
-                    // Pad C if needed (for beta != 0) - NO TRANSPOSE
+                    // Pad C if needed (for non-zero beta) - NO TRANSPOSE
                     if (cNeedsPad)
                     {
                         if (timingEnabled) { sw!.Restart(); }
                         cTemp = AllocateBuffer((int)cSize);
                         if (timingEnabled) { allocTime += sw!.ElapsedTicks; sw.Restart(); }
-                        if (beta != 0.0f)
+                        if (!IsEffectivelyZero(beta))
                         {
                             // C is M×N row-major. Reinterpreted as column-major, it's N×M.
                             // We need mCeiled×nCeiled.
@@ -1412,7 +1476,7 @@ private bool TryExecuteClBlastBaselineGemm(
                         var total = allocTime + packATime + packBTime + packCTime + gemmTime + unpackCTime;
                         double flops = 2.0 * M * N * K;
                         double gflops = flops / (total / ticksPerMs) / 1e6;
-                        Console.WriteLine($"[TIMING-SWAP {M}x{N}x{K}] Alloc={allocTime / ticksPerMs:F2}ms PackA={packATime / ticksPerMs:F2}ms PackB={packBTime / ticksPerMs:F2}ms GEMM={gemmTime / ticksPerMs:F2}ms UnpackC={unpackCTime / ticksPerMs:F2}ms Total={total / ticksPerMs:F2}ms ({gflops:F0} GFLOPS)");
+                        Trace.WriteLine($"[TIMING-SWAP {M}x{N}x{K}] Alloc={allocTime / ticksPerMs:F2}ms PackA={packATime / ticksPerMs:F2}ms PackB={packBTime / ticksPerMs:F2}ms GEMM={gemmTime / ticksPerMs:F2}ms UnpackC={unpackCTime / ticksPerMs:F2}ms Total={total / ticksPerMs:F2}ms ({gflops:F0} GFLOPS)");
                     }
 
                     return true;
@@ -1490,7 +1554,7 @@ private bool TryExecuteClBlastBaselineGemm(
                         cTemp = AllocateBuffer((int)cSize);
                         cBuf = cTemp;
                         if (timingEnabled) { allocTime += sw!.ElapsedTicks; sw.Restart(); }
-                        if (beta != 0.0f)
+                        if (!IsEffectivelyZero(beta))
                         {
                             ClBlastCopyMatrix(C, cBuf, cOne, cTwo, cOne, 0, cOneI, cTwoI, cOneI, 0, true);
                             if (timingEnabled) { Synchronize(); packCTime = sw!.ElapsedTicks; }
@@ -1513,7 +1577,7 @@ private bool TryExecuteClBlastBaselineGemm(
                     {
                         double ticksPerMs = System.Diagnostics.Stopwatch.Frequency / 1000.0;
                         var total = allocTime + packATime + packBTime + packCTime + gemmTime + unpackCTime;
-                        Console.WriteLine($"[TIMING {M}x{N}x{K}] Alloc={allocTime / ticksPerMs:F2}ms PackA={packATime / ticksPerMs:F2}ms PackB={packBTime / ticksPerMs:F2}ms GEMM={gemmTime / ticksPerMs:F2}ms UnpackC={unpackCTime / ticksPerMs:F2}ms Total={total / ticksPerMs:F2}ms");
+                        Trace.WriteLine($"[TIMING {M}x{N}x{K}] Alloc={allocTime / ticksPerMs:F2}ms PackA={packATime / ticksPerMs:F2}ms PackB={packBTime / ticksPerMs:F2}ms GEMM={gemmTime / ticksPerMs:F2}ms UnpackC={unpackCTime / ticksPerMs:F2}ms Total={total / ticksPerMs:F2}ms");
                     }
 
                     return true;
@@ -1588,71 +1652,6 @@ private void CopySubmatrix(IGpuBuffer src, IGpuBuffer dst, int rows, int cols, i
             kernel.Execute2D(rows, cols, localSizeX, localSizeY);
         }
 
-        private bool TryExecutePackedDynamicGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K, float alpha, float beta, GemmConfig config)
-        {
-            if (IsClBlastBaselineKernel(config))
-                return TryExecuteClBlastBaselineGemm(A, B, C, M, N, K, alpha, beta, config);
-
-            int kReg = config.KReg > 0 ? config.KReg : 1;
-            int kUnit = config.TileK * kReg;
-            if (config.TileM <= 0 || config.TileN <= 0 || kUnit <= 0)
-                return false;
-
-            int mPad = CeilDiv(M, config.TileM) * config.TileM;
-            int nPad = CeilDiv(N, config.TileN) * config.TileN;
-            int kPad = CeilDiv(K, kUnit) * kUnit;
-
-            bool useColumnMajorA = config.UseColumnMajorA || IsClBlastBaselineKernel0(config);
-            bool useColumnMajorC = IsClBlastBaselineKernel0(config);
-            bool needsPadding = mPad != M || nPad != N || kPad != K;
-            if (!needsPadding && !useColumnMajorA && !useColumnMajorC)
-                return TryExecuteDynamicGemm(A, B, C, M, N, K, alpha, beta, config);
-
-            long aSize = (long)mPad * kPad;
-            long bSize = (long)kPad * nPad;
-            long cSize = (long)mPad * nPad;
-            if (aSize > int.MaxValue || bSize > int.MaxValue || cSize > int.MaxValue)
-                return false;
-
-            if (EnableTuningDiagnostics)
-            {
-                Console.WriteLine($"[GEMM] Packed GEMM: {M}x{N}x{K} -> {mPad}x{nPad}x{kPad}");
-            }
-
-            using var aPad = AllocateBuffer(mPad * kPad);
-            using var bPad = AllocateBuffer(kPad * nPad);
-            using var cPad = AllocateBuffer(mPad * nPad);
-
-            if (useColumnMajorA)
-                PadCopyTransposeMatrix(A, aPad, M, K, mPad, kPad);
-            else
-                PadCopyMatrix(A, aPad, M, K, mPad, kPad);
-            PadCopyMatrix(B, bPad, K, N, kPad, nPad);
-            if (useColumnMajorC)
-            {
-                if (beta != 0.0f)
-                    PadCopyTransposeMatrix(C, cPad, M, N, mPad, nPad);
-                else
-                    PadCopyMatrix(C, cPad, 0, 0, mPad, nPad);
-            }
-            else
-            {
-                if (beta != 0.0f)
-                    PadCopyMatrix(C, cPad, M, N, mPad, nPad);
-                else
-                    PadCopyMatrix(C, cPad, 0, 0, mPad, nPad);
-            }
-
-            if (!TryExecuteDynamicGemm(aPad, bPad, cPad, mPad, nPad, kPad, alpha, beta, config))
-                return false;
-
-            if (useColumnMajorC)
-                PadCopyFromColumnMajorMatrix(cPad, C, mPad, nPad, M, N);
-            else
-                CopySubmatrix(cPad, C, M, N, nPad, N);
-            return true;
-        }
-
         private static bool IsClBlastBaselineKernel0(GemmConfig config)
         {
             return !string.IsNullOrWhiteSpace(config.KernelName) &&
@@ -1670,6 +1669,11 @@ private static bool IsClBlastBaselineKernel(GemmConfig config)
             return IsClBlastBaselineKernel0(config) || IsClBlastBaselineKernel1(config);
         }
 
+        private static bool IsEffectivelyZero(float value)
+        {
+            return !float.IsNaN(value) && !(value > 0.0f || value < 0.0f);
+        }
+
         private bool TryExecuteDynamicGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K, float alpha, float beta, GemmConfig config)
         {
             if (_dynamicGemm == null)
@@ -1680,7 +1684,7 @@ private bool TryExecuteDynamicGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int
             {
                 if (EnableTuningDiagnostics)
                 {
-                    Console.WriteLine($"[GEMM] Dynamic config invalid: {validationError}");
+                    Trace.WriteLine($"[GEMM] Dynamic config invalid: {validationError}");
                 }
                 return false;
             }
@@ -1696,97 +1700,104 @@ private bool TryExecuteDynamicGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int
             {
                 if (EnableTuningDiagnostics)
                 {
-                    Console.WriteLine($"[GEMM] Dynamic kernel failed ({config.KernelName}): {ex.Message}");
+                    Trace.WriteLine($"[GEMM] Dynamic kernel failed ({config.KernelName}): {ex.Message}");
                 }
                 return false;
             }
         }
 
+        /// <summary>
+        /// Executes a General Matrix Multiplication (GEMM) operation: C = alpha * A * B + beta * C.
+        /// </summary>
+        /// <param name="A">The first input matrix.</param>
+        /// <param name="B">The second input matrix.</param>
+        /// <param name="C">The output matrix (also used as input if beta != 0).</param>
+        /// <param name="M">The number of rows in matrix A and C.</param>
+        /// <param name="N">The number of columns in matrix B and C.</param>
+        /// <param name="K">The number of columns in A and rows in B.</param>
+        /// <param name="alpha">Scaling factor for the product A * B.</param>
+        /// <param name="beta">Scaling factor for matrix C.</param>
+        /// <remarks>
+        /// <para><b>For Beginners:</b> This is the workhorse of AI. It multiplies two grids of numbers together. 
+        /// This method automatically picks the best way to do this based on your GPU's speed and the size 
+        /// of the problem. It prioritizes "dynamic" kernels which are custom-built for your specific hardware.</para>
+        /// </remarks>
         public void Gemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K, float alpha = 1.0f, float beta = 0.0f)
         {
             if (_context == null)
                 throw new InvalidOperationException("OpenCL context not available");
 
-            var bufferA = ((DirectOpenClGpuBuffer)A).Buffer;
-            var bufferB = ((DirectOpenClGpuBuffer)B).Buffer;
-            var bufferC = ((DirectOpenClGpuBuffer)C).Buffer;
-
             bool offlineEnabled = TryGetOfflineTuningMode(out _);
 
-            // DIAGNOSTIC TRACING (always prints for debugging)
-            bool traceEnabled = Environment.GetEnvironmentVariable("AIDOTNET_GEMM_TRACE") == "1";
+            // DIAGNOSTIC TRACING (uses Trace for library-level logging)
+            bool traceEnabled = GetEnvBool("AIDOTNET_GEMM_TRACE");
+            bool enableDynamic = GetEnvBool("AIDOTNET_GEMM_ENABLE_DYNAMIC", defaultValue: false);
 
-            if (!offlineEnabled && TryGetClBlastBaselineConfig(out var baselineConfig))
+            if (enableDynamic && !offlineEnabled && TryGetClBlastBaselineConfig(out var baselineConfig))
             {
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying CLBlast baseline (TileM={baselineConfig.TileM}, TileN={baselineConfig.TileN}, TileK={baselineConfig.TileK})");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying baseline config (TileM={baselineConfig.TileM}, TileN={baselineConfig.TileN}, TileK={baselineConfig.TileK})");
 
-                if (TryExecuteClBlastBaselineGemm(A, B, C, M, N, K, alpha, beta, baselineConfig))
+                if (TryExecutePackedDynamicGemm(A, B, C, M, N, K, alpha, beta, baselineConfig))
                 {
                     if (traceEnabled)
-                        Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: CLBlast baseline executed");
+                        Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: Baseline config executed");
                     return;
                 }
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: CLBlast baseline FAILED");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: Baseline config FAILED");
             }
-            else
+            else if (enableDynamic && traceEnabled)
             {
-                if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SKIP: CLBlast baseline not available (offline={offlineEnabled})");
+                Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SKIP: CLBlast baseline not available (offline={offlineEnabled})");
             }
 
-            if (_dynamicGemm != null && M >= 128 && N >= 128 && K >= 64 &&
+            if (enableDynamic && _dynamicGemm != null && M >= 128 && N >= 128 && K >= 64 &&
                 TryGetTunedConfig(M, N, K, out var tunedConfig))
             {
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying dynamic GEMM");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] Trying dynamic GEMM");
                 if (TryExecutePackedDynamicGemm(A, B, C, M, N, K, alpha, beta, tunedConfig))
                 {
                     if (traceEnabled)
-                        Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: Dynamic GEMM executed");
+                        Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] SUCCESS: Dynamic GEMM executed");
                     return;
                 }
                 if (traceEnabled)
-                    Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: Dynamic GEMM FAILED");
+                    Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: Dynamic GEMM FAILED");
             }
 
-            // FALLBACK: Using our own kernels (NOT CLBlast identical!)
+            // FALLBACK: Using our own kernels (tried when dynamic kernels are disabled or fail)
             if (traceEnabled)
-                Console.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: Using built-in kernel (NOT CLBlast!)");
+                Trace.WriteLine($"[GEMM-TRACE {M}x{N}x{K}] FALLBACK: Using built-in kernels");
 
             // Choose kernel based on matrix size
-            // Use optimized kernel for matrices >= 128 in any dimension
             if (M >= 128 && N >= 128 && K >= 64)
             {
-                // Large matrix - use CLBlast-style register-blocked kernel
-                // Kernel uses 16x16 work group (256 threads), each computes 4x4 outputs = 64x64 tile
-                _logger?.LogWarning("[GEMM {M}x{N}x{K}] FALLBACK kernel: gemm_double_buffered (expected ~30% slower than CLBlast)", M, N, K);
-                var kernel = _kernelCache["gemm_double_buffered"];
+                bool forceUnsafe = GetEnvBool("AIDOTNET_GEMM_UNSAFE");
+                bool safeKernel = !forceUnsafe && M <= 512 && N <= 512 && K <= 512;
+                if (GetEnvBool("AIDOTNET_GEMM_SAFE"))
+                    safeKernel = true;
 
-                kernel.SetArg(0, bufferA.Handle);
-                kernel.SetArg(1, bufferB.Handle);
-                kernel.SetArg(2, bufferC.Handle);
-                kernel.SetArg(3, M);
-                kernel.SetArg(4, N);
-                kernel.SetArg(5, K);
-                kernel.SetArg(6, alpha);
-                kernel.SetArg(7, beta);
-
-                // CRITICAL: Correct global work size calculation for tiled kernel
-                // Each work group processes a 64x64 output tile
-                // Work group size is 16x16 (256 threads)
-                // Each thread computes 4x4 outputs
-                int numTilesM = (M + GemmKernel.TILE_M - 1) / GemmKernel.TILE_M;  // Number of 64-row tiles
-                int numTilesN = (N + GemmKernel.TILE_N - 1) / GemmKernel.TILE_N;  // Number of 64-col tiles
-                int globalSizeX = numTilesM * GemmKernel.WG_SIZE_M;  // 16 threads per tile in M
-                int globalSizeY = numTilesN * GemmKernel.WG_SIZE_N;  // 16 threads per tile in N
-
-                kernel.Execute2D(globalSizeX, globalSizeY, GemmKernel.WG_SIZE_M, GemmKernel.WG_SIZE_N);
+                if (safeKernel)
+                {
+                    if (EnableTuningDiagnostics || traceEnabled)
+                        _logger?.LogDebug("[GEMM {M}x{N}x{K}] Kernel: gemm_medium_tile (safe mode)", M, N, K);
+                    GemmMediumTile(A, B, C, M, N, K, alpha, beta);
+                }
+                else
+                {
+                    if (EnableTuningDiagnostics || traceEnabled)
+                        _logger?.LogDebug("[GEMM {M}x{N}x{K}] Kernel: gemm_double_buffered (fallback)", M, N, K);
+                    GemmDoubleBuffered(A, B, C, M, N, K, alpha, beta);
+                }
             }
             else
             {
                 // Small matrix - use simple kernel (one thread per output)
+                var bufferA = ((DirectOpenClGpuBuffer)A).Buffer;
+                var bufferB = ((DirectOpenClGpuBuffer)B).Buffer;
+                var bufferC = ((DirectOpenClGpuBuffer)C).Buffer;
                 var kernel = _kernelCache["gemm_small"];
 
                 kernel.SetArg(0, bufferA.Handle);
@@ -1798,17 +1809,66 @@ public void Gemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K,
                 kernel.SetArg(6, alpha);
                 kernel.SetArg(7, beta);
 
-                // Simple kernel: one thread per output element
                 var (localSizeX, localSizeY) = CalculateOptimalWorkGroupSize(M, N);
                 kernel.Execute2D(M, N, localSizeX, localSizeY);
             }
+        }
+
+        private bool TryExecutePackedDynamicGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K, float alpha, float beta, GemmConfig config)
+        {
+            int kReg = config.KReg > 0 ? config.KReg : 1;
+            int kUnit = config.TileK * kReg;
+            if (config.TileM <= 0 || config.TileN <= 0 || kUnit <= 0)
+                return false;
+
+            int mPad = CeilDiv(M, config.TileM) * config.TileM;
+            int nPad = CeilDiv(N, config.TileN) * config.TileN;
+            int kPad = CeilDiv(K, kUnit) * kUnit;
+
+            bool useColumnMajorA = config.UseColumnMajorA;
+            bool needsPadding = mPad != M || nPad != N || kPad != K;
+            
+            if (!needsPadding && !useColumnMajorA)
+                return TryExecuteDynamicGemm(A, B, C, M, N, K, alpha, beta, config);
+
+            long aSize = (long)mPad * kPad;
+            long bSize = (long)kPad * nPad;
+            long cSize = (long)mPad * nPad;
+            if (aSize > int.MaxValue || bSize > int.MaxValue || cSize > int.MaxValue)
+                return false;
+
+            if (EnableTuningDiagnostics)
+            {
+                Trace.WriteLine($"[GEMM] Packed GEMM: {M}x{N}x{K} -> {mPad}x{nPad}x{kPad}");
+            }
+
+            using var aPad = AllocateBuffer((int)aSize);
+            using var bPad = AllocateBuffer((int)bSize);
+            using var cPad = AllocateBuffer((int)cSize);
+
+            if (useColumnMajorA)
+                PadCopyTransposeMatrix(A, aPad, M, K, mPad, kPad);
+            else
+                PadCopyMatrix(A, aPad, M, K, mPad, kPad);
 
-            // Note: Synchronization removed from inner loop - caller should sync when needed
+            PadCopyMatrix(B, bPad, K, N, kPad, nPad);
+
+            if (!IsEffectivelyZero(beta))
+                PadCopyMatrix(C, cPad, M, N, mPad, nPad);
+            else
+                ZeroBuffer(cPad, (int)cSize);
+
+            if (!TryExecuteDynamicGemm(aPad, bPad, cPad, mPad, nPad, kPad, alpha, beta, config))
+                return false;
+
+            CopySubmatrix(cPad, C, M, N, nPad, N);
+
+            return true;
         }
 
         public IGpuBuffer MatMul(IGpuBuffer A, IGpuBuffer B, int M, int N, int K)
         {
-            Console.WriteLine($"[OpenClBackend.MatMul] Called: {M}x{N}x{K}");
+            Trace.WriteLine($"[OpenClBackend.MatMul] Called: {M}x{N}x{K}");
             var C = AllocateBuffer(M * N);
             Gemm(A, B, C, M, N, K, 1.0f, 0.0f);
             // Sync only when returning buffer that might be immediately read
@@ -1824,6 +1884,14 @@ public void GemmClblastRdna1(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, in
             ExecuteGemmKernel("gemm_clblast_rdna1", A, B, C, M, N, K, alpha, beta, 64, 8);
         }
 
+        /// <summary>
+        /// Double-buffered GEMM kernel used as the built-in fallback for larger matrices.
+        /// </summary>
+        public void GemmDoubleBuffered(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, int N, int K, float alpha = 1.0f, float beta = 0.0f)
+        {
+            ExecuteGemmKernel("gemm_double_buffered", A, B, C, M, N, K, alpha, beta, 64, 16);
+        }
+
         /// <summary>
         /// Medium tile GEMM kernel.
         /// </summary>
@@ -3506,65 +3574,65 @@ void WriteColored(string message, ConsoleColor color)
             {
                 if (Console.IsOutputRedirected)
                 {
-                    Console.WriteLine(message);
+                    Trace.WriteLine(message);
                     return;
                 }
 
                 var previous = Console.ForegroundColor;
                 Console.ForegroundColor = color;
-                Console.WriteLine(message);
+                Trace.WriteLine(message);
                 Console.ForegroundColor = previous;
             }
 
-            Console.WriteLine();
-            Console.WriteLine("=== OpenCL GEMM Diagnostics ===");
-            Console.WriteLine($"Matrix dimensions: M={diagnostics.M}, N={diagnostics.N}, K={diagnostics.K}");
-            Console.WriteLine($"Kernel: {diagnostics.KernelName}");
-            Console.WriteLine($"Work configuration: Global({diagnostics.GlobalSizeX}x{diagnostics.GlobalSizeY}), Local({diagnostics.LocalSizeX}x{diagnostics.LocalSizeY})");
-            Console.WriteLine($"Work items launched: {diagnostics.WorkItemsLaunched:N0}");
-            Console.WriteLine($"Work groups launched: {diagnostics.WorkGroupsLaunched:N0}");
-            Console.WriteLine();
+            Trace.WriteLine("");
+            Trace.WriteLine("=== OpenCL GEMM Diagnostics ===");
+            Trace.WriteLine($"Matrix dimensions: M={diagnostics.M}, N={diagnostics.N}, K={diagnostics.K}");
+            Trace.WriteLine($"Kernel: {diagnostics.KernelName}");
+            Trace.WriteLine($"Work configuration: Global({diagnostics.GlobalSizeX}x{diagnostics.GlobalSizeY}), Local({diagnostics.LocalSizeX}x{diagnostics.LocalSizeY})");
+            Trace.WriteLine($"Work items launched: {diagnostics.WorkItemsLaunched:N0}");
+            Trace.WriteLine($"Work groups launched: {diagnostics.WorkGroupsLaunched:N0}");
+            Trace.WriteLine("");
 
             if (diagnostics.IsProfilingAvailable && diagnostics.KernelExecutionNs > 0)
             {
-                Console.WriteLine("--- GPU Timing (from OpenCL events) ---");
-                Console.WriteLine($"Queue to Submit: {diagnostics.QueueToSubmitNs / 1e6:F3} ms");
-                Console.WriteLine($"Submit to Start (launch overhead): {diagnostics.SubmitToStartNs / 1e6:F3} ms");
-                Console.WriteLine($"Kernel Execution: {diagnostics.KernelExecutionNs / 1e6:F3} ms");
-                Console.WriteLine($"Total GPU Time: {diagnostics.TotalGpuTimeNs / 1e6:F3} ms");
+                Trace.WriteLine("--- GPU Timing (from OpenCL events) ---");
+                Trace.WriteLine($"Queue to Submit: {diagnostics.QueueToSubmitNs / 1e6:F3} ms");
+                Trace.WriteLine($"Submit to Start (launch overhead): {diagnostics.SubmitToStartNs / 1e6:F3} ms");
+                Trace.WriteLine($"Kernel Execution: {diagnostics.KernelExecutionNs / 1e6:F3} ms");
+                Trace.WriteLine($"Total GPU Time: {diagnostics.TotalGpuTimeNs / 1e6:F3} ms");
             }
             else if (!string.IsNullOrEmpty(diagnostics.ProfilingError))
             {
-                Console.WriteLine($"Profiling error: {diagnostics.ProfilingError}");
+                Trace.WriteLine($"Profiling error: {diagnostics.ProfilingError}");
             }
 
-            Console.WriteLine($"Wall clock time: {diagnostics.WallClockMs:F3} ms");
-            Console.WriteLine();
+            Trace.WriteLine($"Wall clock time: {diagnostics.WallClockMs:F3} ms");
+            Trace.WriteLine("");
 
-            Console.WriteLine("--- Performance Metrics ---");
-            Console.WriteLine($"FLOPS required: {diagnostics.FlopsRequired:N0} ({diagnostics.FlopsRequired / 1e9:F2} GFLOP)");
-            Console.WriteLine($"Bytes transferred: {diagnostics.BytesTransferred:N0} ({diagnostics.BytesTransferred / 1e6:F2} MB)");
-            Console.WriteLine($"Arithmetic intensity: {diagnostics.ArithmeticIntensity:F2} FLOP/byte");
-            Console.WriteLine($"Achieved GFLOPS: {diagnostics.AchievedGflops:F2}");
-            Console.WriteLine($"Achieved bandwidth: {diagnostics.AchievedBandwidthGBps:F2} GB/s");
-            Console.WriteLine($"Compute efficiency: {diagnostics.ComputeEfficiency:F1}% of theoretical peak");
-            Console.WriteLine();
+            Trace.WriteLine("--- Performance Metrics ---");
+            Trace.WriteLine($"FLOPS required: {diagnostics.FlopsRequired:N0} ({diagnostics.FlopsRequired / 1e9:F2} GFLOP)");
+            Trace.WriteLine($"Bytes transferred: {diagnostics.BytesTransferred:N0} ({diagnostics.BytesTransferred / 1e6:F2} MB)");
+            Trace.WriteLine($"Arithmetic intensity: {diagnostics.ArithmeticIntensity:F2} FLOP/byte");
+            Trace.WriteLine($"Achieved GFLOPS: {diagnostics.AchievedGflops:F2}");
+            Trace.WriteLine($"Achieved bandwidth: {diagnostics.AchievedBandwidthGBps:F2} GB/s");
+            Trace.WriteLine($"Compute efficiency: {diagnostics.ComputeEfficiency:F1}% of theoretical peak");
+            Trace.WriteLine("");
 
-            Console.WriteLine("--- Bottleneck Analysis ---");
+            Trace.WriteLine("--- Bottleneck Analysis ---");
             if (diagnostics.SubmitToStartNs > diagnostics.KernelExecutionNs * 0.5 && diagnostics.KernelExecutionNs > 0)
             {
                 WriteColored("WARNING: High launch overhead detected (>50% of kernel time)", ConsoleColor.Yellow);
-                Console.WriteLine("  -> Consider batching multiple small operations");
+                Trace.WriteLine("  -> Consider batching multiple small operations");
             }
             if (diagnostics.IsLikelyMemoryBound)
             {
                 WriteColored("LIKELY MEMORY BOUND: Achieved GFLOPS limited by memory bandwidth", ConsoleColor.Yellow);
-                Console.WriteLine("  -> Consider using data tiling, caching, or reducing data movement");
+                Trace.WriteLine("  -> Consider using data tiling, caching, or reducing data movement");
             }
             else if (diagnostics.ComputeEfficiency < 50)
             {
                 WriteColored("LIKELY COMPUTE BOUND with low efficiency:", ConsoleColor.Red);
-                Console.WriteLine("  -> Check for bank conflicts, divergent warps, or suboptimal work group size");
+                Trace.WriteLine("  -> Check for bank conflicts, divergent warps, or suboptimal work group size");
             }
             else if (diagnostics.ComputeEfficiency < 80)
             {
@@ -3574,7 +3642,7 @@ void WriteColored(string message, ConsoleColor color)
             {
                 WriteColored("GOOD EFFICIENCY: Kernel is well-optimized", ConsoleColor.Green);
             }
-            Console.WriteLine();
+            Trace.WriteLine("");
         }
 
         /// <summary>
@@ -3584,27 +3652,27 @@ public void RunGemmBenchmark(int[] sizes, int warmupIterations = 3, int benchmar
         {
             if (_context == null)
             {
-                Console.WriteLine("OpenCL context not available");
+                Trace.WriteLine("OpenCL context not available");
                 return;
             }
 
             var deviceInfo = GetDeviceInfo();
-            Console.WriteLine("=== OpenCL GEMM Benchmark ===");
-            Console.WriteLine($"Device: {deviceInfo.DeviceName}");
-            Console.WriteLine($"Vendor: {deviceInfo.DeviceVendor}");
-            Console.WriteLine($"Compute Units: {deviceInfo.ComputeUnits}");
-            Console.WriteLine($"Clock: {deviceInfo.ClockFrequencyMHz} MHz");
-            Console.WriteLine($"Theoretical Peak: {deviceInfo.TheoreticalPeakGflops:F0} GFLOPS");
-            Console.WriteLine($"Profiling enabled: {IsProfilingEnabled}");
-            Console.WriteLine();
+            Trace.WriteLine("=== OpenCL GEMM Benchmark ===");
+            Trace.WriteLine($"Device: {deviceInfo.DeviceName}");
+            Trace.WriteLine($"Vendor: {deviceInfo.DeviceVendor}");
+            Trace.WriteLine($"Compute Units: {deviceInfo.ComputeUnits}");
+            Trace.WriteLine($"Clock: {deviceInfo.ClockFrequencyMHz} MHz");
+            Trace.WriteLine($"Theoretical Peak: {deviceInfo.TheoreticalPeakGflops:F0} GFLOPS");
+            Trace.WriteLine($"Profiling enabled: {IsProfilingEnabled}");
+            Trace.WriteLine("");
 
             int sizeIndex = 0;
             foreach (int size in sizes)
             {
                 sizeIndex++;
                 int M = size, N = size, K = size;
-                Console.WriteLine($"[Progress] {sizeIndex}/{sizes.Length} size {size}x{size}x{size}");
-                Console.WriteLine($"--- Matrix size: {size}x{size}x{size} ---");
+                Trace.WriteLine($"[Progress] {sizeIndex}/{sizes.Length} size {size}x{size}x{size}");
+                Trace.WriteLine($"--- Matrix size: {size}x{size}x{size} ---");
 
                 // Allocate buffers
                 var dataA = new float[M * K];
@@ -3657,16 +3725,16 @@ public void RunGemmBenchmark(int[] sizes, int warmupIterations = 3, int benchmar
                 avgBandwidth /= benchmarkIterations;
                 avgLaunchOverhead /= benchmarkIterations;
 
-                Console.WriteLine($"  Kernel: {allDiagnostics[0].KernelName}");
-                Console.WriteLine($"  Avg kernel time: {avgKernelTimeMs:F3} ms");
+                Trace.WriteLine($"  Kernel: {allDiagnostics[0].KernelName}");
+                Trace.WriteLine($"  Avg kernel time: {avgKernelTimeMs:F3} ms");
                 if (avgLaunchOverhead > 0)
                 {
-                    Console.WriteLine($"  Avg launch overhead: {avgLaunchOverhead:F3} ms");
+                    Trace.WriteLine($"  Avg launch overhead: {avgLaunchOverhead:F3} ms");
                 }
-                Console.WriteLine($"  Avg GFLOPS: {avgGflops:F2}");
-                Console.WriteLine($"  Avg bandwidth: {avgBandwidth:F2} GB/s");
-                Console.WriteLine($"  Efficiency: {avgGflops / deviceInfo.TheoreticalPeakGflops * 100:F1}%");
-                Console.WriteLine();
+                Trace.WriteLine($"  Avg GFLOPS: {avgGflops:F2}");
+                Trace.WriteLine($"  Avg bandwidth: {avgBandwidth:F2} GB/s");
+                Trace.WriteLine($"  Efficiency: {avgGflops / deviceInfo.TheoreticalPeakGflops * 100:F1}%");
+                Trace.WriteLine("");
             }
         }
 
@@ -3714,13 +3782,13 @@ public TuningResult[] RunBayesianGemmOptimization(int M, int N, int K, int maxTr
             var capabilities = GpuCapabilities.Detect(ComputeUnits, GlobalMemoryBytes, (int)LocalMemoryBytes,
                 (int)_maxWorkGroupSize, DeviceVendor, DeviceName, _context.Extensions);
 
-            Console.WriteLine("=== Bayesian GEMM Optimization ===");
-            Console.WriteLine($"Matrix: {M}x{N}x{K}, Device: {DeviceName}, Max trials: {maxTrials}");
+            Trace.WriteLine("=== Bayesian GEMM Optimization ===");
+            Trace.WriteLine($"Matrix: {M}x{N}x{K}, Device: {DeviceName}, Max trials: {maxTrials}");
 
             // Print GPU capabilities if diagnostics enabled
             if (EnableTuningDiagnostics)
             {
-                Console.WriteLine("[GPU Capabilities]");
+                Trace.WriteLine("[GPU Capabilities]");
                 Console.Write(capabilities.GetDiagnosticString());
             }
 
@@ -3761,7 +3829,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                     benchmarkFailures++;
                     if (EnableTuningDiagnostics)
                     {
-                        Console.WriteLine($"  [Validation] {config.KernelName}: {validationError}");
+                        Trace.WriteLine($"  [Validation] {config.KernelName}: {validationError}");
                     }
                     database.MarkAsTested(M, N, K, config, 0);
                     return double.NaN;
@@ -3774,7 +3842,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                     {
                         if (EnableTuningDiagnostics)
                         {
-                            Console.WriteLine($"  [Cache] {config.KernelName}: {cachedGflops.Value:F2} GFLOPS");
+                            Trace.WriteLine($"  [Cache] {config.KernelName}: {cachedGflops.Value:F2} GFLOPS");
                         }
 
                         return ops / (cachedGflops.Value * 1e6);
@@ -3808,12 +3876,12 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                 catch (Exception ex)
                 {
                     benchmarkFailures++;
-                    Console.WriteLine($"  Config {config} failed: {ex.Message}");
+                    Trace.WriteLine($"  Config {config} failed: {ex.Message}");
 
                     // Print kernel stats on failure
                     if (EnableTuningDiagnostics && _dynamicGemm != null)
                     {
-                        Console.WriteLine($"  [DynamicGemm Stats] {_dynamicGemm.GetDiagnosticStats()}");
+                        Trace.WriteLine($"  [DynamicGemm Stats] {_dynamicGemm.GetDiagnosticStats()}");
                     }
 
                     database.MarkAsTested(M, N, K, config, 0);
@@ -3833,8 +3901,8 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
             {
                 var cachedConfig = cachedEntry.Value.Config;
                 databaseGflops = cachedEntry.Value.GFlops;  // Use stored GFLOPS as baseline
-                Console.WriteLine($"Using cached configuration: {cachedConfig}");
-                Console.WriteLine($"Database best: {databaseGflops:F2} GFLOPS (threshold to beat)");
+                Trace.WriteLine($"Using cached configuration: {cachedConfig}");
+                Trace.WriteLine($"Database best: {databaseGflops:F2} GFLOPS (threshold to beat)");
 
                 // Re-benchmark to validate config works and add to result set
                 var cachedTimeMs = BenchmarkConfigNoCache(cachedConfig);
@@ -3848,7 +3916,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                         GFlops = revalidatedGflops,
                         IsValid = true
                     };
-                    Console.WriteLine($"Revalidated: {revalidatedGflops:F2} GFLOPS");
+                    Trace.WriteLine($"Revalidated: {revalidatedGflops:F2} GFLOPS");
                 }
             }
 
@@ -3859,10 +3927,10 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
             // Print final statistics
             if (EnableTuningDiagnostics)
             {
-                Console.WriteLine($"\n[Benchmark Stats] Attempts: {benchmarkAttempts}, Failures: {benchmarkFailures}");
+                Trace.WriteLine($"\n[Benchmark Stats] Attempts: {benchmarkAttempts}, Failures: {benchmarkFailures}");
                 if (_dynamicGemm != null)
                 {
-                    Console.WriteLine($"[DynamicGemm Stats] {_dynamicGemm.GetDiagnosticStats()}");
+                    Trace.WriteLine($"[DynamicGemm Stats] {_dynamicGemm.GetDiagnosticStats()}");
                 }
             }
 
@@ -3879,18 +3947,18 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
             if (allResults.Count > 0 && allResults[0].IsValid)
             {
                 var best = allResults[0];
-                Console.WriteLine($"Best: {best.Config} - {best.GFlops:F2} GFLOPS");
+                Trace.WriteLine($"Best: {best.Config} - {best.GFlops:F2} GFLOPS");
 
                 // Only update database if we found something better than the DATABASE best
                 // Note: We compare against databaseGflops (historical best), NOT re-benchmarked value
                 if (best.GFlops > databaseGflops)
                 {
-                    Console.WriteLine($"NEW GLOBAL BEST! {best.GFlops:F2} > {databaseGflops:F2} GFLOPS (previous best)");
+                    Trace.WriteLine($"NEW GLOBAL BEST! {best.GFlops:F2} > {databaseGflops:F2} GFLOPS (previous best)");
                     database.StoreResult(M, N, K, best.Config, best.GFlops);
                 }
                 else
                 {
-                    Console.WriteLine($"No improvement: {best.GFlops:F2} <= {databaseGflops:F2} GFLOPS (database best)");
+                    Trace.WriteLine($"No improvement: {best.GFlops:F2} <= {databaseGflops:F2} GFLOPS (database best)");
                 }
             }
 
@@ -3911,8 +3979,8 @@ public TuningResult[] RunExhaustiveGemmOptimization(int M, int N, int K, int war
             var capabilities = GpuCapabilities.Detect(ComputeUnits, GlobalMemoryBytes, (int)LocalMemoryBytes,
                 (int)_maxWorkGroupSize, DeviceVendor, DeviceName, _context.Extensions);
 
-            Console.WriteLine("=== EXHAUSTIVE GEMM Optimization (CLBlast-style) ===");
-            Console.WriteLine($"Matrix: {M}x{N}x{K}, Device: {DeviceName}");
+            Trace.WriteLine("=== EXHAUSTIVE GEMM Optimization (CLBlast-style) ===");
+            Trace.WriteLine($"Matrix: {M}x{N}x{K}, Device: {DeviceName}");
 
             var dataA = new float[M * K];
             var dataB = new float[K * N];
@@ -3945,7 +4013,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                 {
                     if (EnableTuningDiagnostics)
                     {
-                        Console.WriteLine($"  [Validation] {config.KernelName}: {validationError}");
+                        Trace.WriteLine($"  [Validation] {config.KernelName}: {validationError}");
                     }
                     database.MarkAsTested(M, N, K, config, 0);
                     return double.NaN;
@@ -3958,7 +4026,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                     {
                         if (EnableTuningDiagnostics)
                         {
-                            Console.WriteLine($"  [Cache] {config.KernelName}: {cachedGflops.Value:F2} GFLOPS");
+                            Trace.WriteLine($"  [Cache] {config.KernelName}: {cachedGflops.Value:F2} GFLOPS");
                         }
 
                         return ops / (cachedGflops.Value * 1e6);
@@ -3989,7 +4057,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
                 }
                 catch (Exception ex)
                 {
-                    Console.WriteLine($"  Config {config} failed: {ex.Message}");
+                    Trace.WriteLine($"  Config {config} failed: {ex.Message}");
                     database.MarkAsTested(M, N, K, config, 0);
                     return double.NaN;
                 }
@@ -4003,7 +4071,7 @@ double BenchmarkConfig(GemmConfig config, bool allowCached)
             if (results.Length > 0 && results[0].IsValid)
             {
                 var best = results[0];
-                Console.WriteLine($"EXHAUSTIVE Best: {best.Config} - {best.GFlops:F2} GFLOPS");
+                Trace.WriteLine($"EXHAUSTIVE Best: {best.Config} - {best.GFlops:F2} GFLOPS");
 
                 database.StoreResult(M, N, K, best.Config, best.GFlops);
             }
@@ -4048,7 +4116,7 @@ public double GemmWithClBlast(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, int M, i
 
             if (status != ClBlastNative.StatusCode.Success)
             {
-                Console.WriteLine($"CLBlast SGEMM failed with status: {status}");
+                Trace.WriteLine($"CLBlast SGEMM failed with status: {status}");
                 return -1.0;
             }
 
@@ -4293,7 +4361,7 @@ public RooflineAnalysis ProfileGemm(IGpuBuffer A, IGpuBuffer B, IGpuBuffer C, in
             // Auto-print if profiling is enabled
             if (GetEnvBool(GemmProfileEnvVar))
             {
-                Console.WriteLine(result);
+                Trace.WriteLine(result);
             }
 
             return result;
@@ -4397,7 +4465,7 @@ public string CompareWithClBlast(int M, int N, int K, int warmupRuns = 3, int be
         /// </summary>
         public static void PrintDiagnosticHelp()
         {
-            Console.WriteLine(@"
+            Trace.WriteLine(@"
 === AiDotNet GPU Diagnostic Environment Variables ===
 
 TIMING & TRACING:
@@ -4419,6 +4487,10 @@ GPU SPECS (for accurate roofline analysis):
 
 DEBUG:
   AIDOTNET_FORCE_DIRECT=1     Force XgemmDirect path (skip indirect path)
+  AIDOTNET_GEMM_ENABLE_DYNAMIC=1  Enable dynamic GEMM kernels (default: off)
+  AIDOTNET_GEMM_SAFE=1        Use safe GEMM fallback kernel for correctness
+  AIDOTNET_GEMM_UNSAFE=1      Force gemm_double_buffered fallback kernel
+  AIDOTNET_GEMM_VALIDATE=1    Validate GEMM output for NaN/Inf and fall back to CPU
 
 KERNEL VARIANTS (A/B testing):
   AIDOTNET_GEMM_VARIANT=0     Original CLBlast baseline
@@ -8673,9 +8745,20 @@ private void CopyBuffer(IGpuBuffer src, IGpuBuffer dst, int size)
 
         private void ZeroBuffer(IGpuBuffer buffer, int size)
         {
+            if (size <= 0)
+                return;
+            if (_context != null && _kernelCache.TryGetValue("zero_buffer", out var kernel))
+            {
+                var buf = ((DirectOpenClGpuBuffer)buffer).Buffer;
+                kernel.SetArg(0, buf.Handle);
+                kernel.SetArg(1, size);
+                kernel.Execute1D(size, Math.Min(256, size));
+                return;
+            }
+
             var data = new float[size];
-            var buf = ((DirectOpenClGpuBuffer)buffer).Buffer;
-            buf.CopyFromHost(data);
+            var fallback = ((DirectOpenClGpuBuffer)buffer).Buffer;
+            fallback.CopyFromHost(data);
         }
 
         private void ScaleBuffer(IGpuBuffer buffer, float scale, int size)
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClNativeBindings.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClNativeBindings.cs
index ac5a5ac966..4d6daca358 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClNativeBindings.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/OpenCL/OpenClNativeBindings.cs
@@ -3,6 +3,7 @@
 // Works on ALL .NET versions including .NET Framework 4.6.2.
 
 using System;
+using System.Diagnostics;
 using System.Runtime.InteropServices;
 
 namespace AiDotNet.Tensors.Engines.DirectGpu.OpenCL
@@ -15,6 +16,7 @@ internal static class OpenClNativeBindings
     {
         // OpenCL library name varies by platform
         private const string OpenClLibrary = "OpenCL";
+        public static bool EnableDiagnostics { get; set; }
 
         #region Error Codes
 
@@ -354,12 +356,12 @@ public static extern int GetKernelWorkGroupInfo(
         #region Helper Methods
 
         /// <summary>
-        /// Prints diagnostic information about OpenCL DLL search paths.
+        /// Prints diagnostic information about OpenCL DLL search paths.        
         /// Call this to debug DLL loading issues.
         /// </summary>
         public static void PrintDllSearchDiagnostics()
         {
-            Console.WriteLine("[OpenCL DLL Diagnostics] Searching for OpenCL.dll...");
+            LogDiagnostic("[OpenCL DLL Diagnostics] Searching for OpenCL.dll...");
 
             // Check common Windows locations for OpenCL.dll
             var searchPaths = new[]
@@ -372,12 +374,12 @@ public static void PrintDllSearchDiagnostics()
             foreach (var path in searchPaths)
             {
                 bool exists = System.IO.File.Exists(path);
-                Console.WriteLine($"[OpenCL DLL Diagnostics] {path} - {(exists ? "FOUND" : "not found")}");
+                LogDiagnostic($"[OpenCL DLL Diagnostics] {path} - {(exists ? "FOUND" : "not found")}");
             }
 
             // Check PATH environment variable
-            var pathEnv = Environment.GetEnvironmentVariable("PATH") ?? "";
-            Console.WriteLine("[OpenCL DLL Diagnostics] Checking PATH directories...");
+            var pathEnv = Environment.GetEnvironmentVariable("PATH") ?? "";     
+            LogDiagnostic("[OpenCL DLL Diagnostics] Checking PATH directories...");
 
             var pathDirs = pathEnv.Split(';');
             bool foundInPath = false;
@@ -389,7 +391,7 @@ public static void PrintDllSearchDiagnostics()
                     var openClPath = System.IO.Path.Combine(dir.Trim(), "OpenCL.dll");
                     if (System.IO.File.Exists(openClPath))
                     {
-                        Console.WriteLine($"[OpenCL DLL Diagnostics] Found in PATH: {openClPath}");
+                        LogDiagnostic($"[OpenCL DLL Diagnostics] Found in PATH: {openClPath}");
                         foundInPath = true;
                     }
                 }
@@ -401,13 +403,13 @@ public static void PrintDllSearchDiagnostics()
 
             if (!foundInPath)
             {
-                Console.WriteLine("[OpenCL DLL Diagnostics] OpenCL.dll NOT found in any PATH directory");
+                LogDiagnostic("[OpenCL DLL Diagnostics] OpenCL.dll NOT found in any PATH directory");
             }
 
-            Console.WriteLine("[OpenCL DLL Diagnostics] If OpenCL.dll is missing, install GPU drivers:");
-            Console.WriteLine("[OpenCL DLL Diagnostics]   - NVIDIA: Install GeForce or CUDA drivers");
-            Console.WriteLine("[OpenCL DLL Diagnostics]   - AMD: Install Adrenalin drivers");
-            Console.WriteLine("[OpenCL DLL Diagnostics]   - Intel: Install Intel Graphics drivers or OpenCL Runtime");
+            LogDiagnostic("[OpenCL DLL Diagnostics] If OpenCL.dll is missing, install GPU drivers:");
+            LogDiagnostic("[OpenCL DLL Diagnostics]   - NVIDIA: Install GeForce or CUDA drivers");
+            LogDiagnostic("[OpenCL DLL Diagnostics]   - AMD: Install Adrenalin drivers");
+            LogDiagnostic("[OpenCL DLL Diagnostics]   - Intel: Install Intel Graphics drivers or OpenCL Runtime");
         }
 
         /// <summary>
@@ -589,25 +591,31 @@ public static bool IsAvailable
             {
                 try
                 {
-                    int err = GetPlatformIDs(0, null, out uint numPlatforms);
-                    bool available = err == CL_SUCCESS && numPlatforms > 0;
-                    Console.WriteLine($"[OpenCL Diagnostics] GetPlatformIDs returned error code: {err}, platforms found: {numPlatforms}, available: {available}");
+                    int err = GetPlatformIDs(0, null, out uint numPlatforms); // lgtm[cs/call-to-unmanaged-code]
+                    bool available = err == CL_SUCCESS && numPlatforms > 0;     
+                    LogDiagnostic($"[OpenCL Diagnostics] GetPlatformIDs returned error code: {err}, platforms found: {numPlatforms}, available: {available}");
                     return available;
                 }
                 catch (DllNotFoundException ex)
                 {
-                    Console.WriteLine($"[OpenCL Diagnostics] DllNotFoundException: {ex.Message}");
+                    LogDiagnostic($"[OpenCL Diagnostics] DllNotFoundException: {ex.Message}");
                     PrintDllSearchDiagnostics();
                     return false;
                 }
                 catch (Exception ex)
                 {
-                    Console.WriteLine($"[OpenCL Diagnostics] Exception during OpenCL availability check: {ex.GetType().Name}: {ex.Message}");
+                    LogDiagnostic($"[OpenCL Diagnostics] Exception during OpenCL availability check: {ex.GetType().Name}: {ex.Message}");
                     return false;
                 }
             }
         }
 
+        private static void LogDiagnostic(string message)
+        {
+            if (EnableDiagnostics)
+                Trace.WriteLine(message);
+        }
+
         /// <summary>
         /// Gets a ulong profiling info value from an event.
         /// Timestamps are in nanoseconds.
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpu/Profiling/GemmProfiler.cs b/src/AiDotNet.Tensors/Engines/DirectGpu/Profiling/GemmProfiler.cs
index e42eb79850..6b6f932f32 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpu/Profiling/GemmProfiler.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpu/Profiling/GemmProfiler.cs
@@ -135,15 +135,15 @@ public ProfileResult RunFullProfile()
 
         if (_config.Verbose)
         {
-            Console.WriteLine("=".PadRight(80, '='));
-            Console.WriteLine("GEMM PROFILING SESSION");
-            Console.WriteLine("=".PadRight(80, '='));
-            Console.WriteLine($"Device: {_backend.DeviceName}");
-            Console.WriteLine($"Architecture: {_arch.Name}");
-            Console.WriteLine($"Peak: {_roofline.PeakGflops:F0} GFLOPS, {_roofline.PeakBandwidthGBs:F0} GB/s");
-            Console.WriteLine($"Ridge Point: {_roofline.RidgePoint:F1} FLOPS/byte");
-            Console.WriteLine($"Sizes: {string.Join(", ", _config.Sizes)}");
-            Console.WriteLine("-".PadRight(80, '-'));
+            Trace.WriteLine("=".PadRight(80, '='));
+            Trace.WriteLine("GEMM PROFILING SESSION");
+            Trace.WriteLine("=".PadRight(80, '='));
+            Trace.WriteLine($"Device: {_backend.DeviceName}");
+            Trace.WriteLine($"Architecture: {_arch.Name}");
+            Trace.WriteLine($"Peak: {_roofline.PeakGflops:F0} GFLOPS, {_roofline.PeakBandwidthGBs:F0} GB/s");
+            Trace.WriteLine($"Ridge Point: {_roofline.RidgePoint:F1} FLOPS/byte");
+            Trace.WriteLine($"Sizes: {string.Join(", ", _config.Sizes)}");
+            Trace.WriteLine("-".PadRight(80, '-'));
         }
 
         foreach (var size in _config.Sizes)
@@ -155,14 +155,14 @@ public ProfileResult RunFullProfile()
 
                 if (_config.Verbose)
                 {
-                    Console.WriteLine(entry.GetSummary());
+                    Trace.WriteLine(entry.GetSummary());
                 }
             }
             catch (Exception ex)
             {
                 if (_config.Verbose)
                 {
-                    Console.WriteLine($"{size}x{size}x{size}: ERROR - {ex.Message}");
+                    Trace.WriteLine($"{size}x{size}x{size}: ERROR - {ex.Message}");
                 }
             }
         }
@@ -183,10 +183,10 @@ public ProfileResult RunFullProfile()
 
         if (_config.Verbose)
         {
-            Console.WriteLine("-".PadRight(80, '-'));
-            Console.WriteLine($"Profiling complete in {result.ProfileDurationSeconds:F1}s");
-            Console.WriteLine($"Best: {result.BestGflops:F0} GFLOPS ({result.BestEfficiencyPercent:F1}% efficiency)");
-            Console.WriteLine("=".PadRight(80, '='));
+            Trace.WriteLine("-".PadRight(80, '-'));
+            Trace.WriteLine($"Profiling complete in {result.ProfileDurationSeconds:F1}s");
+            Trace.WriteLine($"Best: {result.BestGflops:F0} GFLOPS ({result.BestEfficiencyPercent:F1}% efficiency)");
+            Trace.WriteLine("=".PadRight(80, '='));
         }
 
         // Export if configured
@@ -194,14 +194,14 @@ public ProfileResult RunFullProfile()
         {
             ProfileExporter.ToJsonFile(result, _config.JsonOutputPath);
             if (_config.Verbose)
-                Console.WriteLine($"JSON exported to: {_config.JsonOutputPath}");
+                Trace.WriteLine($"JSON exported to: {_config.JsonOutputPath}");
         }
 
         if (_config.CsvOutputPath is { Length: > 0 })
         {
             ProfileExporter.ToCsvFile(result, _config.CsvOutputPath);
             if (_config.Verbose)
-                Console.WriteLine($"CSV exported to: {_config.CsvOutputPath}");
+                Trace.WriteLine($"CSV exported to: {_config.CsvOutputPath}");
         }
 
         return result;
@@ -313,9 +313,9 @@ public ProfileResult RunRectangularProfile(int[] mSizes, int[] nSizes, int[] kSi
 
         if (_config.Verbose)
         {
-            Console.WriteLine("=".PadRight(80, '='));
-            Console.WriteLine("RECTANGULAR GEMM PROFILING");
-            Console.WriteLine("=".PadRight(80, '='));
+            Trace.WriteLine("=".PadRight(80, '='));
+            Trace.WriteLine("RECTANGULAR GEMM PROFILING");
+            Trace.WriteLine("=".PadRight(80, '='));
         }
 
         foreach (var m in mSizes)
@@ -331,14 +331,14 @@ public ProfileResult RunRectangularProfile(int[] mSizes, int[] nSizes, int[] kSi
 
                         if (_config.Verbose)
                         {
-                            Console.WriteLine(entry.GetSummary());
+                            Trace.WriteLine(entry.GetSummary());
                         }
                     }
                     catch (Exception ex)
                     {
                         if (_config.Verbose)
                         {
-                            Console.WriteLine($"{m}x{n}x{k}: ERROR - {ex.Message}");
+                            Trace.WriteLine($"{m}x{n}x{k}: ERROR - {ex.Message}");
                         }
                     }
                 }
diff --git a/src/AiDotNet.Tensors/Engines/DirectGpuTensorEngine.cs b/src/AiDotNet.Tensors/Engines/DirectGpuTensorEngine.cs
index e0188d7740..21560e8409 100644
--- a/src/AiDotNet.Tensors/Engines/DirectGpuTensorEngine.cs
+++ b/src/AiDotNet.Tensors/Engines/DirectGpuTensorEngine.cs
@@ -808,9 +808,7 @@ Matrix<T> IEngine.MatrixMultiply<T>(Matrix<T> a, Matrix<T> b)
             if (resultData == null)
                 return base.MatrixMultiply(a, b);
 
-            var result = new Matrix<T>(a.Rows, b.Columns);
-            resultData.AsSpan().CopyTo(result.AsWritableSpan());
-            return result;
+            return new Matrix<T>(a.Rows, b.Columns, resultData);
         }
         catch
         {
diff --git a/src/AiDotNet.Tensors/Engines/Engine.cs b/src/AiDotNet.Tensors/Engines/Engine.cs
index 1a7ca9524a..ab774dcdb7 100644
--- a/src/AiDotNet.Tensors/Engines/Engine.cs
+++ b/src/AiDotNet.Tensors/Engines/Engine.cs
@@ -2,6 +2,7 @@
 // Automatic engine selection based on hardware capabilities
 
 using System;
+using System.Diagnostics;
 using System.Numerics;
 using System.Runtime.Intrinsics;
 using System.Text;
@@ -167,22 +168,22 @@ private static IEngine CreateOptimalEngine()
                 if (directGpu != null && directGpu.IsAvailable)
                 {
                     var directEngine = new DirectGpuTensorEngine(directGpu);
-                    Console.WriteLine($"[Engine] DirectGpu available: {directGpu.BackendName} ({directGpu.DeviceName})");
-                    Console.WriteLine($"[Engine] DirectGpu: {directGpu.ComputeUnits} CUs, {directGpu.GlobalMemoryGB:F1}GB VRAM");
+                    Trace.WriteLine($"[Engine] DirectGpu available: {directGpu.BackendName} ({directGpu.DeviceName})");
+                    Trace.WriteLine($"[Engine] DirectGpu: {directGpu.ComputeUnits} CUs, {directGpu.GlobalMemoryGB:F1}GB VRAM");
                     return directEngine;
                 }
 
-                Console.WriteLine("[Engine] DirectGpu: Not available (no compatible GPU backends found)");
+                Trace.WriteLine("[Engine] DirectGpu: Not available (no compatible GPU backends found)");
             }
             catch (Exception ex)
             {
-                Console.WriteLine($"[Engine] DirectGpu initialization failed: {ex.Message}");
+                Trace.WriteLine($"[Engine] DirectGpu initialization failed: {ex.Message}");
             }
 
             // Fallback to CPU (always available)
             var cpuEngine = new CpuEngine();
-            Console.WriteLine($"[Engine] Auto-selected: {cpuEngine.Name}");
-            Console.WriteLine($"[Engine] No GPU detected - operations will use CPU with SIMD acceleration");
+            Trace.WriteLine($"[Engine] Auto-selected: {cpuEngine.Name}");
+            Trace.WriteLine($"[Engine] No GPU detected - operations will use CPU with SIMD acceleration");
             return cpuEngine;
         }
 
diff --git a/src/AiDotNet.Tensors/Helpers/TensorPrimitivesHelper.cs b/src/AiDotNet.Tensors/Helpers/TensorPrimitivesHelper.cs
index 7cb93a350b..b9c5bb40e1 100644
--- a/src/AiDotNet.Tensors/Helpers/TensorPrimitivesHelper.cs
+++ b/src/AiDotNet.Tensors/Helpers/TensorPrimitivesHelper.cs
@@ -348,5 +348,16 @@ public static T CosineSimilarity(Vector<T> a, Vector<T> b)
         return NumOps.CosineSimilarity(aArray, bArray);
     }
 
+    /// <summary>
+    /// Checks if any element in the vector is NaN or Infinity.
+    /// </summary>
+    /// <param name="x">The source vector.</param>
+    /// <param name="badIndex">The index of the first non-finite value found, or -1 if all are finite.</param>
+    /// <returns>True if any element is non-finite, false otherwise.</returns>
+    public static bool IsAnyNonFinite(Vector<T> x, out int badIndex)
+    {
+        return NumOps.IsAnyNonFinite(x.AsSpan(), out badIndex);
+    }
+
     #endregion
 }
diff --git a/src/AiDotNet.Tensors/Interfaces/IVectorizedOperations.cs b/src/AiDotNet.Tensors/Interfaces/IVectorizedOperations.cs
index 75321489a2..90efb8fd62 100644
--- a/src/AiDotNet.Tensors/Interfaces/IVectorizedOperations.cs
+++ b/src/AiDotNet.Tensors/Interfaces/IVectorizedOperations.cs
@@ -408,6 +408,22 @@ public interface IVectorizedOperations<T>
     /// <param name="destination">The destination span for values of type T.</param>
     void FromHalfSpan(ReadOnlySpan<Half> source, Span<T> destination);
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (not NaN or Infinity).
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">The index of the first non-finite value found, or -1 if all are finite.</param>
+    /// <returns>True if all elements are finite, false otherwise.</returns>
+    bool AllFinite(ReadOnlySpan<T> x, out int badIndex);
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">When this method returns true, contains the index of the first non-finite value found; otherwise, -1.</param>
+    /// <returns>True if any element is non-finite; otherwise, false.</returns>
+    bool IsAnyNonFinite(ReadOnlySpan<T> x, out int badIndex);
+
     #region Vectorized Activation Functions
 
     /// <summary>
diff --git a/src/AiDotNet.Tensors/LinearAlgebra/Matrix.cs b/src/AiDotNet.Tensors/LinearAlgebra/Matrix.cs
index 93e0168d96..f92ffec436 100644
--- a/src/AiDotNet.Tensors/LinearAlgebra/Matrix.cs
+++ b/src/AiDotNet.Tensors/LinearAlgebra/Matrix.cs
@@ -25,6 +25,13 @@ public Matrix(int rows, int columns) : base(rows, columns)
     {
     }
 
+    /// <summary>
+    /// Initializes a new matrix by wrapping an existing array.
+    /// </summary>
+    internal Matrix(int rows, int cols, T[] data) : base(rows, cols, data)
+    {
+    }
+
     /// <summary>
     /// Initializes a new matrix from a collection of collections, where each inner collection represents a row.
     /// </summary>
diff --git a/src/AiDotNet.Tensors/LinearAlgebra/MatrixBase.cs b/src/AiDotNet.Tensors/LinearAlgebra/MatrixBase.cs
index 141484f867..7e48c69d97 100644
--- a/src/AiDotNet.Tensors/LinearAlgebra/MatrixBase.cs
+++ b/src/AiDotNet.Tensors/LinearAlgebra/MatrixBase.cs
@@ -64,6 +64,18 @@ protected MatrixBase(int rows, int cols)
         this._data = new T[rows * cols];
     }
 
+    /// <summary>
+    /// Internal constructor to wrap an existing array.
+    /// </summary>
+    internal MatrixBase(int rows, int cols, T[] data)
+    {
+        _rows = rows;
+        _cols = cols;
+        _data = data ?? throw new ArgumentNullException(nameof(data));
+        if (data.Length < rows * cols)
+            throw new ArgumentException("Data array is too small for the specified dimensions.");
+    }
+
     /// <summary>
     /// Creates a matrix from a collection of row values.
     /// </summary>
diff --git a/src/AiDotNet.Tensors/NumericOperations/ByteOperations.cs b/src/AiDotNet.Tensors/NumericOperations/ByteOperations.cs
index 2eff0a9f5d..3113a1490f 100644
--- a/src/AiDotNet.Tensors/NumericOperations/ByteOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/ByteOperations.cs
@@ -671,6 +671,26 @@ public byte SignOrZero(byte value)
     /// <returns>The value as a double.</returns>
     public double ToDouble(byte value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// Byte values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<byte> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// Byte values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<byte> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/ComplexOperations.cs b/src/AiDotNet.Tensors/NumericOperations/ComplexOperations.cs
index 6f4a6de801..46b1f3c086 100644
--- a/src/AiDotNet.Tensors/NumericOperations/ComplexOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/ComplexOperations.cs
@@ -960,6 +960,33 @@ public double ToDouble(Complex<T> value)
         return _ops.ToDouble(value.Real);
     }
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<Complex<T>> x, out int badIndex)
+    {
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (_ops.IsNaN(x[i].Real) || _ops.IsInfinity(x[i].Real) ||
+                _ops.IsNaN(x[i].Imaginary) || _ops.IsInfinity(x[i].Imaginary))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<Complex<T>> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/DecimalOperations.cs b/src/AiDotNet.Tensors/NumericOperations/DecimalOperations.cs
index 1c4de29555..74e992e1a8 100644
--- a/src/AiDotNet.Tensors/NumericOperations/DecimalOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/DecimalOperations.cs
@@ -691,6 +691,26 @@ public decimal SignOrZero(decimal value)
     /// </summary>
     public double ToDouble(decimal value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// Decimals are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<decimal> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity. 
+    /// Decimals are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<decimal> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/DoubleOperations.cs b/src/AiDotNet.Tensors/NumericOperations/DoubleOperations.cs
index 34dba12d03..3a743389bc 100644
--- a/src/AiDotNet.Tensors/NumericOperations/DoubleOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/DoubleOperations.cs
@@ -801,6 +801,39 @@ public double SignOrZero(double value)
     /// <inheritdoc/>
     public bool SupportsGpuAcceleration => true;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">The index of the first non-finite value found, or -1 if all are finite.</param>
+    /// <returns>True if all elements are finite, false otherwise.</returns>
+    public bool AllFinite(ReadOnlySpan<double> x, out int badIndex)
+    {
+        // Fallback or find the exact index of the non-finite value
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (double.IsNaN(x[i]) || double.IsInfinity(x[i]))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">When this method returns true, contains the index of the first non-finite value found; otherwise, -1.</param>
+    /// <returns>True if any element is non-finite; otherwise, false.</returns>
+    public bool IsAnyNonFinite(ReadOnlySpan<double> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     #region IVectorizedOperations<double> Implementation - SIMD via TensorPrimitivesCore
 
     /// <summary>
diff --git a/src/AiDotNet.Tensors/NumericOperations/FloatOperations.cs b/src/AiDotNet.Tensors/NumericOperations/FloatOperations.cs
index 03974ad5c4..61d9e38647 100644
--- a/src/AiDotNet.Tensors/NumericOperations/FloatOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/FloatOperations.cs
@@ -865,6 +865,39 @@ public float SignOrZero(float value)
     /// </summary>
     public bool SupportsGpuAcceleration => true;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">The index of the first non-finite value found, or -1 if all are finite.</param>
+    /// <returns>True if all elements are finite, false otherwise.</returns>
+    public bool AllFinite(ReadOnlySpan<float> x, out int badIndex)
+    {
+        // Fallback or find the exact index of the non-finite value
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (float.IsNaN(x[i]) || float.IsInfinity(x[i]))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">When this method returns true, contains the index of the first non-finite value found; otherwise, -1.</param>
+    /// <returns>True if any element is non-finite; otherwise, false.</returns>
+    public bool IsAnyNonFinite(ReadOnlySpan<float> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     #region IVectorizedOperations<float> Implementation - SIMD via TensorPrimitives
 
     private static readonly FloatOperations _instance = new();
diff --git a/src/AiDotNet.Tensors/NumericOperations/HalfOperations.cs b/src/AiDotNet.Tensors/NumericOperations/HalfOperations.cs
index 1ca4cfcd42..4ee7792b2e 100644
--- a/src/AiDotNet.Tensors/NumericOperations/HalfOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/HalfOperations.cs
@@ -243,6 +243,36 @@ public Half SignOrZero(Half value)
     /// </remarks>
     public double ToDouble(Half value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<Half> x, out int badIndex)
+    {
+        // Fallback or find the exact index of the non-finite value
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (Half.IsNaN(x[i]) || Half.IsInfinity(x[i]))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    /// <param name="x">The source span.</param>
+    /// <param name="badIndex">When this method returns true, contains the index of the first non-finite value found; otherwise, -1.</param>
+    /// <returns>True if any element is non-finite; otherwise, false.</returns>
+    public bool IsAnyNonFinite(ReadOnlySpan<Half> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => true;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/Int32Operations.cs b/src/AiDotNet.Tensors/NumericOperations/Int32Operations.cs
index 471cec6dff..30d85f8ab4 100644
--- a/src/AiDotNet.Tensors/NumericOperations/Int32Operations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/Int32Operations.cs
@@ -740,8 +740,29 @@ public int SignOrZero(int value)
     /// <returns>The value as a double.</returns>
     public double ToDouble(int value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// Int32 values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<int> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// Int32 values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<int> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
-    public bool SupportsCpuAcceleration => true;
+    public bool SupportsCpuAcceleration => false;
+
 
     /// <inheritdoc/>
     public bool SupportsGpuAcceleration => true;
diff --git a/src/AiDotNet.Tensors/NumericOperations/Int64Operations.cs b/src/AiDotNet.Tensors/NumericOperations/Int64Operations.cs
index 6bdb642240..b4ffbf7e98 100644
--- a/src/AiDotNet.Tensors/NumericOperations/Int64Operations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/Int64Operations.cs
@@ -787,6 +787,26 @@ public long SignOrZero(long value)
     /// <returns>The value as a double.</returns>
     public double ToDouble(long value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// Int64 values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<long> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// Int64 values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<long> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => true;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/MultivectorOperations.cs b/src/AiDotNet.Tensors/NumericOperations/MultivectorOperations.cs
index 2480118642..474df52e97 100644
--- a/src/AiDotNet.Tensors/NumericOperations/MultivectorOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/MultivectorOperations.cs
@@ -315,6 +315,32 @@ public double ToDouble(Multivector<T> value)
         return _ops.ToDouble(value.Scalar);
     }
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<Multivector<T>> x, out int badIndex)
+    {
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (IsNaN(x[i]) || IsInfinity(x[i]))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<Multivector<T>> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     public bool SupportsCpuAcceleration => false;
 
     public bool SupportsGpuAcceleration => false;
diff --git a/src/AiDotNet.Tensors/NumericOperations/OctonionOperations.cs b/src/AiDotNet.Tensors/NumericOperations/OctonionOperations.cs
index f0fd95adfc..16538e6c30 100644
--- a/src/AiDotNet.Tensors/NumericOperations/OctonionOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/OctonionOperations.cs
@@ -315,6 +315,32 @@ public double ToDouble(Octonion<T> value)
         return _ops.ToDouble(value.Scalar);
     }
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<Octonion<T>> x, out int badIndex)
+    {
+        for (int i = 0; i < x.Length; i++)
+        {
+            if (IsNaN(x[i]) || IsInfinity(x[i]))
+            {
+                badIndex = i;
+                return false;
+            }
+        }
+
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<Octonion<T>> x, out int badIndex)
+    {
+        return !AllFinite(x, out badIndex);
+    }
+
     public bool SupportsCpuAcceleration => false;
 
     public bool SupportsGpuAcceleration => false;
diff --git a/src/AiDotNet.Tensors/NumericOperations/SByteOperations.cs b/src/AiDotNet.Tensors/NumericOperations/SByteOperations.cs
index 4e8b7ed756..852073b407 100644
--- a/src/AiDotNet.Tensors/NumericOperations/SByteOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/SByteOperations.cs
@@ -738,6 +738,26 @@ public class SByteOperations : INumericOperations<sbyte>
     /// </summary>
     public double ToDouble(sbyte value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// SByte values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<sbyte> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// SByte values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<sbyte> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/ShortOperations.cs b/src/AiDotNet.Tensors/NumericOperations/ShortOperations.cs
index 93bc453c3a..79b65b1fa0 100644
--- a/src/AiDotNet.Tensors/NumericOperations/ShortOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/ShortOperations.cs
@@ -699,6 +699,26 @@ public short SignOrZero(short value)
     /// </summary>
     public double ToDouble(short value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// Short values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<short> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// Short values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<short> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/UInt16Operations.cs b/src/AiDotNet.Tensors/NumericOperations/UInt16Operations.cs
index 4d95eb99f9..51cf5ce004 100644
--- a/src/AiDotNet.Tensors/NumericOperations/UInt16Operations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/UInt16Operations.cs
@@ -703,6 +703,26 @@ public class UInt16Operations : INumericOperations<ushort>
     /// </summary>
     public double ToDouble(ushort value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// UInt16 values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<ushort> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// UInt16 values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<ushort> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/UInt32Operations.cs b/src/AiDotNet.Tensors/NumericOperations/UInt32Operations.cs
index 33a6cf416b..de49061e76 100644
--- a/src/AiDotNet.Tensors/NumericOperations/UInt32Operations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/UInt32Operations.cs
@@ -712,6 +712,26 @@ public class UInt32Operations : INumericOperations<uint>
     /// </summary>
     public double ToDouble(uint value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// UInt32 values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<uint> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// UInt32 values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<uint> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/UInt64Operations.cs b/src/AiDotNet.Tensors/NumericOperations/UInt64Operations.cs
index f25e5035c1..edb88f4859 100644
--- a/src/AiDotNet.Tensors/NumericOperations/UInt64Operations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/UInt64Operations.cs
@@ -761,6 +761,26 @@ public ulong FromHalf(Half value)
     /// </summary>
     public double ToDouble(ulong value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// UInt64 values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<ulong> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// UInt64 values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<ulong> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/src/AiDotNet.Tensors/NumericOperations/UIntOperations.cs b/src/AiDotNet.Tensors/NumericOperations/UIntOperations.cs
index 12f2a9b669..8a18bafe48 100644
--- a/src/AiDotNet.Tensors/NumericOperations/UIntOperations.cs
+++ b/src/AiDotNet.Tensors/NumericOperations/UIntOperations.cs
@@ -722,6 +722,26 @@ public uint SignOrZero(uint value)
     /// </summary>
     public double ToDouble(uint value) => (double)value;
 
+    /// <summary>
+    /// Checks if all elements in the span are finite (neither NaN nor Infinity).
+    /// UInt values are always finite.
+    /// </summary>
+    public bool AllFinite(ReadOnlySpan<uint> x, out int badIndex)
+    {
+        badIndex = -1;
+        return true;
+    }
+
+    /// <summary>
+    /// Checks if any element in the span is NaN or Infinity.
+    /// UInt values are always finite.
+    /// </summary>
+    public bool IsAnyNonFinite(ReadOnlySpan<uint> x, out int badIndex)
+    {
+        badIndex = -1;
+        return false;
+    }
+
     /// <inheritdoc/>
     public bool SupportsCpuAcceleration => false;
 
diff --git a/tests/AiDotNet.Tensors.Benchmarks/CpuMatMulDiagnostics.cs b/tests/AiDotNet.Tensors.Benchmarks/CpuMatMulDiagnostics.cs
new file mode 100644
index 0000000000..e8b5a50316
--- /dev/null
+++ b/tests/AiDotNet.Tensors.Benchmarks/CpuMatMulDiagnostics.cs
@@ -0,0 +1,97 @@
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using AiDotNet.Tensors.Engines;
+using AiDotNet.Tensors.Helpers;
+using AiDotNet.Tensors.LinearAlgebra;
+using AiDotNet.Tensors.Operators;
+using AiDotNet.Tensors.Benchmarks.Helpers;
+
+namespace AiDotNet.Tensors.Benchmarks;
+
+/// <summary>
+/// Provides diagnostic tools and performance benchmarks for CPU-based matrix multiplication.
+/// </summary>
+/// <remarks>
+/// <para><b>For Beginners:</b> Matrix multiplication is one of the most critical operations in AI. 
+/// This class helps us measure how fast the computer's "brain" (the CPU) can perform these 
+/// calculations and ensures the specialized optimizations are working as intended.</para>
+/// </remarks>
+public static class CpuMatMulDiagnostics
+{
+    /// <summary>
+    /// Runs a comprehensive suite of CPU matrix multiplication benchmarks.
+    /// </summary>
+    /// <remarks>
+    /// <para><b>For Beginners:</b> This method executes various tests on different matrix sizes. 
+    /// It provides a "warmup" phase to let the CPU get ready, and then performs multiple 
+    /// iterations to get an accurate average speed, measured in GFLOPS (Billions of operations per second).</para>
+    /// </remarks>
+    public static void Run()
+    {
+        Console.WriteLine("=== CPU Matrix Multiply Diagnostics ===");
+        Console.WriteLine($"ProcessorCount: {Environment.ProcessorCount}");
+        Console.WriteLine(TensorPrimitivesCore.GetHardwareAccelerationInfo());
+        Console.WriteLine();
+
+        var engine = new CpuEngine();
+        int[] sizes = GetSizes();
+        const int warmupIterations = 1;
+
+        foreach (int size in sizes)
+        {
+            Console.WriteLine($"{size}x{size} MatMul");
+            Console.WriteLine(new string('-', 48));
+
+            var random = RandomHelper.CreateSeededRandom(42);
+            var a = BenchmarkHelper.CreateRandomMatrix(size, size, random);
+            var b = BenchmarkHelper.CreateRandomMatrix(size, size, random);
+
+            for (int i = 0; i < warmupIterations; i++)
+            {
+                _ = engine.MatrixMultiply(a, b);
+            }
+
+            int benchmarkIterations = size >= 2048 ? 2 : 3;
+            double totalMs = 0;
+            double checksum = 0;
+            for (int i = 0; i < benchmarkIterations; i++)
+            {
+                var sw = Stopwatch.StartNew();
+                var result = engine.MatrixMultiply(a, b);
+                sw.Stop();
+
+                totalMs += sw.Elapsed.TotalMilliseconds;
+                checksum += result[0, 0];
+                Console.WriteLine($"  Iteration {i + 1}: {sw.Elapsed.TotalMilliseconds:F2}ms");
+            }
+
+            double avgMs = totalMs / benchmarkIterations;
+            double seconds = avgMs / 1000.0;
+            double flops = 2.0 * size * size * size;
+            double gflops = flops / seconds / 1_000_000_000.0;
+            Console.WriteLine($"  Average: {avgMs:F2}ms  Throughput: {gflops:F2} GFLOPS  Checksum: {checksum:F3}");
+            Console.WriteLine();
+        }
+    }
+
+    private static int[] GetSizes()
+    {
+        string? sizesEnv = Environment.GetEnvironmentVariable("AIDOTNET_CPU_MATMUL_SIZES");
+        if (!string.IsNullOrWhiteSpace(sizesEnv))
+        {
+            var tokens = sizesEnv.Split(new[] { ',', ';', ' ', '|' }, StringSplitOptions.RemoveEmptyEntries);
+            var sizes = new List<int>();
+            foreach (var token in tokens)
+            {
+                if (int.TryParse(token, out int size) && size > 0)
+                    sizes.Add(size);
+            }
+
+            if (sizes.Count > 0)
+                return sizes.ToArray();
+        }
+
+        return new[] { 256, 512, 1024, 2048 };
+    }
+}
diff --git a/tests/AiDotNet.Tensors.Benchmarks/GpuMatMulDiagnostics.cs b/tests/AiDotNet.Tensors.Benchmarks/GpuMatMulDiagnostics.cs
new file mode 100644
index 0000000000..f1a9387297
--- /dev/null
+++ b/tests/AiDotNet.Tensors.Benchmarks/GpuMatMulDiagnostics.cs
@@ -0,0 +1,125 @@
+using System;
+using System.Diagnostics;
+using AiDotNet.Tensors.Engines;
+using AiDotNet.Tensors.Helpers;
+using AiDotNet.Tensors.Interfaces;
+using AiDotNet.Tensors.LinearAlgebra;
+using AiDotNet.Tensors.Benchmarks.Helpers;
+
+namespace AiDotNet.Tensors.Benchmarks;
+
+/// <summary>
+/// Provides diagnostic tools and performance benchmarks for GPU-based matrix multiplication.
+/// </summary>
+/// <remarks>
+/// <para><b>For Beginners:</b> GPUs are incredibly fast at math, but they can sometimes be 
+/// tricky to set up. This class helps us make sure the GPU is giving the exact same 
+/// answers as the CPU, and measures exactly how much faster it is (the speedup).</para>
+/// </remarks>
+public static class GpuMatMulDiagnostics
+{
+    /// <summary>
+    /// Runs a series of GPU correctness and performance tests.
+    /// </summary>
+    /// <remarks>
+    /// <para><b>For Beginners:</b> This method runs the actual tests. It first checks if the 
+    /// GPU is even available, then it tests various math problem sizes to ensure accuracy 
+    /// and measure throughput in GFLOPS (billions of operations per second).</para>
+    /// </remarks>
+    public static void Run()
+    {
+        Console.WriteLine("=== GPU Matrix Multiply Diagnostics ===");
+
+        if (!AiDotNetEngine.AutoDetectAndConfigureGpu())
+        {
+            Console.WriteLine("GPU not available; skipping GPU diagnostics.");
+            return;
+        }
+
+        var gpuEngine = AiDotNetEngine.Current;
+        var cpuEngine = new CpuEngine();
+        Console.WriteLine($"Engine: {gpuEngine.Name}");
+        Console.WriteLine();
+
+        var random = RandomHelper.CreateSeededRandom(123);
+        int[] correctnessSizes = { 128, 256, 512, 1024 };
+
+        foreach (int size in correctnessSizes)
+        {
+            var a = BenchmarkHelper.CreateRandomMatrix(size, size, random);
+            var b = BenchmarkHelper.CreateRandomMatrix(size, size, random);
+
+            Console.WriteLine($"Correctness check (safe): {size}x{size}");
+            SetKernelMode(forceSafe: true, forceUnsafe: false);
+            var safeStats = CompareInternal(cpuEngine, gpuEngine, a, b);
+            PrintStats(safeStats);
+
+            if (size <= 512)
+            {
+                Console.WriteLine($"Correctness check (unsafe): {size}x{size}");
+                SetKernelMode(forceSafe: false, forceUnsafe: true);
+                var unsafeStats = CompareInternal(cpuEngine, gpuEngine, a, b);
+                PrintStats(unsafeStats);
+            }
+
+            Console.WriteLine();
+        }
+
+        SetKernelMode(forceSafe: false, forceUnsafe: false);
+
+        int[] perfSizes = { 1024, 2048 };
+        foreach (int perfSize in perfSizes)
+        {
+            var perfA = BenchmarkHelper.CreateRandomMatrix(perfSize, perfSize, random);
+            var perfB = BenchmarkHelper.CreateRandomMatrix(perfSize, perfSize, random);
+
+            Console.WriteLine($"GPU performance: {perfSize}x{perfSize}");
+            
+            // Warmup pass
+            _ = gpuEngine.MatrixMultiply(perfA, perfB);
+
+            int perfIterations = perfSize >= 2048 ? 2 : 3;
+            double totalMs = 0;
+            for (int i = 0; i < perfIterations; i++)
+            {
+                var sw = Stopwatch.StartNew();
+                _ = gpuEngine.MatrixMultiply(perfA, perfB);
+                sw.Stop();
+                totalMs += sw.Elapsed.TotalMilliseconds;
+                Console.WriteLine($"  Iteration {i + 1}: {sw.Elapsed.TotalMilliseconds:F2}ms");
+            }
+
+            double avgMs = totalMs / perfIterations;
+            double seconds = avgMs / 1000.0;
+            double flops = 2.0 * perfSize * perfSize * perfSize;
+            double gflops = flops / seconds / 1_000_000_000.0;
+            Console.WriteLine($"  Average: {avgMs:F2}ms  Throughput: {gflops:F2} GFLOPS");
+            Console.WriteLine();
+        }
+    }
+
+    private static (double maxError, double avgError, int nonFiniteCount) CompareInternal(
+        CpuEngine cpuEngine,
+        IEngine gpuEngine,
+        Matrix<float> a,
+        Matrix<float> b)
+    {
+        var cpuResult = cpuEngine.MatrixMultiply(a, b);
+        var gpuResult = gpuEngine.MatrixMultiply(a, b);
+
+        return BenchmarkHelper.Compare(cpuResult, gpuResult);
+    }
+
+    private static void PrintStats((double maxError, double avgError, int nonFiniteCount) stats)
+    {
+        Console.WriteLine($"  Max error: {stats.maxError:E3}");
+        Console.WriteLine($"  Avg error: {stats.avgError:E3}");
+        Console.WriteLine($"  Non-finite GPU values: {stats.nonFiniteCount}");
+    }
+
+    private static void SetKernelMode(bool forceSafe, bool forceUnsafe)
+    {
+        Environment.SetEnvironmentVariable("AIDOTNET_GEMM_SAFE", forceSafe ? "1" : "0");
+        Environment.SetEnvironmentVariable("AIDOTNET_GEMM_UNSAFE", forceUnsafe ? "1" : "0");
+    }
+}
\ No newline at end of file
diff --git a/tests/AiDotNet.Tensors.Benchmarks/Helpers/BenchmarkHelper.cs b/tests/AiDotNet.Tensors.Benchmarks/Helpers/BenchmarkHelper.cs
new file mode 100644
index 0000000000..dfa9dff5d1
--- /dev/null
+++ b/tests/AiDotNet.Tensors.Benchmarks/Helpers/BenchmarkHelper.cs
@@ -0,0 +1,83 @@
+using System;
+using AiDotNet.Tensors.LinearAlgebra;
+
+namespace AiDotNet.Tensors.Benchmarks.Helpers;
+
+/// <summary>
+/// Provides shared utility methods for tensor and matrix benchmarks.
+/// </summary>
+/// <remarks>
+/// <para><b>For Beginners:</b> This helper class contains common tasks used in multiple benchmarks, 
+/// such as creating random data and comparing results. Moving this logic here makes the individual 
+/// benchmarks cleaner and easier to maintain.</para>
+/// </remarks>
+public static class BenchmarkHelper
+{
+    /// <summary>
+    /// Creates a matrix of the specified size filled with random float values between -1 and 1.
+    /// </summary>
+    /// <param name="rows">The number of rows in the matrix.</param>
+    /// <param name="cols">The number of columns in the matrix.</param>
+    /// <param name="random">The random number generator to use.</param>
+    /// <returns>A new matrix populated with random values.</returns>
+    /// <remarks>
+    /// <para><b>For Beginners:</b> When testing math operations, we need data to work with. 
+    /// This method fills a grid (matrix) with random numbers so we can test if the GPU 
+    /// and CPU give the same answers even with unpredictable data.</para>
+    /// </remarks>
+    public static Matrix<float> CreateRandomMatrix(int rows, int cols, Random random)
+    {
+        var matrix = new Matrix<float>(rows, cols);
+        var data = matrix.AsWritableSpan();
+        for (int i = 0; i < data.Length; i++)
+        {
+            data[i] = (float)((random.NextDouble() * 2.0) - 1.0);
+        }
+
+        return matrix;
+    }
+
+    /// <summary>
+    /// Compares two result matrices and calculates error statistics.
+    /// </summary>
+    /// <param name="reference">The authoritative reference matrix (usually from the CPU).</param>
+    /// <param name="actual">The matrix to validate (usually from the GPU).</param>
+    /// <returns>A tuple containing maximum error, average error, and count of non-finite values.</returns>
+    /// <remarks>
+    /// <para><b>For Beginners:</b> Computers can sometimes make small mistakes in rounding, or big mistakes 
+    /// if a component isn't working right. This method checks every single number in two results 
+    /// to see how different they are. It helps us guarantee that our fast GPU math is just as 
+    /// accurate as our reliable CPU math.</para>
+    /// </remarks>
+    public static (double maxError, double avgError, int nonFiniteCount) Compare(
+        Matrix<float> reference,
+        Matrix<float> actual)
+    {
+        var refSpan = reference.AsSpan();
+        var actSpan = actual.AsSpan();
+
+        double maxError = 0;
+        double sumError = 0;
+        int count = 0;
+        int nonFiniteCount = 0;
+
+        for (int i = 0; i < actSpan.Length; i++)
+        {
+            float actVal = actSpan[i];
+            if (float.IsNaN(actVal) || float.IsInfinity(actVal))
+            {
+                nonFiniteCount++;
+                continue;
+            }
+
+            double error = Math.Abs(refSpan[i] - actVal);
+            sumError += error;
+            if (error > maxError)
+                maxError = error;
+            count++;
+        }
+
+        double avgError = count > 0 ? sumError / count : double.NaN;
+        return (maxError, avgError, nonFiniteCount);
+    }
+}
diff --git a/tests/AiDotNet.Tensors.Benchmarks/Program.cs b/tests/AiDotNet.Tensors.Benchmarks/Program.cs
index 8b823e32d7..b55af1aced 100644
--- a/tests/AiDotNet.Tensors.Benchmarks/Program.cs
+++ b/tests/AiDotNet.Tensors.Benchmarks/Program.cs
@@ -29,6 +29,18 @@ static void Main(string[] args)
             return;
         }
 
+        if (args[0] == "--cpu-matmul")
+        {
+            CpuMatMulDiagnostics.Run();
+            return;
+        }
+
+        if (args[0] == "--gpu-matmul")
+        {
+            GpuMatMulDiagnostics.Run();
+            return;
+        }
+
 #if !NET462
         // Run cuBLAS vs DirectGpu GEMM benchmark
         if (args[0] == "--cublas")
@@ -63,6 +75,8 @@ static void Main(string[] args)
         Console.WriteLine("  --quick    : Run quick performance validation (default)");
         Console.WriteLine("  --full     : Run full BenchmarkDotNet suite (trigonometric)");
         Console.WriteLine("  --linalg   : Run linear algebra benchmarks vs MathNet.Numerics");
+        Console.WriteLine("  --cpu-matmul : Run CPU matmul diagnostics (timed iterations)");
+        Console.WriteLine("  --gpu-matmul : Run GPU matmul diagnostics (correctness + perf)");
 #if !NET462
         Console.WriteLine("  --cublas   : Run cuBLAS vs DirectGpu GEMM benchmark");
         Console.WriteLine("  --opencl   : Run OpenCL GEMM benchmark (AMD/Intel GPUs)");
diff --git a/tests/AiDotNet.Tests/DirectGpuTests.cs b/tests/AiDotNet.Tests/DirectGpuTests.cs
index 9b7a4dcb1d..84ed9b8a3a 100644
--- a/tests/AiDotNet.Tests/DirectGpuTests.cs
+++ b/tests/AiDotNet.Tests/DirectGpuTests.cs
@@ -649,6 +649,7 @@ public void DirectGpuEngine_NewKernels_Correctness()
         // Test each new kernel against the reference
         var newKernels = new (string Name, Action<IGpuBuffer, IGpuBuffer, IGpuBuffer, int, int, int> Execute)[]
         {
+            ("gemm_double_buffered", (a, b, c, m, n, k) => backend.GemmDoubleBuffered(a, b, c, m, n, k)),
             ("gemm_kreg4", (a, b, c, m, n, k) => backend.GemmKreg4(a, b, c, m, n, k)),
             ("gemm_prefetch", (a, b, c, m, n, k) => backend.GemmPrefetch(a, b, c, m, n, k)),
             ("gemm_wide_vec", (a, b, c, m, n, k) => backend.GemmWideVec(a, b, c, m, n, k)),