Merge pull request #30297 from sanjayk9:dnn-matmul-nobias-gemm

dnn: skip MatMul's zero-fill when there is no bias - #30297

Without a bias, MatMul pre-zeroed the whole output with a single-threaded memset and then still passed beta = 1 into the gemm, forcing every kernel to read and multiply-add against that zeroed buffer. Passing beta = 0 instead removes the memset and lets MLAS/fastGemmThin write the product directly. Measured 31-50% faster on a small, write-heavy, no-bias shape (512x32x512, 8 threads); `MatMul.matmul/41` and the new `MatMul.matmul_dynamicB/41` cover the constant-B and dynamic-B paths.

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [x] I agree to contribute to the project under Apache 2 License.
- [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [x] The PR is proposed to the proper branch
- [ ] There is a reference to the original bug report and related work
- [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [ ] The feature is well documented and sample code can be built with the project CMake

modules/dnn/perf/perf_gemm.cpp

@@ -90,6 +90,9 @@ static const GemmParam_t test_matmul_configs[] = {
     { {1, 1, 3072}, {1, 3072, 768} },
     { {1, 2, 768}, {1, 768, 768} },
     { {1, 4, 768}, {1, 768, 768} },
+
+    // write-heavy, small K, no bias: the output store dominates
+    { {512, 32}, {32, 512} },
 };
 
 struct GemmParamId
@@ -356,6 +359,56 @@ PERF_TEST_P_(MatMul, matmul)
     SANITY_CHECK_NOTHING();
 }
 
+// Same configs as matmul, but with B as a second network input instead of a constant
+// blob - exercises the blobs.empty() dispatch (MLAS / fastGemmBatch "truly dynamic B"),
+// which matmul above never reaches since it always packs B as a blob.
+PERF_TEST_P_(MatMul, matmul_dynamicB)
+{
+    int test_id = (int)get<0>(GetParam());
+    ASSERT_GE(test_id, 0); ASSERT_LT(test_id, MatMulParamId::MATMUL_LAST);
+    const GemmParam_t& params = test_matmul_configs[test_id];
+    auto a_shape = params.a_shape;
+    auto b_shape = params.b_shape;
+    auto trans_a = params.trans_a;
+    auto trans_b = params.trans_b;
+
+    Backend backend_id = get<0>(get<1>(GetParam()));
+    Target target_id = get<1>(get<1>(GetParam()));
+
+    Mat A(a_shape, CV_32F);
+    randu(A, -1.0f, 1.0f);
+    Mat B(b_shape, CV_32F);
+    randu(B, -1.0f, 1.0f);
+
+    LayerParams lp;
+    lp.type = "MatMul";
+    lp.name = "testLayer";
+    lp.set("transA", trans_a);
+    lp.set("transB", trans_b);
+
+    Net net;
+    int id = net.addLayerToPrev(lp.name, lp.type, lp);
+    net.connect(0, 1, id, 1);
+    net.setPreferableBackend(backend_id);
+    net.setPreferableTarget(target_id);
+
+    // warmup
+    {
+        std::vector<std::string> input_names{"A", "B"};
+        net.setInputsNames(input_names);
+        net.setInput(A, input_names[0]);
+        net.setInput(B, input_names[1]);
+        Mat out = net.forward();
+    }
+
+    TEST_CYCLE()
+    {
+        Mat res = net.forward();
+    }
+
+    SANITY_CHECK_NOTHING();
+}
+
 PERF_TEST_P_(MatMul, innerproduct)
 {
     int test_id = (int)get<0>(GetParam());

modules/dnn/src/cuda4dnn/primitives/matmul_broadcast.hpp

@@ -25,7 +25,7 @@ namespace cv { namespace dnn { namespace cuda4dnn {
     public:
         using wrapper_type = GetCUDABackendWrapperType<T>;
 
-        MatMulBroadcastOp(csl::Stream stream_, csl::cublas::Handle handle, const Mat &B, const Mat &bias, bool _transA, bool _transB,
+        MatMulBroadcastOp(csl::Stream stream_, csl::cublas::Handle handle, const Mat &B, const Mat &bias, float beta_, bool _transA, bool _transB,
                  const std::vector<size_t> &A_offsets_, const std::vector<size_t> &B_offsets_, std::vector<size_t> &C_offsets_,
                  size_t batch_)
             : stream(std::move(stream_)), cublasHandle(std::move(handle)), A_offsets(A_offsets_), B_offsets(B_offsets_), C_offsets(C_offsets_), batch(batch_)
@@ -42,6 +42,7 @@ namespace cv { namespace dnn { namespace cuda4dnn {
 
             transA = _transA;
             transB = _transB;
+            beta = beta_;
         }
 
         void forward(
@@ -68,14 +69,20 @@ namespace cv { namespace dnn { namespace cuda4dnn {
             // add bias if exists
             if (!bias_tensor.empty() || inputs.size() >= 3) {
                 csl::TensorView<T> bias;
+                float bias_scale = 1.f;  // a constant bias comes pre-scaled by beta
                 if (bias_tensor.empty()) {
                     auto bias_wrapper = inputs[2].dynamicCast<wrapper_type>();
                     bias = bias_wrapper->getView();
+                    bias_scale = beta;
                 } else {
                     bias = csl::TensorView<T>(bias_tensor);
                 }
 
-                kernels::eltwise_sum_2<T>(stream, output, output, bias);
+                if (bias_scale == 1.f)
+                    kernels::eltwise_sum_2<T>(stream, output, output, bias);
+                else
+                    kernels::eltwise_sum_coeff_2<T>(stream, output, T(1.f), output,
+                                                    T(bias_scale), bias);
             }
         }
 
@@ -85,6 +92,7 @@ namespace cv { namespace dnn { namespace cuda4dnn {
         csl::Tensor<T> input_B_tensor;
         csl::Tensor<T> bias_tensor;
         bool transA, transB;
+        float beta;
 
         std::vector<size_t> A_offsets;
         std::vector<size_t> B_offsets;

modules/dnn/src/layers/cpu_kernels/fast_gemm.hpp

@@ -131,6 +131,9 @@ struct MatMulHelper {
             B_steps[i] = B_steps[i + 1] * B_broadcast_shape[i + 1];
         }
         size_t t, idx;
+        // M*K or N*K can be 0 (empty A/B); row indices are 0 and unread for K == 0 -
+        // matmul_layer.cpp's forward_ocl() (the only reader) bails out before using them.
+        const size_t A_step = (size_t)M * K, B_step = (size_t)N * K;
         for (size_t i = 0; i < batch; i++) {
             C_offsets[i] = i * C_step;
             C_rows[i] = i;
@@ -146,8 +149,8 @@ struct MatMulHelper {
             }
             A_offsets[i] = A_offset;
             B_offsets[i] = B_offset;
-            A_rows[i] = A_offset / (M * K);
-            B_rows[i] = B_offset / (N * K);
+            A_rows[i] = A_step > 0 ? A_offset / A_step : 0;
+            B_rows[i] = B_step > 0 ? B_offset / B_step : 0;
         }
     }
 
@@ -220,4 +223,4 @@ void pagedAttnAVGemm(
 
 }} // cv::dnn
 
-#endif // OPENCV_DNN_FAST_GEMM_HPP
+#endif // OPENCV_DNN_FAST_GEMM_HPP

modules/dnn/src/layers/matmul_layer.cpp

@@ -234,7 +234,7 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
             B_mat = &promoted_B;
         }
 
-        if (B_mat && B_mat->data != last_packed_input_B_data) {
+        if (B_mat && helper.N > 0 && helper.K > 0 && B_mat->data != last_packed_input_B_data) {
             packed_input_B.clear();
             packed_input_B.shrink_to_fit();
             thin_packed_B.clear();
@@ -264,8 +264,10 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                 auto *broadcast_bias_ptr = broadcast_bias.ptr<float>();
 
                 const auto *bias = bias_mat.ptr<const float>();
+                // Stored unscaled: forward() applies beta once (gemm beta or explicitly), and the
+                // OpenCL/nGraph/CUDA paths add scaledBroadcastBias() instead.
                 if (bias_mat.total() == 1) { // [], [1], [1, ...]
-                    float b = (*bias) * beta;
+                    float b = *bias;
                     for (size_t i = 0; i < broadcast_bias.total(); i++) {
                         broadcast_bias_ptr[i] = b;
                     }
@@ -275,7 +277,7 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                     for (size_t i = 0; i < loops; i++) {
                         size_t step = i * inner_size;
                         for (size_t j = 0; j < inner_size; j++) {
-                            broadcast_bias_ptr[step + j] = beta * bias[j];
+                            broadcast_bias_ptr[step + j] = bias[j];
                         }
                     }
                 } else {
@@ -340,14 +342,19 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
 
         const auto *a = A.ptr<const float>();
         auto *y = Y.ptr<float>();
+        // Without a bias the gemm overwrites Y (beta = 0), so Y needs no zero-fill.
+        const bool has_bias = (inputs.size() + blobs.size()) >= 3;
+        const float gemm_beta = has_bias ? beta : 0.f;
         // add bias if existed
-        if ((inputs.size() + blobs.size()) >= 3) {
+        if (has_bias) {
             const auto &shape_Y = shape(Y);
+            // Stored unscaled here; beta is applied exactly once below, either by the gemm's
+            // own beta (M, N, K all > 0) or explicitly when the gemm is skipped.
             if (blobs.empty()) { // bias from input
                 const auto &bias_mat = inputs.back();
                 const auto *bias = bias_mat.ptr<const float>();
                 if (bias_mat.total() == 1) { // [], [1], [1, ...]
-                    float b = (*bias) * beta;
+                    float b = *bias;
                     for (size_t i = 0; i < Y.total(); i++) {
                         y[i] = b;
                     }
@@ -358,7 +365,7 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                         for (int i = r.start; i < r.end; i++) {
                             const size_t output_offset = i * inner_size;
                             for (size_t j = 0; j < inner_size; j++) {
-                                y[output_offset + j] = beta * bias[j];
+                                y[output_offset + j] = bias[j];
                             }
                         }
                     }, double(batches * inner_size * (1 / 1024.0)));
@@ -369,8 +376,18 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                 const auto *bias = broadcast_bias.ptr<const float>();
                 std::memcpy(y, bias, total(shape_Y) * sizeof(float));
             }
-        } else {
-            std::memset(y, 0, Y.total() * sizeof(float));
+        }
+        // Any of M, N, K == 0 means the product contributes nothing, and finalize()
+        // skips packing B whenever N or K is 0 - never reach the gemm dispatch (and its own
+        // beta) below, so apply beta to the bias prefill here instead.
+        if (helper.M == 0 || helper.N == 0 || helper.K == 0) {
+            if (!has_bias)
+                std::memset(y, 0, Y.total() * sizeof(float));
+            else if (beta != 1.f)
+                for (size_t i = 0; i < Y.total(); i++)
+                    y[i] *= beta;
+            fusion.run(Y);
+            return;
         }
 
         if (blobs.empty()) {
@@ -385,29 +402,29 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                 done = mlasSgemmBatch(helper.batch,
                                       helper.A_offsets.data(), helper.B_offsets.data(), helper.C_offsets.data(),
                                       trans_a, trans_b, helper.M, helper.N, helper.K,
-                                      alpha, a, lda_mem, b, ldb_mem, beta, y, helper.ldc);
+                                      alpha, a, lda_mem, b, ldb_mem, gemm_beta, y, helper.ldc);
             }
             if (!done) {
                 fastGemmBatch(helper.batch, helper.A_offsets.data(), helper.B_offsets.data(), helper.C_offsets.data(),
                               helper.M, helper.N, helper.K, alpha, a, helper.lda0, helper.lda1,
-                              b, helper.ldb0, helper.ldb1, beta, y, helper.ldc, opt);
+                              b, helper.ldb0, helper.ldb1, gemm_beta, y, helper.ldc, opt);
             }
         } else if (!thin_packed_B.empty()) {
             fastGemmThin(helper.M, helper.N, helper.K, alpha,
                          a, helper.lda0, helper.lda1,
-                         thin_packed_B.data(), beta,
+                         thin_packed_B.data(), gemm_beta,
                          y, helper.ldc, opt.multi_thread);
         } else if (!packed_input_B.empty()) {
             fastGemmBatch(helper.batch, helper.A_offsets.data(), helper.packed_B_offsets.data(), helper.C_offsets.data(),
                           helper.M, helper.N, helper.K, alpha, a, helper.lda0, helper.lda1,
-                          packed_input_B.data(), beta, y, helper.ldc, opt);
+                          packed_input_B.data(), gemm_beta, y, helper.ldc, opt);
         } else {
             // truly dynamic B (changes every call — no packing cache available)
             const auto &B = inputs[1];
             const auto *b = B.ptr<const float>();
             fastGemmBatch(helper.batch, helper.A_offsets.data(), helper.B_offsets.data(), helper.C_offsets.data(),
                           helper.M, helper.N, helper.K, alpha, a, helper.lda0, helper.lda1,
-                          b, helper.ldb0, helper.ldb1, beta, y, helper.ldc, opt);
+                          b, helper.ldb0, helper.ldb1, gemm_beta, y, helper.ldc, opt);
         }
         fusion.run(Y);
     }
@@ -540,6 +557,10 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
         if (inputs.size() >= 3) {
             return false;
         }
+        // batch_A/batch_B below divide by M*K and N*K; the CPU path handles empty dims
+        if (helper.M == 0 || helper.N == 0 || helper.K == 0) {
+            return false;
+        }
 
         const auto &input_A = inputs[0];
         auto &output = outputs[0];
@@ -556,7 +577,7 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
             }
             if ((inputs.size() + blobs.size() >= 3)) {
                 if (bias_umat.empty()) {
-                    broadcast_bias.copyTo(bias_umat);
+                    scaledBroadcastBias().copyTo(bias_umat);
                 }
             } else {
                 if (bias_umat.empty()) {
@@ -629,6 +650,10 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
             result = std::make_shared<ov::op::v0::MatMul>(input_A_node, input_B_node, trans_a, trans_b);
             if (nodes.size() >= 3) {
                 bias = nodes[2].dynamicCast<InfEngineNgraphNode>()->node;
+                if (beta != 1.f) {
+                    auto coeff = std::make_shared<ov::op::v0::Constant>(ov::element::f32, ov::Shape{}, &beta);
+                    bias = std::make_shared<ov::op::v1::Multiply>(bias, coeff, ov::op::AutoBroadcastType::NUMPY);
+                }
                 result = std::make_shared<ov::op::v1::Add>(result, bias);
             }
         } else {
@@ -636,8 +661,9 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
             auto input_B_node = std::make_shared<ov::op::v0::Constant>(ov::element::f32, input_B_shape, blobs[0].data);
             result = std::make_shared<ov::op::v0::MatMul>(input_A_node, input_B_node, trans_a, trans_b);
             if ((nodes.size() + blobs.size()) >= 3) {
-                const auto bias_shape = shape(broadcast_bias);
-                bias = std::make_shared<ov::op::v0::Constant>(ov::element::f32, std::vector<size_t>(bias_shape.begin(), bias_shape.end()), broadcast_bias.data);
+                const Mat bias_mat = scaledBroadcastBias();
+                const auto bias_shape = shape(bias_mat);
+                bias = std::make_shared<ov::op::v0::Constant>(ov::element::f32, std::vector<size_t>(bias_shape.begin(), bias_shape.end()), bias_mat.data);
                 result = std::make_shared<ov::op::v1::Add>(result, bias);
             }
         }
@@ -678,13 +704,13 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
         if (!blobs.empty()) {
             input_B = blobs.front();
             if (blobs.size() >= 2) {
-                bias = broadcast_bias;
+                bias = scaledBroadcastBias();
             }
         }
 
         CV_CheckFalse(helper.empty(), "DNN/MatMul/CUDA: MatMulHelper is not initialized");
 
-        return make_cuda_node<cuda4dnn::MatMulBroadcastOp>(preferableTarget, std::move(context->stream), std::move(context->cublas_handle), input_B, bias, trans_a, trans_b, helper.A_offsets, helper.B_offsets, helper.C_offsets, helper.batch);
+        return make_cuda_node<cuda4dnn::MatMulBroadcastOp>(preferableTarget, std::move(context->stream), std::move(context->cublas_handle), input_B, bias, beta, trans_a, trans_b, helper.A_offsets, helper.B_offsets, helper.C_offsets, helper.batch);
     }
 #endif // HAVE_CUDA
 
@@ -717,6 +743,20 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
                 auto input_bias_wrapper = inputs[2].dynamicCast<CannBackendWrapper>();
                 auto input_bias_desc = input_bias_wrapper->getTensorDesc();
                 auto input_bias_node = nodes[2].dynamicCast<CannBackendNode>()->getOp();
+                if (beta != 1.f) { // BatchMatMulV2 adds the bias unscaled
+                    Mat beta_mat(1, 1, CV_32F, Scalar(beta));
+                    auto const_beta_node = std::make_shared<CannConstOp>(
+                        beta_mat.data, beta_mat.type(), std::vector<int>{1},
+                        cv::format("%s_beta", name.c_str()));
+                    auto mul = std::make_shared<ge::op::Mul>(
+                        cv::format("%s_bias_beta", name.c_str()));
+                    mul->set_input_x1_by_name(*input_bias_node, "y");
+                    mul->update_input_desc_x1(*input_bias_desc);
+                    mul->set_input_x2_by_name(*(const_beta_node->getOp()), "y");
+                    mul->update_input_desc_x2(*(const_beta_node->getTensorDesc()));
+                    mul->update_output_desc_y(*input_bias_desc);
+                    input_bias_node = mul;
+                }
                 op->set_input_bias_by_name(*input_bias_node, "y");
                 op->update_input_desc_bias(*input_bias_desc);
             }
@@ -726,7 +766,8 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
             op->set_input_x2_by_name(*(const_B_node->getOp()), "y");
             op->update_input_desc_x2(*(const_B_node->getTensorDesc()));
             if ((inputs.size() + blobs.size()) >= 3) { // does not support broadcast bias
-                auto bias_mat = blobs.back();
+                Mat bias_mat;
+                blobs.back().convertTo(bias_mat, -1, beta);  // BatchMatMulV2 adds the bias unscaled
                 auto bias_shape = shape(bias_mat);
 
                 // reshape if 1d
@@ -748,6 +789,14 @@ class MatMulLayerImpl CV_FINAL : public MatMulLayer {
 #endif // HAVE_CANN
 
  private:
+    Mat scaledBroadcastBias() const {
+        if (beta == 1.f)
+            return broadcast_bias;
+        Mat scaled;
+        broadcast_bias.convertTo(scaled, -1, beta);
+        return scaled;
+    }
+
     int real_ndims_C;
 
     std::vector<float> packed_input_B;

modules/dnn/test/test_layers.cpp

@@ -4017,4 +4017,168 @@ TEST(Test_MatMul, PackedBroadcastMultipleSlices)
     normAssert(outputs[0], expected, "MatMul packed broadcast multi-slice mismatch", 1e-4, 1e-4);
 }
 
+// M, N or K == 0 used to divide by zero or read out of bounds.
+TEST(Test_MatMul, EmptyDimsDoNotCrash)
+{
+    const int dims[][3] = {  // {M, K, N}
+        {0, 4, 5},  // empty A (no rows)
+        {3, 4, 0},  // empty B (no columns)
+        {3, 0, 4},  // empty contraction dim
+    };
+    for (const auto &d : dims)
+    {
+        const int M = d[0], K = d[1], N = d[2];
+        for (bool constB : {false, true})
+        {
+            Mat A(M, K, CV_32F);
+            Mat B(K, N, CV_32F);
+
+            LayerParams lp;
+            lp.type = "MatMul";
+            lp.name = "matmul_empty_dims";
+            if (constB)
+                lp.blobs.push_back(B);
+
+            Ptr<Layer> layer = LayerFactory::createLayerInstance(lp.type, lp);
+            ASSERT_TRUE(layer);
+            std::vector<Mat> inputs = constB ? std::vector<Mat>{A} : std::vector<Mat>{A, B};
+            std::vector<Mat> outputs;
+            runLayer(layer, inputs, outputs);
+            ASSERT_EQ(outputs.size(), (size_t)1);
+            ASSERT_EQ(shape(outputs[0]), MatShape({M, N}));
+
+            Mat expected = Mat::zeros(M, N, CV_32F);
+            normAssert(outputs[0], expected, "MatMul empty-dims output mismatch", 1e-4, 1e-4);
+        }
+    }
+}
+
+// K == 0 skips the gemm entirely, so beta must be applied to the bias explicitly instead of
+// relying on the (never-called) gemm's own beta. Rank-2 bias, both as a constant blob
+// (finalize()'s broadcast_bias path) and as a dynamic input (forward()'s own broadcast path).
+TEST(Test_MatMul, EmptyKDimAppliesBetaToBias)
+{
+    const int M = 3, N = 4;
+    Mat bias(M, N, CV_32F);
+    randu(bias, -1.f, 1.f);
+    const float beta = 2.0f;
+
+    for (bool constB : {false, true})
+    {
+        Mat A(M, 0, CV_32F), B(0, N, CV_32F);
+
+        LayerParams lp;
+        lp.type = "MatMul";
+        lp.name = "matmul_empty_k_bias";
+        lp.set("beta", beta);
+        lp.set("real_ndims_C", 2);
+        if (constB) {
+            lp.blobs.push_back(B);
+            lp.blobs.push_back(bias);
+        }
+
+        Ptr<Layer> layer = LayerFactory::createLayerInstance(lp.type, lp);
+        ASSERT_TRUE(layer);
+        std::vector<Mat> inputs = constB ? std::vector<Mat>{A} : std::vector<Mat>{A, B, bias};
+        std::vector<Mat> outputs;
+        runLayer(layer, inputs, outputs);
+
+        Mat expected = beta * bias;
+        normAssert(outputs[0], expected, "MatMul empty-K bias/beta mismatch", 1e-4, 1e-4);
+    }
+}
+
+// Backends add the bias themselves, so each must apply beta for every bias rank, both when
+// B and the bias are constant blobs and when they are network inputs.
+typedef DNNTestLayer Test_MatMul_BiasBeta;
+TEST_P(Test_MatMul_BiasBeta, Accuracy)
+{
+    const int M = 4, K = 8, N = 6;
+    const float beta = 2.f;
+    Mat A(M, K, CV_32F), B(K, N, CV_32F);
+    randu(A, -1.f, 1.f);
+    randu(B, -1.f, 1.f);
+
+    Mat biasScalar(1, 1, CV_32F, Scalar(0.75f));
+    Mat biasRow(1, N, CV_32F), biasFull(M, N, CV_32F);
+    randu(biasRow, -1.f, 1.f);
+    randu(biasFull, -1.f, 1.f);
+
+    struct BiasCase { Mat bias, broadcasted; int ndimsC; };
+    const BiasCase cases[] = {
+        { biasScalar, Mat(M, N, CV_32F, Scalar(0.75f)), 2 },
+        { biasRow.reshape(1, std::vector<int>{N}), repeat(biasRow, M, 1), 1 },
+        { biasFull, biasFull, 2 },
+    };
+    for (bool constInputs : {true, false})
+    for (const BiasCase& c : cases)
+    {
+        LayerParams lp;
+        lp.type = "MatMul";
+        lp.name = "matmul_bias_beta";
+        lp.set("beta", beta);
+        lp.set("real_ndims_C", c.ndimsC);
+
+        Net net;
+        if (constInputs)
+        {
+            lp.blobs.push_back(B);
+            lp.blobs.push_back(c.bias);
+            net.addLayerToPrev(lp.name, lp.type, lp);
+            net.setInput(A);
+        }
+        else
+        {
+            int id = net.addLayerToPrev(lp.name, lp.type, lp);
+            net.connect(0, 1, id, 1);
+            net.connect(0, 2, id, 2);
+            net.setInputsNames({"A", "B", "bias"});
+            net.setInput(A, "A");
+            net.setInput(B, "B");
+            net.setInput(c.bias, "bias");
+        }
+        net.setPreferableBackend(backend);
+        net.setPreferableTarget(target);
+        Mat out = net.forward();
+
+        Mat ref = A * B + beta * c.broadcasted;
+        normAssert(out.reshape(1, M), ref, constInputs ? "constant B and bias" : "B and bias as inputs",
+                   default_l1, default_lInf);
+    }
+}
+INSTANTIATE_TEST_CASE_P(/**/, Test_MatMul_BiasBeta, dnnBackendsAndTargets());
+
+// Without a bias the gemm runs with beta = 0 and must never read Y: poison it with NaN.
+// M = 2 with constant B takes fastGemmThin, M = 37 the packed fastGemmBatch path; dynamic
+// B (constB = false) takes the MLAS/fastGemmBatch "truly dynamic B" path in both cases.
+TEST(Test_MatMul, NoBiasIgnoresOutputContents)
+{
+    const int K = 53, N = 67;
+    for (int M : {2, 37})
+    for (bool constB : {false, true})
+    {
+        Mat A(M, K, CV_32F), B(K, N, CV_32F);
+        randu(A, -1.f, 1.f);
+        randu(B, -1.f, 1.f);
+
+        LayerParams lp;
+        lp.type = "MatMul";
+        lp.name = "matmul_no_bias";
+        if (constB)
+            lp.blobs.push_back(B);
+        Ptr<Layer> layer = LayerFactory::createLayerInstance(lp.type, lp);
+        ASSERT_TRUE(layer);
+
+        std::vector<Mat> inputs = constB ? std::vector<Mat>{A} : std::vector<Mat>{A, B};
+        std::vector<Mat> outputs(1, Mat(M, N, CV_32F, Scalar(std::numeric_limits<float>::quiet_NaN())));
+        std::vector<Mat> internals;
+        layer->finalize(inputs, outputs);
+        layer->forward(inputs, outputs, internals);
+
+        Mat expected;
+        gemm(A, B, 1., noArray(), 0., expected);
+        normAssert(outputs[0], expected, "MatMul no-bias output mismatch", 1e-4, 1e-4);
+    }
+}
+
 }} // namespace