diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8acc4098d6a..612b04ceac3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -33,6 +33,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 
 * Fixed `flags[key]` of `dpnp.ndarray` and `dpnp.tensor.usm_ndarray` to raise `KeyError` for an unknown key [#3141](https://github.com/IntelPython/dpnp/pull/3141)
+* Released the GIL around the oneMKL VM calls in the VM extension, which deadlocked on an in-order queue of a CPU device [#3156](https://github.com/IntelPython/dpnp/pull/3156)
 
 ### Security
 
diff --git a/dpnp/libtensor/extensions/vm/abs.cpp b/dpnp/libtensor/extensions/vm/abs.cpp
index 2c2ade041fb..51e47e23eee 100644
--- a/dpnp/libtensor/extensions/vm/abs.cpp
+++ b/dpnp/libtensor/extensions/vm/abs.cpp
@@ -94,6 +94,8 @@ static sycl::event abs_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::abs(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/acos.cpp b/dpnp/libtensor/extensions/vm/acos.cpp
index 5242a01bbf3..24cdbc83d65 100644
--- a/dpnp/libtensor/extensions/vm/acos.cpp
+++ b/dpnp/libtensor/extensions/vm/acos.cpp
@@ -94,6 +94,8 @@ static sycl::event acos_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::acos(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/acosh.cpp b/dpnp/libtensor/extensions/vm/acosh.cpp
index 3518a3ec94e..8e68b0f4508 100644
--- a/dpnp/libtensor/extensions/vm/acosh.cpp
+++ b/dpnp/libtensor/extensions/vm/acosh.cpp
@@ -94,6 +94,8 @@ static sycl::event acosh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::acosh(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/add.cpp b/dpnp/libtensor/extensions/vm/add.cpp
index 49c0a2faf11..2f418b73153 100644
--- a/dpnp/libtensor/extensions/vm/add.cpp
+++ b/dpnp/libtensor/extensions/vm/add.cpp
@@ -113,6 +113,8 @@ static sycl::event add_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::add(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/arg.cpp b/dpnp/libtensor/extensions/vm/arg.cpp
index 75ab640dbf8..42d54dd57b3 100644
--- a/dpnp/libtensor/extensions/vm/arg.cpp
+++ b/dpnp/libtensor/extensions/vm/arg.cpp
@@ -92,6 +92,8 @@ static sycl::event arg_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::arg(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/asin.cpp b/dpnp/libtensor/extensions/vm/asin.cpp
index 3ad3f42c097..a3e211ba3f9 100644
--- a/dpnp/libtensor/extensions/vm/asin.cpp
+++ b/dpnp/libtensor/extensions/vm/asin.cpp
@@ -94,6 +94,8 @@ static sycl::event asin_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::asin(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/asinh.cpp b/dpnp/libtensor/extensions/vm/asinh.cpp
index 7c3788ebebb..4b62f40bc6f 100644
--- a/dpnp/libtensor/extensions/vm/asinh.cpp
+++ b/dpnp/libtensor/extensions/vm/asinh.cpp
@@ -94,6 +94,8 @@ static sycl::event asinh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::asinh(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/atan.cpp b/dpnp/libtensor/extensions/vm/atan.cpp
index b8df2c0390d..cc3c8940cee 100644
--- a/dpnp/libtensor/extensions/vm/atan.cpp
+++ b/dpnp/libtensor/extensions/vm/atan.cpp
@@ -94,6 +94,8 @@ static sycl::event atan_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::atan(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/atan2.cpp b/dpnp/libtensor/extensions/vm/atan2.cpp
index 643668613fb..672f761634c 100644
--- a/dpnp/libtensor/extensions/vm/atan2.cpp
+++ b/dpnp/libtensor/extensions/vm/atan2.cpp
@@ -102,6 +102,8 @@ static sycl::event atan2_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::atan2(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/atanh.cpp b/dpnp/libtensor/extensions/vm/atanh.cpp
index 2a71dfcad88..8aebc141f7f 100644
--- a/dpnp/libtensor/extensions/vm/atanh.cpp
+++ b/dpnp/libtensor/extensions/vm/atanh.cpp
@@ -94,6 +94,8 @@ static sycl::event atanh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::atanh(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/cbrt.cpp b/dpnp/libtensor/extensions/vm/cbrt.cpp
index 3a37a8f3050..da7a49c14e7 100644
--- a/dpnp/libtensor/extensions/vm/cbrt.cpp
+++ b/dpnp/libtensor/extensions/vm/cbrt.cpp
@@ -91,6 +91,8 @@ static sycl::event cbrt_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::cbrt(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/ceil.cpp b/dpnp/libtensor/extensions/vm/ceil.cpp
index 88698445a34..1bb59632dd4 100644
--- a/dpnp/libtensor/extensions/vm/ceil.cpp
+++ b/dpnp/libtensor/extensions/vm/ceil.cpp
@@ -91,6 +91,8 @@ static sycl::event ceil_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::ceil(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/common.hpp b/dpnp/libtensor/extensions/vm/common.hpp
index a350013ff24..90fc5f00fed 100644
--- a/dpnp/libtensor/extensions/vm/common.hpp
+++ b/dpnp/libtensor/extensions/vm/common.hpp
@@ -50,7 +50,14 @@ namespace ext_ns = ext::common;
 namespace py = pybind11;
 namespace td_ns = dpnp::tensor::type_dispatch;
 
-namespace dpnp::extensions::vm::py_internal
+namespace dpnp::extensions::vm
+{
+// Releases the GIL around a oneMKL VM call. On a CPU device VM waits on its
+// `depends` and computes on the calling thread, so holding the GIL can
+// deadlock with host tasks that need it (e.g. on an in-order queue).
+using vm_gil_release = py::gil_scoped_release;
+
+namespace py_internal
 {
 template <typename output_typesT, typename contig_dispatchT>
 bool need_to_call_unary_ufunc(sycl::queue &exec_q,
@@ -487,4 +494,5 @@ bool need_to_call_binary_ufunc(sycl::queue &exec_q,
                                     __name__##_ContigFactory>(                 \
             contig_dispatch_vector);                                           \
     };
-} // namespace dpnp::extensions::vm::py_internal
+} // namespace py_internal
+} // namespace dpnp::extensions::vm
diff --git a/dpnp/libtensor/extensions/vm/conj.cpp b/dpnp/libtensor/extensions/vm/conj.cpp
index 3a8e38e6799..1ffe557fbf1 100644
--- a/dpnp/libtensor/extensions/vm/conj.cpp
+++ b/dpnp/libtensor/extensions/vm/conj.cpp
@@ -92,6 +92,8 @@ static sycl::event conj_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::conj(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/copysign.cpp b/dpnp/libtensor/extensions/vm/copysign.cpp
index ec5e7d4e23a..dd77b22a2cf 100644
--- a/dpnp/libtensor/extensions/vm/copysign.cpp
+++ b/dpnp/libtensor/extensions/vm/copysign.cpp
@@ -102,6 +102,8 @@ static sycl::event copysign_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::copysign(
         exec_q,
         n, // number of elements to be calculated
diff --git a/dpnp/libtensor/extensions/vm/cos.cpp b/dpnp/libtensor/extensions/vm/cos.cpp
index 98653c1be01..308b3b572ab 100644
--- a/dpnp/libtensor/extensions/vm/cos.cpp
+++ b/dpnp/libtensor/extensions/vm/cos.cpp
@@ -94,6 +94,8 @@ static sycl::event cos_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::cos(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/cosh.cpp b/dpnp/libtensor/extensions/vm/cosh.cpp
index d655b2dcc68..3e02a9b9179 100644
--- a/dpnp/libtensor/extensions/vm/cosh.cpp
+++ b/dpnp/libtensor/extensions/vm/cosh.cpp
@@ -94,6 +94,8 @@ static sycl::event cosh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::cosh(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/div.cpp b/dpnp/libtensor/extensions/vm/div.cpp
index 68de80960d6..bb9d0f339ff 100644
--- a/dpnp/libtensor/extensions/vm/div.cpp
+++ b/dpnp/libtensor/extensions/vm/div.cpp
@@ -113,6 +113,8 @@ static sycl::event div_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::div(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/erf_funcs.cpp b/dpnp/libtensor/extensions/vm/erf_funcs.cpp
index 38796ae15ff..eeabc33a34b 100644
--- a/dpnp/libtensor/extensions/vm/erf_funcs.cpp
+++ b/dpnp/libtensor/extensions/vm/erf_funcs.cpp
@@ -106,6 +106,8 @@ using ew_cmn_ns::unary_contig_impl_fn_ptr_t;
         using resTy = typename OutputType<T>::value_type;                      \
         resTy *y = reinterpret_cast<resTy *>(out_y);                           \
                                                                                \
+        vm_gil_release lock{};                                                 \
+                                                                               \
         return mkl_vm::__name__(                                               \
             exec_q, n, /* number of elements to be calculated*/                \
             a,         /* pointer `a` containing input vector of size n*/      \
diff --git a/dpnp/libtensor/extensions/vm/exp.cpp b/dpnp/libtensor/extensions/vm/exp.cpp
index bbcefdc460b..d9ec198808d 100644
--- a/dpnp/libtensor/extensions/vm/exp.cpp
+++ b/dpnp/libtensor/extensions/vm/exp.cpp
@@ -94,6 +94,8 @@ static sycl::event exp_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::exp(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/exp2.cpp b/dpnp/libtensor/extensions/vm/exp2.cpp
index d25e9b89bcd..933b29903ff 100644
--- a/dpnp/libtensor/extensions/vm/exp2.cpp
+++ b/dpnp/libtensor/extensions/vm/exp2.cpp
@@ -91,6 +91,8 @@ static sycl::event exp2_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::exp2(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/expm1.cpp b/dpnp/libtensor/extensions/vm/expm1.cpp
index 6d62800bc73..ecf5332c443 100644
--- a/dpnp/libtensor/extensions/vm/expm1.cpp
+++ b/dpnp/libtensor/extensions/vm/expm1.cpp
@@ -91,6 +91,8 @@ static sycl::event expm1_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::expm1(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/floor.cpp b/dpnp/libtensor/extensions/vm/floor.cpp
index 91350e69588..f5e1f224ace 100644
--- a/dpnp/libtensor/extensions/vm/floor.cpp
+++ b/dpnp/libtensor/extensions/vm/floor.cpp
@@ -91,6 +91,8 @@ static sycl::event floor_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::floor(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/fmax.cpp b/dpnp/libtensor/extensions/vm/fmax.cpp
index 71000e1392b..e2009107f7e 100644
--- a/dpnp/libtensor/extensions/vm/fmax.cpp
+++ b/dpnp/libtensor/extensions/vm/fmax.cpp
@@ -102,6 +102,8 @@ static sycl::event fmax_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::fmax(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/fmin.cpp b/dpnp/libtensor/extensions/vm/fmin.cpp
index 1a8f331125d..050b0d44383 100644
--- a/dpnp/libtensor/extensions/vm/fmin.cpp
+++ b/dpnp/libtensor/extensions/vm/fmin.cpp
@@ -102,6 +102,8 @@ static sycl::event fmin_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::fmin(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/fmod.cpp b/dpnp/libtensor/extensions/vm/fmod.cpp
index 56481afa5ef..43b0710dbd4 100644
--- a/dpnp/libtensor/extensions/vm/fmod.cpp
+++ b/dpnp/libtensor/extensions/vm/fmod.cpp
@@ -102,6 +102,8 @@ static sycl::event fmod_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::fmod(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/hypot.cpp b/dpnp/libtensor/extensions/vm/hypot.cpp
index 88195916d32..2a63fbeaabd 100644
--- a/dpnp/libtensor/extensions/vm/hypot.cpp
+++ b/dpnp/libtensor/extensions/vm/hypot.cpp
@@ -102,6 +102,8 @@ static sycl::event hypot_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::hypot(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/i0.cpp b/dpnp/libtensor/extensions/vm/i0.cpp
index 68d3d866064..8a5d9e1acec 100644
--- a/dpnp/libtensor/extensions/vm/i0.cpp
+++ b/dpnp/libtensor/extensions/vm/i0.cpp
@@ -91,6 +91,8 @@ static sycl::event i0_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::i0(exec_q,
                       n, // number of elements to be calculated
                       a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/inv.cpp b/dpnp/libtensor/extensions/vm/inv.cpp
index 5575be3cb42..3c2a0c61395 100644
--- a/dpnp/libtensor/extensions/vm/inv.cpp
+++ b/dpnp/libtensor/extensions/vm/inv.cpp
@@ -91,6 +91,8 @@ static sycl::event inv_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::inv(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/ln.cpp b/dpnp/libtensor/extensions/vm/ln.cpp
index 8ce83e2c476..f87ea66291f 100644
--- a/dpnp/libtensor/extensions/vm/ln.cpp
+++ b/dpnp/libtensor/extensions/vm/ln.cpp
@@ -94,6 +94,8 @@ static sycl::event ln_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::ln(exec_q,
                       n, // number of elements to be calculated
                       a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/log10.cpp b/dpnp/libtensor/extensions/vm/log10.cpp
index b175c336a8c..42e2d945f39 100644
--- a/dpnp/libtensor/extensions/vm/log10.cpp
+++ b/dpnp/libtensor/extensions/vm/log10.cpp
@@ -94,6 +94,8 @@ static sycl::event log10_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::log10(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/log1p.cpp b/dpnp/libtensor/extensions/vm/log1p.cpp
index 336ddb8cbc8..1f46d9eab2c 100644
--- a/dpnp/libtensor/extensions/vm/log1p.cpp
+++ b/dpnp/libtensor/extensions/vm/log1p.cpp
@@ -91,6 +91,8 @@ static sycl::event log1p_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::log1p(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/log2.cpp b/dpnp/libtensor/extensions/vm/log2.cpp
index c17f4e19ab1..c269f5b50fb 100644
--- a/dpnp/libtensor/extensions/vm/log2.cpp
+++ b/dpnp/libtensor/extensions/vm/log2.cpp
@@ -91,6 +91,8 @@ static sycl::event log2_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::log2(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/modf.cpp b/dpnp/libtensor/extensions/vm/modf.cpp
index 68b6f945910..387e2acdbbc 100644
--- a/dpnp/libtensor/extensions/vm/modf.cpp
+++ b/dpnp/libtensor/extensions/vm/modf.cpp
@@ -100,6 +100,8 @@ static sycl::event modf_contig_impl(sycl::queue &exec_q,
     fractT *y = reinterpret_cast<fractT *>(out_y);
     intT *z = reinterpret_cast<intT *>(out_z);
 
+    vm_gil_release lock{};
+
     return mkl_vm::modf(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/mul.cpp b/dpnp/libtensor/extensions/vm/mul.cpp
index a25f6794da4..224c9f1f070 100644
--- a/dpnp/libtensor/extensions/vm/mul.cpp
+++ b/dpnp/libtensor/extensions/vm/mul.cpp
@@ -113,6 +113,8 @@ static sycl::event mul_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::mul(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/nextafter.cpp b/dpnp/libtensor/extensions/vm/nextafter.cpp
index 32e4df5895d..cc64e33690b 100644
--- a/dpnp/libtensor/extensions/vm/nextafter.cpp
+++ b/dpnp/libtensor/extensions/vm/nextafter.cpp
@@ -103,6 +103,8 @@ static sycl::event
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::nextafter(
         exec_q,
         n, // number of elements to be calculated
diff --git a/dpnp/libtensor/extensions/vm/pow.cpp b/dpnp/libtensor/extensions/vm/pow.cpp
index 7aed34fe5fd..005458c03ce 100644
--- a/dpnp/libtensor/extensions/vm/pow.cpp
+++ b/dpnp/libtensor/extensions/vm/pow.cpp
@@ -113,6 +113,8 @@ static sycl::event pow_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::pow(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/rint.cpp b/dpnp/libtensor/extensions/vm/rint.cpp
index 5723065c53a..8410b3e8017 100644
--- a/dpnp/libtensor/extensions/vm/rint.cpp
+++ b/dpnp/libtensor/extensions/vm/rint.cpp
@@ -91,6 +91,8 @@ static sycl::event rint_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::rint(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/sin.cpp b/dpnp/libtensor/extensions/vm/sin.cpp
index 9dec89ad9d6..1f7d909e675 100644
--- a/dpnp/libtensor/extensions/vm/sin.cpp
+++ b/dpnp/libtensor/extensions/vm/sin.cpp
@@ -94,6 +94,8 @@ static sycl::event sin_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::sin(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/sinh.cpp b/dpnp/libtensor/extensions/vm/sinh.cpp
index ab3d77b9add..bc082b573eb 100644
--- a/dpnp/libtensor/extensions/vm/sinh.cpp
+++ b/dpnp/libtensor/extensions/vm/sinh.cpp
@@ -94,6 +94,8 @@ static sycl::event sinh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::sinh(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/sqr.cpp b/dpnp/libtensor/extensions/vm/sqr.cpp
index 80e30762a0b..b711a4d0795 100644
--- a/dpnp/libtensor/extensions/vm/sqr.cpp
+++ b/dpnp/libtensor/extensions/vm/sqr.cpp
@@ -91,6 +91,8 @@ static sycl::event sqr_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::sqr(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/sqrt.cpp b/dpnp/libtensor/extensions/vm/sqrt.cpp
index 2a1bb7f7a99..e388788fe3b 100644
--- a/dpnp/libtensor/extensions/vm/sqrt.cpp
+++ b/dpnp/libtensor/extensions/vm/sqrt.cpp
@@ -94,6 +94,8 @@ static sycl::event sqrt_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::sqrt(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/sub.cpp b/dpnp/libtensor/extensions/vm/sub.cpp
index 022ec5bfefa..63ee87a7183 100644
--- a/dpnp/libtensor/extensions/vm/sub.cpp
+++ b/dpnp/libtensor/extensions/vm/sub.cpp
@@ -113,6 +113,8 @@ static sycl::event sub_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T1, T2>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::sub(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing 1st input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/tan.cpp b/dpnp/libtensor/extensions/vm/tan.cpp
index a4ff4842c4a..df15e73d9e0 100644
--- a/dpnp/libtensor/extensions/vm/tan.cpp
+++ b/dpnp/libtensor/extensions/vm/tan.cpp
@@ -94,6 +94,8 @@ static sycl::event tan_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::tan(exec_q,
                        n, // number of elements to be calculated
                        a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/tanh.cpp b/dpnp/libtensor/extensions/vm/tanh.cpp
index 9734a2fc718..9cd25a51e53 100644
--- a/dpnp/libtensor/extensions/vm/tanh.cpp
+++ b/dpnp/libtensor/extensions/vm/tanh.cpp
@@ -94,6 +94,8 @@ static sycl::event tanh_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::tanh(exec_q,
                         n, // number of elements to be calculated
                         a, // pointer `a` containing input vector of size n
diff --git a/dpnp/libtensor/extensions/vm/trunc.cpp b/dpnp/libtensor/extensions/vm/trunc.cpp
index 57a749916af..2e572e90d9a 100644
--- a/dpnp/libtensor/extensions/vm/trunc.cpp
+++ b/dpnp/libtensor/extensions/vm/trunc.cpp
@@ -91,6 +91,8 @@ static sycl::event trunc_contig_impl(sycl::queue &exec_q,
     using resTy = typename OutputType<T>::value_type;
     resTy *y = reinterpret_cast<resTy *>(out_y);
 
+    vm_gil_release lock{};
+
     return mkl_vm::trunc(exec_q,
                          n, // number of elements to be calculated
                          a, // pointer `a` containing input vector of size n
diff --git a/dpnp/tests/test_vm_gil.py b/dpnp/tests/test_vm_gil.py
new file mode 100644
index 00000000000..8dccdbbd26b
--- /dev/null
+++ b/dpnp/tests/test_vm_gil.py
@@ -0,0 +1,68 @@
+"""oneMKL VM calls in the VM extension must release the GIL.
+
+On a CPU device oneMKL VM waits on the dependent events and computes on the
+calling thread. On an in-order queue those events are ordered after the host
+tasks dpnp submits to keep arguments alive, and these host tasks acquire the
+GIL, so a VM call holding the GIL never returns.
+"""
+
+import subprocess
+import sys
+
+import dpctl
+import pytest
+
+import dpnp.libtensor._vm_impl as vmi
+
+# Prints "novm" if oneMKL VM is not used for the subtraction, or "ok".
+_VM_AFTER_HOST_TASKS = """
+import dpctl, dpnp
+import dpnp.libtensor._vm_impl as vmi
+q = dpctl.SyclQueue("cpu", property="in_order")
+xp = dpnp.linspace(10.0, 0.0, 5, sycl_queue=q)
+a, b = xp[-1:], xp[:1]
+out = dpnp.empty_like(a)
+if not vmi._mkl_sub_to_call(q, a.get_array(), b.get_array(), out.get_array()):
+    print("novm")
+    raise SystemExit
+for _ in range(100):
+    xs = xp[dpnp.argsort(xp)]  # submits host tasks keeping arguments alive
+    xs[-1:] - xs[:1]
+q.wait()
+print("ok")
+"""
+
+_TIMEOUT = 120  # seconds, a run takes about one
+
+
+def _cpu_device_with_fp64():
+    for dev in dpctl.get_devices(device_type="cpu"):
+        if dev.has_aspect_fp64:
+            return dev
+    return None
+
+
+@pytest.mark.skipif(not vmi._is_available(), reason="requires oneMKL VM")
+@pytest.mark.skipif(
+    _cpu_device_with_fp64() is None,
+    reason="requires a CPU device with fp64 support",
+)
+def test_in_order_queue_after_host_tasks():
+    try:
+        res = subprocess.run(
+            [sys.executable, "-c", _VM_AFTER_HOST_TASKS],
+            capture_output=True,
+            text=True,
+            timeout=_TIMEOUT,
+        )
+    except subprocess.TimeoutExpired:
+        pytest.fail(
+            f"no progress in {_TIMEOUT}s: a oneMKL VM call holds the GIL "
+            "while waiting on events behind host tasks. It needs "
+            "py::gil_scoped_release."
+        )
+    assert res.returncode == 0, res.stderr
+    out = res.stdout.strip()
+    if out == "novm":
+        pytest.skip("oneMKL VM is not used for the subtraction")
+    assert out == "ok", out