diff --git a/cpp/include/raft/matrix/detail/select_k-inl.cuh b/cpp/include/raft/matrix/detail/select_k-inl.cuh
index bbdbcb886f..b0dddb0a3d 100644
--- a/cpp/include/raft/matrix/detail/select_k-inl.cuh
+++ b/cpp/include/raft/matrix/detail/select_k-inl.cuh
@@ -299,6 +299,17 @@ void select_k(raft::resources const& handle,
       return detail::select::warpsort::
         select_k_impl<T, IdxT, detail::select::warpsort::warp_sort_distributed_ext, row_layout>(
           handle, in_val, in_idx, batch_size, len, k, out_val, out_idx, select_min);
+    case SelectAlgo::kWarpDistributedShmStable:
+      // Only the keys that the IVF searches select on, so other instantiations do not grow.
+      if constexpr (std::is_same_v<T, float> || std::is_same_v<T, half>) {
+        return detail::select::warpsort::select_k_impl<
+          T,
+          IdxT,
+          detail::select::warpsort::warp_sort_distributed_ext_stable,
+          row_layout>(handle, in_val, in_idx, batch_size, len, k, out_val, out_idx, select_min);
+      } else {
+        RAFT_FAIL("kWarpDistributedShmStable supports float and half keys only");
+      }
     case SelectAlgo::kWarpAuto:
       return detail::select::warpsort::select_k<T, IdxT, row_layout>(
         handle, in_val, in_idx, batch_size, len, k, out_val, out_idx, select_min);
diff --git a/cpp/include/raft/matrix/detail/select_warpsort.cuh b/cpp/include/raft/matrix/detail/select_warpsort.cuh
index 3341156cea..51b6a2bc20 100644
--- a/cpp/include/raft/matrix/detail/select_warpsort.cuh
+++ b/cpp/include/raft/matrix/detail/select_warpsort.cuh
@@ -23,6 +23,7 @@
 #include <rmm/resource_ref.hpp>
 
 #include <cub/util_type.cuh>  // cub::Traits
+#include <cuda/std/variant>   // cuda::std::monostate
 #include <cuda/stream>
 
 #include <algorithm>
@@ -109,6 +110,16 @@ _RAFT_DEVICE _RAFT_FORCEINLINE auto is_ordered(T left, T right) -> bool
   if constexpr (!Ascending) { return left > right; }
 }
 
+/** `is_ordered` with an equal-key tie broken on the index, as `swap_needed` does. */
+template <bool Ascending, typename T, typename IdxT>
+_RAFT_DEVICE _RAFT_FORCEINLINE auto is_ordered(T left, T right, IdxT left_idx, IdxT right_idx)
+  -> bool
+{
+  if (left != right) { return is_ordered<Ascending>(left, right); }
+  if constexpr (Ascending) { return left_idx < right_idx; }
+  if constexpr (!Ascending) { return left_idx > right_idx; }
+}
+
 }  // namespace
 
 /**
@@ -126,8 +137,12 @@ _RAFT_DEVICE _RAFT_FORCEINLINE auto is_ordered(T left, T right) -> bool
  * @tparam IdxT
  *   the type of payload (normally, indices of elements), i.e.
  *   the content sorted alongside the keys.
+ * @tparam TieBreak
+ *   order an equal-key tie on the index rather than on arrival, so that the selected elements
+ *   do not depend on the input order. Without it the code is the same as before the parameter
+ *   existed.
  */
