diff --git a/include/onnxruntime/core/session/onnxruntime_c_api.h b/include/onnxruntime/core/session/onnxruntime_c_api.h
index 3c6df85c04..1a6d50bde9 100644
--- a/include/onnxruntime/core/session/onnxruntime_c_api.h
+++ b/include/onnxruntime/core/session/onnxruntime_c_api.h
@@ -61,6 +61,8 @@ extern "C" {
 #define _Check_return_
 #define _Outptr_result_maybenull_
 #define _In_reads_(X)
+#define _Inout_updates_(X)
+#define _Out_writes_(X)
 #define _Inout_updates_all_(X)
 #define _Out_writes_bytes_all_(X)
 #define _Out_writes_all_(X)
diff --git a/onnxruntime/contrib_ops/cpu/quantization/quant_gemm.cc b/onnxruntime/contrib_ops/cpu/quantization/quant_gemm.cc
index fcf774ee19..ff8ad09082 100644
--- a/onnxruntime/contrib_ops/cpu/quantization/quant_gemm.cc
+++ b/onnxruntime/contrib_ops/cpu/quantization/quant_gemm.cc
@@ -2,6 +2,7 @@
 // Licensed under the MIT License.
 
 #include "core/common/safeint.h"
+#include "core/common/narrow.h"
 #include "core/providers/cpu/math/gemm_base.h"
 #include "core/providers/cpu/math/gemm_helper.h"
 #include "core/providers/cpu/quantization/matmul_integer_base.h"
@@ -28,9 +29,9 @@ class QGemm : protected GemmBase, public MatMulIntegerBase {
     if (!helper.State().IsOK())
       return helper.State();
 
-    size_t M = SafeInt<size_t>(helper.M());
-    size_t N = SafeInt<size_t>(helper.N());
-    size_t K = SafeInt<size_t>(helper.K());
+    ptrdiff_t M = helper.M();
+    ptrdiff_t N = helper.N();
+    ptrdiff_t K = helper.K();
 
     // validate scales and zero points
     const auto* a_zp = context->Input<Tensor>(IN_A_ZERO_POINT);
@@ -39,7 +40,7 @@ class QGemm : protected GemmBase, public MatMulIntegerBase {
     const auto* a_scale = context->Input<Tensor>(IN_A_SCALE);
     const auto* b_scale = context->Input<Tensor>(IN_B_SCALE);
     const auto* y_scale = context->Input<Tensor>(IN_Y_SCALE);
-    CheckInputs(a_zp, b_zp, y_zp, a_scale, b_scale, y_scale, helper);
+    ORT_RETURN_IF_ERROR(CheckInputs(a_zp, b_zp, y_zp, a_scale, b_scale, y_scale, helper));
 
     AllocatorPtr allocator;
     ORT_RETURN_IF_ERROR(context->GetTempSpaceAllocator(&allocator));
@@ -66,7 +67,7 @@ class QGemm : protected GemmBase, public MatMulIntegerBase {
       }
     }
 
-    auto y = context->Output(OUT_Y, {SafeInt<int64_t>(M), SafeInt<int64_t>(N)});
+    auto y = context->Output(OUT_Y, {M, N});
     if (M == 0 || N == 0) return Status::OK();
 
     // prepare output buffer of GEMM
@@ -78,14 +79,16 @@ class QGemm : protected GemmBase, public MatMulIntegerBase {
       gemm_output_buffer.emplace(DataTypeImpl::GetType<int32_t>(), outputshape, allocator);
       gemm_output_data = gemm_output_buffer->MutableData<int32_t>();
     } else {
+      // y_scale is NULL. Then y_zp must also be NULL. In this case Y's type must be int32_t. This is checked by the
+      // OP schema type. So the following cast is safe.
       gemm_output_data = static_cast<int32_t*>(y->MutableDataRaw());
     }
 
     if (c != nullptr) {
-      GemmBroadcastBias(M, N, 1.f, c->Data<int32_t>(), &(c->Shape()), gemm_output_data);
+      GemmBroadcastBias<int32_t>(M, N, 1, c->Data<int32_t>(), &(c->Shape()), gemm_output_data);
     }
 
-    MLAS_GEMM_QUANT_SHAPE_PARAMS gemm_shape{M, N, K, a_is_signed, b_is_signed, c != nullptr};
+    MLAS_GEMM_QUANT_SHAPE_PARAMS gemm_shape{narrow<size_t>(M), narrow<size_t>(N), narrow<size_t>(K), a_is_signed, b_is_signed, c != nullptr};
     MLAS_GEMM_QUANT_DATA_PARAMS gemm_param;
 
     gemm_param.A = a_data;
@@ -137,29 +140,30 @@ class QGemm : protected GemmBase, public MatMulIntegerBase {
     OUT_Y = 0
   };
 
-  static void CheckInputs(const Tensor* a_zp, const Tensor* b_zp, const Tensor* y_zp,
-                          const Tensor* a_scale, const Tensor* b_scale, const Tensor* y_scale, const GemmHelper& helper) {
-    ORT_ENFORCE(IsScalarOr1ElementVector(a_scale),
-                "QGemm : scale of input a must be a scalar or 1D tensor of size 1");
-    ORT_ENFORCE(IsScalarOr1ElementVector(a_zp),
-                "QGemm : zero point of input a must be a scalar or 1D tensor of size 1");
+  static Status CheckInputs(const Tensor* a_zp, const Tensor* b_zp, const Tensor* y_zp,
+                            const Tensor* a_scale, const Tensor* b_scale, const Tensor* y_scale, const GemmHelper& helper) {
+    ORT_RETURN_IF_NOT(IsScalarOr1ElementVector(a_scale),
+                      "QGemm : scale of input a must be a scalar or 1D tensor of size 1");
+    ORT_RETURN_IF_NOT(IsScalarOr1ElementVector(a_zp),
+                      "QGemm : zero point of input a must be a scalar or 1D tensor of size 1");
 
     const auto& b_zp_shape = b_zp->Shape();
     const auto& b_scale_shape = b_scale->Shape();
-    ORT_ENFORCE(b_zp_shape.NumDimensions() == 0 ||
-                    (b_zp_shape.NumDimensions() == 1 && (b_zp_shape[0] == 1 || b_zp_shape[0] == helper.N())),
-                "QGemm : zero point of input b must be a scalar or 1D tensor of size 1 or N");
-    ORT_ENFORCE(b_scale_shape.NumDimensions() == 0 ||
-                    (b_scale_shape.NumDimensions() == 1 && (b_scale_shape[0] == 1 || b_scale_shape[0] == helper.N())),
-                "QGemm : scale of input b must be a scalar or 1D tensor of size 1 or N");
-    ORT_ENFORCE(b_scale_shape.NumDimensions() == b_zp_shape.NumDimensions() &&
-                    (b_scale_shape.NumDimensions() == 0 || (b_scale_shape[0] == b_zp_shape[0])),
-                "QGemm : zero point and scale of input b should have same shape size");
+    ORT_RETURN_IF_NOT(b_zp_shape.NumDimensions() == 0 ||
+                          (b_zp_shape.NumDimensions() == 1 && (b_zp_shape[0] == 1 || b_zp_shape[0] == helper.N())),
+                      "QGemm : zero point of input b must be a scalar or 1D tensor of size 1 or N");
+    ORT_RETURN_IF_NOT(b_scale_shape.NumDimensions() == 0 ||
+                          (b_scale_shape.NumDimensions() == 1 && (b_scale_shape[0] == 1 || b_scale_shape[0] == helper.N())),
+                      "QGemm : scale of input b must be a scalar or 1D tensor of size 1 or N");
+    ORT_RETURN_IF_NOT(b_scale_shape.NumDimensions() == b_zp_shape.NumDimensions() &&
+                          (b_scale_shape.NumDimensions() == 0 || (b_scale_shape[0] == b_zp_shape[0])),
+                      "QGemm : zero point and scale of input b should have same shape size");
 
-    ORT_ENFORCE(y_zp == nullptr || IsScalarOr1ElementVector(y_zp),
-                "QGemm : zero point of y must be null or a scalar or 1D tensor of size 1");
-    ORT_ENFORCE(y_scale == nullptr || IsScalarOr1ElementVector(y_scale),
-                "QGemm : scale of y must be null or a scalar or 1D tensor of size 1");
+    ORT_RETURN_IF_NOT(y_zp == nullptr || IsScalarOr1ElementVector(y_zp),
+                      "QGemm : zero point of y must be null or a scalar or 1D tensor of size 1");
+    ORT_RETURN_IF_NOT(y_scale == nullptr || IsScalarOr1ElementVector(y_scale),
+                      "QGemm : scale of y must be null or a scalar or 1D tensor of size 1");
+    return Status::OK();
   }
 
   std::vector<float> ComputeOutputScale(const Tensor* a_scale, const Tensor* b_scale, const Tensor* y_scale) const {
diff --git a/onnxruntime/core/providers/cpu/cpu_execution_provider.cc b/onnxruntime/core/providers/cpu/cpu_execution_provider.cc
index 7bf322f825..2b9d69d7f1 100644
--- a/onnxruntime/core/providers/cpu/cpu_execution_provider.cc
+++ b/onnxruntime/core/providers/cpu/cpu_execution_provider.cc
@@ -162,6 +162,9 @@ class ONNX_OPERATOR_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, Aco
 class ONNX_OPERATOR_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, Atan);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, 8, float, Gemm);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, 8, double, Gemm);
+#ifdef MLAS_F16VEC_INTRINSICS_SUPPORTED
+class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, 8, MLFloat16, Gemm);
+#endif
 class ONNX_OPERATOR_VERSIONED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 1, 10, Hardmax);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 1, 10, float, LogSoftmax);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 1, 10, double, LogSoftmax);
