From 138d8e5f33d053f2dfb40a158a0d0008b42f0853 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 27 Oct 2026 12:47:37 +0100 Subject: [PATCH 1/4] Add minmax implementation --- CHANGELOG.md | 1 + doc/reference/math.rst | 2 + dpnp/__init__.py | 2 + dpnp/dpnp_iface_statistics.py | 123 +++ dpnp/tensor/CMakeLists.txt | 1 + dpnp/tensor/_reduction.py | 133 ++++ .../libtensor/include/kernels/reductions.hpp | 707 ++++++++++++++++++ .../libtensor/source/reductions/minmax.cpp | 469 ++++++++++++ .../libtensor/source/reductions/minmax.hpp | 46 ++ .../source/reductions/reduction_common.cpp | 2 + dpnp/tests/test_statistics.py | 232 ++++++ dpnp/tests/test_sycl_queue.py | 1 + dpnp/tests/test_usm_type.py | 1 + 13 files changed, 1720 insertions(+) create mode 100644 dpnp/tensor/libtensor/source/reductions/minmax.cpp create mode 100644 dpnp/tensor/libtensor/source/reductions/minmax.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index e1ddf0416861..520fea0102c6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added * Added support for boolean input in `dpnp.positive` and the `+` operator of `dpnp.ndarray` [#3128](https://github.com/IntelPython/dpnp/pull/3128) +* Added implementation of `dpnp.minmax`, computing the minimum and the maximum of an array in a single pass by a new fused reduction kernel [#XXXX](https://github.com/IntelPython/dpnp/pull/XXXX) ### Changed diff --git a/doc/reference/math.rst b/doc/reference/math.rst index fd4af791767a..d26b5cc91b35 100644 --- a/doc/reference/math.rst +++ b/doc/reference/math.rst @@ -208,6 +208,8 @@ Extrema finding fmin nanmin + minmax + Miscellaneous ------------- diff --git a/dpnp/__init__.py b/dpnp/__init__.py index 68a683d8095d..0f726406841b 100644 --- a/dpnp/__init__.py +++ b/dpnp/__init__.py @@ -547,6 +547,7 @@ mean, median, min, + minmax, ptp, std, var, @@ -1016,6 +1017,7 @@ "mean", "median", "min", + "minmax", "nanmax", "nanmean", "nanmedian", diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 114072407ca8..8a5baec22f82 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -50,6 +50,7 @@ import dpnp.backend.extensions.statistics._statistics_impl as statistics_ext import dpnp.tensor as dpt import dpnp.tensor._tensor_elementwise_impl as ti +import dpnp.tensor._tensor_impl as dti from .dpnp_utils import get_usm_allocations from .dpnp_utils.dpnp_utils_common import ( @@ -59,6 +60,7 @@ from .dpnp_utils.dpnp_utils_reduction import dpnp_wrap_reduction_call from .dpnp_utils.dpnp_utils_statistics import dpnp_cov, dpnp_median from .tensor._numpy_helper import normalize_axis_index +from .tensor._reduction import _minmax as dpt_minmax def _count_reduce_items(arr, axis, where=True): @@ -1409,6 +1411,127 @@ def min(a, axis=None, out=None, keepdims=False, initial=None, where=True): ) +def minmax(a, axis=None, out=None, keepdims=False, initial=None, where=True): + """ + Return the minimum and maximum of an array or along an axis. + + This is equivalent to ``(dpnp.min(a, ...), dpnp.max(a, ...))`` but computes + both the minimum and the maximum in a single pass over `a`. + + For full documentation refer to :obj:`numpy.minmax`. + + Parameters + ---------- + a : {dpnp.ndarray, usm_ndarray} + Input array. + axis : {None, int or tuple of ints}, optional + Axis or axes along which to operate. By default, flattened input is + used. If this is a tuple of integers, the minimum and the maximum are + selected over multiple axes, instead of a single axis or all the axes + as before. + + Default: ``None``. + out : {None, tuple of None, dpnp.ndarray or usm_ndarray}, optional + A tuple ``(min, max)`` of two arrays in which to place the result. + Each array must be of the same shape and buffer length as the expected + output. A ``None`` entry means that a new array is allocated for that + result. + + Default: ``None``. + keepdims : {None, bool}, optional + If this is set to ``True``, the axes which are reduced are left in the + result as dimensions with size one. With this option, the result will + broadcast correctly against the input array. + + Default: ``False``. + + Returns + ------- + out : tuple of dpnp.ndarray + A tuple ``(min, max)`` holding the minimum and maximum of `a`. If + `axis` is ``None``, the results are zero-dimensional arrays. If `axis` + is an integer, the results are arrays of dimension ``a.ndim - 1``. If + `axis` is a tuple, the results are arrays of dimension + ``a.ndim - len(axis)``. + + Limitations + ----------- + Parameters `where`, and `initial` are only supported with their default + values. Otherwise ``NotImplementedError`` exception will be raised. + + See Also + -------- + :obj:`dpnp.min` : Return the minimum of an array. + :obj:`dpnp.max` : Return the maximum of an array. + :obj:`dpnp.ptp` : Range of values (maximum - minimum) along an axis. + :obj:`dpnp.nanmin` : The minimum value of an array along a given axis, + ignores NaNs. + :obj:`dpnp.nanmax` : The maximum value of an array along a given axis, + ignores NaNs. + + Notes + ----- + NaN values are propagated, that is if at least one item is NaN, the + corresponding output value will be NaN as well. To ignore NaN values + use :obj:`dpnp.nanmin` and :obj:`dpnp.nanmax`. + + Examples + -------- + >>> import dpnp as np + >>> a = np.arange(4).reshape((2, 2)) + >>> a + array([[0, 1], + [2, 3]]) + >>> np.minmax(a) # min and max of the flattened array + (array(0), array(3)) + + >>> np.minmax(a, axis=0) # along the first axis + (array([0, 1]), array([2, 3])) + >>> np.minmax(a, axis=1) # along the second axis + (array([0, 2]), array([1, 3])) + + >>> b = np.arange(5, dtype=float) + >>> b[2] = np.nan + >>> np.minmax(b) + (array(nan), array(nan)) + + """ + + dpnp.check_limitations(initial=initial, where=where) + usm_a = dpnp.get_usm_ndarray(a) + + if out is None: + out = (None, None) + elif not isinstance(out, tuple): + raise TypeError("'out' must be a tuple of arrays") + elif len(out) != 2: + raise ValueError("'out' tuple must have exactly two entries") + + if all(res is not None for res in out): + dpnp.check_supported_arrays_type(*out) + # pylint: disable=protected-access + if dti._array_overlap(*map(dpnp.get_usm_ndarray, out)): + raise ValueError("Output arrays cannot overlap") + + usm_out = [] + for res in out: + if res is not None: + dpnp.check_supported_arrays_type(res) + # the result is computed in the data type of `a` + if res.dtype != a.dtype: + res = dpnp.astype(res, a.dtype, copy=False) + res = dpnp.get_usm_ndarray(res) + usm_out.append(res) + + usm_res = dpt_minmax( + usm_a, axis=axis, keepdims=keepdims, out=tuple(usm_out) + ) + return tuple( + dpnp.get_result_array(x, res, casting="unsafe") + for x, res in zip(usm_res, out) + ) + + def ptp( a, /, diff --git a/dpnp/tensor/CMakeLists.txt b/dpnp/tensor/CMakeLists.txt index 4ec0bada6104..55956a1ac380 100644 --- a/dpnp/tensor/CMakeLists.txt +++ b/dpnp/tensor/CMakeLists.txt @@ -194,6 +194,7 @@ set(_reduction_sources ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/logsumexp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/max.cpp ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/min.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/minmax.cpp ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/prod.cpp ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/reduce_hypot.cpp ${CMAKE_CURRENT_SOURCE_DIR}/libtensor/source/reductions/sum.cpp diff --git a/dpnp/tensor/_reduction.py b/dpnp/tensor/_reduction.py index 782fc2b0b442..89c1c6bd5bd4 100644 --- a/dpnp/tensor/_reduction.py +++ b/dpnp/tensor/_reduction.py @@ -138,6 +138,139 @@ def _comparison_over_axis(x, axis, keepdims, out, _reduction_fn): return out +def _minmax(x, /, *, axis=None, keepdims=False, out=None): + """ + Calculates the minimum and the maximum of the input array ``x`` in a + single pass. + + Args: + x (usm_ndarray): + input array. + axis (Optional[int, Tuple[int, ...]]): + axis or axes along which the extrema must be computed. + If ``None``, they are computed over the entire array. + Default: ``None``. + keepdims (Optional[bool]): + if ``True``, the reduced axes (dimensions) are included in the + results as singleton dimensions. Default: ``False``. + out (Optional[Tuple[Optional[usm_ndarray], Optional[usm_ndarray]]]): + the arrays into which the minima and the maxima are written. + Each array must have the expected shape and the data type of + ``x``. A ``None`` entry means that a new array is returned. + Default: ``None``. + + Returns: + Tuple[usm_ndarray, usm_ndarray]: + the minima and the maxima of ``x``, with its data type. + """ + if not isinstance(x, dpt.usm_ndarray): + raise TypeError(f"Expected dpnp.tensor.usm_ndarray, got {type(x)}") + + nd = x.ndim + if axis is None: + axis = tuple(range(nd)) + perm = list(axis) + x_tmp = x + else: + if not isinstance(axis, (tuple, list)): + axis = (axis,) + axis = normalize_axis_tuple(axis, nd, "axis") + perm = [i for i in range(nd) if i not in axis] + list(axis) + x_tmp = dpt.permute_dims(x, perm) + red_nd = len(axis) + if any([x_tmp.shape[i] == 0 for i in range(-red_nd, 0)]): + raise ValueError("reduction cannot be performed over zero-size axes") + res_shape = x_tmp.shape[: nd - red_nd] + exec_q = x.sycl_queue + res_dt = x.dtype + res_usm_type = x.usm_type + + if keepdims: + final_res_shape = tuple( + x.shape[d] if d not in axis else 1 for d in range(nd) + ) + else: + final_res_shape = res_shape + + if out is None: + out = (None, None) + orig_out = list(out) + out = list(out) + for i, res in enumerate(out): + if res is None: + out[i] = dpt.empty( + res_shape, + dtype=res_dt, + usm_type=res_usm_type, + sycl_queue=exec_q, + ) + continue + + if not isinstance(res, dpt.usm_ndarray): + raise TypeError( + f"output array must be of usm_ndarray type, got {type(res)}" + ) + if not res.flags.writable: + raise ValueError("provided `out` array is read-only") + if not res.shape == final_res_shape: + raise ValueError( + "The shape of input and output arrays are inconsistent. " + f"Expected output shape is {final_res_shape}, got {res.shape}" + ) + if res_dt != res.dtype: + raise ValueError( + f"Output array of type {res_dt} is needed, got {res.dtype}" + ) + if dpt.get_execution_queue((exec_q, res.sycl_queue)) is None: + raise dpt.ExecutionPlacementError( + "Input and output allocation queues are not compatible" + ) + if keepdims: + res = dpt.squeeze(res, axis=axis) + orig_out[i] = res + if ti._array_overlap(x, res): + res = dpt.empty_like(res) + out[i] = res + + _manager = SequentialOrderManager[exec_q] + dep_evs = _manager.submitted_events + if red_nd == 0: + for res in out: + ht_e_cpy, cpy_e = ti._copy_usm_ndarray_into_usm_ndarray( + src=x_tmp, dst=res, sycl_queue=exec_q, depends=dep_evs + ) + _manager.add_event_pair(ht_e_cpy, cpy_e) + red_evs = _manager.submitted_events + else: + hev, red_ev = tri._minmax_over_axis( + src=x_tmp, + trailing_dims_to_reduce=red_nd, + dst_min=out[0], + dst_max=out[1], + sycl_queue=exec_q, + depends=dep_evs, + ) + _manager.add_event_pair(hev, red_ev) + red_evs = [red_ev] + + for i, res in enumerate(out): + if not (orig_out[i] is None or orig_out[i] is res): + ht_e_cpy2, cpy2_e = ti._copy_usm_ndarray_into_usm_ndarray( + src=res, dst=orig_out[i], sycl_queue=exec_q, depends=red_evs + ) + _manager.add_event_pair(ht_e_cpy2, cpy2_e) + out[i] = orig_out[i] + + if keepdims: + res_shape = res_shape + (1,) * red_nd + inv_perm = sorted(range(nd), key=lambda d: perm[d]) + out = [ + dpt.permute_dims(dpt.reshape(res, res_shape), inv_perm) + for res in out + ] + return tuple(out) + + def _reduction_over_axis( x, axis, diff --git a/dpnp/tensor/libtensor/include/kernels/reductions.hpp b/dpnp/tensor/libtensor/include/kernels/reductions.hpp index 04a988ceba95..aa9e4b2a30d3 100644 --- a/dpnp/tensor/libtensor/include/kernels/reductions.hpp +++ b/dpnp/tensor/libtensor/include/kernels/reductions.hpp @@ -3315,4 +3315,711 @@ sycl::event search_axis0_over_group_temps_contig_impl( } } +/* ============ Fused minimum and maximum reduction ============ */ + +template +struct SequentialMinMaxReduction +{ +private: + const argT *inp_ = nullptr; + argT *out_min_ = nullptr; + argT *out_max_ = nullptr; + argT min_identity_; + argT max_identity_; + InputOutputIterIndexerT inp_out_iter_indexer_; + InputRedIndexerT inp_reduced_dims_indexer_; + std::size_t reduction_max_gid_ = 0; + +public: + SequentialMinMaxReduction( + const argT *inp, + argT *res_min, + argT *res_max, + const argT &min_identity_val, + const argT &max_identity_val, + const InputOutputIterIndexerT &arg_res_iter_indexer, + const InputRedIndexerT &arg_reduced_dims_indexer, + std::size_t reduction_size) + : inp_(inp), out_min_(res_min), out_max_(res_max), + min_identity_(min_identity_val), max_identity_(max_identity_val), + inp_out_iter_indexer_(arg_res_iter_indexer), + inp_reduced_dims_indexer_(arg_reduced_dims_indexer), + reduction_max_gid_(reduction_size) + { + } + + void operator()(sycl::id<1> id) const + { + auto const &inp_out_iter_offsets_ = inp_out_iter_indexer_(id[0]); + const ssize_t &inp_iter_offset = + inp_out_iter_offsets_.get_first_offset(); + const ssize_t &out_min_iter_offset = + inp_out_iter_offsets_.get_second_offset(); + const ssize_t &out_max_iter_offset = + inp_out_iter_offsets_.get_third_offset(); + + const MinOpT min_op{}; + const MaxOpT max_op{}; + + argT red_min(min_identity_); + argT red_max(max_identity_); + for (std::size_t m = 0; m < reduction_max_gid_; ++m) { + const ssize_t inp_reduction_offset = inp_reduced_dims_indexer_(m); + const ssize_t inp_offset = inp_iter_offset + inp_reduction_offset; + + const argT val = normalize_bool(inp_[inp_offset]); + red_min = min_op(red_min, val); + red_max = max_op(red_max, val); + } + + out_min_[out_min_iter_offset] = red_min; + out_max_[out_max_iter_offset] = red_max; + } +}; + +template +using is_real_floating = std::bool_constant || + std::is_same_v>; + +/* @brief Reduces both extrema over a work-group with SYCL group algorithms */ +template +struct MinMaxGroupReducer +{ + template + void operator()(const GroupT &wg, argT &red_min, argT &red_max) const + { + if constexpr (std::is_same_v) { + red_min = sycl::all_of_group(wg, red_min); + red_max = sycl::any_of_group(wg, red_max); + } + else if constexpr (is_real_floating::value) { + // a NaN reaches both extrema of a work-item + if (sycl::any_of_group(wg, std::isnan(red_min))) { + // propagate the NaN of the first work-item having one + const std::uint32_t lid = wg.get_local_linear_id(); + const std::uint32_t wgs = wg.get_local_linear_range(); + const std::uint32_t nan_lid = + sycl::reduce_over_group(wg, std::isnan(red_min) ? lid : wgs, + sycl::minimum()); + red_min = sycl::group_broadcast(wg, red_min, nan_lid); + red_max = sycl::group_broadcast(wg, red_max, nan_lid); + } + else { + red_min = + sycl::reduce_over_group(wg, red_min, sycl::minimum()); + red_max = + sycl::reduce_over_group(wg, red_max, sycl::maximum()); + } + } + else { + red_min = sycl::reduce_over_group(wg, red_min, MinOpT()); + red_max = sycl::reduce_over_group(wg, red_max, MaxOpT()); + } + } +}; + +/* @brief Whether MinMaxGroupReducer supports the type */ +template +struct can_use_minmax_group_reducer +{ + static constexpr bool value = + std::is_same_v || + (is_real_floating::value && + can_use_reduce_over_group, argT>::value && + can_use_reduce_over_group, argT>::value) || + (can_use_reduce_over_group::value && + can_use_reduce_over_group::value); +}; + +/* @brief Reduces both extrema over a work-group with + * custom_reduce_over_group */ +template +struct CustomMinMaxGroupReducer +{ +private: + SlmT local_mem_; + +public: + CustomMinMaxGroupReducer(const SlmT &local_mem) : local_mem_(local_mem) {} + + template + void operator()(const GroupT &wg, argT &red_min, argT &red_max) const + { + red_min = + su_ns::custom_reduce_over_group(wg, local_mem_, red_min, MinOpT()); + // local memory is reused for the maximum + sycl::group_barrier(wg, sycl::memory_scope::work_group); + red_max = + su_ns::custom_reduce_over_group(wg, local_mem_, red_max, MaxOpT()); + } +}; + +/* + When First is true, the input is read once to update both extrema. + Otherwise, the partial minima and maxima of a previous pass are reduced. +*/ +template +struct MinMaxReductionOverGroupFunctor +{ +private: + const argT *inp_min_ = nullptr; + const argT *inp_max_ = nullptr; + argT *out_min_ = nullptr; + argT *out_max_ = nullptr; + argT min_identity_; + argT max_identity_; + InputOutputIterIndexerT inp_out_iter_indexer_; + InputRedIndexerT inp_reduced_dims_indexer_; + GroupReducerT group_reducer_; + std::size_t reduction_max_gid_ = 0; + std::size_t iter_gws_ = 1; + std::size_t reductions_per_wi = 16; + +public: + MinMaxReductionOverGroupFunctor( + const argT *inp_min, + const argT *inp_max, + argT *res_min, + argT *res_max, + const argT &min_identity_val, + const argT &max_identity_val, + const InputOutputIterIndexerT &arg_res_iter_indexer, + const InputRedIndexerT &arg_reduced_dims_indexer, + const GroupReducerT &group_reducer, + std::size_t reduction_size, + std::size_t iteration_size, + std::size_t reduction_size_per_wi) + : inp_min_(inp_min), inp_max_(inp_max), out_min_(res_min), + out_max_(res_max), min_identity_(min_identity_val), + max_identity_(max_identity_val), + inp_out_iter_indexer_(arg_res_iter_indexer), + inp_reduced_dims_indexer_(arg_reduced_dims_indexer), + group_reducer_(group_reducer), reduction_max_gid_(reduction_size), + iter_gws_(iteration_size), reductions_per_wi(reduction_size_per_wi) + { + } + + void operator()(sycl::nd_item<1> it) const + { + const std::size_t reduction_lid = it.get_local_id(0); + const std::size_t wg = + it.get_local_range(0); // 0 <= reduction_lid < wg + + const std::size_t iter_gid = it.get_group(0) % iter_gws_; + const std::size_t reduction_batch_id = it.get_group(0) / iter_gws_; + const std::size_t n_reduction_groups = + it.get_group_range(0) / iter_gws_; + + // work-items operates over input with indices + // inp_data_id = reduction_batch_id * wg * reductions_per_wi + m * wg + // + reduction_lid + // for 0 <= m < reductions_per_wi + + const auto &inp_out_iter_offsets_ = inp_out_iter_indexer_(iter_gid); + const auto &inp_iter_offset = inp_out_iter_offsets_.get_first_offset(); + const auto &out_min_iter_offset = + inp_out_iter_offsets_.get_second_offset(); + const auto &out_max_iter_offset = + inp_out_iter_offsets_.get_third_offset(); + + const MinOpT min_op{}; + const MaxOpT max_op{}; + + argT red_min(min_identity_); + argT red_max(max_identity_); + std::size_t arg_reduce_gid0 = + reduction_lid + reduction_batch_id * wg * reductions_per_wi; + for (std::size_t m = 0; m < reductions_per_wi; ++m) { + std::size_t arg_reduce_gid = arg_reduce_gid0 + m * wg; + + if (arg_reduce_gid < reduction_max_gid_) { + auto inp_reduction_offset = + inp_reduced_dims_indexer_(arg_reduce_gid); + auto inp_offset = inp_iter_offset + inp_reduction_offset; + + if constexpr (First) { + const argT val = normalize_bool(inp_min_[inp_offset]); + red_min = min_op(red_min, val); + red_max = max_op(red_max, val); + } + else { + red_min = min_op(red_min, inp_min_[inp_offset]); + red_max = max_op(red_max, inp_max_[inp_offset]); + } + } + } + + auto work_group = it.get_group(); + group_reducer_(work_group, red_min, red_max); + + if (work_group.leader()) { + // each group writes to a different memory location + out_min_[out_min_iter_offset * n_reduction_groups + + reduction_batch_id] = red_min; + out_max_[out_max_iter_offset * n_reduction_groups + + reduction_batch_id] = red_max; + } + } +}; + +template +class minmax_seq_krn; + +template +class minmax_over_group_krn; + +template +sycl::event + submit_minmax_reduction(sycl::queue &exec_q, + const argTy *arg_min, + const argTy *arg_max, + argTy *res_min, + argTy *res_max, + std::size_t wg, + std::size_t iter_nelems, + std::size_t reduction_nelems, + std::size_t reductions_per_wi, + std::size_t reduction_groups, + const InputOutputIterIndexerT &in_out_iter_indexer, + const ReductionIndexerT &reduction_indexer, + const std::vector &depends) +{ + static constexpr argTy min_identity_val = + su_ns::Identity::value; + static constexpr argTy max_identity_val = + su_ns::Identity::value; + + sycl::event red_ev = exec_q.submit([&](sycl::handler &cgh) { + cgh.depends_on(depends); + + auto globalRange = sycl::range<1>{iter_nelems * reduction_groups * wg}; + auto localRange = sycl::range<1>{wg}; + auto ndRange = sycl::nd_range<1>(globalRange, localRange); + + if constexpr (can_use_minmax_group_reducer::value) { + using GroupReducerT = MinMaxGroupReducer; + using KernelName = class minmax_over_group_krn< + argTy, MinOpT, MaxOpT, InputOutputIterIndexerT, + ReductionIndexerT, GroupReducerT, First>; + + cgh.parallel_for( + ndRange, + MinMaxReductionOverGroupFunctor< + argTy, MinOpT, MaxOpT, InputOutputIterIndexerT, + ReductionIndexerT, GroupReducerT, First>( + arg_min, arg_max, res_min, res_max, min_identity_val, + max_identity_val, in_out_iter_indexer, reduction_indexer, + GroupReducerT{}, reduction_nelems, iter_nelems, + reductions_per_wi)); + } + else { + using SlmT = sycl::local_accessor; + SlmT local_memory = SlmT(localRange, cgh); + + using GroupReducerT = + CustomMinMaxGroupReducer; + using KernelName = class minmax_over_group_krn< + argTy, MinOpT, MaxOpT, InputOutputIterIndexerT, + ReductionIndexerT, GroupReducerT, First>; + + cgh.parallel_for( + ndRange, + MinMaxReductionOverGroupFunctor< + argTy, MinOpT, MaxOpT, InputOutputIterIndexerT, + ReductionIndexerT, GroupReducerT, First>( + arg_min, arg_max, res_min, res_max, min_identity_val, + max_identity_val, in_out_iter_indexer, reduction_indexer, + GroupReducerT{local_memory}, reduction_nelems, iter_nelems, + reductions_per_wi)); + } + }); + return red_ev; +} + +/* + Computes the minimum and the maximum in one pass over the input: + `inp_iter_indexer` maps an iteration index to the offset of its reduction + in the input, `res_min_iter_indexer` and `res_max_iter_indexer` map it to + the offsets of its results, and `reduction_indexer` maps a reduction index + to an input offset relative to the iteration offset. +*/ +template +sycl::event + minmax_over_group_temps_impl(sycl::queue &exec_q, + std::size_t iter_nelems, + std::size_t reduction_nelems, + const argTy *arg_tp, + argTy *res_min_tp, + argTy *res_max_tp, + const InputIterIndexerT &inp_iter_indexer, + const ResIterIndexerT &res_min_iter_indexer, + const ResIterIndexerT &res_max_iter_indexer, + const ReductionIndexerT &reduction_indexer, + const std::vector &depends) +{ + using dpnp::tensor::offset_utils::NoOpIndexer; + using dpnp::tensor::offset_utils::Strided1DIndexer; + using dpnp::tensor::offset_utils::ThreeOffsets_CombinedIndexer; + + using ResInputOutputIterIndexerT = + ThreeOffsets_CombinedIndexer; + const ResInputOutputIterIndexerT res_in_out_iter_indexer{ + inp_iter_indexer, res_min_iter_indexer, res_max_iter_indexer}; + + const sycl::device &d = exec_q.get_device(); + const auto &sg_sizes = d.get_info(); + std::size_t wg = choose_workgroup_size<4>(reduction_nelems, sg_sizes); + + if (reduction_nelems < wg) { + static constexpr argTy min_identity_val = + su_ns::Identity::value; + static constexpr argTy max_identity_val = + su_ns::Identity::value; + + sycl::event comp_ev = exec_q.submit([&](sycl::handler &cgh) { + cgh.depends_on(depends); + + using KernelName = class minmax_seq_krn; + + cgh.parallel_for( + sycl::range<1>(iter_nelems), + SequentialMinMaxReduction( + arg_tp, res_min_tp, res_max_tp, min_identity_val, + max_identity_val, res_in_out_iter_indexer, + reduction_indexer, reduction_nelems)); + }); + + return comp_ev; + } + + static constexpr std::size_t preferred_reductions_per_wi = 8; + // prevents running out of resources on CPU + std::size_t max_wg = reduction_detail::get_work_group_size(d); + + if (reduction_nelems <= preferred_reductions_per_wi * max_wg) { + // Perform reduction using one 1 work-group per iteration, + // can output directly to res + + if (iter_nelems == 1) { + // increase GPU occupancy + wg = max_wg; + } + std::size_t reductions_per_wi = + std::max(1, (reduction_nelems + wg - 1) / wg); + + std::size_t reduction_groups = + (reduction_nelems + reductions_per_wi * wg - 1) / + (reductions_per_wi * wg); + assert(reduction_groups == 1); + + sycl::event comp_ev = + submit_minmax_reduction( + exec_q, arg_tp, arg_tp, res_min_tp, res_max_tp, wg, iter_nelems, + reduction_nelems, reductions_per_wi, reduction_groups, + res_in_out_iter_indexer, reduction_indexer, depends); + + return comp_ev; + } + + // more than one work-groups is needed, requires temporaries + std::size_t reduction_groups = + (reduction_nelems + preferred_reductions_per_wi * wg - 1) / + (preferred_reductions_per_wi * wg); + assert(reduction_groups > 1); + + std::size_t second_iter_reduction_groups_ = + (reduction_groups + preferred_reductions_per_wi * wg - 1) / + (preferred_reductions_per_wi * wg); + + // partial minima are followed by partial maxima + const std::size_t tmp_alloc_size = + iter_nelems * (reduction_groups + second_iter_reduction_groups_); + auto tmp_owner = dpnp::tensor::alloc_utils::smart_malloc_device( + 2 * tmp_alloc_size, exec_q); + + argTy *min_tmp = tmp_owner.get(); + argTy *min_tmp2 = min_tmp + reduction_groups * iter_nelems; + argTy *max_tmp = min_tmp + tmp_alloc_size; + argTy *max_tmp2 = max_tmp + reduction_groups * iter_nelems; + + sycl::event first_reduction_ev; + { + using InputOutputIterIndexerT = + ThreeOffsets_CombinedIndexer; + + static constexpr NoOpIndexer noop_tmp_indexer{}; + const InputOutputIterIndexerT in_out_iter_indexer{ + inp_iter_indexer, noop_tmp_indexer, noop_tmp_indexer}; + + first_reduction_ev = submit_minmax_reduction( + exec_q, arg_tp, arg_tp, min_tmp, max_tmp, wg, iter_nelems, + reduction_nelems, preferred_reductions_per_wi, reduction_groups, + in_out_iter_indexer, reduction_indexer, depends); + } + + std::size_t remaining_reduction_nelems = reduction_groups; + sycl::event dependent_ev = first_reduction_ev; + + using TmpReductionIndexerT = NoOpIndexer; + static constexpr TmpReductionIndexerT tmp_reduction_indexer{}; + + while (remaining_reduction_nelems > preferred_reductions_per_wi * max_wg) { + std::size_t reduction_groups_ = (remaining_reduction_nelems + + preferred_reductions_per_wi * wg - 1) / + (preferred_reductions_per_wi * wg); + assert(reduction_groups_ > 1); + + // keep reducing + using InputOutputIterIndexerT = + ThreeOffsets_CombinedIndexer; + + static constexpr NoOpIndexer noop_tmp_indexer{}; + const InputOutputIterIndexerT in_out_iter_indexer{ + Strided1DIndexer{/* size */ iter_nelems, + /* step */ remaining_reduction_nelems}, + noop_tmp_indexer, noop_tmp_indexer}; + + sycl::event partial_reduction_ev = + submit_minmax_reduction( + exec_q, min_tmp, max_tmp, min_tmp2, max_tmp2, wg, iter_nelems, + remaining_reduction_nelems, preferred_reductions_per_wi, + reduction_groups_, in_out_iter_indexer, tmp_reduction_indexer, + {dependent_ev}); + + remaining_reduction_nelems = reduction_groups_; + std::swap(min_tmp, min_tmp2); + std::swap(max_tmp, max_tmp2); + dependent_ev = std::move(partial_reduction_ev); + } + + // final reduction to res + using InputOutputIterIndexerT = + ThreeOffsets_CombinedIndexer; + + const InputOutputIterIndexerT in_out_iter_indexer{ + Strided1DIndexer{/* size */ iter_nelems, + /* step */ remaining_reduction_nelems}, + res_min_iter_indexer, res_max_iter_indexer}; + + wg = max_wg; + std::size_t reductions_per_wi = + std::max(1, (remaining_reduction_nelems + wg - 1) / wg); + + reduction_groups = + (remaining_reduction_nelems + reductions_per_wi * wg - 1) / + (reductions_per_wi * wg); + assert(reduction_groups == 1); + + sycl::event final_reduction_ev = + submit_minmax_reduction( + exec_q, min_tmp, max_tmp, res_min_tp, res_max_tp, wg, iter_nelems, + remaining_reduction_nelems, reductions_per_wi, reduction_groups, + in_out_iter_indexer, tmp_reduction_indexer, {dependent_ev}); + + sycl::event cleanup_host_task_event = + dpnp::tensor::alloc_utils::async_smart_free( + exec_q, {final_reduction_ev}, tmp_owner); + + // FIXME: do not return host-task event + // Instead collect all host-tasks to a list + + return cleanup_host_task_event; +} + +typedef sycl::event (*minmax_strided_impl_fn_ptr)( + sycl::queue &, + std::size_t, + std::size_t, + const char *, + char *, + char *, + int, + const ssize_t *, + ssize_t, + ssize_t, + ssize_t, + int, + const ssize_t *, + ssize_t, + const std::vector &); + +template +sycl::event minmax_over_group_temps_strided_impl( + sycl::queue &exec_q, + std::size_t iter_nelems, // number of reductions (num. of rows in a + // matrix when reducing over rows) + std::size_t reduction_nelems, // size of each reduction (length of rows, + // i.e. number of columns) + const char *arg_cp, + char *res_min_cp, + char *res_max_cp, + int iter_nd, + // shape, input strides, minimum strides, maximum strides + const ssize_t *iter_shape_and_strides, + ssize_t iter_arg_offset, + ssize_t iter_res_min_offset, + ssize_t iter_res_max_offset, + int red_nd, + const ssize_t *reduction_shape_stride, + ssize_t reduction_arg_offset, + const std::vector &depends) +{ + using dpnp::tensor::offset_utils::StridedIndexer; + using dpnp::tensor::offset_utils::UnpackedStridedIndexer; + + const argTy *arg_tp = reinterpret_cast(arg_cp); + argTy *res_min_tp = reinterpret_cast(res_min_cp); + argTy *res_max_tp = reinterpret_cast(res_max_cp); + + // Only 2*iter_nd entries describing shape and strides of iterated + // dimensions of input array from iter_shape_and_strides are going + // to be accessed by inp_iter_indexer + const StridedIndexer inp_iter_indexer{iter_nd, iter_arg_offset, + iter_shape_and_strides}; + const UnpackedStridedIndexer res_min_iter_indexer{ + iter_nd, iter_res_min_offset, + /* shape */ iter_shape_and_strides, + /* strides */ iter_shape_and_strides + 2 * iter_nd}; + const UnpackedStridedIndexer res_max_iter_indexer{ + iter_nd, iter_res_max_offset, + /* shape */ iter_shape_and_strides, + /* strides */ iter_shape_and_strides + 3 * iter_nd}; + const StridedIndexer reduction_indexer{red_nd, reduction_arg_offset, + reduction_shape_stride}; + + return minmax_over_group_temps_impl( + exec_q, iter_nelems, reduction_nelems, arg_tp, res_min_tp, res_max_tp, + inp_iter_indexer, res_min_iter_indexer, res_max_iter_indexer, + reduction_indexer, depends); +} + +typedef sycl::event (*minmax_contig_impl_fn_ptr)( + sycl::queue &, + std::size_t, + std::size_t, + const char *, + char *, + char *, + ssize_t, + ssize_t, + ssize_t, + ssize_t, + const std::vector &); + +template +sycl::event minmax_axis1_over_group_temps_contig_impl( + sycl::queue &exec_q, + std::size_t iter_nelems, // number of reductions (num. of rows in a + // matrix when reducing over rows) + std::size_t reduction_nelems, // size of each reduction (length of rows, + // i.e. number of columns) + const char *arg_cp, + char *res_min_cp, + char *res_max_cp, + ssize_t iter_arg_offset, + ssize_t iter_res_min_offset, + ssize_t iter_res_max_offset, + ssize_t reduction_arg_offset, + const std::vector &depends) +{ + using dpnp::tensor::offset_utils::NoOpIndexer; + using dpnp::tensor::offset_utils::Strided1DIndexer; + + const argTy *arg_tp = reinterpret_cast(arg_cp) + + iter_arg_offset + reduction_arg_offset; + argTy *res_min_tp = + reinterpret_cast(res_min_cp) + iter_res_min_offset; + argTy *res_max_tp = + reinterpret_cast(res_max_cp) + iter_res_max_offset; + + const Strided1DIndexer rows_indexer{/* size */ iter_nelems, + /* step */ reduction_nelems}; + static constexpr NoOpIndexer result_indexer{}; + static constexpr NoOpIndexer reduction_indexer{}; + + return minmax_over_group_temps_impl( + exec_q, iter_nelems, reduction_nelems, arg_tp, res_min_tp, res_max_tp, + rows_indexer, result_indexer, result_indexer, reduction_indexer, + depends); +} + +template +sycl::event minmax_axis0_over_group_temps_contig_impl( + sycl::queue &exec_q, + std::size_t iter_nelems, // number of reductions (num. of cols in a + // matrix when reducing over cols) + std::size_t reduction_nelems, // size of each reduction (length of cols, + // i.e. number of rows) + const char *arg_cp, + char *res_min_cp, + char *res_max_cp, + ssize_t iter_arg_offset, + ssize_t iter_res_min_offset, + ssize_t iter_res_max_offset, + ssize_t reduction_arg_offset, + const std::vector &depends) +{ + using dpnp::tensor::offset_utils::NoOpIndexer; + using dpnp::tensor::offset_utils::Strided1DIndexer; + + const argTy *arg_tp = reinterpret_cast(arg_cp) + + iter_arg_offset + reduction_arg_offset; + argTy *res_min_tp = + reinterpret_cast(res_min_cp) + iter_res_min_offset; + argTy *res_max_tp = + reinterpret_cast(res_max_cp) + iter_res_max_offset; + + static constexpr NoOpIndexer columns_indexer{}; + static constexpr NoOpIndexer result_indexer{}; + const Strided1DIndexer reduction_indexer{/* size */ reduction_nelems, + /* step */ iter_nelems}; + + return minmax_over_group_temps_impl( + exec_q, iter_nelems, reduction_nelems, arg_tp, res_min_tp, res_max_tp, + columns_indexer, result_indexer, result_indexer, reduction_indexer, + depends); +} + } // namespace dpnp::tensor::kernels diff --git a/dpnp/tensor/libtensor/source/reductions/minmax.cpp b/dpnp/tensor/libtensor/source/reductions/minmax.cpp new file mode 100644 index 000000000000..672d82bd5a5a --- /dev/null +++ b/dpnp/tensor/libtensor/source/reductions/minmax.cpp @@ -0,0 +1,469 @@ +//***************************************************************************** +// Copyright (c) 2026, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// - Neither the name of the copyright holder nor the names of its contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +// THE POSSIBILITY OF SUCH DAMAGE. +//***************************************************************************** +// +//===---------------------------------------------------------------------===// +/// +/// \file +/// This file defines functions of dpnp.tensor._tensor_reductions_impl +/// extension. +//===---------------------------------------------------------------------===// + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "dpnp4pybind11.hpp" +#include +#include + +#include "kernels/reductions.hpp" +#include "simplify_iteration_space.hpp" +#include "utils/memory_overlap.hpp" +#include "utils/offset_utils.hpp" +#include "utils/output_validation.hpp" +#include "utils/sycl_alloc_utils.hpp" +#include "utils/sycl_utils.hpp" +#include "utils/type_dispatch.hpp" +#include "utils/type_dispatch_building.hpp" + +#include "minmax.hpp" + +namespace dpnp::tensor::py_internal +{ + +namespace py = pybind11; +namespace td_ns = dpnp::tensor::type_dispatch; +namespace su_ns = dpnp::tensor::sycl_utils; + +namespace impl +{ + +using dpnp::tensor::kernels::minmax_strided_impl_fn_ptr; +static minmax_strided_impl_fn_ptr + minmax_over_axis_strided_temps_dispatch_vector[td_ns::num_types]; + +using dpnp::tensor::kernels::minmax_contig_impl_fn_ptr; +static minmax_contig_impl_fn_ptr + minmax_over_axis1_contig_temps_dispatch_vector[td_ns::num_types]; +static minmax_contig_impl_fn_ptr + minmax_over_axis0_contig_temps_dispatch_vector[td_ns::num_types]; + +/* @brief Operators used by min and max reductions with temporaries */ +template +using MinOpT = + std::conditional_t && !std::is_same_v, + sycl::minimum, + su_ns::Minimum>; + +template +using MaxOpT = + std::conditional_t && !std::is_same_v, + sycl::maximum, + su_ns::Maximum>; + +template +struct MinMaxOverAxisTempsStridedFactory +{ + fnT get() const + { + return dpnp::tensor::kernels::minmax_over_group_temps_strided_impl< + srcTy, MinOpT, MaxOpT>; + } +}; + +template +struct MinMaxOverAxis1TempsContigFactory +{ + fnT get() const + { + return dpnp::tensor::kernels::minmax_axis1_over_group_temps_contig_impl< + srcTy, MinOpT, MaxOpT>; + } +}; + +template +struct MinMaxOverAxis0TempsContigFactory +{ + fnT get() const + { + return dpnp::tensor::kernels::minmax_axis0_over_group_temps_contig_impl< + srcTy, MinOpT, MaxOpT>; + } +}; + +void populate_minmax_over_axis_dispatch_vectors(void) +{ + using td_ns::DispatchVectorBuilder; + + DispatchVectorBuilder + dvb1; + dvb1.populate_dispatch_vector( + minmax_over_axis_strided_temps_dispatch_vector); + + DispatchVectorBuilder + dvb2; + dvb2.populate_dispatch_vector( + minmax_over_axis1_contig_temps_dispatch_vector); + + DispatchVectorBuilder + dvb3; + dvb3.populate_dispatch_vector( + minmax_over_axis0_contig_temps_dispatch_vector); +} + +} // namespace impl + +/*! @brief Computes the minimum and the maximum over the trailing axes in a + * single pass */ +std::pair + py_minmax_over_axis(const dpnp::tensor::usm_ndarray &src, + int trailing_dims_to_reduce, // comp over this many + // trailing indexes + const dpnp::tensor::usm_ndarray &dst_min, + const dpnp::tensor::usm_ndarray &dst_max, + sycl::queue &exec_q, + const std::vector &depends) +{ + int src_nd = src.get_ndim(); + int iteration_nd = src_nd - trailing_dims_to_reduce; + if (trailing_dims_to_reduce <= 0 || iteration_nd < 0) { + throw py::value_error("Trailing_dim_to_reduce must be positive, but no " + "greater than rank of the array being reduced"); + } + + int dst_nd = dst_min.get_ndim(); + if (dst_nd != iteration_nd || dst_max.get_ndim() != iteration_nd) { + throw py::value_error("Destination array rank does not match input " + "array rank and number of reduced dimensions"); + } + + const py::ssize_t *src_shape_ptr = src.get_shape_raw(); + const py::ssize_t *dst_min_shape_ptr = dst_min.get_shape_raw(); + const py::ssize_t *dst_max_shape_ptr = dst_max.get_shape_raw(); + + bool same_shapes = true; + for (int i = 0; same_shapes && (i < dst_nd); ++i) { + same_shapes = same_shapes && + (src_shape_ptr[i] == dst_min_shape_ptr[i]) && + (src_shape_ptr[i] == dst_max_shape_ptr[i]); + } + + if (!same_shapes) { + throw py::value_error("Destination shape does not match unreduced " + "dimensions of the input shape"); + } + + if (!dpnp::utils::queues_are_compatible(exec_q, {src, dst_min, dst_max})) { + throw py::value_error( + "Execution queue is not compatible with allocation queues"); + } + + dpnp::tensor::validation::CheckWritable::throw_if_not_writable(dst_min); + dpnp::tensor::validation::CheckWritable::throw_if_not_writable(dst_max); + + std::size_t dst_nelems = dst_min.get_size(); + + if (dst_nelems == 0) { + return std::make_pair(sycl::event(), sycl::event()); + } + + std::size_t reduction_nelems(1); + for (int i = dst_nd; i < src_nd; ++i) { + reduction_nelems *= static_cast(src_shape_ptr[i]); + } + + if (reduction_nelems == 0) { + throw py::value_error( + "reduction cannot be performed over zero-size axes"); + } + + // check that the arrays do not overlap + auto const &overlap = dpnp::tensor::overlap::MemoryOverlap(); + if (overlap(src, dst_min) || overlap(src, dst_max) || + overlap(dst_min, dst_max)) { + throw py::value_error("Arrays index overlapping segments of memory"); + } + + dpnp::tensor::validation::AmpleMemory::throw_if_not_ample(dst_min, + dst_nelems); + dpnp::tensor::validation::AmpleMemory::throw_if_not_ample(dst_max, + dst_nelems); + + int src_typenum = src.get_typenum(); + if (src_typenum != dst_min.get_typenum() || + src_typenum != dst_max.get_typenum()) { + throw py::value_error("Destination arrays must have the same data " + "type as the input array"); + } + + const auto &array_types = td_ns::usm_ndarray_types(); + int src_typeid = array_types.typenum_to_lookup_id(src_typenum); + + using impl::minmax_over_axis0_contig_temps_dispatch_vector; + using impl::minmax_over_axis1_contig_temps_dispatch_vector; + using impl::minmax_over_axis_strided_temps_dispatch_vector; + + // handle special case when both reduction and iteration are 1D contiguous + bool is_src_c_contig = src.is_c_contiguous(); + bool is_src_f_contig = src.is_f_contiguous(); + bool is_dst_c_contig = + dst_min.is_c_contiguous() && dst_max.is_c_contiguous(); + bool is_dst_f_contig = + dst_min.is_f_contiguous() && dst_max.is_f_contiguous(); + + static constexpr py::ssize_t zero_offset = 0; + + if ((is_src_c_contig && is_dst_c_contig) || + (is_src_f_contig && dst_nelems == 1)) { + auto fn = minmax_over_axis1_contig_temps_dispatch_vector[src_typeid]; + if (fn != nullptr) { + std::size_t iter_nelems = dst_nelems; + + sycl::event minmax_over_axis_contig_ev = + fn(exec_q, iter_nelems, reduction_nelems, src.get_data(), + dst_min.get_data(), dst_max.get_data(), + zero_offset, // iteration_src_offset + zero_offset, // iteration_dst_min_offset + zero_offset, // iteration_dst_max_offset + zero_offset, // reduction_src_offset + depends); + + sycl::event keep_args_event = dpnp::utils::keep_args_alive( + exec_q, {src, dst_min, dst_max}, {minmax_over_axis_contig_ev}); + + return std::make_pair(keep_args_event, minmax_over_axis_contig_ev); + } + } + else if (is_src_f_contig && + ((is_dst_c_contig && dst_nd == 1) || is_dst_f_contig)) { + auto fn = minmax_over_axis0_contig_temps_dispatch_vector[src_typeid]; + if (fn != nullptr) { + std::size_t iter_nelems = dst_nelems; + + sycl::event minmax_over_axis_contig_ev = + fn(exec_q, iter_nelems, reduction_nelems, src.get_data(), + dst_min.get_data(), dst_max.get_data(), + zero_offset, // iteration_src_offset + zero_offset, // iteration_dst_min_offset + zero_offset, // iteration_dst_max_offset + zero_offset, // reduction_src_offset + depends); + + sycl::event keep_args_event = dpnp::utils::keep_args_alive( + exec_q, {src, dst_min, dst_max}, {minmax_over_axis_contig_ev}); + + return std::make_pair(keep_args_event, minmax_over_axis_contig_ev); + } + } + + auto const &src_shape_vecs = src.get_shape_vector(); + auto const &src_strides_vecs = src.get_strides_vector(); + auto const &dst_min_strides_vecs = dst_min.get_strides_vector(); + auto const &dst_max_strides_vecs = dst_max.get_strides_vector(); + + int reduction_nd = trailing_dims_to_reduce; + const py::ssize_t *reduction_shape_ptr = src_shape_ptr + dst_nd; + using shT = std::vector; + shT reduction_src_strides(std::begin(src_strides_vecs) + dst_nd, + std::end(src_strides_vecs)); + + shT simplified_reduction_shape; + shT simplified_reduction_src_strides; + py::ssize_t reduction_src_offset(0); + + simplify_iteration_space_1( + reduction_nd, reduction_shape_ptr, reduction_src_strides, + // output + simplified_reduction_shape, simplified_reduction_src_strides, + reduction_src_offset); + + const py::ssize_t *iteration_shape_ptr = src_shape_ptr; + + shT iteration_src_strides(std::begin(src_strides_vecs), + std::begin(src_strides_vecs) + iteration_nd); + + shT simplified_iteration_shape; + shT simplified_iteration_src_strides; + shT simplified_iteration_dst_min_strides; + shT simplified_iteration_dst_max_strides; + py::ssize_t iteration_src_offset(0); + py::ssize_t iteration_dst_min_offset(0); + py::ssize_t iteration_dst_max_offset(0); + + if (iteration_nd == 0) { + if (dst_nelems != 1) { + throw std::runtime_error("iteration_nd == 0, but dst_nelems != 1"); + } + iteration_nd = 1; + simplified_iteration_shape.push_back(1); + simplified_iteration_src_strides.push_back(0); + simplified_iteration_dst_min_strides.push_back(0); + simplified_iteration_dst_max_strides.push_back(0); + } + else { + simplify_iteration_space_3( + iteration_nd, iteration_shape_ptr, iteration_src_strides, + dst_min_strides_vecs, dst_max_strides_vecs, + // output + simplified_iteration_shape, simplified_iteration_src_strides, + simplified_iteration_dst_min_strides, + simplified_iteration_dst_max_strides, iteration_src_offset, + iteration_dst_min_offset, iteration_dst_max_offset); + } + + if ((reduction_nd == 1) && (iteration_nd == 1)) { + bool mat_reduce_over_axis1 = false; + bool mat_reduce_over_axis0 = false; + bool array_reduce_all_elems = false; + std::size_t iter_nelems = dst_nelems; + + bool unit_dst_strides = + (simplified_iteration_dst_min_strides[0] == 1) && + (simplified_iteration_dst_max_strides[0] == 1); + + if (simplified_reduction_src_strides[0] == 1) { + array_reduce_all_elems = (simplified_iteration_shape[0] == 1); + mat_reduce_over_axis1 = + unit_dst_strides && + (static_cast( + simplified_iteration_src_strides[0]) == reduction_nelems); + } + else if (static_cast( + simplified_reduction_src_strides[0]) == iter_nelems) { + mat_reduce_over_axis0 = + unit_dst_strides && (simplified_iteration_src_strides[0] == 1); + } + + if (mat_reduce_over_axis1 || array_reduce_all_elems) { + auto fn = + minmax_over_axis1_contig_temps_dispatch_vector[src_typeid]; + if (fn != nullptr) { + sycl::event minmax_over_axis1_contig_ev = + fn(exec_q, iter_nelems, reduction_nelems, src.get_data(), + dst_min.get_data(), dst_max.get_data(), + iteration_src_offset, iteration_dst_min_offset, + iteration_dst_max_offset, reduction_src_offset, depends); + + sycl::event keep_args_event = dpnp::utils::keep_args_alive( + exec_q, {src, dst_min, dst_max}, + {minmax_over_axis1_contig_ev}); + + return std::make_pair(keep_args_event, + minmax_over_axis1_contig_ev); + } + } + else if (mat_reduce_over_axis0) { + auto fn = + minmax_over_axis0_contig_temps_dispatch_vector[src_typeid]; + if (fn != nullptr) { + sycl::event minmax_over_axis0_contig_ev = + fn(exec_q, iter_nelems, reduction_nelems, src.get_data(), + dst_min.get_data(), dst_max.get_data(), + iteration_src_offset, iteration_dst_min_offset, + iteration_dst_max_offset, reduction_src_offset, depends); + + sycl::event keep_args_event = dpnp::utils::keep_args_alive( + exec_q, {src, dst_min, dst_max}, + {minmax_over_axis0_contig_ev}); + + return std::make_pair(keep_args_event, + minmax_over_axis0_contig_ev); + } + } + } + + auto fn = minmax_over_axis_strided_temps_dispatch_vector[src_typeid]; + if (fn == nullptr) { + throw std::runtime_error("Datatypes are not supported"); + } + + std::vector host_task_events{}; + using dpnp::tensor::offset_utils::device_allocate_and_pack; + auto arrays_metainfo_packing_triple_ = + device_allocate_and_pack( + exec_q, host_task_events, + // iteration metadata + simplified_iteration_shape, simplified_iteration_src_strides, + simplified_iteration_dst_min_strides, + simplified_iteration_dst_max_strides, + // reduction metadata + simplified_reduction_shape, simplified_reduction_src_strides); + auto tmp_owner = std::move(std::get<0>(arrays_metainfo_packing_triple_)); + const auto ©_metadata_ev = std::get<2>(arrays_metainfo_packing_triple_); + const py::ssize_t *temp_allocation_ptr = tmp_owner.get(); + + const py::ssize_t *iter_shape_and_strides = temp_allocation_ptr; + const py::ssize_t *reduction_shape_stride = + temp_allocation_ptr + 4 * simplified_iteration_shape.size(); + + std::vector all_deps; + all_deps.reserve(depends.size() + 1); + all_deps.resize(depends.size()); + std::copy(depends.begin(), depends.end(), all_deps.begin()); + all_deps.push_back(copy_metadata_ev); + + auto minmax_ev = fn(exec_q, dst_nelems, reduction_nelems, src.get_data(), + dst_min.get_data(), dst_max.get_data(), iteration_nd, + iter_shape_and_strides, iteration_src_offset, + iteration_dst_min_offset, iteration_dst_max_offset, + reduction_nd, // number dimensions being reduced + reduction_shape_stride, reduction_src_offset, all_deps); + + sycl::event temp_cleanup_ev = dpnp::tensor::alloc_utils::async_smart_free( + exec_q, {minmax_ev}, tmp_owner); + host_task_events.push_back(temp_cleanup_ev); + + sycl::event keep_args_event = dpnp::utils::keep_args_alive( + exec_q, {src, dst_min, dst_max}, host_task_events); + + return std::make_pair(keep_args_event, minmax_ev); +} + +void init_minmax(py::module_ m) +{ + impl::populate_minmax_over_axis_dispatch_vectors(); + + m.def("_minmax_over_axis", &py_minmax_over_axis, + "Computes the minimum and the maximum of `src` over its trailing " + "`trailing_dims_to_reduce` axes in a single pass, storing them to " + "`dst_min` and `dst_max`.", + py::arg("src"), py::arg("trailing_dims_to_reduce"), + py::arg("dst_min"), py::arg("dst_max"), py::arg("sycl_queue"), + py::arg("depends") = py::list()); +} + +} // namespace dpnp::tensor::py_internal diff --git a/dpnp/tensor/libtensor/source/reductions/minmax.hpp b/dpnp/tensor/libtensor/source/reductions/minmax.hpp new file mode 100644 index 000000000000..41ae9673b77c --- /dev/null +++ b/dpnp/tensor/libtensor/source/reductions/minmax.hpp @@ -0,0 +1,46 @@ +//***************************************************************************** +// Copyright (c) 2026, Intel Corporation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// - Neither the name of the copyright holder nor the names of its contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +// THE POSSIBILITY OF SUCH DAMAGE. +//***************************************************************************** +// +//===---------------------------------------------------------------------===// +/// +/// \file +/// This file defines functions of dpnp.tensor._tensor_reductions_impl +/// extension. +//===---------------------------------------------------------------------===// + +#pragma once +#include + +namespace py = pybind11; + +namespace dpnp::tensor::py_internal +{ + +extern void init_minmax(py::module_ m); + +} // namespace dpnp::tensor::py_internal diff --git a/dpnp/tensor/libtensor/source/reductions/reduction_common.cpp b/dpnp/tensor/libtensor/source/reductions/reduction_common.cpp index b99cfce24d41..22218f2b2f91 100644 --- a/dpnp/tensor/libtensor/source/reductions/reduction_common.cpp +++ b/dpnp/tensor/libtensor/source/reductions/reduction_common.cpp @@ -42,6 +42,7 @@ #include "logsumexp.hpp" #include "max.hpp" #include "min.hpp" +#include "minmax.hpp" #include "prod.hpp" #include "reduce_hypot.hpp" #include "sum.hpp" @@ -61,6 +62,7 @@ void init_reduction_functions(py::module_ m) init_logsumexp(m); init_max(m); init_min(m); + init_minmax(m); init_prod(m); init_reduce_hypot(m); init_sum(m); diff --git a/dpnp/tests/test_statistics.py b/dpnp/tests/test_statistics.py index 3d5387d3de39..f9393a9241f7 100644 --- a/dpnp/tests/test_statistics.py +++ b/dpnp/tests/test_statistics.py @@ -853,6 +853,238 @@ def test_error(self, func): getattr(dpnp, func)(ia, initial=6) +class TestMinMax: + @pytest.mark.parametrize("axis", [None, 0, 1, -1, 2, -2, (1, 2), (0, -2)]) + @pytest.mark.parametrize("keepdims", [False, True]) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) + def test_basic(self, axis, keepdims, dtype): + a = generate_random_numpy_array((4, 4, 6, 8), dtype=dtype) + ia = dpnp.array(a) + + result = dpnp.minmax(ia, axis=axis, keepdims=keepdims) + assert isinstance(result, tuple) and len(result) == 2 + expected = numpy.min(a, axis=axis, keepdims=keepdims) + assert_dtype_allclose(result[0], expected) + expected = numpy.max(a, axis=axis, keepdims=keepdims) + assert_dtype_allclose(result[1], expected) + + # sequential, one work-group and several work-groups per reduction + @pytest.mark.parametrize("n", [7, 1000, 100_000]) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) + def test_size(self, n, dtype): + a = generate_random_numpy_array((3, n), dtype=dtype) + ia = dpnp.array(a) + + for x, ix in [(a, ia), (a[:, ::-2], ia[:, ::-2])]: + for axis in [None, 0, 1]: + result = dpnp.minmax(ix, axis=axis) + assert_dtype_allclose(result[0], numpy.min(x, axis=axis)) + assert_dtype_allclose(result[1], numpy.max(x, axis=axis)) + + # several passes over partial results + @pytest.mark.parametrize("dtype", [dpnp.int8, dpnp.float32]) + def test_large(self, dtype): + a = generate_random_numpy_array(6_000_000, dtype=dtype) + a[123] = -100 + a[-5] = 100 + ia = dpnp.array(a) + + result = dpnp.minmax(ia) + assert_dtype_allclose(result[0], numpy.min(a)) + assert_dtype_allclose(result[1], numpy.max(a)) + + @pytest.mark.parametrize("order", ["C", "F"]) + @pytest.mark.parametrize("axis", [None, 0, 1, 2, (0, 2), (1, 2)]) + def test_layout(self, order, axis): + a = generate_random_numpy_array((7, 300, 6), order=order) + ia = dpnp.array(a, order=order) + + for x, ix in [(a, ia), (a[::2, ::-3], ia[::2, ::-3])]: + result = dpnp.minmax(ix, axis=axis) + assert_dtype_allclose(result[0], numpy.min(x, axis=axis)) + assert_dtype_allclose(result[1], numpy.max(x, axis=axis)) + + @pytest.mark.parametrize("keepdims", [False, True]) + def test_out_strided(self, keepdims): + a = generate_random_numpy_array((5, 600, 7), dtype=dpnp.int32) + ia = dpnp.array(a) + exp_min = numpy.min(a, axis=1, keepdims=keepdims) + exp_max = numpy.max(a, axis=1, keepdims=keepdims) + + # both outputs with different strides + out_min = dpnp.empty((5, 14), dtype=ia.dtype)[:, ::2] + out_max = dpnp.empty((7, 5), dtype=ia.dtype).T + if keepdims: + out_min, out_max = out_min[:, None], out_max[:, None] + + result = dpnp.minmax( + ia, axis=1, keepdims=keepdims, out=(out_min, out_max) + ) + assert result[0] is out_min and result[1] is out_max + assert_array_equal(result[0], exp_min) + assert_array_equal(result[1], exp_max) + + @with_requires("numpy>=2.6.0.dev0") + @pytest.mark.parametrize("axis", [None, 0, 1, (0, 1)]) + @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) + def test_numpy(self, axis, dtype): + a = generate_random_numpy_array((3, 4), dtype=dtype) + ia = dpnp.array(a) + + result = dpnp.minmax(ia, axis) + expected = numpy.minmax(a, axis) + assert_dtype_allclose(result[0], expected[0]) + assert_dtype_allclose(result[1], expected[1]) + + @pytest.mark.parametrize("axis", [None, 0, 1]) + @pytest.mark.parametrize("shape", [(4, 5), (3, 1000), (3, 100_000)]) + @pytest.mark.parametrize("dtype", get_float_complex_dtypes()) + def test_nan(self, axis, shape, dtype): + a = generate_random_numpy_array(shape, dtype=dtype) + a[0, 1] = numpy.inf + a[0, -2] = -numpy.inf + a[1, 2] = a[-1, -1] = numpy.nan + ia = dpnp.array(a) + + result = dpnp.minmax(ia, axis=axis) + assert_dtype_allclose(result[0], numpy.min(a, axis=axis)) + assert_dtype_allclose(result[1], numpy.max(a, axis=axis)) + + def test_0d(self): + ia = dpnp.array(3) + result = dpnp.minmax(ia) + assert_array_equal(result[0], 3) + assert_array_equal(result[1], 3) + + def test_usm_ndarray(self): + a = numpy.array([[3, 1, 2], [5, 0, 4]]) + ia = dpt.asarray(a) + + res_min, res_max = dpnp.minmax(ia, axis=1) + assert isinstance(res_min, dpnp.ndarray) + assert isinstance(res_max, dpnp.ndarray) + assert_array_equal(res_min, numpy.min(a, axis=1)) + assert_array_equal(res_max, numpy.max(a, axis=1)) + + def test_out(self): + a = numpy.arange(12, dtype=numpy.float32).reshape((2, 2, 3)) + ia = dpnp.array(a) + exp_min, exp_max = numpy.min(a, axis=0), numpy.max(a, axis=0) + + # out is a tuple of dpnp_array + out = (dpnp.empty_like(ia[0]), dpnp.empty_like(ia[0])) + result = dpnp.minmax(ia, axis=0, out=out) + assert result[0] is out[0] and result[1] is out[1] + assert_allclose(result[0], exp_min) + assert_allclose(result[1], exp_max) + + # out is a tuple of usm_ndarray + dpt_out = dpt.empty((2, 3), dtype=ia.dtype) + out = (dpt_out, dpt.empty_like(dpt_out)) + result = dpnp.minmax(ia, axis=0, out=out) + assert result[0].get_array() is out[0] + assert result[1].get_array() is out[1] + assert_allclose(result[0], exp_min) + assert_allclose(result[1], exp_max) + + # out has None entry + out_max = dpnp.empty_like(ia[0]) + result = dpnp.minmax(ia, axis=0, out=(None, out_max)) + assert result[1] is out_max + assert_allclose(result[0], exp_min) + assert_allclose(result[1], exp_max) + + def test_out_dtype(self): + a = numpy.arange(12, dtype=numpy.float32).reshape((3, 4)) - 2.5 + ia = dpnp.array(a) + + out = (dpnp.empty(4, dtype=dpnp.int32), dpnp.empty(4, dtype="i2")) + result = dpnp.minmax(ia, axis=0, out=out) + assert result[0] is out[0] and result[1] is out[1] + + expected = numpy.min(a, axis=0, out=numpy.empty(4, dtype=numpy.int32)) + assert_dtype_allclose(result[0], expected) + expected = numpy.max(a, axis=0, out=numpy.empty(4, dtype="i2")) + assert_dtype_allclose(result[1], expected) + + @pytest.mark.parametrize("idx", [0, 1]) + def test_out_overlap_input(self, idx): + # the 1st row holds the maximum of some columns + a = numpy.array([[5, 1, 9, 3], [1, 7, 0, 2], [2, 3, 4, 8]]) + ia = dpnp.array(a) + exp_min, exp_max = numpy.min(a, axis=0), numpy.max(a, axis=0) + + out = [None, None] + out[idx] = ia[0] + result = dpnp.minmax(ia, axis=0, out=tuple(out)) + assert result[idx] is out[idx] + assert_array_equal(result[0], exp_min) + assert_array_equal(result[1], exp_max) + + def test_out_error(self): + ia = dpnp.arange(12).reshape((3, 4)) + out = dpnp.empty(4, dtype=ia.dtype) + + # out is not a tuple + with pytest.raises(TypeError, match="must be a tuple"): + dpnp.minmax(ia, axis=0, out=out) + with pytest.raises(TypeError, match="must be a tuple"): + dpnp.minmax(ia, axis=0, out=[out, dpnp.empty_like(out)]) + + # out has incorrect number of entries + with pytest.raises(ValueError, match="exactly two entries"): + dpnp.minmax(ia, axis=0, out=(out,)) + with pytest.raises(ValueError, match="exactly two entries"): + dpnp.minmax(ia, axis=0, out=(out, None, None)) + + # out entry is numpy array + with pytest.raises(TypeError): + dpnp.minmax(ia, axis=0, out=(numpy.empty(4), None)) + with pytest.raises(TypeError): + dpnp.minmax(ia, axis=0, out=(None, numpy.empty(4))) + + # out entry has incorrect shape + with pytest.raises(ValueError): + dpnp.minmax(ia, axis=0, out=(dpnp.empty(3), None)) + with pytest.raises(ValueError): + dpnp.minmax(ia, axis=0, out=(None, dpnp.empty(3))) + + # out entries overlap + with pytest.raises(ValueError, match="cannot overlap"): + dpnp.minmax(ia, axis=0, out=(out, out)) + buf = dpnp.empty(6, dtype=ia.dtype) + with pytest.raises(ValueError, match="cannot overlap"): + dpnp.minmax(ia, axis=0, out=(buf[:4], buf[2:])) + + @pytest.mark.parametrize("axis", [None, 0]) + def test_empty_error(self, axis): + ia = dpnp.empty((0, 3)) + with pytest.raises(ValueError): + dpnp.minmax(ia, axis=axis) + + def test_empty(self): + a = numpy.empty((0, 3), dtype=numpy.float32) + ia = dpnp.array(a) + + result = dpnp.minmax(ia, axis=1) + assert_dtype_allclose(result[0], numpy.min(a, axis=1)) + assert_dtype_allclose(result[1], numpy.max(a, axis=1)) + + def test_error(self): + ia = dpnp.arange(5) + # where is not supported + with pytest.raises(NotImplementedError): + dpnp.minmax(ia, where=False) + + # initial is not supported + with pytest.raises(NotImplementedError): + dpnp.minmax(ia, initial=6) + + # input is numpy array + with pytest.raises(TypeError): + dpnp.minmax(numpy.arange(5)) + + class TestMean: @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) @pytest.mark.parametrize("axis", [None, 0, 1, (0, 1)]) diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py index e0f882f2175d..0a4fc10c3f31 100644 --- a/dpnp/tests/test_sycl_queue.py +++ b/dpnp/tests/test_sycl_queue.py @@ -378,6 +378,7 @@ def test_1in_1out(func, data, device): "func, data", [ pytest.param("frexp", numpy.arange(9)), + pytest.param("minmax", [1.0, 2.0, 4.0, 7.0]), pytest.param("modf", [0, 3.5]), ], ) diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py index 567078c3eff2..8292c10daecc 100644 --- a/dpnp/tests/test_usm_type.py +++ b/dpnp/tests/test_usm_type.py @@ -660,6 +660,7 @@ def test_1in_1out(func, data, usm_type): "func, data", [ pytest.param("frexp", numpy.arange(9)), + pytest.param("minmax", [1.0, 2.0, 4.0, 7.0]), pytest.param("modf", [0, 3.5]), ], ) From 01f4f9ea4d1af33843e4c75100b4e6130e27aa18 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Fri, 30 Oct 2026 19:03:20 +0100 Subject: [PATCH 2/4] Extend dpnp.minmax tests and drop an unused shape vector Cover float16, NaN at the first, middle and last position (also a complex NaN in the imaginary part only), integer limits of all integer types, multi-pass reductions of several rows, axis=() with out, and read-only and other-queue out errors. Link dpnp.minmax from dpnp.min and dpnp.max. --- dpnp/dpnp_iface_statistics.py | 2 + .../libtensor/source/reductions/minmax.cpp | 1 - dpnp/tests/test_statistics.py | 82 +++++++++++++++++-- 3 files changed, 75 insertions(+), 10 deletions(-) diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index 8a5baec22f82..c39a9240732f 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -1096,6 +1096,7 @@ def max(a, axis=None, out=None, keepdims=False, initial=None, where=True): See Also -------- :obj:`dpnp.min` : Return the minimum of an array. + :obj:`dpnp.minmax` : Return the minimum and maximum of an array. :obj:`dpnp.maximum` : Element-wise maximum of two arrays, propagates NaNs. :obj:`dpnp.fmax` : Element-wise maximum of two arrays, ignores NaNs. :obj:`dpnp.nanmax` : The maximum value of an array along a given axis, @@ -1371,6 +1372,7 @@ def min(a, axis=None, out=None, keepdims=False, initial=None, where=True): See Also -------- :obj:`dpnp.max` : Return the maximum of an array. + :obj:`dpnp.minmax` : Return the minimum and maximum of an array. :obj:`dpnp.minimum` : Element-wise minimum of two arrays, propagates NaNs. :obj:`dpnp.fmin` : Element-wise minimum of two arrays, ignores NaNs. :obj:`dpnp.nanmin` : The minimum value of an array along a given axis, diff --git a/dpnp/tensor/libtensor/source/reductions/minmax.cpp b/dpnp/tensor/libtensor/source/reductions/minmax.cpp index 672d82bd5a5a..855d89028386 100644 --- a/dpnp/tensor/libtensor/source/reductions/minmax.cpp +++ b/dpnp/tensor/libtensor/source/reductions/minmax.cpp @@ -290,7 +290,6 @@ std::pair } } - auto const &src_shape_vecs = src.get_shape_vector(); auto const &src_strides_vecs = src.get_strides_vector(); auto const &dst_min_strides_vecs = dst_min.get_strides_vector(); auto const &dst_max_strides_vecs = dst_max.get_strides_vector(); diff --git a/dpnp/tests/test_statistics.py b/dpnp/tests/test_statistics.py index f9393a9241f7..0fedf7359e45 100644 --- a/dpnp/tests/test_statistics.py +++ b/dpnp/tests/test_statistics.py @@ -10,6 +10,7 @@ import dpnp import dpnp.tensor as dpt +from dpnp.exceptions import ExecutionPlacementError from .helper import ( assert_dtype_allclose, @@ -19,6 +20,7 @@ get_complex_dtypes, get_float_complex_dtypes, get_float_dtypes, + get_integer_dtypes, has_support_aspect64, numpy_version, ) @@ -856,7 +858,9 @@ def test_error(self, func): class TestMinMax: @pytest.mark.parametrize("axis", [None, 0, 1, -1, 2, -2, (1, 2), (0, -2)]) @pytest.mark.parametrize("keepdims", [False, True]) - @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) + @pytest.mark.parametrize( + "dtype", get_all_dtypes(no_none=True, no_float16=False) + ) def test_basic(self, axis, keepdims, dtype): a = generate_random_numpy_array((4, 4, 6, 8), dtype=dtype) ia = dpnp.array(a) @@ -870,7 +874,9 @@ def test_basic(self, axis, keepdims, dtype): # sequential, one work-group and several work-groups per reduction @pytest.mark.parametrize("n", [7, 1000, 100_000]) - @pytest.mark.parametrize("dtype", get_all_dtypes(no_none=True)) + @pytest.mark.parametrize( + "dtype", get_all_dtypes(no_none=True, no_float16=False) + ) def test_size(self, n, dtype): a = generate_random_numpy_array((3, n), dtype=dtype) ia = dpnp.array(a) @@ -882,16 +888,17 @@ def test_size(self, n, dtype): assert_dtype_allclose(result[1], numpy.max(x, axis=axis)) # several passes over partial results + @pytest.mark.parametrize("axis", [None, 1]) @pytest.mark.parametrize("dtype", [dpnp.int8, dpnp.float32]) - def test_large(self, dtype): - a = generate_random_numpy_array(6_000_000, dtype=dtype) - a[123] = -100 - a[-5] = 100 + def test_large(self, axis, dtype): + a = generate_random_numpy_array((2, 6_000_000), dtype=dtype) + a[0, 123] = -100 + a[1, -5] = 100 ia = dpnp.array(a) - result = dpnp.minmax(ia) - assert_dtype_allclose(result[0], numpy.min(a)) - assert_dtype_allclose(result[1], numpy.max(a)) + result = dpnp.minmax(ia, axis=axis) + assert_dtype_allclose(result[0], numpy.min(a, axis=axis)) + assert_dtype_allclose(result[1], numpy.max(a, axis=axis)) @pytest.mark.parametrize("order", ["C", "F"]) @pytest.mark.parametrize("axis", [None, 0, 1, 2, (0, 2), (1, 2)]) @@ -950,6 +957,35 @@ def test_nan(self, axis, shape, dtype): assert_dtype_allclose(result[0], numpy.min(a, axis=axis)) assert_dtype_allclose(result[1], numpy.max(a, axis=axis)) + @pytest.mark.parametrize("n", [7, 1000, 100_000]) + @pytest.mark.parametrize( + "dtype", get_float_complex_dtypes(no_float16=False) + ) + def test_nan_position(self, n, dtype): + a = generate_random_numpy_array(n, dtype=dtype) + nans = [numpy.nan] + if dpnp.issubdtype(dtype, dpnp.complexfloating): + nans.append(complex(0, numpy.nan)) + + for pos in [0, n // 2, -1]: + for nan in nans: + x = a.copy() + x[pos] = nan + result = dpnp.minmax(dpnp.array(x)) + assert_dtype_allclose(result[0], numpy.min(x)) + assert_dtype_allclose(result[1], numpy.max(x)) + + # all items equal to the identity of the other reduction + @pytest.mark.parametrize("n", [7, 1000, 100_000]) + @pytest.mark.parametrize("dtype", get_integer_dtypes(all_int_types=True)) + def test_integer_limits(self, n, dtype): + info = numpy.iinfo(dtype) + for val in [info.min, info.max]: + a = numpy.full(n, val, dtype=dtype) + result = dpnp.minmax(dpnp.array(a)) + assert_array_equal(result[0], numpy.min(a)) + assert_array_equal(result[1], numpy.max(a)) + def test_0d(self): ia = dpnp.array(3) result = dpnp.minmax(ia) @@ -1007,6 +1043,23 @@ def test_out_dtype(self): expected = numpy.max(a, axis=0, out=numpy.empty(4, dtype="i2")) assert_dtype_allclose(result[1], expected) + def test_out_no_axis(self): + a = numpy.arange(6, dtype=numpy.float32).reshape((2, 3)) + ia = dpnp.array(a) + expected = numpy.min(a, axis=()), numpy.max(a, axis=()) + + out = (dpnp.empty_like(ia), dpnp.empty_like(ia)) + result = dpnp.minmax(ia, axis=(), out=out) + assert result[0] is out[0] and result[1] is out[1] + assert_array_equal(result[0], expected[0]) + assert_array_equal(result[1], expected[1]) + + # out overlaps the input + result = dpnp.minmax(ia, axis=(), out=(ia, None)) + assert result[0] is ia + assert_array_equal(result[0], expected[0]) + assert_array_equal(result[1], expected[1]) + @pytest.mark.parametrize("idx", [0, 1]) def test_out_overlap_input(self, idx): # the 1st row holds the maximum of some columns @@ -1049,6 +1102,17 @@ def test_out_error(self): with pytest.raises(ValueError): dpnp.minmax(ia, axis=0, out=(None, dpnp.empty(3))) + # out entry is read-only + out_read_only = dpnp.empty_like(out) + out_read_only.flags.writable = False + with pytest.raises(ValueError, match="read-only"): + dpnp.minmax(ia, axis=0, out=(None, out_read_only)) + + # out entry is allocated on another queue + out_q = dpnp.empty_like(out, sycl_queue=dpctl.SyclQueue()) + with pytest.raises((ValueError, ExecutionPlacementError)): + dpnp.minmax(ia, axis=0, out=(out_q, None)) + # out entries overlap with pytest.raises(ValueError, match="cannot overlap"): dpnp.minmax(ia, axis=0, out=(out, out)) From 2100f682bf4b03f244b238ac499842314cccfb84 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Fri, 30 Oct 2026 19:04:38 +0100 Subject: [PATCH 3/4] Set the PR link in the changelog entry --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 520fea0102c6..fc87d782c59b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added * Added support for boolean input in `dpnp.positive` and the `+` operator of `dpnp.ndarray` [#3128](https://github.com/IntelPython/dpnp/pull/3128) -* Added implementation of `dpnp.minmax`, computing the minimum and the maximum of an array in a single pass by a new fused reduction kernel [#XXXX](https://github.com/IntelPython/dpnp/pull/XXXX) +* Added implementation of `dpnp.minmax`, computing the minimum and the maximum of an array in a single pass by a new fused reduction kernel [#3135](https://github.com/IntelPython/dpnp/pull/3135) ### Changed From d7c57b0f093fa6843bec3a2efeccaa9a8ccc7ddb Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Fri, 30 Oct 2026 19:29:46 +0100 Subject: [PATCH 4/4] Document the out differences of dpnp.minmax from NumPy --- dpnp/dpnp_iface_statistics.py | 8 ++++++++ dpnp/tests/test_statistics.py | 16 ++++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py index c39a9240732f..0d297ef82dca 100644 --- a/dpnp/dpnp_iface_statistics.py +++ b/dpnp/dpnp_iface_statistics.py @@ -1461,6 +1461,14 @@ def minmax(a, axis=None, out=None, keepdims=False, initial=None, where=True): Parameters `where`, and `initial` are only supported with their default values. Otherwise ``NotImplementedError`` exception will be raised. + If an array in `out` has a data type other than that of `a`, the result is + computed in the data type of `a` and then cast, while NumPy casts `a` to + the data type of `out` first. The results differ when the cast does not + preserve the order of values, e.g. for a narrowing integer cast. + + The arrays in `out` must not overlap in memory. Interleaved views, such as + ``(r[:, 0], r[:, 1])``, are rejected as well, while NumPy accepts them. + See Also -------- :obj:`dpnp.min` : Return the minimum of an array. diff --git a/dpnp/tests/test_statistics.py b/dpnp/tests/test_statistics.py index 0fedf7359e45..18f2b27bb5be 100644 --- a/dpnp/tests/test_statistics.py +++ b/dpnp/tests/test_statistics.py @@ -1043,6 +1043,17 @@ def test_out_dtype(self): expected = numpy.max(a, axis=0, out=numpy.empty(4, dtype="i2")) assert_dtype_allclose(result[1], expected) + def test_out_dtype_narrowing(self): + a = numpy.array([300, 100]) + ia = dpnp.array(a) + + out = (dpnp.empty((), dtype=dpnp.int8), dpnp.empty((), dtype=dpnp.int8)) + result = dpnp.minmax(ia, out=out) + + # unlike NumPy, the result is cast after the reduction, not the input + assert_array_equal(result[0], numpy.min(a).astype(numpy.int8)) + assert_array_equal(result[1], numpy.max(a).astype(numpy.int8)) + def test_out_no_axis(self): a = numpy.arange(6, dtype=numpy.float32).reshape((2, 3)) ia = dpnp.array(a) @@ -1120,6 +1131,11 @@ def test_out_error(self): with pytest.raises(ValueError, match="cannot overlap"): dpnp.minmax(ia, axis=0, out=(buf[:4], buf[2:])) + # unlike NumPy, interleaved out entries are rejected as overlapping + res = dpnp.empty((4, 2), dtype=ia.dtype) + with pytest.raises(ValueError, match="cannot overlap"): + dpnp.minmax(ia, axis=0, out=(res[:, 0], res[:, 1])) + @pytest.mark.parametrize("axis", [None, 0]) def test_empty_error(self, axis): ia = dpnp.empty((0, 3))