diff --git a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_small_m.wgsl.template b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_k_parallel.wgsl.template
similarity index 82%
rename from onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_small_m.wgsl.template
rename to onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_k_parallel.wgsl.template
index a53e62921b427..2bed945685d16 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_small_m.wgsl.template
+++ b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_k_parallel.wgsl.template
@@ -1,6 +1,15 @@
 // Copyright (c) Microsoft Corporation. All rights reserved.
 // Licensed under the MIT License.
 
+// K-parallel DP4A projection: each workgroup computes one input row and
+// tile_size output columns. Its 128 threads process four columns concurrently,
+// using 32 K lanes per column and 32 elements per lane in each K iteration.
+// The group stages 1024 activations and their scales in workgroup memory,
+// reusing them across output columns, including additional four-column groups.
+// Each lane accumulates scaled integer dot products; the 32 lane partials are
+// reduced through inter_results to produce each output value after processing K.
+// This exposes parallelism along K when small M or N yields few output tiles.
+
 #param acc_f32
 #param tile_size
 #param tile_size_k_vec
@@ -15,35 +24,13 @@
 
 #use .getByOffset .setByOffset
 
-
 #include "quantization/dp4a_matmul_common.wgsl.template"
 
-// This algorithm works to compute dot product of k in parallel, by processing k at each step amongst tile_size_k_vec threads,
-// and utilizing the remaining threads in the workgroup to process additional rows of b in parallel (such that the values in shared memory for A can be reused).
-// For each load of k, the tile_size_k_vec threads also reload B tile_size/num_concurrent_b_rows times to compute partial dot products of other B rows
-// in order to complete all tile_size b rows in this workgroup and also reusing the loaded in register values of a.
-
-// 1. Each workgroup handles tile_size_k_vec * k_vectorization_in_b (32) columns and num_concurrent_b_rows of matrix B at a time,
-// iterating over the columns to compute a partial dot product.
-// 2. Uses vec4 vectorization where each K represents 32 elements of matrix B
-
-// 1. Workgroup Responsibility:
-//    - Processes one row of matrix A
-//    - Handles tile_size rows of matrix B
-//
-// 2. Computation Process:
-//    - Reads [tile_size][tile_size_k_vec] block of B data at a time
-//    - Each thread within workgroup computes dot products of 32 A*B elements since each K represents 32 elements of matrix B
-//    - Stores intermediate results in shared memory (inter_results)
-//    - Iterates through columns accumulating results in inter_results
-//    - Performs final reduction sum in inter_results for output
-// sub_tile_count is the number of concurrent b rows processed by the workgroup.
-
 const double_tile_size_k_vec = 2 * tile_size_k_vec;
 
-// Precision of the cross-tile accumulator, selected by the
-// "enableMatmulFp32Accumulation" provider option. SDP8AI itself is an exact integer
-// dot product; it is the accumulation over K that can overflow the f16 max (~65504).
+// Narrow outputs use FP32 cross-tile accumulation to limit reduction error;
+// other outputs follow "enableMatmulFp32Accumulation". SDP8AI itself is an exact
+// integer dot product; accumulation over K can overflow the f16 max (~65504).
 #if acc_f32
 alias acc_element_t = f32;
 #else
@@ -64,7 +51,7 @@ fn loadSHMA(batch: u32, a_global: u32, kidx_v: u32, col: u32)
     }
 
     tile_A[col] = a.getByOffset(batch*uniforms.M*uniforms.K16+a_global*uniforms.K16+k_offset);