@@ -340,6 +343,9 @@ class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOn
 class ONNX_OPERATOR_VERSIONED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 10, Flatten);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 10, float, Gemm);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 10, double, Gemm);
+#ifdef MLAS_F16VEC_INTRINSICS_SUPPORTED
+class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 10, MLFloat16, Gemm);
+#endif
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 12, float, MatMul);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 12, double, MatMul);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 12, int32_t, MatMul);
@@ -481,6 +487,9 @@ class ONNX_OPERATOR_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, Sp
 class ONNX_OPERATOR_VERSIONED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, ScatterND);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, float, Gemm);
 class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, double, Gemm);
+#ifdef MLAS_F16VEC_INTRINSICS_SUPPORTED
+class ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, MLFloat16, Gemm);
+#endif
 class ONNX_OPERATOR_VERSIONED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, GatherElements);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, uint8_t, BitShift);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, uint32_t, BitShift);
@@ -577,6 +586,9 @@ class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain,
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, string, Expand);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, float, Gemm);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, double, Gemm);
+#ifdef MLAS_F16VEC_INTRINSICS_SUPPORTED
+class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, MLFloat16, Gemm);
+#endif
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, float, MatMul);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, double, MatMul);
 class ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, int32_t, MatMul);
@@ -2414,6 +2426,12 @@ Status RegisterFp16Kernels(KernelRegistry& kernel_registry) {
       BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 14, MLFloat16, Relu)>,
       BuildKernelCreateInfo<ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 6, 15, MLFloat16, LeakyRelu)>,
       BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 16, MLFloat16, LeakyRelu)>,
+      BuildKernelCreateInfo<ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 7, 8, MLFloat16, Gemm)>,
+      BuildKernelCreateInfo<ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 9, 10,
+                                                                            MLFloat16, Gemm)>,
+      BuildKernelCreateInfo<ONNX_OPERATOR_VERSIONED_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 11, 12, MLFloat16, Gemm)>,
+
+      BuildKernelCreateInfo<ONNX_OPERATOR_TYPED_KERNEL_CLASS_NAME(kCpuExecutionProvider, kOnnxDomain, 13, MLFloat16, Gemm)>,
   };
 
   for (auto& function_table_entry : function_table) {
diff --git a/onnxruntime/core/providers/cpu/math/gemm.cc b/onnxruntime/core/providers/cpu/math/gemm.cc
index 389f8bad4b..5a886cce9d 100644
--- a/onnxruntime/core/providers/cpu/math/gemm.cc
+++ b/onnxruntime/core/providers/cpu/math/gemm.cc
@@ -1,6 +1,7 @@
 // Copyright (c) Microsoft Corporation. All rights reserved.
 // Licensed under the MIT License.
 
+#include <onnxruntime_config.h>
 #include "core/providers/cpu/math/gemm.h"
 #include "core/common/narrow.h"
 #include "core/common/safeint.h"
@@ -25,6 +26,13 @@ ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
     double,
     KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<double>()),
     Gemm<double>);
+ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
+    Gemm,
+    7,
+    8,
+    MLFloat16,
+    KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<MLFloat16>()),
+    Gemm<MLFloat16>);
 
 // opset 9 added support for additional types (int32, uint32, int64, uint64), however we haven't enabled those yet.
 ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
@@ -41,6 +49,13 @@ ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
     double,
     KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<double>()),
     Gemm<double>);
+ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
+    Gemm,
+    9,
+    10,
+    MLFloat16,
+    KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<MLFloat16>()),
+    Gemm<MLFloat16>);
 
 // opset 11 made bias input 'C' optional
 ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
@@ -57,6 +72,13 @@ ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
     double,
     KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<double>()),
     Gemm<double>);
+ONNX_CPU_OPERATOR_VERSIONED_TYPED_KERNEL(
+    Gemm,
+    11,
+    12,
+    MLFloat16,
+    KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<MLFloat16>()),
+    Gemm<MLFloat16>);
 
 // opset 13 Adds BFloat16 support but we are not supporting it yet
 ONNX_CPU_OPERATOR_TYPED_KERNEL(
@@ -71,6 +93,12 @@ ONNX_CPU_OPERATOR_TYPED_KERNEL(
     double,
     KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<double>()),
     Gemm<double>);
+ONNX_CPU_OPERATOR_TYPED_KERNEL(
+    Gemm,
+    13,
+    MLFloat16,
+    KernelDefBuilder().TypeConstraint("T", DataTypeImpl::GetTensorType<MLFloat16>()),
+    Gemm<MLFloat16>);
 
 bool GemmPackBFp32(AllocatorPtr& alloc,
                    const Tensor& tensor_b,
@@ -112,10 +140,10 @@ bool GemmPackBFp32(AllocatorPtr& alloc,
 
 template <typename T>
 void Gemm<T>::ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
-                          int64_t M, int64_t N, int64_t K,
-                          float alpha,
+                          ptrdiff_t M, ptrdiff_t N, ptrdiff_t K,
+                          T alpha,
                           const T* a_data, const T* b_data,
-                          float beta,
+                          T beta,
                           const T* c_data, const TensorShape* c_shape,
                           T* y_data,
                           concurrency::ThreadPool* thread_pool) {
@@ -127,7 +155,7 @@ void Gemm<T>::ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
   GemmBroadcastBias(M, N, beta, c_data, c_shape, y_data);
 
   math::Gemm<T>(trans_a, trans_b,
-                narrow<ptrdiff_t>(M), narrow<ptrdiff_t>(N), narrow<ptrdiff_t>(K),
+                M, N, K,
                 alpha,
                 a_data,
                 b_data,
@@ -138,8 +166,69 @@ void Gemm<T>::ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
                 thread_pool);
 }
 
+template <>
+void Gemm<MLFloat16>::ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
+                                  ptrdiff_t M, ptrdiff_t N, ptrdiff_t K,
+                                  MLFloat16 alpha,
+                                  const MLFloat16* a_data, const MLFloat16* b_data,
+                                  MLFloat16 beta,
+                                  const MLFloat16* c_data, const TensorShape* c_shape,
+                                  MLFloat16* y_data,
+                                  concurrency::ThreadPool* thread_pool) {
+  // if input is empty tensor, return directly as nothing need to be calculated.
+  if (M == 0 || N == 0)
+    return;
+
+#if defined(__GNUC__) && defined(HAS_CLASS_MEMACCESS)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wclass-memaccess"
+#endif
+  // MLFloat16's constructor is explicit, so here we need to use memset
+  if (c_data == nullptr)
+    memset(&beta, 0, sizeof(MLFloat16));
+#if defined(__GNUC__) && defined(HAS_CLASS_MEMACCESS)
+#pragma GCC diagnostic pop
+#endif
+#ifdef MLAS_F16VEC_INTRINSICS_SUPPORTED
+  bool support_mlas = false;
+  if (c_shape == nullptr) {
+    support_mlas = true;
+  } else if (c_shape->NumDimensions() == 1 && (*c_shape)[0] == N) {
+    support_mlas = true;
+  } else if (c_shape->NumDimensions() == 2 && (((*c_shape)[0] == 1 && (*c_shape)[1] == N) || ((*c_shape)[0] == N && (*c_shape)[1] == 1))) {
+    support_mlas = true;
+  }
+  if (trans_a == CblasNoTrans && trans_b == CblasNoTrans && support_mlas && alpha.ToFloat() == 1.0 && beta.ToFloat() == 1.0) {
+    MLAS_HALF_GEMM_DATA_PARAMS data;
+    data.A = a_data;
+    data.lda = K;
+    data.B = b_data;
+    data.ldb = N;
+    data.C = y_data;
+    data.ldc = N;
+    if (c_shape != nullptr) {
+      data.Bias = c_data;
+    }
+    MlasHalfGemmBatch(M, N, K, 1, &data, thread_pool);
+    return;
+  }
+#endif
+  // Fallback to Eigen
+  // Broadcast the bias as needed if bias is given
+  GemmBroadcastBias(M, N, beta, c_data, c_shape, y_data);
+#if defined(__GNUC__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wstrict-aliasing"
+#endif
+  math::Gemm<Eigen::half>(trans_a, trans_b, M, N, K, *reinterpret_cast<Eigen::half*>(&alpha),
+                          reinterpret_cast<const Eigen::half*>(a_data), reinterpret_cast<const Eigen::half*>(b_data), *reinterpret_cast<Eigen::half*>(&beta), reinterpret_cast<Eigen::half*>(y_data), thread_pool);
+#if defined(__GNUC__)
+#pragma GCC diagnostic pop
+#endif
+}
+
 template void Gemm<float>::ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
-                                       int64_t M, int64_t N, int64_t K,
+                                       ptrdiff_t M, ptrdiff_t N, ptrdiff_t K,
                                        float alpha,
                                        const float* a_data, const float* b_data,
                                        float beta,
@@ -196,16 +285,15 @@ Status Gemm<float>::UseSharedPrePackedBuffers(std::vector<BufferUniquePtr>& prep
 }
 
 template <typename T>
