diff --git a/CHANGELOG.md b/CHANGELOG.md index cfdf45ce920..d9fc960c2e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -82,6 +82,8 @@ This release is compatible with NumPy 2.5. * Fixed `dpnp.tensor.top_k` aborting for `k=0` by returning empty result arrays without launching a zero-sized kernel [#3022](https://github.com/IntelPython/dpnp/pull/3022) * Fixed comparison functions (`dpnp.equal`, `dpnp.not_equal`, `dpnp.less`, `dpnp.less_equal`, `dpnp.greater`, `dpnp.greater_equal`) and `dpnp.divide` raising `OverflowError` when comparing an integer array against a Python integer scalar outside the array dtype's range [#3017](https://github.com/IntelPython/dpnp/pull/3017) * Fixed a crash in boolean-mask advanced indexing (`dpnp.ndarray` get/set item) when the selection is empty (e.g. a scalar `False` index that injects a length-0 axis) [#3019](https://github.com/IntelPython/dpnp/pull/3019) +* Released the GIL before the remaining blocking OneMKL BLAS and LAPACK calls to prevent host tasks contention, completing the work started in [#2850](https://github.com/IntelPython/dpnp/pull/2850) [#3027](https://github.com/IntelPython/dpnp/pull/3027) + ### Security diff --git a/dpnp/backend/extensions/blas/dot.hpp b/dpnp/backend/extensions/blas/dot.hpp index 300a5348c53..755dd046d2e 100644 --- a/dpnp/backend/extensions/blas/dot.hpp +++ b/dpnp/backend/extensions/blas/dot.hpp @@ -30,12 +30,15 @@ #include +#include + #include "dot_common.hpp" namespace dpnp::extensions::blas { namespace mkl_blas = oneapi::mkl::blas; namespace type_utils = dpnp::tensor::type_utils; +namespace py = pybind11; template static sycl::event dot_impl(sycl::queue &exec_q, @@ -58,6 +61,10 @@ static sycl::event dot_impl(sycl::queue &exec_q, sycl::event dot_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + dot_event = mkl_blas::column_major::dot(exec_q, n, // size of the input vectors x, // Pointer to vector x. diff --git a/dpnp/backend/extensions/blas/dotc.hpp b/dpnp/backend/extensions/blas/dotc.hpp index 535f79156fb..051c85b7b16 100644 --- a/dpnp/backend/extensions/blas/dotc.hpp +++ b/dpnp/backend/extensions/blas/dotc.hpp @@ -30,12 +30,15 @@ #include +#include + #include "dot_common.hpp" namespace dpnp::extensions::blas { namespace mkl_blas = oneapi::mkl::blas; namespace type_utils = dpnp::tensor::type_utils; +namespace py = pybind11; template static sycl::event dotc_impl(sycl::queue &exec_q, @@ -58,6 +61,10 @@ static sycl::event dotc_impl(sycl::queue &exec_q, sycl::event dotc_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + dotc_event = mkl_blas::column_major::dotc(exec_q, n, // size of the input vectors diff --git a/dpnp/backend/extensions/blas/dotu.hpp b/dpnp/backend/extensions/blas/dotu.hpp index aafe9721339..f517d89bf21 100644 --- a/dpnp/backend/extensions/blas/dotu.hpp +++ b/dpnp/backend/extensions/blas/dotu.hpp @@ -30,12 +30,15 @@ #include +#include + #include "dot_common.hpp" namespace dpnp::extensions::blas { namespace mkl_blas = oneapi::mkl::blas; namespace type_utils = dpnp::tensor::type_utils; +namespace py = pybind11; template static sycl::event dotu_impl(sycl::queue &exec_q, @@ -58,6 +61,10 @@ static sycl::event dotu_impl(sycl::queue &exec_q, sycl::event dotu_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + dotu_event = mkl_blas::column_major::dotu(exec_q, n, // size of the input vectors diff --git a/dpnp/backend/extensions/blas/gemm.cpp b/dpnp/backend/extensions/blas/gemm.cpp index bf91f5cfb0b..b7b39c3bc66 100644 --- a/dpnp/backend/extensions/blas/gemm.cpp +++ b/dpnp/backend/extensions/blas/gemm.cpp @@ -113,6 +113,10 @@ static sycl::event gemm_impl(sycl::queue &exec_q, c, ldc, deps); } }; + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + gemm_event = gemm_func( exec_q, transA, // Defines the transpose operation for matrix A: diff --git a/dpnp/backend/extensions/blas/gemm_batch.cpp b/dpnp/backend/extensions/blas/gemm_batch.cpp index f1bc9f30b5d..d4e09162a6e 100644 --- a/dpnp/backend/extensions/blas/gemm_batch.cpp +++ b/dpnp/backend/extensions/blas/gemm_batch.cpp @@ -125,6 +125,10 @@ static sycl::event gemm_batch_impl(sycl::queue &exec_q, strideb, beta, c, ldc, stridec, batch_size, deps); } }; + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + gemm_batch_event = gemm_batch_func( exec_q, transA, // Defines the transpose operation for matrix A: diff --git a/dpnp/backend/extensions/lapack/geqrf.cpp b/dpnp/backend/extensions/lapack/geqrf.cpp index 083c408ea3d..03e9c5060fd 100644 --- a/dpnp/backend/extensions/lapack/geqrf.cpp +++ b/dpnp/backend/extensions/lapack/geqrf.cpp @@ -86,6 +86,10 @@ static sycl::event geqrf_impl(sycl::queue &exec_q, sycl::event geqrf_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); geqrf_event = mkl_lapack::geqrf( diff --git a/dpnp/backend/extensions/lapack/geqrf_batch.cpp b/dpnp/backend/extensions/lapack/geqrf_batch.cpp index 6e62d4dde33..739e2a82b54 100644 --- a/dpnp/backend/extensions/lapack/geqrf_batch.cpp +++ b/dpnp/backend/extensions/lapack/geqrf_batch.cpp @@ -94,6 +94,10 @@ static sycl::event geqrf_batch_impl(sycl::queue &exec_q, sycl::event geqrf_batch_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); geqrf_batch_event = mkl_lapack::geqrf_batch( diff --git a/dpnp/backend/extensions/lapack/gesv.cpp b/dpnp/backend/extensions/lapack/gesv.cpp index 072efd169c8..3bab2e5b35a 100644 --- a/dpnp/backend/extensions/lapack/gesv.cpp +++ b/dpnp/backend/extensions/lapack/gesv.cpp @@ -112,6 +112,10 @@ static sycl::event gesv_impl(sycl::queue &exec_q, #if defined(USE_ONEMATH) sycl::event getrf_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + getrf_event = mkl_lapack::getrf( exec_q, n, // The order of the square matrix A (0 ≤ n). @@ -169,6 +173,10 @@ static sycl::event gesv_impl(sycl::queue &exec_q, } #else try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + comp_event = mkl_lapack::gesv( exec_q, n, // The order of the square matrix A diff --git a/dpnp/backend/extensions/lapack/gesvd.cpp b/dpnp/backend/extensions/lapack/gesvd.cpp index ebc6000671b..02b38897714 100644 --- a/dpnp/backend/extensions/lapack/gesvd.cpp +++ b/dpnp/backend/extensions/lapack/gesvd.cpp @@ -99,6 +99,10 @@ static sycl::event gesvd_impl(sycl::queue &exec_q, sycl::event gesvd_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + gesvd_event = mkl_lapack::gesvd( exec_q, jobu, // Character specifying how to compute the matrix U: diff --git a/dpnp/backend/extensions/lapack/getrf.cpp b/dpnp/backend/extensions/lapack/getrf.cpp index 42021639d63..c31213b3862 100644 --- a/dpnp/backend/extensions/lapack/getrf.cpp +++ b/dpnp/backend/extensions/lapack/getrf.cpp @@ -87,6 +87,10 @@ static sycl::event getrf_impl(sycl::queue &exec_q, sycl::event getrf_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); getrf_event = mkl_lapack::getrf( diff --git a/dpnp/backend/extensions/lapack/getrf_batch.cpp b/dpnp/backend/extensions/lapack/getrf_batch.cpp index b4e30dae57d..91779455737 100644 --- a/dpnp/backend/extensions/lapack/getrf_batch.cpp +++ b/dpnp/backend/extensions/lapack/getrf_batch.cpp @@ -96,6 +96,10 @@ static sycl::event getrf_batch_impl(sycl::queue &exec_q, sycl::event getrf_batch_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); getrf_batch_event = mkl_lapack::getrf_batch( diff --git a/dpnp/backend/extensions/lapack/heevd.cpp b/dpnp/backend/extensions/lapack/heevd.cpp index 002a7ad2ef8..86ba9776588 100644 --- a/dpnp/backend/extensions/lapack/heevd.cpp +++ b/dpnp/backend/extensions/lapack/heevd.cpp @@ -42,6 +42,7 @@ namespace dpnp::extensions::lapack { +namespace py = pybind11; namespace mkl_lapack = oneapi::mkl::lapack; namespace type_utils = dpnp::tensor::type_utils; @@ -73,6 +74,10 @@ static sycl::event heevd_impl(sycl::queue &exec_q, sycl::event heevd_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + heevd_event = mkl_lapack::heevd( exec_q, jobz, // 'jobz == job::vec' means eigenvalues and eigenvectors are diff --git a/dpnp/backend/extensions/lapack/orgqr_batch.cpp b/dpnp/backend/extensions/lapack/orgqr_batch.cpp index 840ef82a09b..0d996bef2b8 100644 --- a/dpnp/backend/extensions/lapack/orgqr_batch.cpp +++ b/dpnp/backend/extensions/lapack/orgqr_batch.cpp @@ -96,6 +96,10 @@ static sycl::event orgqr_batch_impl(sycl::queue &exec_q, sycl::event orgqr_batch_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); orgqr_batch_event = mkl_lapack::orgqr_batch( diff --git a/dpnp/backend/extensions/lapack/potrf.cpp b/dpnp/backend/extensions/lapack/potrf.cpp index d7edc2d244f..a14ce223458 100644 --- a/dpnp/backend/extensions/lapack/potrf.cpp +++ b/dpnp/backend/extensions/lapack/potrf.cpp @@ -83,6 +83,10 @@ static sycl::event potrf_impl(sycl::queue &exec_q, sycl::event potrf_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); potrf_event = mkl_lapack::potrf( diff --git a/dpnp/backend/extensions/lapack/potrf_batch.cpp b/dpnp/backend/extensions/lapack/potrf_batch.cpp index 58bdad46293..8193b32d4f0 100644 --- a/dpnp/backend/extensions/lapack/potrf_batch.cpp +++ b/dpnp/backend/extensions/lapack/potrf_batch.cpp @@ -90,6 +90,10 @@ static sycl::event potrf_batch_impl(sycl::queue &exec_q, sycl::event potrf_batch_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); potrf_batch_event = mkl_lapack::potrf_batch( diff --git a/dpnp/backend/extensions/lapack/syevd.cpp b/dpnp/backend/extensions/lapack/syevd.cpp index 7280e415f27..f2be6214f60 100644 --- a/dpnp/backend/extensions/lapack/syevd.cpp +++ b/dpnp/backend/extensions/lapack/syevd.cpp @@ -42,6 +42,7 @@ namespace dpnp::extensions::lapack { +namespace py = pybind11; namespace mkl_lapack = oneapi::mkl::lapack; namespace type_utils = dpnp::tensor::type_utils; @@ -73,6 +74,10 @@ static sycl::event syevd_impl(sycl::queue &exec_q, sycl::event syevd_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + syevd_event = mkl_lapack::syevd( exec_q, jobz, // 'jobz == job::vec' means eigenvalues and eigenvectors are diff --git a/dpnp/backend/extensions/lapack/ungqr.cpp b/dpnp/backend/extensions/lapack/ungqr.cpp index 7bea568f048..42b14059176 100644 --- a/dpnp/backend/extensions/lapack/ungqr.cpp +++ b/dpnp/backend/extensions/lapack/ungqr.cpp @@ -87,6 +87,10 @@ static sycl::event ungqr_impl(sycl::queue &exec_q, sycl::event ungqr_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); ungqr_event = mkl_lapack::ungqr( diff --git a/dpnp/backend/extensions/lapack/ungqr_batch.cpp b/dpnp/backend/extensions/lapack/ungqr_batch.cpp index 6452a30e2a8..1ac7220064b 100644 --- a/dpnp/backend/extensions/lapack/ungqr_batch.cpp +++ b/dpnp/backend/extensions/lapack/ungqr_batch.cpp @@ -96,6 +96,10 @@ static sycl::event ungqr_batch_impl(sycl::queue &exec_q, sycl::event ungqr_batch_event; try { + // Release GIL to avoid serialization of host task submissions + // to the same queue in OneMKL + py::gil_scoped_release lock{}; + scratchpad = sycl::malloc_device(scratchpad_size, exec_q); ungqr_batch_event = mkl_lapack::ungqr_batch( diff --git a/dpnp/tests/test_blas_lapack_gil.py b/dpnp/tests/test_blas_lapack_gil.py new file mode 100644 index 00000000000..9aec868d079 --- /dev/null +++ b/dpnp/tests/test_blas_lapack_gil.py @@ -0,0 +1,169 @@ +"""Blocking oneMKL calls in the BLAS/LAPACK extensions must release the GIL. + +Progress of a competing thread is compared against ``SyclQueue.wait()``, which +is ``nogil`` and therefore the best rate achievable on the machine. +""" + +import sys +import threading +import time + +import dpctl +import pytest + +import dpnp + +from .helper import has_support_aspect64, is_gpu_device + +# Smaller sizes stop discriminating: the calls either stay asynchronous or +# block too briefly for a stable measurement. +_SIZE = 2048 + +_MIN_RATIO = 0.10 + +_BACKLOG = 2 # queued matmuls, so the measured call has something to wait on + +_TRIALS = 5 # samples averaged per measurement + + +class _Ticker: + """Counts how often a competing Python thread gets scheduled.""" + + def __enter__(self): + self.ticks = 0 + self._stop = False + self._thread = threading.Thread(target=self._run, daemon=True) + self._thread.start() + return self + + def _run(self): + while not self._stop: + self.ticks += 1 + time.sleep(0) + + def __exit__(self, *exc): + self._stop = True + self._thread.join(timeout=5) + return False + + +@pytest.mark.skipif(not has_support_aspect64(), reason="requires fp64 support") +class TestBlockingCallsReleaseGil: + @pytest.fixture(autouse=True) + def _switch_interval(self): + # Stop CPython handing the GIL over on its own timer. + previous = sys.getswitchinterval() + sys.setswitchinterval(0.001) + yield + sys.setswitchinterval(previous) + + @pytest.fixture + def mats(self): + spd = dpnp.eye(_SIZE, dtype="f8") * float(_SIZE) + return spd, dpnp.eye(_SIZE, dtype="f8") + 0.1 + + def _assert_releases_gil(self, name, spd, fn): + queue = spd.sycl_queue + + def rate(ticker, func): + total_ticks = 0 + total_ms = 0.0 + for _ in range(_TRIALS): + for _ in range(_BACKLOG): + dpnp.matmul(spd, spd) + ticker.ticks = 0 + start = time.perf_counter() + func() + total_ms += 1000 * (time.perf_counter() - start) + total_ticks += ticker.ticks + queue.wait() + return total_ticks / max(total_ms, 1e-3) + + fn() # warm up JIT + queue.wait() + + with _Ticker() as ticker: + measured = rate(ticker, fn) + reference = rate(ticker, queue.wait) + + assert reference > 0, "reference measurement produced no ticks" + ratio = measured / reference + assert ratio >= _MIN_RATIO, ( + f"{name} holds the GIL while blocking: {measured:.2f} ticks/ms vs " + f"{reference:.2f} for nogil queue.wait() (ratio {ratio:.3f}, need " + f">= {_MIN_RATIO}). The oneMKL call needs py::gil_scoped_release." + ) + + @pytest.mark.slow + def test_potrf(self, mats): + spd, _ = mats + self._assert_releases_gil( + "potrf", spd, lambda: dpnp.linalg.cholesky(spd) + ) + + @pytest.mark.slow + def test_getrf(self, mats): + spd, gen = mats + self._assert_releases_gil("getrf", spd, lambda: dpnp.linalg.det(gen)) + + @pytest.mark.slow + def test_syevd(self, mats): + spd, _ = mats + self._assert_releases_gil("syevd", spd, lambda: dpnp.linalg.eigh(spd)) + + @pytest.mark.slow + def test_gesv(self, mats): + spd, gen = mats + rhs = dpnp.ones(_SIZE, dtype="f8") + self._assert_releases_gil( + "gesv", spd, lambda: dpnp.linalg.solve(gen, rhs) + ) + + +# GPU only: on a CPU device oneMKL already saturates every core, so two threads +# contend for the same hardware and the ratio stays ~1.0 either way. +@pytest.mark.slow +@pytest.mark.skipif(not is_gpu_device(), reason="requires a GPU device") +@pytest.mark.skipif(not has_support_aspect64(), reason="requires fp64 support") +def test_multithreaded_linalg_overlaps(): + """Two threads on two queues must run oneMKL work concurrently.""" + max_ratio = 0.85 + + device = dpctl.select_default_device() + queues = [dpctl.SyclQueue(device, property="in_order") for _ in range(2)] + mats = [ + dpnp.eye(_SIZE, dtype="f8", sycl_queue=q) * float(_SIZE) for q in queues + ] + + def work(idx): + for _ in range(2): + dpnp.linalg.eigh(mats[idx]) + queues[idx].wait() + + for idx in range(2): # warm up JIT + work(idx) + + previous = sys.getswitchinterval() + sys.setswitchinterval(0.001) + try: + start = time.perf_counter() + work(0) + work(1) + serial = time.perf_counter() - start + + start = time.perf_counter() + threads = [threading.Thread(target=work, args=(i,)) for i in range(2)] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + parallel = time.perf_counter() - start + finally: + sys.setswitchinterval(previous) + + ratio = parallel / serial + assert ratio <= max_ratio, ( + f"two threads running eigh did not overlap: {parallel * 1000:.1f} ms " + f"parallel vs {serial * 1000:.1f} ms serial (ratio {ratio:.3f}, need " + f"<= {max_ratio}). The oneMKL call is holding the GIL." + )