-template <int Capacity, bool Ascending, typename T, typename IdxT>
+template <int Capacity, bool Ascending, typename T, typename IdxT, bool TieBreak = false>
 class warp_sort {
   static_assert(is_a_power_of_two(Capacity));
   static_assert(std::is_default_constructible_v<IdxT>);
@@ -138,6 +153,13 @@ class warp_sort {
    *  i.e. `Ascending ? upper_bound<T>() : lower_bound<T>()`.
    */
   static constexpr T kDummy = Ascending ? upper_bound<T>() : lower_bound<T>();
+  /**
+   * The index of an empty slot. With `TieBreak` it loses every tie, so an empty slot never
+   * outranks a real element whose key equals `kDummy`.
+   */
+  static constexpr IdxT kDummyIdx =
+    TieBreak ? (Ascending ? upper_bound<IdxT>() : lower_bound<IdxT>()) : IdxT{};
+  static constexpr bool kTieBreak = TieBreak;
   /** Width of the subwarp. */
   static constexpr int kWarpWidth = std::min<int>(Capacity, WarpSize);
   /** The number of elements to select. */
@@ -157,7 +179,7 @@ class warp_sort {
 #pragma unroll
     for (int i = 0; i < kMaxArrLen; i++) {
       val_arr_[i] = kDummy;
-      idx_arr_[i] = IdxT{};
+      idx_arr_[i] = kDummyIdx;
     }
   }
 
@@ -188,7 +210,13 @@ class warp_sort {
       for (int i = kMaxArrLen - 1; i >= 0; --i, idx += kWarpWidth) {
         if (idx < k) {
           T t = in[idx];
-          if (is_ordered<Ascending>(t, val_arr_[i])) {
+          bool take;
+          if constexpr (TieBreak) {
+            take = is_ordered<Ascending>(t, val_arr_[i], in_idx[idx], idx_arr_[i]);
+          } else {
+            take = is_ordered<Ascending>(t, val_arr_[i]);
+          }
+          if (take) {
             val_arr_[i] = t;
             idx_arr_[i] = in_idx[idx];
           }
@@ -196,7 +224,7 @@ class warp_sort {
       }
     }
     if (kWarpWidth < WarpSize || do_merge) {
-      util::bitonic<kMaxArrLen>(Ascending, kWarpWidth).merge(val_arr_, idx_arr_);
+      util::bitonic<kMaxArrLen, TieBreak>(Ascending, kWarpWidth).merge(val_arr_, idx_arr_);
     }
   }
 
@@ -259,12 +287,19 @@ class warp_sort {
     for (int i = std::min(kMaxArrLen, PerThreadSizeIn); i > 0; i--) {
       T& key  = val_arr_[kMaxArrLen - i];
       T other = keys_in[PerThreadSizeIn - i];
-      if (is_ordered<Ascending>(other, key)) {
+      bool take;
+      if constexpr (TieBreak) {
+        take =
+          is_ordered<Ascending>(other, key, ids_in[PerThreadSizeIn - i], idx_arr_[kMaxArrLen - i]);
+      } else {
+        take = is_ordered<Ascending>(other, key);
+      }
+      if (take) {
         key                      = other;
         idx_arr_[kMaxArrLen - i] = ids_in[PerThreadSizeIn - i];
       }
     }
-    util::bitonic<kMaxArrLen>(Ascending, kWarpWidth).merge(val_arr_, idx_arr_);
+    util::bitonic<kMaxArrLen, TieBreak>(Ascending, kWarpWidth).merge(val_arr_, idx_arr_);
   }
 };
 
@@ -480,12 +515,13 @@ class warp_sort_distributed : public warp_sort<Capacity, Ascending, T, IdxT> {
  * The same as `warp_sort_distributed`, but keeps the temporary value and index buffers
  * in the given external pointers (normally, a shared memory pointer should be passed in).
  */
-template <int Capacity, bool Ascending, typename T, typename IdxT>
-class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT> {
+template <int Capacity, bool Ascending, typename T, typename IdxT, bool TieBreak = false>
+class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT, TieBreak> {
  public:
-  using warp_sort<Capacity, Ascending, T, IdxT>::kDummy;
-  using warp_sort<Capacity, Ascending, T, IdxT>::kWarpWidth;
-  using warp_sort<Capacity, Ascending, T, IdxT>::k;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::kDummy;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::kDummyIdx;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::kWarpWidth;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::k;
 
   constexpr static auto mem_required(uint32_t block_size) -> size_t
   {
@@ -493,13 +529,17 @@ class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT>
   }
 
   _RAFT_DEVICE warp_sort_distributed_ext(int k, T* val_buf, IdxT* idx_buf, T limit = kDummy)
-    : warp_sort<Capacity, Ascending, T, IdxT>(k),
+    : warp_sort<Capacity, Ascending, T, IdxT, TieBreak>(k),
       val_buf_(val_buf),
       idx_buf_(idx_buf),
       buf_len_(0),
       k_th_(limit)
   {
     val_buf_[laneId()] = kDummy;
+    if constexpr (TieBreak) {
+      k_th_idx_          = kDummyIdx;
+      idx_buf_[laneId()] = kDummyIdx;
+    }
   }
 
   _RAFT_DEVICE static auto init_blockwide(int k, uint8_t* shmem, T limit = kDummy)
@@ -516,12 +556,18 @@ class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT>
     auto warp_offset = Pow2<WarpSize>::roundDown(threadIdx.x);
     val_buf += warp_offset;
     idx_buf += warp_offset;
-    return warp_sort_distributed_ext<Capacity, Ascending, T, IdxT>{k, val_buf, idx_buf, limit};
+    return warp_sort_distributed_ext<Capacity, Ascending, T, IdxT, TieBreak>{
+      k, val_buf, idx_buf, limit};
   }
 
   _RAFT_DEVICE void add(T val, IdxT idx)
   {
-    bool do_add = is_ordered<Ascending>(val, k_th_);
+    bool do_add;
+    if constexpr (TieBreak) {
+      do_add = is_ordered<Ascending>(val, k_th_, idx, k_th_idx_);
+    } else {
+      do_add = is_ordered<Ascending>(val, k_th_);
+    }
     // mask tells which lanes in the warp have valid items to be added
     uint32_t mask = ballot(do_add);
     if (mask == 0) { return; }
@@ -563,6 +609,7 @@ class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT>
     // NB on using srcLane: it's ok if it is outside the warp size / width;
     //                      the modulo op will be done inside the __shfl_sync.
     k_th_ = shfl(val_arr_[kMaxArrLen - 1], k - 1, kWarpWidth);
+    if constexpr (TieBreak) { k_th_idx_ = shfl(idx_arr_[kMaxArrLen - 1], k - 1, kWarpWidth); }
   }
 
   _RAFT_DEVICE _RAFT_FORCEINLINE void merge_buf_()
@@ -571,20 +618,43 @@ class warp_sort_distributed_ext : public warp_sort<Capacity, Ascending, T, IdxT>
     T buf_val          = val_buf_[laneId()];
     IdxT buf_idx       = idx_buf_[laneId()];
     val_buf_[laneId()] = kDummy;
-    util::bitonic<1>(!Ascending, kWarpWidth).sort(buf_val, buf_idx);
+    if constexpr (TieBreak) { idx_buf_[laneId()] = kDummyIdx; }
+    util::bitonic<1, TieBreak>(!Ascending, kWarpWidth).sort(buf_val, buf_idx);
     this->template merge_in<1>(&buf_val, &buf_idx);
     set_k_th_();  // contains warp sync
   }
 
-  using warp_sort<Capacity, Ascending, T, IdxT>::kMaxArrLen;
-  using warp_sort<Capacity, Ascending, T, IdxT>::val_arr_;
-  using warp_sort<Capacity, Ascending, T, IdxT>::idx_arr_;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::kMaxArrLen;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::val_arr_;
+  using warp_sort<Capacity, Ascending, T, IdxT, TieBreak>::idx_arr_;
 
   T* val_buf_;
   IdxT* idx_buf_;
   uint32_t buf_len_;  // 0 <= buf_len_ < WarpSize
 
   T k_th_;
+  [[no_unique_address]] std::conditional_t<TieBreak, IdxT, cuda::std::monostate> k_th_idx_;
+};
+
+/**
+ * `warp_sort_distributed_ext` that orders an equal-key tie on the index, so the selected elements
+ * do not depend on the input order. Selected by `SelectAlgo::kWarpDistributedShmStable`.
+ */
+template <int Capacity, bool Ascending, typename T, typename IdxT>
+class warp_sort_distributed_ext_stable
+  : public warp_sort_distributed_ext<Capacity, Ascending, T, IdxT, true> {
+  using base_t = warp_sort_distributed_ext<Capacity, Ascending, T, IdxT, true>;
+
+ public:
+  using base_t::base_t;
+
+  _RAFT_DEVICE static auto init_blockwide(int k, uint8_t* shmem, T limit = base_t::kDummy)
+  {
+    return warp_sort_distributed_ext_stable{base_t::init_blockwide(k, shmem, limit)};
+  }
+
+ private:
+  _RAFT_DEVICE explicit warp_sort_distributed_ext_stable(const base_t& queue) : base_t(queue) {}
 };
 
 /**
@@ -786,6 +856,7 @@ __launch_bounds__(256) RAFT_KERNEL block_kernel(const T* in,
   //   2. Fewer expensive template instantiations
   using bits_t       = typename cub::Traits<T>::UnsignedBits;
   using bq_t         = block_sort<WarpSortClass, Capacity, Ascending, bits_t, IdxT>;
+  using warp_sort_t  = WarpSortClass<Capacity, Ascending, bits_t, IdxT>;
   uint8_t* warp_smem = bq_t::queue_t::mem_required(blockDim.x) > 0 ? smem_buf_bytes : nullptr;
   bq_t queue(k, warp_smem);
 
@@ -795,8 +866,10 @@ __launch_bounds__(256) RAFT_KERNEL block_kernel(const T* in,
   for (IdxT i = threadIdx.x + blockIdx.x * blockDim.x; i < per_thread_lim; i += stride) {
     // Twiddle the input value to ensure proper comparison of floating-point and signed int values
     queue.add(i < len ? cub::Traits<T>::TwiddleIn(__ldcs(reinterpret_cast<const bits_t*>(in) + i))
-                      : WarpSortClass<Capacity, Ascending, bits_t, IdxT>::kDummy,
-              (i < len && in_idx != nullptr) ? __ldcs(in_idx + i) : i);
+                      : warp_sort_t::kDummy,
+              (i < len && in_idx != nullptr)         ? __ldcs(in_idx + i)
+              : (warp_sort_t::kTieBreak && i >= len) ? warp_sort_t::kDummyIdx
+                                                     : i);
   }
 
   // * write out the result
@@ -974,6 +1047,12 @@ struct LaunchThreshold<warp_sort_distributed_ext> {
   static constexpr int len_factor_for_single_block = 32;
 };
 
+template <>
+struct LaunchThreshold<warp_sort_distributed_ext_stable> {
+  static constexpr int len_factor_for_multi_block  = 2;
+  static constexpr int len_factor_for_single_block = 32;
+};
+
 template <>
 struct LaunchThreshold<warp_sort_immediate> {
   static constexpr int len_factor_for_choosing     = 4;
diff --git a/cpp/include/raft/matrix/select_k_types.hpp b/cpp/include/raft/matrix/select_k_types.hpp
index c6f42468e4..8c239bb531 100644
--- a/cpp/include/raft/matrix/select_k_types.hpp
+++ b/cpp/include/raft/matrix/select_k_types.hpp
@@ -1,5 +1,5 @@
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION.
+ * SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 
@@ -67,6 +67,11 @@ enum class SelectAlgo : uint8_t {
    * in the given external pointers (normally, a shared memory pointer should be passed in).
    */
   kWarpDistributedShm = 8,
+  /**
+   * The same as `kWarpDistributedShm`, but orders an equal-key tie by the payload index instead of
+   * by lane arrival, so a repeated selection over the same input returns the same elements.
+   */
+  kWarpDistributedShmStable = 9,
 };
 
 inline auto operator<<(std::ostream& os, const SelectAlgo& algo) -> std::ostream&
@@ -74,6 +79,8 @@ inline auto operator<<(std::ostream& os, const SelectAlgo& algo) -> std::ostream
   auto underlying_value = static_cast<std::underlying_type<SelectAlgo>::type>(algo);
 
   switch (algo) {
+    case SelectAlgo::kWarpDistributedShmStable:
+      return os << "kWarpDistributedShmStable=" << underlying_value;
     case SelectAlgo::kAuto: return os << "kAuto=" << underlying_value;
     case SelectAlgo::kRadix8bits: return os << "kRadix8bits=" << underlying_value;
     case SelectAlgo::kRadix11bits: return os << "kRadix11bits=" << underlying_value;
diff --git a/cpp/include/raft/util/bitonic_sort.cuh b/cpp/include/raft/util/bitonic_sort.cuh
index 318b47e0dd..23aca23bde 100644
--- a/cpp/include/raft/util/bitonic_sort.cuh
+++ b/cpp/include/raft/util/bitonic_sort.cuh
@@ -1,5 +1,5 @@
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION.
+ * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 
@@ -27,6 +27,22 @@ _RAFT_DEVICE _RAFT_FORCEINLINE void conditional_assign(bool cond, T& ptr, T x)
   if (cond) { ptr = x; }
 }
 
+/** First element of a payload pack. k-selection always passes the index first. */
+template <typename P0, typename... Rest>
+_RAFT_DEVICE _RAFT_FORCEINLINE P0 first_payload(P0 head, Rest...)
+{
+  return head;
+}
+
+/** Whether `a` must move after `b`, breaking an equal-key tie on the index. */
+template <typename KeyT, typename IdxT>
+_RAFT_DEVICE _RAFT_FORCEINLINE bool swap_needed(bool ascending, KeyT a, KeyT b, IdxT ia, IdxT ib)
+{
+  if (a != b) { return ascending ? (a > b) : (a < b); }
+  // Must flip with `ascending`, or paired lanes disagree on a tie and duplicate an element.
+  return ascending ? (ia > ib) : (ia < ib);
+}
+
 }  // namespace
 
 /**
@@ -67,9 +83,13 @@ _RAFT_DEVICE _RAFT_FORCEINLINE void conditional_assign(bool cond, T& ptr, T x)
  *   number of elements processed in each thread;
  *   i.e. the total data size is `Size * warp_width`.
  *   Must be power-of-two.
+ * @tparam TieBreak
+ *   order an equal-key tie on the first payload, which k-selection uses for the index, so
+ *   that a repeated sort of the same input gives the same order. Without it the code is the
+ *   same as before the parameter existed.
  *
  */
-template <int Size = 1>
+template <int Size = 1, bool TieBreak = false>
 class bitonic {
   static_assert(isPo2(Size));
 
@@ -115,7 +135,7 @@ class bitonic {
   _RAFT_DEVICE _RAFT_FORCEINLINE void merge(KeyT* __restrict__ keys,
                                             PayloadTs* __restrict__... payloads) const
   {
-    return bitonic<Size>::merge_impl(ascending_, warp_width_, keys, payloads...);
+    return bitonic<Size, TieBreak>::merge_impl(ascending_, warp_width_, keys, payloads...);
   }
 
   /**
@@ -134,7 +154,7 @@ class bitonic {
   _RAFT_DEVICE _RAFT_FORCEINLINE void sort(KeyT* __restrict__ keys,
                                            PayloadTs* __restrict__... payloads) const
   {
-    return bitonic<Size>::sort_impl(ascending_, warp_width_, keys, payloads...);
+    return bitonic<Size, TieBreak>::sort_impl(ascending_, warp_width_, keys, payloads...);
   }
 
   /**
@@ -173,7 +193,7 @@ class bitonic {
   const int warp_width_;
   const bool ascending_;
 
-  template <int AnotherSize>
+  template <int AnotherSize, bool AnotherTieBreak>
   friend class bitonic;
 
   template <typename KeyT, typename... PayloadTs>
@@ -192,7 +212,17 @@ class bitonic {
           const int other_i = i + stride;
           KeyT& key         = keys[i];
           KeyT& other       = keys[other_i];
-          if (ascending ? key > other : key < other) {
+          bool do_swap;
+          if constexpr (TieBreak && sizeof...(PayloadTs) > 0) {
+            do_swap = swap_needed(ascending,
+                                  key,
+                                  other,
+                                  first_payload(payloads[i]...),
+                                  first_payload(payloads[other_i]...));
+          } else {
+            do_swap = ascending ? key > other : key < other;
+          }
+          if (do_swap) {
             swap(key, other);
             (swap(payloads[i], payloads[other_i]), ...);
           }
@@ -206,7 +236,16 @@ class bitonic {
       for (int stride = (warp_width >> 1); stride > 0; stride >>= 1) {
         const bool is_second = lane & stride;
         const KeyT other     = shfl_xor(key, stride, warp_width);
-        const bool do_assign = (ascending != is_second) ? key > other : key < other;
+        const bool asc       = (ascending != is_second);
+        bool do_assign;
+        if constexpr (TieBreak && sizeof...(PayloadTs) > 0) {
+          // The shuffle is unconditional: it is a warp-collective operation.
+          const auto my_id    = first_payload(payloads[i]...);
+          const auto other_id = shfl_xor(my_id, stride, warp_width);
+          do_assign           = swap_needed(asc, key, other, my_id, other_id);
+        } else {
+          do_assign = asc ? key > other : key < other;
+        }
 
         conditional_assign(do_assign, key, other);
         // NB: don't put shfl_xor in a conditional; it must be called by all threads in a warp.
@@ -225,14 +264,14 @@ class bitonic {
     if constexpr (Size == 1) {
       const int lane = laneId();
       for (int width = 2; width < warp_width; width <<= 1) {
-        bitonic<1>::merge_impl(lane & width, width, keys, payloads...);
+        bitonic<1, TieBreak>::merge_impl(lane & width, width, keys, payloads...);
       }
     } else {
       constexpr int kSize2 = Size / 2;
-      bitonic<kSize2>::sort_impl(false, warp_width, keys, payloads...);
-      bitonic<kSize2>::sort_impl(true, warp_width, keys + kSize2, (payloads + kSize2)...);
+      bitonic<kSize2, TieBreak>::sort_impl(false, warp_width, keys, payloads...);
+      bitonic<kSize2, TieBreak>::sort_impl(true, warp_width, keys + kSize2, (payloads + kSize2)...);
     }
-    bitonic<Size>::merge_impl(ascending, warp_width, keys, payloads...);
+    bitonic<Size, TieBreak>::merge_impl(ascending, warp_width, keys, payloads...);
   }
 };
 
diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt
index e1f7a1f755..2d12815e1d 100644
--- a/cpp/tests/CMakeLists.txt
+++ b/cpp/tests/CMakeLists.txt
@@ -219,8 +219,8 @@ if(BUILD_TESTS)
   )
 
   ConfigureTest(
-    NAME MATRIX_SELECT_TEST PATH matrix/select_k.cu matrix/select_k_edgecases.cu LIB
-    EXPLICIT_INSTANTIATE_ONLY GPUS 1 PERCENT 100
+    NAME MATRIX_SELECT_TEST PATH matrix/select_k.cu matrix/select_k_edgecases.cu
+    matrix/select_k_stable.cu LIB EXPLICIT_INSTANTIATE_ONLY GPUS 1 PERCENT 100
   )
 
   ConfigureTest(
diff --git a/cpp/tests/matrix/select_k_stable.cu b/cpp/tests/matrix/select_k_stable.cu
new file mode 100644
index 0000000000..b9e802cdde
--- /dev/null
+++ b/cpp/tests/matrix/select_k_stable.cu
@@ -0,0 +1,155 @@
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+
+#include <raft/core/copy.hpp>
+#include <raft/core/device_mdarray.hpp>
+#include <raft/core/device_resources.hpp>
+#include <raft/core/resource/cuda_stream.hpp>
+#include <raft/matrix/select_k.cuh>
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <cstdint>
+#include <cstring>
+#include <numeric>
+#include <optional>
+#include <random>
+#include <vector>
+
+namespace raft::matrix {
+
+namespace {
+
+// The ids that `kWarpDistributedShmStable` must return, sorted: the best keys, and on an equal key
+// the smaller id for select-min and the larger id for select-max.
+auto expected_ids(const std::vector<float>& keys,
+                  const std::vector<int64_t>& ids,
+                  int k,
+                  bool select_min) -> std::vector<int64_t>
+{
+  std::vector<size_t> order(keys.size());
+  std::iota(order.begin(), order.end(), 0);
+  std::sort(order.begin(), order.end(), [&](size_t a, size_t b) {
+    if (keys[a] != keys[b]) { return select_min ? keys[a] < keys[b] : keys[a] > keys[b]; }
+    return select_min ? ids[a] < ids[b] : ids[a] > ids[b];
+  });
+  std::vector<int64_t> out;
+  for (int i = 0; i < k; i++) {
+    out.push_back(ids[order[i]]);
+  }
+  std::sort(out.begin(), out.end());
+  return out;
+}
+
+auto stable_select(const raft::device_resources& handle,
+                   const std::vector<float>& keys,
+                   const std::vector<int64_t>& ids,
+                   int k,
+                   bool select_min) -> std::vector<int64_t>
+{
+  auto stream   = raft::resource::get_cuda_stream(handle);
+  int64_t len   = keys.size();
+  auto in_keys  = raft::make_device_matrix<float, int64_t>(handle, 1, len);
+  auto in_ids   = raft::make_device_matrix<int64_t, int64_t>(handle, 1, len);
+  auto out_keys = raft::make_device_matrix<float, int64_t>(handle, 1, k);
+  auto out_ids  = raft::make_device_matrix<int64_t, int64_t>(handle, 1, k);
+  raft::copy(in_keys.data_handle(), keys.data(), len, stream);
+  raft::copy(in_ids.data_handle(), ids.data(), len, stream);
+  raft::matrix::select_k<float, int64_t>(handle,
+                                         raft::make_const_mdspan(in_keys.view()),
+                                         std::make_optional(raft::make_const_mdspan(in_ids.view())),
+                                         out_keys.view(),
+                                         out_ids.view(),
+                                         select_min,
+                                         /*sorted=*/true,
+                                         SelectAlgo::kWarpDistributedShmStable);
+  std::vector<int64_t> out(k);
+  raft::copy(out.data(), out_ids.data_handle(), k, stream);
+  raft::resource::sync_stream(handle, stream);
+  std::sort(out.begin(), out.end());
+  return out;
+}
+
+// Selects from the same (key, id) pairs in two input orders.
+void check_order_independent(const std::vector<float>& keys, const std::vector<int64_t>& ids, int k)
+{
+  raft::device_resources handle;
+  std::mt19937 rng(1234);
+  for (bool select_min : {true, false}) {
+    auto expected = expected_ids(keys, ids, k, select_min);
+    for (int order = 0; order < 2; order++) {
+      std::vector<size_t> perm(keys.size());
+      std::iota(perm.begin(), perm.end(), 0);
+      std::shuffle(perm.begin(), perm.end(), rng);
+      std::vector<float> perm_keys(keys.size());
+      std::vector<int64_t> perm_ids(keys.size());
+      for (size_t i = 0; i < perm.size(); i++) {
+        perm_keys[i] = keys[perm[i]];
+        perm_ids[i]  = ids[perm[i]];
+      }
+      EXPECT_EQ(stable_select(handle, perm_keys, perm_ids, k, select_min), expected)
+        << "select_min=" << select_min << " order=" << order << " len=" << keys.size()
+        << " k=" << k;
+    }
+  }
+}
+
+auto distinct_ids(size_t n) -> std::vector<int64_t>
+{
+  std::vector<int64_t> ids(n);
+  std::iota(ids.begin(), ids.end(), 0);
+  std::shuffle(ids.begin(), ids.end(), std::mt19937(42));
+  return ids;
+}
+
+auto from_bits(uint32_t bits) -> float
+{
+  float f;
+  std::memcpy(&f, &bits, sizeof(f));
+  return f;
+}
+
+}  // namespace
+
+TEST(SelectKStable, AllKeysEqual)
+{
+  for (size_t len : {4096, 300000}) {
+    for (int k : {10, 100}) {
+      check_order_independent(std::vector<float>(len, 1.0f), distinct_ids(len), k);
+    }
+  }
+}
+
+TEST(SelectKStable, TiesAtTheKthPlace)
+{
+  for (size_t len : {4096, 300000}) {
+    std::vector<float> keys(len);
+    std::mt19937 rng(7);
+    for (auto& key : keys) {
+      key = float(rng() % 4);
+    }
+    check_order_independent(keys, distinct_ids(len), 100);
+  }
+}
+
+// A real key that twiddles to the same bits as the empty-slot key must outrank the empty slots.
+// These NaN bit patterns twiddle to the all-ones (select-min) and all-zeros (select-max) keys.
+TEST(SelectKStable, RealKeysEqualToTheEmptyKeyWin)
+{
+  raft::device_resources handle;
+  std::vector<int64_t> ids = {11, 3, 7, 5, 9};
+  for (bool select_min : {true, false}) {
+    float key = from_bits(select_min ? 0x7FFFFFFFu : 0xFFFFFFFFu);
+    std::vector<float> keys(ids.size(), key);
+    auto got = stable_select(handle, keys, ids, 8, select_min);
+    for (int64_t id : ids) {
+      EXPECT_NE(std::find(got.begin(), got.end(), id), got.end())
+        << "select_min=" << select_min << " missing id " << id;
+    }
+  }
+}
+
+}  // namespace raft::matrix