-void Gemm<T>::ComputeActivation(T* y_data, size_t y_size, concurrency::ThreadPool* thread_pool) const {
+void Gemm<T>::ComputeActivation(_Inout_updates_(y_size) T* y_data, ptrdiff_t y_size, _Inout_opt_ concurrency::ThreadPool* thread_pool) const {
   if (activation_) {
     std::unique_ptr<functors::ElementWiseRangedTransform<T>> f(activation_->Copy());
     f->input = y_data;
     f->output = y_data;
-    std::ptrdiff_t total_len = static_cast<std::ptrdiff_t>(y_size);
     double cost = f->Cost();
     functors::ElementWiseRangedTransform<T>* c(f.get());
     concurrency::ThreadPool::TryParallelFor(
-        thread_pool, total_len,
+        thread_pool, y_size,
         {static_cast<float>(sizeof(T)), static_cast<float>(sizeof(T)), cost},
         [c](std::ptrdiff_t first, std::ptrdiff_t last) { (*c)(first, last); });
   }
@@ -226,9 +314,9 @@ Status Gemm<T>::Compute(OpKernelContext* context) const {
   if (!helper.State().IsOK())
     return helper.State();
 
-  int64_t M = helper.M();
-  int64_t N = helper.N();
-  int64_t K = helper.K();
+  ptrdiff_t M = helper.M();
+  ptrdiff_t N = helper.N();
+  ptrdiff_t K = helper.K();
 
   auto Y = context->Output(0, {M, N});
 
@@ -243,6 +331,48 @@ Status Gemm<T>::Compute(OpKernelContext* context) const {
   ComputeGemm(trans_A_, trans_B_, M, N, K, alpha_, A->Data<T>(), B->Data<T>(), beta_,
               c_data, c_shape, y_data, thread_pool);
 
+  ComputeActivation(y_data, SafeInt<ptrdiff_t>(M) * N, thread_pool);
+
+  return Status::OK();
+}
+
+template <>
+Status Gemm<MLFloat16>::Compute(OpKernelContext* context) const {
+  concurrency::ThreadPool* thread_pool = context->GetOperatorThreadPool();
+
+  const auto* A = context->Input<Tensor>(0);
+  const auto* B = packed_b_ ? nullptr : context->Input<Tensor>(1);
+  const auto* C = context->Input<Tensor>(2);
+
+  // Bias could be missing. Treat as scalar 0 if that is the case.
+  GemmHelper helper(A->Shape(), trans_A_ != CblasNoTrans, B ? B->Shape() : b_shape_, trans_B_ != CblasNoTrans,
+                    C != nullptr ? C->Shape() : TensorShape({}));
+
+  if (!helper.State().IsOK())
+    return helper.State();
+
+  ptrdiff_t M = helper.M();
+  ptrdiff_t N = helper.N();
+  ptrdiff_t K = helper.K();
+
+  auto Y = context->Output(0, {M, N});
+
+  // if input is empty tensor, return as nothing need to be calculated and we've set the shape for the output
+  if (M == 0 || N == 0)
+    return Status::OK();
+
+  MLFloat16* y_data = Y->MutableData<MLFloat16>();
+
+  const MLFloat16* c_data = C != nullptr ? C->Data<MLFloat16>() : nullptr;
+  const TensorShape* c_shape = C != nullptr ? &C->Shape() : nullptr;
+
+  if (B) {
+    ComputeGemm(trans_A_, trans_B_, M, N, K, static_cast<MLFloat16>(alpha_), A->Data<MLFloat16>(), B->Data<MLFloat16>(), static_cast<MLFloat16>(beta_),
+                c_data, c_shape, y_data, thread_pool);
+  } else {
+    ORT_NOT_IMPLEMENTED("Prepacking of B is supported by MLAS half gemm API, but not implemented by this kernel yet");
+  }
+
   ComputeActivation(y_data, SafeInt<size_t>(M) * N, thread_pool);
 
   return Status::OK();
@@ -263,9 +393,9 @@ Status Gemm<float>::Compute(OpKernelContext* context) const {
   if (!helper.State().IsOK())
     return helper.State();
 
-  int64_t M = helper.M();
-  int64_t N = helper.N();
-  int64_t K = helper.K();
+  ptrdiff_t M = helper.M();
+  ptrdiff_t N = helper.N();
+  ptrdiff_t K = helper.K();
 
   auto Y = context->Output(0, {M, N});
 
diff --git a/onnxruntime/core/providers/cpu/math/gemm.h b/onnxruntime/core/providers/cpu/math/gemm.h
index 2cfcc69a2b..9539497325 100644
--- a/onnxruntime/core/providers/cpu/math/gemm.h
+++ b/onnxruntime/core/providers/cpu/math/gemm.h
@@ -29,10 +29,10 @@ class Gemm : protected GemmBase, public OpKernel {
                                    /*out*/ bool& used_shared_buffers) override;
 
   static void ComputeGemm(CBLAS_TRANSPOSE trans_a, CBLAS_TRANSPOSE trans_b,
-                          int64_t M, int64_t N, int64_t K,
-                          float alpha,
+                          ptrdiff_t M, ptrdiff_t N, ptrdiff_t K,
+                          T alpha,
                           const T* a_data, const T* b_data,
-                          float beta,
+                          T beta,
                           const T* c_data, const TensorShape* c_shape,
                           T* y_data,
                           concurrency::ThreadPool* thread_pool);
@@ -44,7 +44,7 @@ class Gemm : protected GemmBase, public OpKernel {
   // For fused gemm + activation
   std::unique_ptr<functors::ElementWiseRangedTransform<T>> activation_;
 
-  void ComputeActivation(T* y_data, size_t y_size, concurrency::ThreadPool* thread_pool) const;
+  void ComputeActivation(_Inout_updates_(y_size) T* y_data, ptrdiff_t y_size, _Inout_opt_ concurrency::ThreadPool* thread_pool) const;
 };
 
 }  // namespace onnxruntime
diff --git a/onnxruntime/core/providers/cpu/math/gemm_helper.h b/onnxruntime/core/providers/cpu/math/gemm_helper.h
index ef3d983fba..f37b00ac2c 100644
--- a/onnxruntime/core/providers/cpu/math/gemm_helper.h
+++ b/onnxruntime/core/providers/cpu/math/gemm_helper.h
@@ -3,7 +3,10 @@
 
 #pragma once
 #include "core/common/common.h"
+#include "core/common/narrow.h"
 #include "core/util/math_cpuonly.h"
+#include "core/framework/tensor_shape.h"
+#include "core/session/onnxruntime_c_api.h"
 
 namespace onnxruntime {
 
@@ -14,20 +17,31 @@ class GemmHelper {
     ORT_ENFORCE(left.NumDimensions() == 2 || left.NumDimensions() == 1);
     ORT_ENFORCE(right.NumDimensions() == 2);
 
+    for (size_t i = 0; i != left.NumDimensions(); ++i) {
+      ORT_ENFORCE(left[i] >= 0);
+      ORT_ENFORCE(left[i] <= std::numeric_limits<ptrdiff_t>::max());
+    }
+
+    for (size_t i = 0; i != right.NumDimensions(); ++i) {
+      ORT_ENFORCE(right[i] >= 0);
+      ORT_ENFORCE(right[i] <= std::numeric_limits<ptrdiff_t>::max());
+    }
+
     if (trans_left) {
-      M_ = left.NumDimensions() == 2 ? left[1] : left[0];
-      K_ = left.NumDimensions() == 2 ? left[0] : 1;
+      M_ = left.NumDimensions() == 2 ? static_cast<ptrdiff_t>(left[1]) : static_cast<ptrdiff_t>(left[0]);
+      K_ = left.NumDimensions() == 2 ? static_cast<ptrdiff_t>(left[0]) : 1;
     } else {
-      M_ = left.NumDimensions() == 2 ? left[0] : 1;
-      K_ = left.NumDimensions() == 2 ? left[1] : left[0];
+      M_ = left.NumDimensions() == 2 ? static_cast<ptrdiff_t>(left[0]) : 1;
+      K_ = left.NumDimensions() == 2 ? static_cast<ptrdiff_t>(left[1])
+                                     : static_cast<ptrdiff_t>(left[0]);
     }
 
     int k_dim;
     if (trans_right) {
-      N_ = right[0];
+      N_ = static_cast<ptrdiff_t>(right[0]);
       k_dim = 1;
     } else {
-      N_ = right[1];
+      N_ = static_cast<ptrdiff_t>(right[1]);
       k_dim = 0;
     }
 
@@ -45,13 +59,13 @@ class GemmHelper {
     ORT_ENFORCE(M_ >= 0 && K_ > 0 && N_ >= 0);
   }
 
-  int64_t M() const { return M_; }
-  int64_t N() const { return N_; }
-  int64_t K() const { return K_; }
+  ptrdiff_t M() const { return M_; }
+  ptrdiff_t N() const { return N_; }
+  ptrdiff_t K() const { return K_; }
   Status State() const { return status_; }
 
  private:
-  bool IsValidBroadcast(const TensorShape& bias_shape, int64_t M, int64_t N) {
+  static bool IsValidBroadcast(const TensorShape& bias_shape, ptrdiff_t M, ptrdiff_t N) {
     // valid shapes are (,) , (1, N) , (M, 1) , (M, N)
     if (bias_shape.NumDimensions() > 2)
       return false;
@@ -66,32 +80,33 @@ class GemmHelper {
   }
 
  private:
-  int64_t M_;
-  int64_t K_;
-  int64_t N_;
+  GemmHelper() = default;
+  ptrdiff_t M_;
+  ptrdiff_t K_;
+  ptrdiff_t N_;
   Status status_;
 };
 
 template <typename T>
-void GemmBroadcastBias(int64_t M, int64_t N, float beta,
-                       const T* c_data, const TensorShape* c_shape,
-                       T* y_data) {
+void GemmBroadcastBias(ptrdiff_t M, ptrdiff_t N, T beta,
+                       _In_opt_ const T* c_data, _In_opt_ const TensorShape* c_shape,
+                       _Out_writes_(M* N) T* y_data) {
   // Broadcast the bias as needed if bias is given
   if (beta != 0 && c_data != nullptr) {
     ORT_ENFORCE(c_shape != nullptr, "c_shape is required if c_data is provided");
-    auto output_mat = EigenMatrixMapRowMajor<T>(y_data, onnxruntime::narrow<size_t>(M), onnxruntime::narrow<size_t>(N));
+    auto output_mat = EigenMatrixMapRowMajor<T>(y_data, M, N);
     if (c_shape->Size() == 1) {
       // C is (), (1,) or (1, 1), set the scalar
       output_mat.setConstant(*c_data);
     } else if (c_shape->NumDimensions() == 1 || (*c_shape)[0] == 1) {
       // C is (N,) or (1, N)
-      output_mat.rowwise() = ConstEigenVectorMap<T>(c_data, onnxruntime::narrow<size_t>(N)).transpose();
+      output_mat.rowwise() = ConstEigenVectorMap<T>(c_data, N).transpose();
     } else if ((*c_shape)[1] == 1) {
       // C is (M, 1)
-      output_mat.colwise() = ConstEigenVectorMap<T>(c_data, onnxruntime::narrow<size_t>(M));
+      output_mat.colwise() = ConstEigenVectorMap<T>(c_data, M);
     } else {
       // C is (M, N), no broadcast needed.
-      output_mat = ConstEigenMatrixMapRowMajor<T>(c_data, onnxruntime::narrow<size_t>(M), onnxruntime::narrow<size_t>(N));
+      output_mat = ConstEigenMatrixMapRowMajor<T>(c_data, M, N);
     }
   }
 }
diff --git a/onnxruntime/core/providers/cpu/ml/linearclassifier.cc b/onnxruntime/core/providers/cpu/ml/linearclassifier.cc
index 02332281ee..c964a98cf2 100644
--- a/onnxruntime/core/providers/cpu/ml/linearclassifier.cc
+++ b/onnxruntime/core/providers/cpu/ml/linearclassifier.cc
@@ -2,6 +2,7 @@
 // Licensed under the MIT License.
 
 #include "core/providers/cpu/ml/linearclassifier.h"
+#include "core/common/narrow.h"
 #include "core/providers/cpu/math/gemm.h"
 
 namespace onnxruntime {
@@ -45,7 +46,7 @@ LinearClassifier::LinearClassifier(const OpKernelInfo& info)
 // intercepts_: [num_targets]
 // scores: X * coefficients_^T + intercepts_: [num_batches, num_targets]
 void LinearClassifier::ComputeImpl(const gsl::span<const float> input,
-                                   int64_t num_batches, int64_t num_features, int64_t num_targets,
+                                   ptrdiff_t num_batches, ptrdiff_t num_features, ptrdiff_t num_targets,
                                    const std::vector<float>& coefficients,
                                    const std::vector<float>& intercepts,
                                    Tensor& labels_output, Tensor& scores_output,
@@ -139,8 +140,10 @@ Status LinearClassifier::Compute(OpKernelContext* ctx) const {
                   "Input shape needs to be at least a single dimension.");
   }
 
-  int64_t num_batches = input_shape.NumDimensions() == 1 ? 1 : input_shape[0];
-  int64_t num_features = input_shape.NumDimensions() == 1 ? input_shape[0] : input_shape[1];
+  ptrdiff_t num_batches = input_shape.NumDimensions() == 1 ? 1 : narrow<ptrdiff_t>(input_shape[0]);
+  ptrdiff_t num_features = input_shape.NumDimensions() == 1 ? narrow<ptrdiff_t>(
+                                                                  input_shape[0])
+                                                            : narrow<ptrdiff_t>(input_shape[1]);
 
   Tensor* Y = ctx->Output(0, {num_batches});
 
diff --git a/onnxruntime/core/providers/cpu/ml/linearclassifier.h b/onnxruntime/core/providers/cpu/ml/linearclassifier.h
index 012e1cec91..e5f88dcb6b 100644
--- a/onnxruntime/core/providers/cpu/ml/linearclassifier.h
+++ b/onnxruntime/core/providers/cpu/ml/linearclassifier.h
@@ -17,7 +17,7 @@ class LinearClassifier final : public OpKernel {
   Status Compute(OpKernelContext* context) const override;
 
  private:
-  void ComputeImpl(const gsl::span<const float> input, int64_t num_batches, int64_t num_features, int64_t num_targets,
+  void ComputeImpl(const gsl::span<const float> input, ptrdiff_t num_batches, ptrdiff_t num_features, ptrdiff_t num_targets,
                    const std::vector<float>& coefficients,
                    const std::vector<float>& intercepts,
                    Tensor& labels_output,
@@ -27,7 +27,7 @@ class LinearClassifier final : public OpKernel {
                    concurrency::ThreadPool* threadpool) const;
 
   int64_t multi_class_;
-  int64_t class_count_;
+  ptrdiff_t class_count_;
   POST_EVAL_TRANSFORM post_transform_;
   bool using_strings_;
   std::vector<float> coefficients_;
diff --git a/onnxruntime/core/providers/cpu/ml/linearregressor.cc b/onnxruntime/core/providers/cpu/ml/linearregressor.cc
index a5dcdfbed8..6ed5545e70 100644
--- a/onnxruntime/core/providers/cpu/ml/linearregressor.cc
+++ b/onnxruntime/core/providers/cpu/ml/linearregressor.cc
@@ -2,6 +2,7 @@
 // Licensed under the MIT License.
 
 #include "core/providers/cpu/ml/linearregressor.h"
+#include "core/common/narrow.h"
 #include "core/providers/cpu/math/gemm.h"
 
 namespace onnxruntime {
@@ -35,7 +36,7 @@ LinearRegressor::LinearRegressor(const OpKernelInfo& info)
 // intercepts_: optional [num_targets].
 // Output: X * coefficients_^T + intercepts_: [num_batches, num_targets]
 template <typename T>
-static Status ComputeImpl(const Tensor& input, int64_t num_batches, int64_t num_features, int64_t num_targets,
+static Status ComputeImpl(const Tensor& input, ptrdiff_t num_batches, ptrdiff_t num_features, ptrdiff_t num_targets,
                           const std::vector<float>& coefficients,
                           const std::vector<float>* intercepts, Tensor& output,
                           POST_EVAL_TRANSFORM post_transform,
@@ -79,8 +80,9 @@ Status LinearRegressor::Compute(OpKernelContext* ctx) const {
                            input_shape.NumDimensions());
   }
 
-  int64_t num_batches = input_shape.NumDimensions() <= 1 ? 1 : input_shape[0];
-  int64_t num_features = input_shape.NumDimensions() <= 1 ? input_shape.Size() : input_shape[1];
+  ptrdiff_t num_batches = input_shape.NumDimensions() <= 1 ? 1 : narrow<ptrdiff_t>(input_shape[0]);
+  ptrdiff_t num_features = input_shape.NumDimensions() <= 1 ? narrow<ptrdiff_t>(input_shape.Size())
+                                                            : narrow<ptrdiff_t>(input_shape[1]);
   Tensor& Y = *ctx->Output(0, {num_batches, num_targets_});
   concurrency::ThreadPool* tp = ctx->GetOperatorThreadPool();
 
@@ -88,7 +90,7 @@ Status LinearRegressor::Compute(OpKernelContext* ctx) const {
 
   switch (element_type) {
     case ONNX_NAMESPACE::TensorProto_DataType_FLOAT: {
-      status = ComputeImpl<float>(X, num_batches, num_features, num_targets_, coefficients_,
+      status = ComputeImpl<float>(X, num_batches, num_features, narrow<ptrdiff_t>(num_targets_), coefficients_,
                                   use_intercepts_ ? &intercepts_ : nullptr,
                                   Y, post_transform_, tp);
 
diff --git a/onnxruntime/core/providers/cpu/ml/svmclassifier.cc b/onnxruntime/core/providers/cpu/ml/svmclassifier.cc
index a9a3f99152..8c356b4c62 100644
--- a/onnxruntime/core/providers/cpu/ml/svmclassifier.cc
+++ b/onnxruntime/core/providers/cpu/ml/svmclassifier.cc
@@ -47,7 +47,7 @@ SVMClassifier::SVMClassifier(const OpKernelInfo& info)
   class_count_ = 0;
   for (size_t i = 0; i < vectors_per_class_.size(); i++) {
     starting_vector_.push_back(vector_count_);
-    vector_count_ += vectors_per_class_[i];
+    vector_count_ += narrow<ptrdiff_t>(vectors_per_class_[i]);
   }
 
   using_strings_ = false;
diff --git a/onnxruntime/core/providers/cpu/ml/svmclassifier.h b/onnxruntime/core/providers/cpu/ml/svmclassifier.h
index b24a06fc2c..e2ba20e08e 100644
--- a/onnxruntime/core/providers/cpu/ml/svmclassifier.h
+++ b/onnxruntime/core/providers/cpu/ml/svmclassifier.h
@@ -32,7 +32,7 @@ class SVMCommon {
 
   template <typename T>
   void batched_kernel_dot(const gsl::span<const T> a, const gsl::span<const T> b,
-                          int64_t m, int64_t n, int64_t k,
+                          ptrdiff_t m, ptrdiff_t n, ptrdiff_t k,
                           float scalar_C,
                           const gsl::span<T> out,
                           concurrency::ThreadPool* threadpool) const {
@@ -115,9 +115,9 @@ class SVMClassifier final : public OpKernel, private SVMCommon {
   Status ComputeImpl(OpKernelContext& ctx, gsl::span<const float> x_data, const TensorShape& x_shape) const;
 
   bool weights_are_all_positive_;
-  int64_t feature_count_;
-  int64_t class_count_;
-  int64_t vector_count_;
+  ptrdiff_t feature_count_;
+  ptrdiff_t class_count_;
+  ptrdiff_t vector_count_;
   bool using_strings_;
   std::vector<int64_t> vectors_per_class_;
   std::vector<int64_t> starting_vector_;
diff --git a/onnxruntime/core/providers/cpu/ml/svmregressor.cc b/onnxruntime/core/providers/cpu/ml/svmregressor.cc
index 496589b08f..68367470a6 100644
--- a/onnxruntime/core/providers/cpu/ml/svmregressor.cc
+++ b/onnxruntime/core/providers/cpu/ml/svmregressor.cc
@@ -16,9 +16,11 @@ template <typename T>
 SVMRegressor<T>::SVMRegressor(const OpKernelInfo& info)
     : OpKernel(info),
       SVMCommon(info),
-      vector_count_(info.GetAttrOrDefault<int64_t>("n_supports", 0)),
       support_vectors_(info.GetAttrsOrDefault<float>("support_vectors")),
       post_transform_(MakeTransform(info.GetAttrOrDefault<std::string>("post_transform", "NONE"))) {
+  int64_t vector_count = 0;
+  ORT_ENFORCE(info.GetAttr<int64_t>("n_supports", &vector_count).IsOK());
+  vector_count_ = narrow<ptrdiff_t>(vector_count);
   ORT_ENFORCE(info.GetAttrs<float>("rho", rho_).IsOK());
   ORT_ENFORCE(info.GetAttrs<float>("coefficients", coefficients_).IsOK());
   ORT_ENFORCE(!coefficients_.empty());
@@ -40,9 +42,9 @@ template <typename T>
 Status SVMRegressor<T>::Compute(OpKernelContext* ctx) const {
   const auto* X = ctx->Input<Tensor>(0);
 
-  int64_t num_features = X->Shape().NumDimensions() == 1 ? X->Shape()[0] : X->Shape()[1];
-  int64_t num_batches = X->Shape().NumDimensions() == 1 ? 1 : X->Shape()[0];
-  ORT_ENFORCE(num_features == feature_count_);
+  ptrdiff_t num_features = X->Shape().NumDimensions() == 1 ? narrow<ptrdiff_t>(X->Shape()[0]) : narrow<ptrdiff_t>(X->Shape()[1]);
+  ptrdiff_t num_batches = X->Shape().NumDimensions() == 1 ? 1 : narrow<ptrdiff_t>(X->Shape()[0]);
+  ORT_RETURN_IF_NOT(num_features == feature_count_ && num_features >= 0 && num_batches >= 0, "Invalid argument");
 
   // X: [num_batches, feature_count_] where features could be coefficients or support vectors
   // coefficients_: [vector_count_]
diff --git a/onnxruntime/core/providers/cpu/ml/svmregressor.h b/onnxruntime/core/providers/cpu/ml/svmregressor.h
index bf43623884..0e72ba5db6 100644
--- a/onnxruntime/core/providers/cpu/ml/svmregressor.h
+++ b/onnxruntime/core/providers/cpu/ml/svmregressor.h
@@ -24,8 +24,8 @@ class SVMRegressor final : public OpKernel, private SVMCommon {
 
  private:
   bool one_class_;
-  int64_t feature_count_;
-  int64_t vector_count_;
+  ptrdiff_t feature_count_;
+  ptrdiff_t vector_count_;
   std::vector<float> rho_;
   std::vector<float> coefficients_;
   std::vector<float> support_vectors_;
diff --git a/onnxruntime/core/providers/rocm/math/gemm.cc b/onnxruntime/core/providers/rocm/math/gemm.cc
index 60705429f0..8252741830 100644
--- a/onnxruntime/core/providers/rocm/math/gemm.cc
+++ b/onnxruntime/core/providers/rocm/math/gemm.cc
@@ -72,9 +72,9 @@ Status Gemm<T>::ComputeInternal(OpKernelContext* ctx) const {
   if (!helper.State().IsOK())
     return helper.State();
 
-  int M = gsl::narrow_cast<int>(helper.M());
-  int N = gsl::narrow_cast<int>(helper.N());
-  int K = gsl::narrow_cast<int>(helper.K());
+  ptrdiff_t M = helper.M();
+  ptrdiff_t N = helper.N();
+  ptrdiff_t K = helper.K();
   auto* Y = ctx->Output(0, {M, N});
   HipT* out_data = reinterpret_cast<HipT*>(Y->MutableData<T>());
 
diff --git a/onnxruntime/core/util/math_cpu.cc b/onnxruntime/core/util/math_cpu.cc
index 06cdb0a646..b4554e0f0d 100644
--- a/onnxruntime/core/util/math_cpu.cc
+++ b/onnxruntime/core/util/math_cpu.cc
@@ -80,6 +80,50 @@ void Gemm<float, ThreadPool>(CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, ptr
   MlasGemm(TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, N, threadpool);
 }
 
+template <>
+void Gemm<Eigen::half, ThreadPool>(CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, ptrdiff_t M,
+                                   ptrdiff_t N, ptrdiff_t K, Eigen::half alpha, const Eigen::half* A, const Eigen::half* B, Eigen::half beta,
+                                   Eigen::half* C, ThreadPool*) {
+  auto C_mat = EigenMatrixMap<Eigen::half>(C, N, M);
+  if (beta == static_cast<Eigen::half>(0)) {
+    C_mat.setZero();
+  } else {
+    C_mat *= beta;
+  }
+  switch (TransA) {
+    case CblasNoTrans: {
+      switch (TransB) {
+        case CblasNoTrans:
+          C_mat.noalias() += alpha * (ConstEigenMatrixMap<Eigen::half>(B, N, K) *
+                                      ConstEigenMatrixMap<Eigen::half>(A, K, M));
+          return;
+        case CblasTrans:
+          C_mat.noalias() += alpha * (ConstEigenMatrixMap<Eigen::half>(B, K, N).transpose() *
+                                      ConstEigenMatrixMap<Eigen::half>(A, K, M));
+          return;
+        default:
+          ORT_THROW("CblasNoTrans Unexpected CBLAS_TRANSPOSE for TransB of ", TransB);
+      }
+    }
+    case CblasTrans: {
+      switch (TransB) {
+        case CblasNoTrans:
+          C_mat.noalias() += alpha * (ConstEigenMatrixMap<Eigen::half>(B, N, K) *
+                                      ConstEigenMatrixMap<Eigen::half>(A, M, K).transpose());
+          return;
+        case CblasTrans:
+          C_mat.noalias() += alpha * (ConstEigenMatrixMap<Eigen::half>(B, K, N).transpose() *
+                                      ConstEigenMatrixMap<Eigen::half>(A, M, K).transpose());
+          return;
+        default:
+          ORT_THROW("CblasTrans Unexpected CBLAS_TRANSPOSE for TransB of ", TransB);
+      }
+    }
+    default:
+      ORT_THROW("Unexpected CBLAS_TRANSPOSE for TransA of ", TransA);
+  }
+}
+
 #ifdef MLAS_SUPPORTS_GEMM_DOUBLE
 template <>
 void Gemm<double, ThreadPool>(CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, ptrdiff_t M,
diff --git a/onnxruntime/test/providers/cpu/math/gemm_test.cc b/onnxruntime/test/providers/cpu/math/gemm_test.cc
index 068fc2a505..36ab867f1b 100644
--- a/onnxruntime/test/providers/cpu/math/gemm_test.cc
+++ b/onnxruntime/test/providers/cpu/math/gemm_test.cc
@@ -2,6 +2,7 @@
 // Licensed under the MIT License.
 
 #include "gtest/gtest.h"
+#include "core/mlas/inc/mlas.h"
 #include "core/framework/run_options.h"
 #include "test/common/cuda_op_test_utils.h"
 #include "test/providers/provider_test_utils.h"
@@ -24,45 +25,7 @@ const constexpr auto run_with_tunable_op = &run_options;
 
 }  // namespace
 
-template <typename T>
-void TestGemmNoTrans() {
-  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
-    OpTester test("Gemm");
-
-    test.AddAttribute("transA", (int64_t)0);
-    test.AddAttribute("transB", (int64_t)0);
-    test.AddAttribute("alpha", 1.0f);
-    test.AddAttribute("beta", 1.0f);
-
-    test.AddInput<T>("A", {2, 4},
-                     {1.0f, 2.0f, 3.0f, 4.0f,
-                      -1.0f, -2.0f, -3.0f, -4.0f});
-    test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f), b_is_initializer);
-    test.AddInput<T>("C", {2, 3}, std::vector<T>(6, 1.0f), c_is_initializer);
-    test.AddOutput<T>("Y", {2, 3},
-                      {11.0f, 11.0f, 11.0f,
-                       -9.0f, -9.0f, -9.0f});
-    test.Config(run_with_tunable_op)
-        .RunWithConfig();
-  };
-
-  run_test(false, false);
-  // NNAPI EP requires weight to be an initializer
-  run_test(true, false);
-  // CoreML EP requires weight and bias both to be initializers
-  run_test(true, true);
-}
-
-TEST(GemmOpTest, GemmNoTrans_float) {
-  TestGemmNoTrans<float>();
-}
-
-TEST(GemmOpTest, GemmNoTrans_double) {
-  TestGemmNoTrans<double>();
-}
-
 // Only CUDA and ROCM kernel has float 16 support
-#if defined(USE_CUDA) || defined(USE_ROCM)
 TEST(GemmOpTest, GemmNoTrans_f16) {
 #ifdef USE_CUDA
   int min_cuda_architecture = 530;
@@ -71,7 +34,7 @@ TEST(GemmOpTest, GemmNoTrans_f16) {
     return;
   }
 #endif
-  OpTester test("Gemm");
+  OpTester test("Gemm", 13);
 
   test.AddAttribute("transA", (int64_t)0);
   test.AddAttribute("transB", (int64_t)0);
@@ -102,7 +65,6 @@ TEST(GemmOpTest, GemmNoTrans_f16) {
       .Config(run_with_tunable_op)
       .RunWithConfig();
 }
-#endif
 
 #if defined(USE_CUDA) || defined(USE_ROCM) || defined(USE_DNNL)
 TEST(GemmOpTest, GemmNoTrans_bfloat16) {
@@ -147,270 +109,6 @@ TEST(GemmOpTest, GemmNoTrans_bfloat16) {
 }
 #endif  // USE_CUDA USE_RCOM USE_DNNL
 
-template <typename T>
-void TestGemmBroadcast() {
-  auto run_test = [](bool b_is_initializer, bool c_is_initializer) {
-    OpTester test("Gemm");
-
-    test.AddAttribute("transA", (int64_t)0);
-    test.AddAttribute("transB", (int64_t)0);
-    test.AddAttribute("alpha", 1.0f);
-    test.AddAttribute("beta", 1.0f);
-
-    test.AddInput<T>("A", {2, 4},
-                     {1.0f, 2.0f, 3.0f, 4.0f,
-                      -1.0f, -2.0f, -3.0f, -4.0f});
-    test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f), b_is_initializer);
-    test.AddInput<T>("C", {3}, std::vector<T>{1.0f, 2.0f, 3.0f}, c_is_initializer);
-    test.AddOutput<T>("Y", {2, 3},
-                      {11.0f, 12.0f, 13.0f,
-                       -9.0f, -8.0f, -7.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#endif
-    test.Config(run_with_tunable_op)
-        .RunWithConfig();
-  };
-
-  run_test(false, false);
-  // NNAPI EP requires weight to be an initializer
-  run_test(true, false);
-  // CoreML EP requires weight and bias both to be initializers
-  run_test(true, true);
-}
-
-TEST(GemmOpTest, GemmBroadcast) {
-  TestGemmBroadcast<float>();
-  TestGemmBroadcast<double>();
-}
-
-template <typename T>
-static void TestGemmTrans() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)1);
-  test.AddAttribute("transB", (int64_t)1);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 1.0f);
-
-  test.AddInput<T>("A", {4, 2},
-                   {1.0f, -1.0f,
-                    2.0f, -2.0f,
-                    3.0f, -3.0f,
-                    4.0f, -4.0f});
-  test.AddInput<T>("B", {3, 4}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {3}, std::vector<T>(3, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {11.0f, 11.0f, 11.0f,
-                     -9.0f, -9.0f, -9.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#endif
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmTrans) {
-  TestGemmTrans<float>();
-  TestGemmTrans<double>();
-}
-
-// NNAPI EP's GEMM only works as A*B', add case only B is transposed
-// Also test NNAPI EP's handling of non-1D bias (C of Gemm)
-template <typename T>
-static void TestGemmTransB() {
-  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
-    OpTester test("Gemm");
-
-    test.AddAttribute("transA", (int64_t)0);
-    test.AddAttribute("transB", (int64_t)1);
-    test.AddAttribute("alpha", 1.0f);
-    test.AddAttribute("beta", 1.0f);
-
-    test.AddInput<T>("A", {2, 4},
-                     {1.0f, 2.0f, 3.0f, 4.0f,
-                      -1.0f, -2.0f, -3.0f, -4.0f});
-    test.AddInput<T>("B", {3, 4}, std::vector<T>(12, 1.0f), b_is_initializer);
-    test.AddInput<T>("C", {1, 3}, std::vector<T>(3, 1.0f), c_is_initializer);
-    test.AddOutput<T>("Y", {2, 3},
-                      {11.0f, 11.0f, 11.0f,
-                       -9.0f, -9.0f, -9.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#endif
-    test.Config(run_with_tunable_op)
-        .RunWithConfig();
-  };
-  run_test(false, false);
-  // CoreML EP requires weight and bias both to be initializers
-  run_test(true, true);
-}
-
-TEST(GemmOpTest, GemmTransB) {
-  TestGemmTransB<float>();
-  TestGemmTransB<double>();
-}
-
-// NNAPI EP's GEMM only works as A*B', add case only B is transposed
-// Also test NNAPI EP's handling of non-1D bias (C of Gemm) which is broadcastable but not valid for NNAPI
-template <typename T>
-static void TestGemmTransB_1() {
-  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
-    OpTester test("Gemm");
-
-    test.AddAttribute("transA", (int64_t)0);
-    test.AddAttribute("transB", (int64_t)1);
-    test.AddAttribute("alpha", 1.0f);
-    test.AddAttribute("beta", 1.0f);
-
-    test.AddInput<T>("A", {2, 4},
-                     {1.0f, 2.0f, 3.0f, 4.0f,
-                      -1.0f, -2.0f, -3.0f, -4.0f});
-    test.AddInput<T>("B", {3, 4}, std::vector<T>(12, 1.0f), b_is_initializer);
-    test.AddInput<T>("C", {2, 1}, std::vector<T>(2, 1.0f), c_is_initializer);
-    test.AddOutput<T>("Y", {2, 3},
-                      {11.0f, 11.0f, 11.0f,
-                       -9.0f, -9.0f, -9.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#endif
-    test.Config(run_with_tunable_op)
-        .RunWithConfig();
-  };
-  run_test(false, false);
-  // CoreML EP requires weight and bias both to be initializers
-  run_test(true, true);
-}
-
-TEST(GemmOpTest, GemmTransB_1) {
-  TestGemmTransB_1<float>();
-  TestGemmTransB_1<double>();
-}
-
-template <typename T>
-void TestGemmAlpha() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 0.5f);
-  test.AddAttribute("beta", 1.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {3}, std::vector<T>(3, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {6.0f, 6.0f, 6.0f,
-                     -4.0f, -4.0f, -4.0f});
-  // test.AddOutput<T>("Y", {2, 3},
-  //                   {5.0f, 5.0f, 5.0f,
-  //                    -5.0f, -5.0f, -5.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#else
-  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
-#endif
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmAlpha) {
-  TestGemmAlpha<float>();
-  TestGemmAlpha<double>();
-}
-
-template <typename T>
-void TestGemmBeta() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 2.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {3}, std::vector<T>(3, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {12.0f, 12.0f, 12.0f,
-                     -8.0f, -8.0f, -8.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#else
-  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
-#endif
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmBeta) {
-  TestGemmBeta<float>();
-  TestGemmBeta<double>();
-}
-
-template <typename T>
-void TestGemmAlphaBeta() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 0.5f);
-  test.AddAttribute("beta", 2.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {3}, std::vector<T>(3, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {7.0f, 7.0f, 7.0f,
-                     -3.0f, -3.0f, -3.0f});
-#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
-  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
-#else
-  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
-#endif
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmAlphaBeta) {
-  TestGemmAlphaBeta<float>();
-  TestGemmAlphaBeta<double>();
-}
-
-template <typename T>
-void TestGemmNaN() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 0.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {2, 3}, std::vector<T>(6, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {10.0f, 10.0f, 10.0f,
-                     -10.0f, -10.0f, -10.0f});
-
-  // TensorRT: Seg fault in parser
-  test.ConfigExcludeEps({kTensorrtExecutionProvider})
-      .Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmNaN) {
-  TestGemmNaN<float>();
-  TestGemmNaN<double>();
-}
 #if defined(USE_DNNL)
 TEST(GemmOpTest, GemmNaN_bfloat16) {
 #ifdef USE_DNNL
@@ -444,32 +142,6 @@ TEST(GemmOpTest, GemmNaN_bfloat16) {
 }
 #endif  //  USE_DNNL
 
-template <typename T>
-void TestGemmScalarBroadcast() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 1.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {1}, std::vector<T>{1.0f});
-  test.AddOutput<T>("Y", {2, 3},
-                    {11.0f, 11.0f, 11.0f,
-                     -9.0f, -9.0f, -9.0f});
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmScalarBroadcast) {
-  TestGemmScalarBroadcast<float>();
-  TestGemmScalarBroadcast<double>();
-}
-
 #if defined(USE_DNNL)
 TEST(GemmOpTest, GemmScalarBroadcast_bfloat16) {
 #ifdef USE_DNNL
@@ -502,31 +174,6 @@ TEST(GemmOpTest, GemmScalarBroadcast_bfloat16) {
 }
 #endif  //  USE_DNNL
 
-template <typename T>
-void TestGemm2DBroadcast_1() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 1.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {2, 1}, std::vector<T>{1.0, 2.0f});
-  test.AddOutput<T>("Y", {2, 3},
-                    {11.0f, 11.0f, 11.0f,
-                     -8.0f, -8.0f, -8.0f});
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, Gemm2DBroadcast_1) {
-  TestGemm2DBroadcast_1<float>();
-  TestGemm2DBroadcast_1<double>();
-}
 #if defined(USE_DNNL)
 TEST(GemmOpTest, Gemm2DBroadcast_1_bfloat16) {
 #ifdef USE_DNNL
@@ -559,33 +206,6 @@ TEST(GemmOpTest, Gemm2DBroadcast_1_bfloat16) {
 }
 #endif  //  USE_DNNL
 
-template <typename T>
-void TestGemm2DBroadcast_2() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 1.0f);
-
-  // Same as GemmBroadcast, but adding the unnecessary second dimension.
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {1, 3}, std::vector<T>{1.0f, 2.0f, 3.0f});
-  test.AddOutput<T>("Y", {2, 3},
-                    {11.0f, 12.0f, 13.0f,
-                     -9.0f, -8.0f, -7.0f});
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, Gemm2DBroadcast_2) {
-  TestGemm2DBroadcast_2<float>();
-  TestGemm2DBroadcast_2<double>();
-}
-
 #if defined(USE_DNNL)
 TEST(GemmOpTest, Gemm2DBroadcast_2_bfloat16) {
 #ifdef USE_DNNL
@@ -618,32 +238,6 @@ TEST(GemmOpTest, Gemm2DBroadcast_2_bfloat16) {
 }
 #endif  //  USE_DNNL
 
-template <typename T>
-void TestGemmFalseBroadcast() {
-  OpTester test("Gemm");
-
-  test.AddAttribute("transA", (int64_t)0);
-  test.AddAttribute("transB", (int64_t)0);
-  test.AddAttribute("alpha", 1.0f);
-  test.AddAttribute("beta", 1.0f);
-
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {2, 3}, std::vector<T>{1.0f, 1.0f, 1.0f, 2.0f, 2.0f, 2.0f});
-  test.AddOutput<T>("Y", {2, 3},
-                    {11.0f, 11.0f, 11.0f,
-                     -8.0f, -8.0f, -8.0f});
-  test.Config(run_with_tunable_op)
-      .RunWithConfig();
-}
-
-TEST(GemmOpTest, GemmFalseBroadcast) {
-  TestGemmFalseBroadcast<float>();
-  TestGemmFalseBroadcast<double>();
-}
-
 #if defined(USE_DNNL)
 TEST(GemmOpTest, GemmFalseBroadcast_2_bfloat16) {
 #ifdef USE_DNNL
@@ -675,34 +269,363 @@ TEST(GemmOpTest, GemmFalseBroadcast_2_bfloat16) {
   test.Run(OpTester::ExpectResult::kExpectSuccess, "", {}, nullptr, &execution_providers);
 }
 #endif  //  USE_DNNL
+
 template <typename T>
-void TestGemmEmptyTensor() {
+class GemmOpTypedTests : public ::testing::Test {
+};
+
+// On CPUs without fp16 instructions the tests will output a warning:
+// "registered execution providers CPUExecutionProvider were unable to run the model"
+// , then they will still pass.
+using GemmOpTypedTestsTypes = ::testing::Types<float, double, MLFloat16>;
+TYPED_TEST_SUITE(GemmOpTypedTests, GemmOpTypedTestsTypes);
+
+TYPED_TEST(GemmOpTypedTests, TestGemmScalarBroadcast) {
   OpTester test("Gemm");
 
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 1.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {1}, std::vector<TypeParam>{static_cast<TypeParam>(1.0f)});
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                             static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f)});
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+TYPED_TEST(GemmOpTypedTests, TestGemm2DBroadcast_2) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 1.0f);
+
+  // Same as GemmBroadcast, but adding the unnecessary second dimension.
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {1, 3}, std::vector<TypeParam>{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f)});
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(12.0f), static_cast<TypeParam>(13.0f),
+                             static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-7.0f)});
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmFalseBroadcast) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 1.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {2, 3}, std::vector<TypeParam>{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(1.0f), static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(2.0f)});
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                             static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f)});
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmBroadcast) {
+  auto run_test = [](bool b_is_initializer, bool c_is_initializer) {
+    OpTester test("Gemm");
+
+    test.AddAttribute("transA", (int64_t)0);
+    test.AddAttribute("transB", (int64_t)0);
+    test.AddAttribute("alpha", 1.0f);
+    test.AddAttribute("beta", 1.0f);
+
+    test.AddInput<TypeParam>("A", {2, 4},
+                             {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                              static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+    test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)), b_is_initializer);
+    test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f)}, c_is_initializer);
+    test.AddOutput<TypeParam>("Y", {2, 3},
+                              {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(12.0f), static_cast<TypeParam>(13.0f),
+                               static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-7.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#endif
+    test.Config(run_with_tunable_op)
+        .RunWithConfig();
+  };
+
+  run_test(false, false);
+  // NNAPI EP requires weight to be an initializer
+  run_test(true, false);
+  // CoreML EP requires weight and bias both to be initializers
+  run_test(true, true);
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmTrans) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)1);
+  test.AddAttribute("transB", (int64_t)1);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 1.0f);
+
+  test.AddInput<TypeParam>("A", {4, 2},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(-1.0f),
+                            static_cast<TypeParam>(2.0f), static_cast<TypeParam>(-2.0f),
+                            static_cast<TypeParam>(3.0f), static_cast<TypeParam>(-3.0f),
+                            static_cast<TypeParam>(4.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {3, 4}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                             static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#endif
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+// NNAPI EP's GEMM only works as A*B', add case only B is transposed
+// Also test NNAPI EP's handling of non-1D bias (C of Gemm)
+TYPED_TEST(GemmOpTypedTests, TestGemmTransB) {
+  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
+    OpTester test("Gemm");
+
+    test.AddAttribute("transA", (int64_t)0);
+    test.AddAttribute("transB", (int64_t)1);
+    test.AddAttribute("alpha", 1.0f);
+    test.AddAttribute("beta", 1.0f);
+
+    test.AddInput<TypeParam>("A", {2, 4},
+                             {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                              static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+    test.AddInput<TypeParam>("B", {3, 4}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)), b_is_initializer);
+    test.AddInput<TypeParam>("C", {1, 3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)), c_is_initializer);
+    test.AddOutput<TypeParam>("Y", {2, 3},
+                              {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                               static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#endif
+    test.Config(run_with_tunable_op)
+        .RunWithConfig();
+  };
+  run_test(false, false);
+  // CoreML EP requires weight and bias both to be initializers
+  run_test(true, true);
+}
+
+// NNAPI EP's GEMM only works as A*B', add case only B is transposed
+// Also test NNAPI EP's handling of non-1D bias (C of Gemm) which is broadcastable but not valid for NNAPI
+TYPED_TEST(GemmOpTypedTests, TestGemmTransB_1) {
+  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
+    OpTester test("Gemm");
+
+    test.AddAttribute("transA", (int64_t)0);
+    test.AddAttribute("transB", (int64_t)1);
+    test.AddAttribute("alpha", 1.0f);
+    test.AddAttribute("beta", 1.0f);
+
+    test.AddInput<TypeParam>("A", {2, 4},
+                             {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                              static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+    test.AddInput<TypeParam>("B", {3, 4}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)), b_is_initializer);
+    test.AddInput<TypeParam>("C", {2, 1}, std::vector<TypeParam>(2, static_cast<TypeParam>(1.0f)), c_is_initializer);
+    test.AddOutput<TypeParam>("Y", {2, 3},
+                              {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                               static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+    test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#endif
+    test.Config(run_with_tunable_op)
+        .RunWithConfig();
+  };
+  run_test(false, false);
+  // CoreML EP requires weight and bias both to be initializers
+  run_test(true, true);
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmAlpha) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 0.5f);
+  test.AddAttribute("beta", 1.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(6.0f), static_cast<TypeParam>(6.0f), static_cast<TypeParam>(6.0f),
+                             static_cast<TypeParam>(-4.0f), static_cast<TypeParam>(-4.0f), static_cast<TypeParam>(-4.0f)});
+  // test.AddOutput<TypeParam>("Y", {2, 3},
+  //                   {5.0f, 5.0f, 5.0f,
+  //                    -5.0f, -5.0f, -5.0f});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#else
+  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
+#endif
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmBeta) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 2.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(12.0f), static_cast<TypeParam>(12.0f), static_cast<TypeParam>(12.0f),
+                             static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#else
+  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
+#endif
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmNaN) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 0.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {2, 3}, std::vector<TypeParam>(6, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(10.0f), static_cast<TypeParam>(10.0f), static_cast<TypeParam>(10.0f),
+                             static_cast<TypeParam>(-10.0f), static_cast<TypeParam>(-10.0f), static_cast<TypeParam>(-10.0f)});
+
+  // TensorRT: Seg fault in parser
+  test.ConfigExcludeEps({kTensorrtExecutionProvider})
+      .Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmAlphaBeta) {
+  OpTester test("Gemm", 13);
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 0.5f);
+  test.AddAttribute("beta", 2.0f);
+
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(7.0f), static_cast<TypeParam>(7.0f), static_cast<TypeParam>(7.0f),
+                             static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-3.0f)});
+#if defined(OPENVINO_CONFIG_GPU_FP16) || defined(OPENVINO_CONFIG_GPU_FP32)
+  test.ConfigExcludeEps({kOpenVINOExecutionProvider});  // OpenVINO: Temporarily disabled due to accuracy issues
+#else
+  test.ConfigExcludeEps({kTensorrtExecutionProvider});  // TensorRT: Seg fault in parser
+#endif
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+// C is 1D
+TYPED_TEST(GemmOpTypedTests, TestGemm2DBroadcast_1) {
+  OpTester test("Gemm");
+
+  test.AddAttribute("transA", (int64_t)0);
+  test.AddAttribute("transB", (int64_t)0);
+  test.AddAttribute("alpha", 1.0f);
+  test.AddAttribute("beta", 1.0f);
+  std::array<TypeParam, 8> a_data{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                                  static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)};
+
+  test.AddInput<TypeParam>("A", {2, 4}, a_data.data(), a_data.size());
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {2, 1}, std::vector<TypeParam>{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f)});
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                             static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f), static_cast<TypeParam>(-8.0f)});
+  test.Config(run_with_tunable_op)
+      .RunWithConfig();
+}
+
+TYPED_TEST(GemmOpTypedTests, TestGemmNoTrans) {
+  auto run_test = [](bool b_is_initializer, bool c_is_initializer = false) {
+    OpTester test("Gemm", 13);
+
+    test.AddAttribute("transA", (int64_t)0);
+    test.AddAttribute("transB", (int64_t)0);
+    test.AddAttribute("alpha", 1.0f);
+    test.AddAttribute("beta", 1.0f);
+
+    std::array<TypeParam, 8> a_data{static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                                    static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)};
+
+    test.AddInput<TypeParam>("A", {2, 4}, a_data.data(), a_data.size());
+    test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)), b_is_initializer);
+    test.AddInput<TypeParam>("C", {2, 3}, std::vector<TypeParam>(6, static_cast<TypeParam>(1.0f)), c_is_initializer);
+    test.AddOutput<TypeParam>("Y", {2, 3},
+                              {static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f), static_cast<TypeParam>(11.0f),
+                               static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f), static_cast<TypeParam>(-9.0f)});
+    test.Config(run_with_tunable_op)
+        .RunWithConfig();
+  };
+
+  run_test(false, false);
+  // NNAPI EP requires weight to be an initializer
+  run_test(true, false);
+  // CoreML EP requires weight and bias both to be initializers
+  run_test(true, true);
+}
+TYPED_TEST(GemmOpTypedTests, GemmEmptyTensor) {
+  OpTester test("Gemm", 13);
+
   test.AddAttribute("transA", static_cast<int64_t>(0));
   test.AddAttribute("transB", static_cast<int64_t>(0));
   test.AddAttribute("alpha", 1.0f);
   test.AddAttribute("beta", 1.0f);
 
-  test.AddInput<T>("A", {0, 4},
-                   {});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddInput<T>("C", {3}, std::vector<T>(3, 1.0f));
-  test.AddOutput<T>("Y", {0, 3},
-                    {});
+  test.AddInput<TypeParam>("A", {0, 4},
+                           {});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddInput<TypeParam>("C", {3}, std::vector<TypeParam>(3, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {0, 3},
+                            {});
   // TensorRT: doesn't support dynamic shape yet
   test.ConfigExcludeEps({kTensorrtExecutionProvider, kDnnlExecutionProvider, kQnnExecutionProvider})
       .Config(run_with_tunable_op)
       .RunWithConfig();
 }
-
-TEST(GemmOpTest, GemmEmptyTensor) {
-  TestGemmEmptyTensor<float>();
-  TestGemmEmptyTensor<double>();
-}
-
-template <typename T>
-static void TestGemmNoBiasOpset11() {
+TYPED_TEST(GemmOpTypedTests, MissingBias) {
   OpTester test("Gemm", 11);
 
   test.AddAttribute("transA", static_cast<int64_t>(0));
@@ -710,13 +633,13 @@ static void TestGemmNoBiasOpset11() {
   test.AddAttribute("alpha", 1.0f);
   test.AddAttribute("beta", 1.0f);
 
-  test.AddInput<T>("A", {2, 4},
-                   {1.0f, 2.0f, 3.0f, 4.0f,
-                    -1.0f, -2.0f, -3.0f, -4.0f});
-  test.AddInput<T>("B", {4, 3}, std::vector<T>(12, 1.0f));
-  test.AddOutput<T>("Y", {2, 3},
-                    {10.0f, 10.0f, 10.0f,
-                     -10.0f, -10.0f, -10.0f});
+  test.AddInput<TypeParam>("A", {2, 4},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f),
+                            static_cast<TypeParam>(-1.0f), static_cast<TypeParam>(-2.0f), static_cast<TypeParam>(-3.0f), static_cast<TypeParam>(-4.0f)});
+  test.AddInput<TypeParam>("B", {4, 3}, std::vector<TypeParam>(12, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 3},
+                            {static_cast<TypeParam>(10.0f), static_cast<TypeParam>(10.0f), static_cast<TypeParam>(10.0f),
+                             static_cast<TypeParam>(-10.0f), static_cast<TypeParam>(-10.0f), static_cast<TypeParam>(-10.0f)});
   // tensorRT don't seem to support missing bias
   std::unordered_set<std::string> excluded_provider_types{kTensorrtExecutionProvider};
   // QNN Linux result diff 0.011714935302734375 exceed the threshold
@@ -728,33 +651,22 @@ static void TestGemmNoBiasOpset11() {
       .RunWithConfig();
 }
 
-TEST(GemmOpTest, GemmNoBiasOpset11) {
-  TestGemmNoBiasOpset11<float>();
-  TestGemmNoBiasOpset11<double>();
-}
-
-template <typename T>
-static void TestGemmWithAlphaOpset11() {
+TYPED_TEST(GemmOpTypedTests, TestGemmWithAlphaOpset11) {
   OpTester test("Gemm", 11);
 
   test.AddAttribute("alpha", 2.0f);
 
-  test.AddInput<T>("A", {2, 2},
-                   {1.0f, 2.0f, 3.0f, 4.0f});
-  test.AddInput<T>("B", {2, 2}, std::vector<T>(4, 1.0f));
-  test.AddOutput<T>("Y", {2, 2},
-                    {6.0f, 6.0f, 14.0f, 14.0f});
+  test.AddInput<TypeParam>("A", {2, 2},
+                           {static_cast<TypeParam>(1.0f), static_cast<TypeParam>(2.0f), static_cast<TypeParam>(3.0f), static_cast<TypeParam>(4.0f)});
+  test.AddInput<TypeParam>("B", {2, 2}, std::vector<TypeParam>(4, static_cast<TypeParam>(1.0f)));
+  test.AddOutput<TypeParam>("Y", {2, 2},
+                            {static_cast<TypeParam>(6.0f), static_cast<TypeParam>(6.0f), static_cast<TypeParam>(14.0f), static_cast<TypeParam>(14.0f)});
   // tensorRT don't seem to support missing bias
   test.ConfigExcludeEps({kTensorrtExecutionProvider})
       .Config(run_with_tunable_op)
       .RunWithConfig();
 }
 
-TEST(GemmOpTest, GemmWithAlphaOpset11) {
-  TestGemmWithAlphaOpset11<float>();
-  TestGemmWithAlphaOpset11<double>();
-}
-
 #ifndef ENABLE_TRAINING
 // Prepacking is disabled in training builds so no need to test the feature in a training build.
 TEST(GemmOpTest, SharedPrepackedWeights) {
diff --git a/onnxruntime/test/providers/cpu/model_tests.cc b/onnxruntime/test/providers/cpu/model_tests.cc
index 26f8cba2fc..6c5d1399f5 100644
--- a/onnxruntime/test/providers/cpu/model_tests.cc
+++ b/onnxruntime/test/providers/cpu/model_tests.cc
@@ -1125,7 +1125,10 @@ TEST_P(ModelTest, Run) {
                                                     ORT_TSTR("SSD")};
     all_disabled_tests.insert(std::begin(x86_disabled_tests), std::end(x86_disabled_tests));
 #endif
-
+    // fp16 models have different outputs with different kinds of hardware. We need to disable all fp16 models
+    all_disabled_tests.insert(ORT_TSTR("fp16_shufflenet"));
+    all_disabled_tests.insert(ORT_TSTR("fp16_inception_v1"));
+    all_disabled_tests.insert(ORT_TSTR("fp16_tiny_yolov2"));
     std::vector<std::basic_string<ORTCHAR_T>> paths;
 #if defined(NDEBUG) || defined(RUN_MODELTEST_IN_DEBUG_MODE)
 #ifdef _WIN32