-    if (col < scale_a_size_in_tile_a)
+    if (col < scale_a_size_in_tile_a && kidx_v/8 + col < uniforms.K16/8)
     {
     // kidx_v - covers 16 values of k in input_a
     scale_A[col] = scales_a.getByOffset(batch*uniforms.M*(uniforms.K/128) + a_global*(uniforms.K/128) + kidx_v/8 + col);
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.cc b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.cc
index d9c5f5839a3bf..89bb45f513b3c 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.cc
+++ b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.cc
@@ -10,13 +10,12 @@ namespace contrib {
 namespace webgpu {
 
 Status DP4AMatMulQuantizeProgram::GenerateShaderCode(ShaderHelper& shader) const {
-  const auto& a = shader.AddInput("input_a", ShaderUsage::UseUniform | ShaderUsage::UseIndicesTypeAlias | ShaderUsage::UseValueTypeAlias | ShaderUsage::UseElementTypeAlias);
+  const auto& a = shader.AddInput("input_a", ShaderUsage::UseUniform | ShaderUsage::UseIndicesTypeAlias |
+                                                 ShaderUsage::UseValueTypeAlias | ShaderUsage::UseElementTypeAlias);
   const auto& output = shader.AddOutput("output", ShaderUsage::UseUniform);
   const auto& scales = shader.AddOutput("scales", ShaderUsage::UseUniform);
-  return WGSL_TEMPLATE_APPLY(shader, "quantization/dp4a_quantize.wgsl.template",
-                             WGSL_TEMPLATE_VARIABLE(a, a),
-                             WGSL_TEMPLATE_VARIABLE(output, output),
-                             WGSL_TEMPLATE_VARIABLE(scales, scales));
+  return WGSL_TEMPLATE_APPLY(shader, "quantization/dp4a_quantize.wgsl.template", WGSL_TEMPLATE_VARIABLE(a, a),
+                             WGSL_TEMPLATE_VARIABLE(output, output), WGSL_TEMPLATE_VARIABLE(scales, scales));
 }
 
 Status DP4AMatMulNBitsProgram::GenerateShaderCode(ShaderHelper& shader) const {
@@ -52,7 +51,7 @@ Status DP4AMatMulNBitsProgram::GenerateShaderCode(ShaderHelper& shader) const {
 }
 
 // scale_A components = 1, b components = 4, output components = 1
-Status DP4AMatMulNBitsSmallMProgram::GenerateShaderCode(ShaderHelper& shader) const {
+Status DP4AMatMulNBitsKParallelProgram::GenerateShaderCode(ShaderHelper& shader) const {
   const auto& a = shader.AddInput("input_a", ShaderUsage::UseUniform);
   const auto& scales_a = shader.AddInput("scales_a", ShaderUsage::UseUniform);
   const auto& b = shader.AddInput("input_b", ShaderUsage::UseUniform);
@@ -72,7 +71,7 @@ Status DP4AMatMulNBitsSmallMProgram::GenerateShaderCode(ShaderHelper& shader) co
   const uint32_t sub_tile_count = WorkgroupSizeX() / tile_size_k_vec_;
   ORT_ENFORCE(tile_size_ % sub_tile_count == 0, "tile_size_ must be divisible by sub_tile_count");
 
-  return WGSL_TEMPLATE_APPLY(shader, "quantization/dp4a_matmul_small_m.wgsl.template",
+  return WGSL_TEMPLATE_APPLY(shader, "quantization/dp4a_matmul_k_parallel.wgsl.template",
                              WGSL_TEMPLATE_PARAMETER(acc_f32, acc_f32_),
                              WGSL_TEMPLATE_PARAMETER(broadcast_a_row, broadcast_a_row_),
                              WGSL_TEMPLATE_PARAMETER(has_bias, has_bias_),
@@ -106,51 +105,74 @@ Status ApplyDP4AMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Tensor
                                   onnxruntime::webgpu::ComputeContext& context,
                                   Tensor* y,
                                   const uint32_t weight_index,
-                                  const Tensor* weight_index_indirect) {
+                                  const Tensor* weight_index_indirect, DP4AQuantizedInput* quantized_input) {
   constexpr uint32_t kVec4Components = 4;
   constexpr uint32_t kVec2Components = 2;
   constexpr uint32_t kU32Components = 4;
 
-  constexpr uint32_t kBlockSizeA = 128;
-  DP4AMatMulQuantizeProgram quantize_program;
-  quantize_program.SetWorkgroupSize(64);
-  uint32_t tile_size = 64 * kVec4Components;
-  quantize_program.SetDispatchGroupSize((batch_count * M * K + tile_size - 1) / tile_size, 1, 1);
-  TensorShape a_quant_shape{batch_count, M, K / kU32Components};
-  Tensor a_quant = context.CreateGPUTensor(DataTypeImpl::GetType<uint32_t>(), a_quant_shape);
-  TensorShapeVector a_scales_dims({batch_count, 1, M, K / kBlockSizeA});
-  Tensor a_scale = context.CreateGPUTensor(a->DataType(), a_scales_dims);
-  quantize_program.AddInputs({{a, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(kVec4Components)}})
-      .AddOutputs({{&a_quant, ProgramTensorMetadataDependency::Rank, a_quant.Shape(), 1},
-                   {&a_scale, ProgramTensorMetadataDependency::Rank, 1}})
-      .AddUniformVariable({batch_count * M * K / kU32Components});
-  ORT_RETURN_IF_ERROR(context.RunProgram(quantize_program));
+  DP4AQuantizedInput local_input;
+  auto& input = quantized_input ? *quantized_input : local_input;
+  if (!input.values) {
+    constexpr uint32_t kBlockSizeA = 128;
+    DP4AMatMulQuantizeProgram quantize_program;
+    quantize_program.SetWorkgroupSize(64);
+    uint32_t tile_size = 64 * kVec4Components;
+    quantize_program.SetDispatchGroupSize((batch_count * M * K + tile_size - 1) / tile_size, 1, 1);
+    TensorShape a_quant_shape{batch_count, M, K / kU32Components};
+    Tensor a_quant = context.CreateGPUTensor(DataTypeImpl::GetType<uint32_t>(), a_quant_shape);
+    TensorShapeVector a_scales_dims({batch_count, 1, M, K / kBlockSizeA});
+    Tensor a_scale = context.CreateGPUTensor(a->DataType(), a_scales_dims);
+    quantize_program.AddInputs({{a, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(kVec4Components)}})
+        .AddOutputs({{&a_quant, ProgramTensorMetadataDependency::Rank, a_quant.Shape(), 1},
+                     {&a_scale, ProgramTensorMetadataDependency::Rank, 1}})
+        .AddUniformVariable({batch_count * M * K / kU32Components});
+    ORT_RETURN_IF_ERROR(context.RunProgram(quantize_program));
+    input.values.emplace(std::move(a_quant));
+    input.scales.emplace(std::move(a_scale));
+  }
+  const auto& a_quant = *input.values;
+  const auto& a_scale = *input.scales;
 
   const bool has_zero_points = zero_points != nullptr;
   const bool has_bias = bias != nullptr;
   const bool has_weight_idx_indirect = weight_index_indirect != nullptr;
   const bool has_weight_idx = weight_index != 0 || has_weight_idx_indirect;
   const bool single_scale_weights = (block_size == K * N);
-  if (has_weight_idx_indirect || M < min_M_for_tile_optimization) {
-    uint32_t tile_size_k_vec = 32;
-    uint32_t tile_size_n = 4;
+  const bool small_n = N <= 32;
+  // Each K-parallel iteration covers 32 lanes * 32 elements = 1024 K values.
+  // For narrow N, conservatively require two iterations to amortize reduction
+  // overhead; shorter K benefits more from the tiled kernel's reuse across M.
+  const bool use_k_parallel = M < min_M_for_tile_optimization || (small_n && K >= 2048);
+  // Indirect expert weights can differ by row, so they also require the row-wise kernel.
+  if (use_k_parallel || has_weight_idx_indirect) {
+    const uint32_t tile_size_k_vec = 32;
+    // Keep more workgroups available when there are only a few input rows.
+    const uint32_t tile_size_n = small_n && M >= min_M_for_tile_optimization ? N : 4;
 
     const uint32_t b_components = (nbits == 2 ? kVec2Components : kVec4Components);
     const bool broadcast_a = dispatch_M > M;
-    const bool acc_f32 = context.EnableMatmulFp32Accumulation();
-    DP4AMatMulNBitsSmallMProgram mul_program{tile_size_k_vec, tile_size_n, nbits, has_zero_points, has_bias, has_weight_idx, has_weight_idx_indirect, single_scale_weights, broadcast_a, acc_f32};
+    // The narrow projection changes the K reduction order. Keep partial sums
+    // and the final reduction in FP32, then cast once when storing the output.
+    const bool acc_f32 = small_n || context.EnableMatmulFp32Accumulation();
+    DP4AMatMulNBitsKParallelProgram mul_program{
+        tile_size_k_vec, tile_size_n, nbits, has_zero_points, has_bias, has_weight_idx,
+        has_weight_idx_indirect, single_scale_weights, broadcast_a, acc_f32};
     uint32_t num_N_tile = (N + tile_size_n - 1) / tile_size_n;
     mul_program.SetWorkgroupSize(128);
     mul_program.SetDispatchGroupSize(batch_count * dispatch_M * num_N_tile);
-    mul_program.AddInputs({{&a_quant, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(kVec4Components)},
-                           {&a_scale, ProgramTensorMetadataDependency::TypeAndRank, 1},
-                           {b, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(b_components * kU32Components)},
-                           {scales, ProgramTensorMetadataDependency::TypeAndRank, 1}})
-        .AddUniformVariables({batch_count, M, N, K, K / 16, K / 32, block_size, num_N_tile, zero_blocks_per_col, weight_index, dispatch_M})
+    mul_program
+        .AddInputs({{&a_quant, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(kVec4Components)},
+                    {&a_scale, ProgramTensorMetadataDependency::TypeAndRank, 1},
+                    {b, ProgramTensorMetadataDependency::TypeAndRank, static_cast<int>(b_components * kU32Components)},
+                    {scales, ProgramTensorMetadataDependency::TypeAndRank, 1}})
+        .AddUniformVariables({batch_count, M, N, K, K / 16, K / 32, block_size, num_N_tile, zero_blocks_per_col,
+                              weight_index, dispatch_M})
         .AddOutput({y, ProgramTensorMetadataDependency::TypeAndRank, 1})
-        .CacheHint(nbits, tile_size_k_vec, tile_size_n, has_zero_points, single_scale_weights, has_bias, has_weight_idx, has_weight_idx_indirect, broadcast_a, acc_f32);
+        .CacheHint(nbits, tile_size_k_vec, tile_size_n, has_zero_points, single_scale_weights, has_bias, has_weight_idx,
+                   has_weight_idx_indirect, broadcast_a, acc_f32);
     if (has_zero_points) {
-      mul_program.AddInput({zero_points, ProgramTensorMetadataDependency::None, {(zero_points->Shape().Size() + 3) / 4}, 4});
+      mul_program.AddInput(
+          {zero_points, ProgramTensorMetadataDependency::None, {(zero_points->Shape().Size() + 3) / 4}, 4});
     }
     if (has_bias) {
       mul_program.AddInput({bias, ProgramTensorMetadataDependency::None});
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.h b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.h
index 66dc5819e37a8..e627fc7d60dcd 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.h
+++ b/onnxruntime/contrib_ops/webgpu/quantization/dp4a_matmul_nbits.h
@@ -4,6 +4,7 @@
 #pragma once
 
 #include <limits>
+#include <optional>
 
 #include "core/providers/webgpu/program.h"
 #include "core/providers/webgpu/webgpu_kernel.h"
@@ -14,6 +15,12 @@ namespace webgpu {
 
 using namespace onnxruntime::webgpu;
 
+// Invocation-local workspace, reusable only by projections of the same input tensor.
+struct DP4AQuantizedInput {
+  std::optional<Tensor> values;
+  std::optional<Tensor> scales;
+};
+
 class DP4AMatMulQuantizeProgram final : public Program<DP4AMatMulQuantizeProgram> {
  public:
   DP4AMatMulQuantizeProgram() : Program{"DP4AMatMulQuantize"} {}
@@ -59,23 +66,23 @@ class DP4AMatMulNBitsProgram final : public Program<DP4AMatMulNBitsProgram> {
   bool acc_f32_;
 };
 
-class DP4AMatMulNBitsSmallMProgram final : public Program<DP4AMatMulNBitsSmallMProgram> {
+class DP4AMatMulNBitsKParallelProgram final : public Program<DP4AMatMulNBitsKParallelProgram> {
  public:
-  DP4AMatMulNBitsSmallMProgram(uint32_t tile_size_k_vec, uint32_t tile_size, uint32_t nbits,
-                               bool has_zero_points, bool has_bias,
-                               bool has_weight_idx, bool has_weight_idx_indirect, bool single_scale_weights,
-                               bool broadcast_a_row = false,
-                               bool acc_f32 = false) : Program{"DP4AMatMulNBitsSmallMProgram"},
-                                                       tile_size_k_vec_(tile_size_k_vec),
-                                                       tile_size_(tile_size),
-                                                       nbits_(nbits),
-                                                       has_bias_(has_bias),
-                                                       has_zero_points_(has_zero_points),
-                                                       has_weight_idx_(has_weight_idx),
-                                                       has_weight_idx_indirect_(has_weight_idx_indirect),
-                                                       single_scale_weights_(single_scale_weights),
-                                                       broadcast_a_row_(broadcast_a_row),
-                                                       acc_f32_(acc_f32) {}
+  DP4AMatMulNBitsKParallelProgram(uint32_t tile_size_k_vec, uint32_t tile_size, uint32_t nbits,
+                                  bool has_zero_points, bool has_bias,
+                                  bool has_weight_idx, bool has_weight_idx_indirect, bool single_scale_weights,
+                                  bool broadcast_a_row = false,
+                                  bool acc_f32 = false) : Program{"DP4AMatMulNBitsKParallelProgram"},
+                                                          tile_size_k_vec_(tile_size_k_vec),
+                                                          tile_size_(tile_size),
+                                                          nbits_(nbits),
+                                                          has_bias_(has_bias),
+                                                          has_zero_points_(has_zero_points),
+                                                          has_weight_idx_(has_weight_idx),
+                                                          has_weight_idx_indirect_(has_weight_idx_indirect),
+                                                          single_scale_weights_(single_scale_weights),
+                                                          broadcast_a_row_(broadcast_a_row),
+                                                          acc_f32_(acc_f32) {}
   Status GenerateShaderCode(ShaderHelper& sh) const override;
   WEBGPU_PROGRAM_DEFINE_UNIFORM_VARIABLES(
       {"batch_count", ProgramUniformVariableDataType::Uint32},
@@ -117,22 +124,18 @@ Status ApplyDP4AMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Tensor
                                   onnxruntime::webgpu::ComputeContext& context,
                                   Tensor* y,
                                   const uint32_t weight_index,
-                                  const Tensor* weight_index_indirect = nullptr);
+                                  const Tensor* weight_index_indirect = nullptr,
+                                  DP4AQuantizedInput* quantized_input = nullptr);
 
 // The optional M / has_weight_idx_indirect / y arguments fold the original
 // dispatch-precondition (DP4A is preferred when M is large enough, or
 // unconditionally on FP32 outputs and Qualcomm GPUs) into the feasibility check
 // so callers don't need a separate wrapper. Defaults make the precondition
 // trivially satisfied for callers that only want the feasibility check.
-bool CanApplyDP4AMatrixMatMulNBits(onnxruntime::webgpu::ComputeContext& context,
-                                   uint64_t accuracy_level,
-                                   uint32_t block_size,
-                                   uint32_t N,
-                                   uint32_t K,
-                                   uint32_t components_k,
+bool CanApplyDP4AMatrixMatMulNBits(onnxruntime::webgpu::ComputeContext& context, uint64_t accuracy_level,
+                                   uint32_t block_size, uint32_t N, uint32_t K, uint32_t components_k,
                                    uint32_t M = std::numeric_limits<uint32_t>::max(),
-                                   bool has_weight_idx_indirect = false,
-                                   const Tensor* y = nullptr);
+                                   bool has_weight_idx_indirect = false, const Tensor* y = nullptr);
 
 }  // namespace webgpu
 }  // namespace contrib
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.cc b/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.cc
index 404f06f858869..70762af0d163e 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.cc
+++ b/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.cc
@@ -66,7 +66,7 @@ Status MatMulNBitsWideTileProgram::GenerateShaderCode(ShaderHelper& shader) cons
                              WGSL_TEMPLATE_VARIABLE(scales, scales));
 }
 
-// Apply similar idea with DP4AMatMulNBitsSmallMProgram algorithm.
+// Uses the same K-parallel reduction strategy as DP4AMatMulNBitsKParallelProgram.
 Status MatMulNBitsProgram::GenerateShaderCode(ShaderHelper& shader) const {
   const auto& a = shader.AddInput("input_a", ShaderUsage::UseValueTypeAlias);
   const auto& b = shader.AddInput("input_b");
@@ -293,7 +293,7 @@ Status ApplyMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
                         Tensor* y,
                         const uint32_t weight_index,
                         const Tensor* weight_index_indirect,
-                        uint32_t override_M) {
+                        uint32_t override_M, DP4AQuantizedInput* quantized_input) {
   TensorShape b_shape({N_op, K_op});
   MatMulComputeHelper helper;
   ORT_RETURN_IF_ERROR(helper.Compute(a->Shape(), b_shape, false, true));
@@ -353,7 +353,7 @@ Status ApplyMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
   // DP4A Q2 path now supports custom zero points via a 1024-entry LUT (4 zero-point sections × 256 byte values).
   if (CanApplyDP4AMatrixMatMulNBits(context, accuracy_level, block_size, N, K, components_a,
                                     M, has_weight_idx_indirect, y)) {
-    return ApplyDP4AMatrixMatMulNBits(a, b, scales, zero_points, bias, batch_count, M, dispatch_M, N, K, block_size, zero_blocks_per_col, kMinMForTileOptimization, static_cast<uint32_t>(nbits), context, y, weight_index, weight_index_indirect);
+    return ApplyDP4AMatrixMatMulNBits(a, b, scales, zero_points, bias, batch_count, M, dispatch_M, N, K, block_size, zero_blocks_per_col, kMinMForTileOptimization, static_cast<uint32_t>(nbits), context, y, weight_index, weight_index_indirect, quantized_input);
   }
 
   // This program is optimized for FP32/FP16 Block32 prefill.
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.h b/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.h
index f8c2f91b78396..e2205f37920dd 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.h
+++ b/onnxruntime/contrib_ops/webgpu/quantization/matmul_nbits.h
@@ -12,6 +12,8 @@ namespace webgpu {
 
 using namespace onnxruntime::webgpu;
 
+struct DP4AQuantizedInput;
+
 class MatMulNBitsWideTileProgram final : public Program<MatMulNBitsWideTileProgram> {
  public:
   MatMulNBitsWideTileProgram(bool has_zero_points, bool has_bias, bool has_weight_idx, bool has_weight_idx_indirect, uint32_t tile_m, uint32_t tile_n, uint32_t nbits, uint32_t subgroup_min_size, bool acc_f32)
@@ -107,7 +109,8 @@ Status ApplyMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
                         int64_t K_op, int64_t N_op, int64_t block_size_op, int64_t accuracy_level, int64_t bits_op,
                         onnxruntime::webgpu::ComputeContext& context, Tensor* y, const uint32_t weight_index = 0,
                         const Tensor* weight_index_indirect = nullptr,
-                        uint32_t override_M = 0);
+                        uint32_t override_M = 0,
+                        DP4AQuantizedInput* quantized_input = nullptr);
 
 }  // namespace webgpu
 }  // namespace contrib
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.cc b/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.cc
index 8f6b9f6763108..2396dbefd7fa8 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.cc
+++ b/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.cc
@@ -175,56 +175,40 @@ class SubgroupMatrixMatMulNBitsTailCopyProgram final : public Program<SubgroupMa
       {"output_offset", ProgramUniformVariableDataType::Uint32});
 };
 
-Status ApplySubgroupMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
-                                      const Tensor* zero_points, const Tensor* bias,
-                                      uint32_t M,
-                                      uint32_t N,
-                                      uint32_t K,
-                                      uint32_t nbits,
-                                      uint32_t zero_blocks_per_col,
-                                      const SubgroupMatrixConfig& config,
-                                      onnxruntime::webgpu::ComputeContext& context,
-                                      Tensor* y,
-                                      const uint32_t weight_index,
-                                      const Tensor* weight_index_indirect) {
-  const bool has_zero_points = zero_points != nullptr;
-  const bool has_bias = bias != nullptr;
-  const bool has_weight_idx_indirect = weight_index_indirect != nullptr;
-  const bool has_weight_idx = weight_index > 0 || has_weight_idx_indirect;
+namespace {
+
+struct MatMulNBitsTiling {
+  uint32_t m;
+  uint32_t n;
+  uint32_t k;
+  uint32_t workgroup_size;
+};
 
-  // Determine tile sizes first (needed for prepack padding).
-  uint32_t tile_m = 32;
-  uint32_t tile_n = 64;
-  uint32_t tile_size_k = 32;
-  uint32_t workgroup_size = 128;
+MatMulNBitsTiling SelectMatMulNBitsTiling(const SubgroupMatrixConfig& config,
+                                          uint32_t M, uint32_t N, uint32_t K, bool has_bias) {
   if (config.Is(8, 16, 16)) {
-    // Default: 64x64 tile size, 256 thread.
-    tile_m = 64;
-    tile_n = 64;
-    workgroup_size = 256;
-    // Bias output uses a scratch tile, so keep it at 64x64 for workgroup memory limits.
     if (!has_bias) {
       if (M >= 2048 && N % 256 == 0) {
-        // Large M (M>=2048): 128x256 tile size, 512 threads.
-        tile_m = 128;
-        tile_n = 256;
-        workgroup_size = 512;
-      } else if (M >= 256 && N % 128 == 0) {
-        // Mid-sized M (256<=M<2048): 128x128 tile size, 512 threads.
-        tile_m = 128;
-        tile_n = 128;
-        workgroup_size = 512;
+        return {128, 256, 32, 512};
+      }
+      if (M >= 256 && N % 128 == 0) {
+        return {128, 128, 32, 512};
       }
     }
-  } else if (config.Is(16, 16, 16)) {
+    return {64, 64, 32, 256};
+  }
+  if (config.Is(16, 16, 16)) {
     if (config.subgroupSize == 64) {
+      uint32_t tile_m;
+      uint32_t tile_n;
+      uint32_t tile_k;
       // A wave64 subgroup computes a 64x64 micro-tile. Small M distributes waves along N;
       // tall tiles share each dequantized B tile across row waves.
       if (M <= 64) {
         tile_m = 64;
         // Avoid a fully padded second column wave when N is exactly 64.
         tile_n = N == 64 ? 64 : 128;
-        tile_size_k = 64;
+        tile_k = 64;
       } else {
         // Use four waves when that does not add M padding over a two-wave tile and M is large
         // enough to amortize the larger workgroup.
@@ -235,53 +219,74 @@ Status ApplySubgroupMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Te
 
         // Larger K blocks reduce barriers. Four-wave tiles retain enough occupancy for BK=128;
         // larger two-wave grids use BK=64, while small two-wave grids favor BK=32.
-        tile_size_k = M < 256 ? 32 : (use_four_waves || M <= 384 ? 128 : 64);
+        tile_k = M < 256 ? 32 : (use_four_waves || M <= 384 ? 128 : 64);
       }
       // The K loop and cooperative-matrix loads do not mask a partial BK tile. K is guaranteed
       // to be divisible by 32, so step down until BK divides it exactly.
-      if (K % tile_size_k != 0) {
-        tile_size_k = tile_size_k == 128 && K % 64 == 0 ? 64 : 32;
+      if (K % tile_k != 0) {
+        tile_k = tile_k == 128 && K % 64 == 0 ? 64 : 32;
       }
-      workgroup_size = (tile_m / 64) * (tile_n / 64) * config.subgroupSize;
-    } else {
-      // Wave32: four subgroups in a 2x2 grid.
-      tile_m = 128;
-      tile_n = 128;
-      workgroup_size = 4 * config.subgroupSize;
+      return {tile_m, tile_n, tile_k, (tile_m / 64) * (tile_n / 64) * config.subgroupSize};
     }
+    // Wave32: four subgroups in a 2x2 grid.
+    return {128, 128, 32, 4 * config.subgroupSize};
+  }
+  return {32, 64, 32, 128};
+}
+
+}  // namespace
+
+Status PrepackMatMulNBitsInput(const Tensor* a, uint32_t M, uint32_t K, uint32_t padded_M,
+                               const SubgroupMatrixConfig& config, onnxruntime::webgpu::ComputeContext& context,
+                               Tensor& packed) {
+  PrepackProgram program{config.M, config.K};
+  program.SetWorkgroupSize(config.subgroupSize);
+  if (context.HasFeature(wgpu::FeatureName::SubgroupSizeControl)) {
+    program.SetSubgroupSize(config.subgroupSize);
   }
+  ORT_ENFORCE(K % config.K == 0 && padded_M % config.M == 0);
+  program.SetDispatchGroupSize(padded_M / config.M, K / config.K, 1);
+  packed = context.CreateGPUTensor(a->DataType(), TensorShape{padded_M, K});
+  program.AddInputs({{a, ProgramTensorMetadataDependency::TypeAndRank, 1}})
+      .AddOutputs({{&packed, ProgramTensorMetadataDependency::Rank, packed.Shape(), 1}})
+      .AddUniformVariables({{M}, {K}})
+      .CacheHint(config.M, config.K);
+  return context.RunProgram(program);
+}
+
+Status ApplySubgroupMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
+                                      const Tensor* zero_points, const Tensor* bias,
+                                      uint32_t M,
+                                      uint32_t N,
+                                      uint32_t K,
+                                      uint32_t nbits,
+                                      uint32_t zero_blocks_per_col,
+                                      const SubgroupMatrixConfig& config,
+                                      onnxruntime::webgpu::ComputeContext& context,
+                                      Tensor* y,
+                                      const uint32_t weight_index,
+                                      const Tensor* weight_index_indirect,
+                                      const Tensor* prepacked_a) {
+  const bool has_zero_points = zero_points != nullptr;
+  const bool has_bias = bias != nullptr;
+  const bool has_weight_idx_indirect = weight_index_indirect != nullptr;
+  const bool has_weight_idx = weight_index > 0 || has_weight_idx_indirect;
+
+  auto [tile_m, tile_n, tile_size_k, workgroup_size] = SelectMatMulNBitsTiling(config, M, N, K, has_bias);
   ORT_ENFORCE(workgroup_size <= context.DeviceLimits().maxComputeWorkgroupSizeX &&
                   workgroup_size <= context.DeviceLimits().maxComputeInvocationsPerWorkgroup,
               "Subgroup matrix MatMulNBits workgroup size exceeds device limits: ", workgroup_size);
 
   // If applicable, layout optimization of input matrix A(MxK) can be used for SubgroupMatrixLoad.
   Tensor a_prepack;
-  if (config.needsPrepack) {
-    const auto m = config.M;
-    const auto k = config.K;
-
-    // Optimize the layout of input matrix A(MxK) for SubgroupMatrixLoad.
-    PrepackProgram prepack_program{m, k};
-    prepack_program.SetWorkgroupSize(config.subgroupSize);
-    if (context.HasFeature(wgpu::FeatureName::SubgroupSizeControl)) {
-      prepack_program.SetSubgroupSize(config.subgroupSize);
-    }
-
-    // Pad M to workgroup tile size so all subgroups read valid prepacked data.
-    const uint32_t padded_M = CeilDiv(M, tile_m) * tile_m;
-    const auto dispatch_group_size_x = padded_M / m;
-    ORT_ENFORCE(K % k == 0, "K must be a multiple of ", k);
-    const auto dispatch_group_size_y = K / k;
-    // Each workgroup will process one subgroup matrix of size m x k.
-    prepack_program.SetDispatchGroupSize(dispatch_group_size_x, dispatch_group_size_y, 1);
-
-    TensorShape a_prepack_shape{padded_M, K};
-    a_prepack = context.CreateGPUTensor(a->DataType(), a_prepack_shape);
-    prepack_program.AddInputs({{a, ProgramTensorMetadataDependency::TypeAndRank, 1}})
-        .AddOutputs({{&a_prepack, ProgramTensorMetadataDependency::Rank, a_prepack.Shape(), 1}})
-        .AddUniformVariables({{M}, {K}})
-        .CacheHint(m, k);
-    ORT_RETURN_IF_ERROR(context.RunProgram(prepack_program));
+  if (prepacked_a != nullptr) {
+    ORT_ENFORCE(config.needsPrepack && prepacked_a->DataType() == a->DataType() &&
+                prepacked_a->Shape().NumDimensions() == 2 && prepacked_a->Shape()[1] == K &&
+                prepacked_a->Shape()[0] >= CeilDiv(M, tile_m) * tile_m);
+    a = prepacked_a;
+  } else if (config.needsPrepack) {
+    ORT_RETURN_IF_ERROR(PrepackMatMulNBitsInput(a, M, K, CeilDiv(M, tile_m) * tile_m,
+                                                config, context, a_prepack));
     a = &a_prepack;
   }
 
diff --git a/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.h b/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.h
index 240c476675e72..a0e037152e4a2 100644
--- a/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.h
+++ b/onnxruntime/contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.h
@@ -58,6 +58,12 @@ class SubgroupMatrixMatMulNBitsProgram final : public Program<SubgroupMatrixMatM
   bool has_tail_buffer_;
 };
 
+// Prepack one invocation's input. padded_M must cover every consumer's row tile.
+Status PrepackMatMulNBitsInput(const Tensor* a, uint32_t M, uint32_t K, uint32_t padded_M,
+                               const SubgroupMatrixConfig& config,
+                               onnxruntime::webgpu::ComputeContext& context, Tensor& packed);
+
+// prepacked_a must use the same input and subgroup layout, with sufficient row padding.
 Status ApplySubgroupMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Tensor* scales,
                                       const Tensor* zero_points, const Tensor* bias,
                                       uint32_t M,
@@ -69,7 +75,8 @@ Status ApplySubgroupMatrixMatMulNBits(const Tensor* a, const Tensor* b, const Te
                                       onnxruntime::webgpu::ComputeContext& context,
                                       Tensor* y,
                                       const uint32_t weight_index,
-                                      const Tensor* weight_index_indirect = nullptr);
+                                      const Tensor* weight_index_indirect = nullptr,
+                                      const Tensor* prepacked_a = nullptr);
 
 bool CanApplySubgroupMatrixMatMulNBits(onnxruntime::webgpu::ComputeContext& context,
                                        uint64_t accuracy_level,
diff --git a/onnxruntime/test/contrib_ops/matmul_2bits_test.cc b/onnxruntime/test/contrib_ops/matmul_2bits_test.cc
index 34544efa0ac2e..465327410bac9 100644
--- a/onnxruntime/test/contrib_ops/matmul_2bits_test.cc
+++ b/onnxruntime/test/contrib_ops/matmul_2bits_test.cc
@@ -61,6 +61,7 @@ struct TestOptions2Bits {
 
   // Run the node on the CUDA EP instead of the default (CPU / WebGPU) set.
   bool use_cuda{false};
+  bool disable_cpu_ep_fallback{false};
 
   std::optional<float> output_abs_error{};
   std::optional<float> output_rel_error{};
@@ -79,7 +80,8 @@ struct TestOptions2Bits {
 }
 
 template <typename T1>
-void RunTest2Bits(const TestOptions2Bits& opts) {
+void RunTest2Bits(const TestOptions2Bits& opts,
+                  std::vector<std::unique_ptr<IExecutionProvider>>&& explicit_eps = {}) {
   SCOPED_TRACE(opts);
   ASSERT_FALSE(opts.has_typed_zero_point && !opts.has_zero_point);
 
@@ -230,6 +232,19 @@ void RunTest2Bits(const TestOptions2Bits& opts) {
     test.SetOutputRelErr("Y", *opts.output_rel_error);
   }
 
+  if (opts.disable_cpu_ep_fallback) {
+    SessionOptions session_options;
+    session_options.use_per_session_threads = false;
+    ASSERT_STATUS_OK(session_options.config_options.AddConfigEntry(kOrtSessionOptionsDisableCPUEPFallback, "1"));
+    test.Config(session_options);
+  }
+
+  if (!explicit_eps.empty()) {
+    test.ConfigEps(std::move(explicit_eps));
+    test.RunWithConfig();
+    return;
+  }
+
   std::vector<std::unique_ptr<IExecutionProvider>> execution_providers;
   if (opts.use_cuda) {
 #ifdef USE_CUDA
@@ -255,6 +270,18 @@ void RunTest2Bits(const TestOptions2Bits& opts) {
     test.ConfigEps(std::move(execution_providers));
     test.RunWithConfig();
   }
+#ifdef USE_WEBGPU
+  if constexpr (std::is_same<T1, MLFloat16>::value) {
+    // FP16 WebGPU cases must execute too; an unused OpTester traps in Debug
+    // and silently omits coverage in Release.
+    execution_providers.push_back(DefaultWebGpuExecutionProvider());
+    SessionOptions session_options;
+    ASSERT_STATUS_OK(session_options.config_options.AddConfigEntry(kOrtSessionOptionsDisableCPUEPFallback, "1"));
+    test.Config(session_options);
+    test.ConfigEps(std::move(execution_providers));
+    test.RunWithConfig();
+  }
+#endif
 }
 
 template <typename AType, int M, int N, int K, int block_size, int accuracy_level = 0>
@@ -1520,6 +1547,28 @@ TEST(MatMul2Bits, Float32_2b_BlkLen32_Accuracy0) {
   TestMatMul2BitsTyped<float, 100, 288, 256, 32, 0>();
 }
 
+#ifdef USE_WEBGPU
+TEST(MatMulNBits, WebGpu_DP4A_KParallel_Int2) {
+  for (bool has_zero_point : {false, true}) {
+    TestOptions2Bits opts{};
+    opts.M = 65;
+    opts.N = 32;
+    opts.K = 2048;
+    opts.accuracy_level = 4;
+    opts.disable_cpu_ep_fallback = true;
+    opts.has_zero_point = has_zero_point;
+    opts.has_bias = true;
+    opts.output_abs_error = 0.1f;
+    opts.output_rel_error = 0.02f;
+    auto ep = DefaultWebGpuExecutionProvider();
+    if (!ep) GTEST_SKIP() << "WebGPU unavailable.";
+    std::vector<std::unique_ptr<IExecutionProvider>> eps;
+    eps.push_back(std::move(ep));
+    RunTest2Bits<float>(opts, std::move(eps));
+  }
+}
+#endif
+
 #if defined(USE_WEBGPU) && !defined(ORT_USE_EP_API_ADAPTERS)
 
 namespace {
diff --git a/onnxruntime/test/contrib_ops/matmul_4bits_test.cc b/onnxruntime/test/contrib_ops/matmul_4bits_test.cc
index 4d7434ef398dd..f9a2df44cc836 100644
--- a/onnxruntime/test/contrib_ops/matmul_4bits_test.cc
+++ b/onnxruntime/test/contrib_ops/matmul_4bits_test.cc
@@ -134,7 +134,7 @@ void RunTest(const TestOptions& opts,
                 "unexpected type for T1");
 
 #ifdef USE_CUDA
-  if (opts.accuracy_level != 0 && !opts.legacy_shape) {
+  if (explicit_eps.empty() && opts.accuracy_level != 0 && !opts.legacy_shape) {
     return;  // CUDA EP does not handle accuracy level, so only test one level to avoid unnecessary tests.
   }
 #endif
@@ -452,6 +452,29 @@ TEST(MatMulNBits, Float32_4b_Accuracy0) {
   TestMatMulNBitsTyped<float, 100, 288, 1234, 16, 0>();
 }
 
+#ifdef USE_WEBGPU
+TEST(MatMulNBits, WebGpu_DP4A_KParallel_Int4) {
+  for (bool has_zero_point : {false, true}) {
+    TestOptions opts{};
+    opts.batch_count = 2;
+    opts.M = 65;
+    opts.N = 32;
+    opts.K = 2048;
+    opts.accuracy_level = 4;
+    opts.disable_cpu_ep_fallback = true;
+    opts.has_zero_point = has_zero_point;
+    opts.has_bias = true;
+    opts.output_abs_error = 0.1f;
+    opts.output_rel_error = 0.02f;
+    auto ep = DefaultWebGpuExecutionProvider();
+    if (!ep) GTEST_SKIP() << "WebGPU unavailable.";
+    std::vector<std::unique_ptr<IExecutionProvider>> eps;
+    eps.push_back(std::move(ep));
+    RunTest<float>(opts, std::move(eps));
+  }
+}
+#endif
+
 TEST(MatMulNBits, Float32_4b_Accuracy1) {
   TestMatMulNBitsTyped<float, 1, 1, 16, 16, 1>();
   TestMatMulNBitsTyped<float, 1, 288, 1024, 128, 1>();
@@ -1268,7 +1291,7 @@ TEST(MatMulNBits, Float16_LargeK_AccumulatorOverflow) {
 //   accuracy_level 4, M = 8 -> dp4a_matmul.wgsl.template (where the adapter supports it)
 // CanApplyDP4AMatrixMatMulNBits also needs subgroups, a non-Apple vendor and M >= 4
 // (kMinMForTileOptimization) for an fp16 output, so the M = 2 case below only reaches
-// dp4a_matmul_small_m on an fp32-output or Qualcomm adapter; elsewhere it lands on the generic
+// dp4a_matmul_k_parallel on an fp32-output or Qualcomm adapter; elsewhere it lands on the generic
 // kernel and is one more shape for it rather than dp4a coverage.
 // With the option on and an adapter that has subgroup matrices, the subgroup-matrix path declines
 // itself (its cooperative-matrix result type is f16 and cannot honour the request) and the dispatch
diff --git a/onnxruntime/test/contrib_ops/matmul_8bits_test.cc b/onnxruntime/test/contrib_ops/matmul_8bits_test.cc
index 65ea7f5bf2ecb..a0b5d18e65bc7 100644
--- a/onnxruntime/test/contrib_ops/matmul_8bits_test.cc
+++ b/onnxruntime/test/contrib_ops/matmul_8bits_test.cc
@@ -3,6 +3,7 @@
 
 #ifndef ORT_MINIMAL_BUILD
 #include <algorithm>
+#include <cmath>
 #include <optional>
 
 #include "gtest/gtest.h"
@@ -57,6 +58,7 @@ struct TestOptions8Bits {
   bool zp_is_typed{false};  // true: zero_points use same type as scales (T); false: uint8
   bool has_g_idx{false};
   bool has_bias{false};
+  bool disable_cpu_ep_fallback{false};
 
   // When set, RunTest8Bits validates cross-session sharing of the pre-packed weights instead of
   // doing a single run. The model is run in two CPU sessions that use the same container.
@@ -77,11 +79,12 @@ struct TestOptions8Bits {
 }
 
 template <typename T1>
-void RunTest8Bits(const TestOptions8Bits& opts) {
+void RunTest8Bits(const TestOptions8Bits& opts,
+                  std::vector<std::unique_ptr<IExecutionProvider>>&& explicit_eps = {}) {
   SCOPED_TRACE(opts);
 
 #ifdef USE_CUDA
-  if (opts.accuracy_level != 0) {
+  if (explicit_eps.empty() && opts.accuracy_level != 0) {
     return;  // CUDA EP does not handle accuracy level, so only test one level to avoid unnecessary tests.
   }
 #endif
@@ -240,6 +243,19 @@ void RunTest8Bits(const TestOptions8Bits& opts) {
     return;
   }
 
+  if (opts.disable_cpu_ep_fallback) {
+    SessionOptions session_options;
+    session_options.use_per_session_threads = false;
+    ASSERT_STATUS_OK(session_options.config_options.AddConfigEntry(kOrtSessionOptionsDisableCPUEPFallback, "1"));
+    test.Config(session_options);
+  }
+
+  if (!explicit_eps.empty()) {
+    test.ConfigEps(std::move(explicit_eps));
+    test.RunWithConfig();
+    return;
+  }
+
   std::vector<std::unique_ptr<IExecutionProvider>> execution_providers;
 #ifdef USE_CUDA
   execution_providers.emplace_back(DefaultCudaExecutionProvider());
@@ -1057,6 +1073,34 @@ TEST(MatMulNBits, BFloat16_Int8_SmallMFallback) {
 }
 #endif
 
+#if defined(USE_WEBGPU)
+TEST(MatMulNBits, WebGpu_DP4A_KParallel_Int8) {
+  // Small-M K tails, both sides of the K-parallel threshold, and both output widths.
+  constexpr int64_t shapes[][3] = {
+      {3, 16, 1920}, {4, 32, 1920}, {4, 32, 2048}, {65, 16, 4096}, {65, 32, 4096}};
+  for (const auto& [m, n, k] : shapes) {
+    for (bool has_zero_point : {false, true}) {
+      TestOptions8Bits opts{};
+      opts.M = m;
+      opts.N = n;
+      opts.K = k;
+      opts.accuracy_level = 4;
+      opts.disable_cpu_ep_fallback = true;
+      opts.has_zero_point = has_zero_point;
+      opts.has_bias = true;
+      // Activation-quantization error for these Gaussian inputs grows with sqrt(K).
+      opts.output_abs_error = 0.05f * std::sqrt(static_cast<float>(k) / 256.0f);
+      opts.output_rel_error = 0.02f;
+      auto ep = DefaultWebGpuExecutionProvider();
+      if (!ep) GTEST_SKIP() << "WebGPU unavailable.";
+      std::vector<std::unique_ptr<IExecutionProvider>> eps;
+      eps.push_back(std::move(ep));
+      RunTest8Bits<float>(opts, std::move(eps));
+    }
+  }
+}
+#endif
+
 #if !defined(USE_CUDA) && !defined(USE_WEBGPU)
 #ifndef ENABLE_TRAINING
 // Pre-packing (and therefore cross-session sharing of pre-packed weights) is disabled in a full
diff --git a/onnxruntime/test/providers/webgpu/contrib_ops/matmul_nbits_test.cc b/onnxruntime/test/providers/webgpu/contrib_ops/matmul_nbits_test.cc
new file mode 100644
index 0000000000000..7c805305f9a9e
--- /dev/null
+++ b/onnxruntime/test/providers/webgpu/contrib_ops/matmul_nbits_test.cc
@@ -0,0 +1,112 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+#if !defined(ORT_MINIMAL_BUILD) && !defined(DISABLE_CONTRIB_OPS)
+
+#include <cstdint>
+#include <string_view>
+#include <utility>
+
+#include "gtest/gtest.h"
+
+#include "core/common/inlined_containers.h"
+#include "core/providers/webgpu/webgpu_context.h"
+#include "core/providers/webgpu/webgpu_provider_options.h"
+#include "core/session/onnxruntime_session_options_config_keys.h"
+#include "test/providers/provider_test_utils.h"
+#include "test/util/include/asserts.h"
+#include "test/util/include/default_providers.h"
+
+namespace onnxruntime {
+namespace test {
+namespace {
+
+void RunFloat16NarrowDP4AAccumulationTest(int64_t bits) {
+  ConfigOptions config_options;
+  ASSERT_STATUS_OK(config_options.AddConfigEntry(webgpu::options::kEnableMatmulFp32Accumulation,
+                                                 webgpu::options::kEnableMatmulFp32Accumulation_OFF));
+  auto provider = WebGpuExecutionProviderWithOptions(config_options);
+  if (!provider) {
+    GTEST_SKIP() << "WebGPU execution provider is unavailable.";
+  }
+  const auto& context = webgpu::WebGpuContextFactory::GetContext(provider->GetDeviceId());
+  if (!context.DeviceHasFeature(wgpu::FeatureName::ShaderF16)) {
+    GTEST_SKIP() << "WebGPU adapter does not support shader-f16.";
+  }
+  // Match DP4A dispatch: Apple adapters deliberately use other MatMulNBits paths.
+  if (!context.DeviceHasFeature(wgpu::FeatureName::Subgroups) ||
+      context.AdapterInfo().vendor == std::string_view{"apple"}) {
+    GTEST_SKIP() << "WebGPU DP4A path is unavailable.";
+  }
+
+  // N=32 excludes subgroup-matrix dispatch. M=4 and K=4096 select the narrow
+  // K-parallel DP4A path, including its unconditional FP32 accumulation.
+  constexpr int64_t M = 4, N = 32, K = 4096, block_size = 32;
+  constexpr int64_t blocks = K / block_size;
+  const int64_t values_per_byte = 8 / bits;
+  const int64_t blob_size = block_size / values_per_byte;
+  const int32_t zero_point = 1 << (bits - 1);
+  InlinedVector<MLFloat16> a(M * K);
+  InlinedVector<uint8_t> packed_weights(N * blocks * blob_size, 0);
+  InlinedVector<MLFloat16> scales(N * blocks, MLFloat16(16.0f));
+  InlinedVector<MLFloat16> expected(M * N);
+
+  // +/-127 quantizes exactly with activation scale 1. Each 32-element block
+  // contributes +/-65024, which fits FP16, but the first two K iterations sum
+  // to +/-130048. Later blocks cancel that sum; one omitted negative weight
+  // leaves exactly +/-2032. FP16 partial sums overflow before cancellation.
+  for (int64_t m = 0; m < M; ++m) {
+    const float sign = m % 2 == 0 ? 1.0f : -1.0f;
+    for (int64_t k = 0; k < K; ++k) {
+      a[m * K + k] = MLFloat16(sign * 127.0f);
+    }
+    for (int64_t n = 0; n < N; ++n) {
+      expected[m * N + n] = MLFloat16(sign * 2032.0f);
+    }
+  }
+  for (int64_t n = 0; n < N; ++n) {
+    for (int64_t k = 0; k < K; ++k) {
+      const int32_t weight = k == K - 1 ? 0 : (k < K / 2 ? 1 : -1);
+      const auto quantized = static_cast<uint8_t>(zero_point + weight);
+      packed_weights[(n * K + k) / values_per_byte] |=
+          static_cast<uint8_t>(quantized << ((k % values_per_byte) * bits));
+    }
+  }
+
+  OpTester test("MatMulNBits", 1, kMSDomain);
+  test.AddAttribute<int64_t>("K", K);
+  test.AddAttribute<int64_t>("N", N);
+  test.AddAttribute<int64_t>("bits", bits);
+  test.AddAttribute<int64_t>("block_size", block_size);
+  test.AddAttribute<int64_t>("accuracy_level", 4);
+  test.AddInput<MLFloat16>("A", {M, K}, a.data(), a.size());
+  test.AddInput<uint8_t>("B", {N, blocks, blob_size}, packed_weights.data(), packed_weights.size(), true);
+  test.AddInput<MLFloat16>("scales", {N, blocks}, scales.data(), scales.size(), true);
+  test.AddOutput<MLFloat16>("Y", {M, N}, expected.data(), expected.size());
+  test.SetOutputAbsErr("Y", 0.0f);
+  test.SetOutputRelErr("Y", 0.0f);
+
+  SessionOptions options;
+  options.use_per_session_threads = false;
+  ASSERT_STATUS_OK(options.config_options.AddConfigEntry(kOrtSessionOptionsDisableCPUEPFallback, "1"));
+  test.Config(options).ConfigEp(std::move(provider)).RunWithConfig();
+}
+
+}  // namespace
+
+TEST(WebGpuMatMulNBitsTest, Float16_Int2_NarrowDP4AAccumulation) {
+  RunFloat16NarrowDP4AAccumulationTest(2);
+}
+
+TEST(WebGpuMatMulNBitsTest, Float16_Int4_NarrowDP4AAccumulation) {
+  RunFloat16NarrowDP4AAccumulationTest(4);
+}
+
+TEST(WebGpuMatMulNBitsTest, Float16_Int8_NarrowDP4AAccumulation) {
+  RunFloat16NarrowDP4AAccumulationTest(8);
+}
+
+}  // namespace test
+}  // namespace onnxruntime
+
+#endif  // !defined(ORT_MINIMAL_BUILD) && !defined(DISABLE_CONTRIB_OPS)
diff --git a/onnxruntime/test/providers/webgpu/matmul_nbits_input_preparation_test.cc b/onnxruntime/test/providers/webgpu/matmul_nbits_input_preparation_test.cc
new file mode 100644
index 0000000000000..d3264071cca14
--- /dev/null
+++ b/onnxruntime/test/providers/webgpu/matmul_nbits_input_preparation_test.cc
@@ -0,0 +1,170 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+#if defined(USE_WEBGPU) && !defined(DISABLE_CONTRIB_OPS) && !defined(ORT_USE_EP_API_ADAPTERS) && !defined(__wasm__) && !defined(USE_EXTERNAL_DAWN)
+
+#include <memory>
+#include <optional>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "gtest/gtest.h"
+#include "contrib_ops/webgpu/quantization/dp4a_matmul_nbits.h"
+#include "contrib_ops/webgpu/quantization/subgroup_matrix_matmul_nbits.h"
+#include "core/framework/ort_value_name_idx_map.h"
+#include "core/graph/model.h"
+#include "core/providers/webgpu/allocator.h"
+#include "core/providers/webgpu/compute_context.h"
+#include "core/providers/webgpu/webgpu_provider_factory_creator.h"
+#include "test/test_environment.h"
+#include "test/util/include/asserts.h"
+
+namespace onnxruntime {
+namespace test {
+namespace {
+
+using namespace webgpu;
+namespace quantization = contrib::webgpu;
+
+class PreparationTestKernel final : public OpKernel {
+ public:
+  explicit PreparationTestKernel(const OpKernelInfo& info) : OpKernel(info) {}
+  Status Compute(OpKernelContext*) const override { return Status::OK(); }
+};
+
+class PreparationTestContext final : public OpKernelContext {
+ public:
+  explicit PreparationTestContext(AllocatorPtr allocator)
+      : OpKernelContext(nullptr, DefaultLoggingManager().DefaultLogger(), nullptr), allocator_(std::move(allocator)) {}
+  Status GetTempSpaceAllocator(AllocatorPtr* output) const override {
+    *output = allocator_;
+    return Status::OK();
+  }
+
+ private:
+  AllocatorPtr allocator_;
+};
+
+class MatMulNBitsInputPreparationTest : public ::testing::Test {
+ protected:
+  void SetUp() override {
+    ep_ = WebGpuProviderFactoryCreator::Create(options_)->CreateProvider();
+    ep_->SetLogger(&DefaultLoggingManager().DefaultLogger());
+    allocator_ = std::make_shared<GpuBufferAllocator>(
+        0, [this]() -> const BufferManager& { return Ep().BufferManager(); },
+        [this]() -> CommandRecordingState& { return Ep().Recording(); }, false);
+  }
+
+  WebGpuExecutionProvider& Ep() { return static_cast<WebGpuExecutionProvider&>(*ep_); }
+
+  template <typename T>
+  Tensor Upload(const TensorShape& shape, std::vector<T> values) {
+    Tensor tensor(DataTypeImpl::GetType<T>(), shape, allocator_);
+    Ep().BufferManager().Upload(Ep().Recording(), values.data(),
+                                reinterpret_cast<WGPUBuffer>(tensor.MutableDataRaw()), tensor.SizeInBytes());
+    return tensor;
+  }
+
+  template <typename T>
+  std::vector<T> Download(Tensor& tensor) {
+    std::vector<T> values(tensor.Shape().Size());
+    Ep().BufferManager().Download(Ep().Recording(), reinterpret_cast<WGPUBuffer>(tensor.MutableDataRaw()),
+                                  values.data(), tensor.SizeInBytes());
+    return values;
+  }
+
+  template <typename T>
+  void Check(bool subgroup_matrix) {
+    auto& device = WebGpuContextFactory::GetContext(0);
+    if constexpr (std::is_same_v<T, MLFloat16>) {
+      if (!device.DeviceHasFeature(wgpu::FeatureName::ShaderF16)) {
+        GTEST_SKIP() << "FP16 is unavailable.";
+      }
+    }
+    Model model("input_preparation", false, DefaultLoggingManager().DefaultLogger());
+    auto& node = model.MainGraph().AddNode("unused", "Identity", "", {}, {});
+    auto kernel_def = KernelDefBuilder().SetName("Identity").Provider(kWebGpuExecutionProvider).SinceVersion(1).Build();
+    const std::unordered_map<int, OrtValue> initializers;
+    const OrtValueNameIdxMap values;
+    const DataTransferManager transfers;
+    const AllocatorMap allocators;
+    OpKernelInfo info(node, *kernel_def, Ep(), initializers, values, transfers, allocators, options_);
+    PreparationTestKernel kernel(info);
+    PreparationTestContext kernel_context(allocator_);
+    webgpu::ComputeContext context(device, Ep(), kernel, kernel_context);
+    constexpr uint32_t M = 65, K = 256;
+    std::optional<SubgroupMatrixConfig> config;
+    if (subgroup_matrix) {
+      if (!quantization::CanApplySubgroupMatrixMatMulNBits(context, 4, 32, 1, 64, K, 4,
+                                                           std::is_same_v<T, MLFloat16>, config, M, false) ||
+          !config->needsPrepack) {
+        GTEST_SKIP() << "A compatible subgroup-matrix layout is unavailable.";
+      }
+    } else if (!quantization::CanApplyDP4AMatrixMatMulNBits(context, 4, 32, 64, K, 4)) {
+      GTEST_SKIP() << "DP4A is unavailable.";
+    }
+    // Each invocation owns a new workspace; two projections reuse it within that invocation.
+    for (int invocation = 0; invocation < 2; ++invocation) {
+      std::vector<T> input(M * K);
+      for (size_t i = 0; i < input.size(); ++i) {
+        input[i] = T(static_cast<float>(static_cast<int>(i % 17) - 8 + invocation) * 0.125f);
+      }
+      auto a = Upload<T>({M, K}, input);
+      quantization::DP4AQuantizedInput workspace;
+      Tensor packed;
+      if (subgroup_matrix) {
+        // M=65 exercises padding; 256 rows cover the supported workgroup row tiles.
+        ASSERT_STATUS_OK(quantization::PrepackMatMulNBitsInput(&a, M, K, 256, *config, context, packed));
+      }
+      const void* quantized_values = nullptr;
+      const void* quantized_scales = nullptr;
+      for (uint32_t projection = 0; projection < 2; ++projection) {
+        const uint32_t N = subgroup_matrix ? 64 : (projection == 0 ? 16 : 64);
+        const uint32_t bits = projection == 0 ? 8 : 4;
+        auto b = Upload<uint8_t>({N, K / 32, 4 * bits},
+                                 std::vector<uint8_t>(N * K * bits / 8, bits == 8 ? 131 : 0x99));
+        auto scales = Upload<T>({N, K / 32}, std::vector<T>(N * K / 32, T(0.125f)));
+        Tensor shared(DataTypeImpl::GetType<T>(), TensorShape{M, N}, allocator_);
+        Tensor fresh(DataTypeImpl::GetType<T>(), TensorShape{M, N}, allocator_);
+        if (subgroup_matrix) {
+          ASSERT_STATUS_OK(quantization::ApplySubgroupMatrixMatMulNBits(
+              &a, &b, &scales, nullptr, nullptr, M, N, K, bits, K / 32, *config, context, &shared,
+              0, nullptr, &packed));
+          ASSERT_STATUS_OK(quantization::ApplySubgroupMatrixMatMulNBits(
+              &a, &b, &scales, nullptr, nullptr, M, N, K, bits, K / 32, *config, context, &fresh, 0));
+        } else {
+          ASSERT_STATUS_OK(quantization::ApplyDP4AMatrixMatMulNBits(
+              &a, &b, &scales, nullptr, nullptr, 1, M, M, N, K, 32, K / 32, 4, bits,
+              context, &shared, 0, nullptr, &workspace));
+          ASSERT_TRUE(workspace.values.has_value());
+          ASSERT_TRUE(workspace.scales.has_value());
+          if (projection == 0) {
+            quantized_values = workspace.values->DataRaw();
+            quantized_scales = workspace.scales->DataRaw();
+          } else {
+            EXPECT_EQ(workspace.values->DataRaw(), quantized_values);
+            EXPECT_EQ(workspace.scales->DataRaw(), quantized_scales);
+          }
+          ASSERT_STATUS_OK(quantization::ApplyDP4AMatrixMatMulNBits(
+              &a, &b, &scales, nullptr, nullptr, 1, M, M, N, K, 32, K / 32, 4, bits, context, &fresh, 0));
+        }
+        EXPECT_EQ(Download<T>(shared), Download<T>(fresh));
+      }
+    }
+  }
+
+ private:
+  ConfigOptions options_;
+  std::unique_ptr<IExecutionProvider> ep_;
+  AllocatorPtr allocator_;
+};
+
+TEST_F(MatMulNBitsInputPreparationTest, SharedDP4AInputFloat32) { Check<float>(false); }
+TEST_F(MatMulNBitsInputPreparationTest, SharedDP4AInputFloat16) { Check<MLFloat16>(false); }
+TEST_F(MatMulNBitsInputPreparationTest, SharedSubgroupMatrixInputFloat16) { Check<MLFloat16>(true); }
+
+}  // namespace
+}  // namespace test
+}  // namespace onnxruntime
+#endif