Parallelize Winograd GEMM stage and fuse bias+activation

* use activation + blas-fused + b-packed matmul
* parallelize that heavy matmul loop when applicable

Author: Human9000-bit

crates/yscv-kernels/src/ops/conv/gemm_conv.rs

@@ -555,6 +555,8 @@ fn winograd_conv2d_nhwc(
     pad_right: usize,
     activation: Activation,
 ) -> Result<Tensor, KernelError> {
+    use rayon::iter::{IntoParallelIterator, ParallelIterator};
+
     let padded_h = in_h + pad_top + pad_bottom;
     let padded_w = in_w + pad_left + pad_right;
     let out_h = padded_h - 2; // (padded_h - 3) / 1 + 1
@@ -573,6 +575,8 @@ fn winograd_conv2d_nhwc(
     let mut output = AlignedVec::<f32>::uninitialized(batch * out_h * out_w * c_out);
 
     for b in 0..batch {
+        use crate::{GemmEpilogue, core::scope_ctx::par_chunks_mut_dispatch};
+
         let in_batch = &input[b * in_h * in_w * c_in..(b + 1) * in_h * in_w * c_in];
 
         // 2. Input transform: for each tile, for each channel, compute B^T * d * B
@@ -611,11 +615,44 @@ fn winograd_conv2d_nhwc(
         //    V[alpha]: [n_tiles, c_in], U[alpha]: [c_in, c_out]
         //    M[alpha]: [n_tiles, c_out]
         let mut m_buf = vec![0.0f32; 16 * n_tiles * c_out];
-        for a in 0..16 {
-            let v_slice = &v[a * n_tiles * c_in..(a + 1) * n_tiles * c_in];
-            let u_slice = &u[a * c_in * c_out..(a + 1) * c_in * c_out];
-            let m_slice = &mut m_buf[a * n_tiles * c_out..(a + 1) * n_tiles * c_out];
-            super::super::matmul::blas_sgemm(v_slice, u_slice, m_slice, n_tiles, c_in, c_out);
+        let epilogue = GemmEpilogue {
+            activation,
+            bias: bias.map(|b| b.as_ptr()),
+            residual: None,
+        };
+        let config = ParallelMatmulConfig::default();
+        if should_parallelize_len(m_buf.len(), config.min_parallel_output_elements, None) {
+            let packed_u: Vec<_> = (0..16)
+                .into_par_iter()
+                .map(|a| {
+                    use crate::pack_b_for_session;
+
+                    let u_slice = &u[a * c_in * c_out..(a + 1) * c_in * c_out];
+                    pack_b_for_session(u_slice, c_in, c_out)
+                })
+                .collect();
+            par_chunks_mut_dispatch(&mut m_buf, n_tiles * c_out, |a, chunk| {
+                use crate::matmul_2d_slices_fused_maybe_packed;
+
+                let v_slice = &v[a * n_tiles * c_in..(a + 1) * n_tiles * c_in];
+                let u_slice = &u[a * c_in * c_out..(a + 1) * c_in * c_out];
+                let packed = Some(packed_u[a].as_ref());
+                matmul_2d_slices_fused_maybe_packed(
+                    v_slice, n_tiles, c_in, u_slice, c_out, chunk, packed, epilogue, config, None,
+                );
+            });
+        } else {
+            for a in 0..16 {
+                use crate::matmul_2d_slices_fused_maybe_packed;
+
+                let v_slice = &v[a * n_tiles * c_in..(a + 1) * n_tiles * c_in];
+                let u_slice = &u[a * c_in * c_out..(a + 1) * c_in * c_out];
+                let m_slice = &mut m_buf[a * n_tiles * c_out..(a + 1) * n_tiles * c_out];
+
+                matmul_2d_slices_fused_maybe_packed(
+                    v_slice, n_tiles, c_in, u_slice, c_out, m_slice, None, epilogue, config, None,
+                );
+            }
         }
 
         // 4. Output transform: A^T * M * A → 2×2 output per tile, with bias + activation
@@ -638,24 +675,16 @@ fn winograd_conv2d_nhwc(
                     let mut otile = [0.0f32; 4];
                     winograd_output_tile(&mt, &mut otile);
 
-                    // Add bias
-                    let bias_val = bias.map_or(0.0, |bd| bd[co]);
+                    // Store input tile
                     for dy in 0..valid_h {
                         for dx in 0..valid_w {
                             let idx = (oy0 + dy) * out_w * c_out + (ox0 + dx) * c_out + co;
-                            out_batch[idx] = otile[dy * 2 + dx] + bias_val;
+                            out_batch[idx] = otile[dy * 2 + dx];
                         }
                     }
                 }
             }
         }
-
-        // Apply activation on the whole batch output
-        match activation {
-            Activation::Silu => silu_slice_inplace(out_batch),
-            Activation::Relu => relu_slice_inplace(out_batch),
-            Activation::None => {}
-        }
     }
 
     Tensor::from_aligned(vec![batch, out_h, out_w, c_out], output).map_err(Into::into)