diff --git a/cpp/include/raft/cache/cache_util.cuh b/cpp/include/raft/cache/cache_util.cuh
index ce8ef9a095..a65227c402 100644
--- a/cpp/include/raft/cache/cache_util.cuh
+++ b/cpp/include/raft/cache/cache_util.cuh
@@ -41,9 +41,9 @@ namespace cache {
  * @param [in] n the number of elements that need to be collected
  * @param [out] out vectors collected from the cache, size [n_vec * n]
  */
-template <typename math_t>
-__global__ void get_vecs(const math_t *cache, int n_vec, const int *cache_idx,
-                         int n, math_t *out) {
+template <typename math_t, typename idx_t, typename int_t>
+__global__ void get_vecs(const math_t *cache, int_t n_vec,
+                         const idx_t *cache_idx, int_t n, math_t *out) {
   int tid = threadIdx.x + blockIdx.x * blockDim.x;
   int row = tid % n_vec;  // row idx
   if (tid < n_vec * n) {
diff --git a/cpp/include/raft/sparse/linalg/degree.cuh b/cpp/include/raft/sparse/linalg/degree.cuh
index 9bd322c90a..ef6a067c39 100644
--- a/cpp/include/raft/sparse/linalg/degree.cuh
+++ b/cpp/include/raft/sparse/linalg/degree.cuh
@@ -43,11 +43,11 @@ namespace linalg {
  * @param nnz the size of the rows array
  * @param results array to place results
  */
-template <int TPB_X = 64>
-__global__ void coo_degree_kernel(const int *rows, int nnz, int *results) {
+template <int TPB_X = 64, typename T = int>
+__global__ void coo_degree_kernel(const T *rows, int nnz, T *results) {
   int row = (blockIdx.x * TPB_X) + threadIdx.x;
   if (row < nnz) {
-    raft::myAtomicAdd(results + rows[row], 1);
+    atomicAdd(results + rows[row], (T)1);
   }
 }
 
@@ -59,8 +59,8 @@ __global__ void coo_degree_kernel(const int *rows, int nnz, int *results) {
  * @param results: output result array
  * @param stream: cuda stream to use
  */
-template <int TPB_X = 64>
-void coo_degree(const int *rows, int nnz, int *results, cudaStream_t stream) {
+template <int TPB_X = 64, typename T = int>
+void coo_degree(const T *rows, int nnz, T *results, cudaStream_t stream) {
   dim3 grid_rc(raft::ceildiv(nnz, TPB_X), 1, 1);
   dim3 blk_rc(TPB_X, 1, 1);
 
diff --git a/cpp/include/raft/sparse/selection/knn.cuh b/cpp/include/raft/sparse/selection/knn.cuh
index 3566939bc4..49573a679d 100644
--- a/cpp/include/raft/sparse/selection/knn.cuh
+++ b/cpp/include/raft/sparse/selection/knn.cuh
@@ -31,8 +31,6 @@
 #include <raft/sparse/coo.cuh>
 #include <raft/sparse/csr.cuh>
 #include <raft/sparse/distance/distance.cuh>
-#include <raft/sparse/selection/selection.cuh>
-
 #include <raft/spatial/knn/knn.hpp>
 
 #include <raft/cudart_utils.h>
@@ -339,8 +337,9 @@ class sparse_knn_t {
     if (metric == raft::distance::DistanceType::InnerProduct) ascending = false;
 
     // kernel to slice first (min) k cols and copy into batched merge buffer
-    select_k(batch_dists, batch_indices, batch_rows, batch_cols, out_dists,
-             out_indices, ascending, n_neighbors, handle.get_stream());
+    raft::spatial::knn::select_k(batch_dists, batch_indices, batch_rows,
+                                 batch_cols, out_dists, out_indices, ascending,
+                                 n_neighbors, handle.get_stream());
   }
 
   void compute_distances(csr_batcher_t<value_idx, value_t> &idx_batcher,
diff --git a/cpp/include/raft/sparse/selection/knn_graph.cuh b/cpp/include/raft/sparse/selection/knn_graph.cuh
index 1cdd66f516..3df1c77081 100644
--- a/cpp/include/raft/sparse/selection/knn_graph.cuh
+++ b/cpp/include/raft/sparse/selection/knn_graph.cuh
@@ -88,11 +88,12 @@ void conv_indices(in_t *inds, out_t *out, size_t size, cudaStream_t stream) {
  * @param[in] n number of observations (columns) in X
  * @param[in] metric distance metric to use when constructing neighborhoods
  * @param[out] out output edge list
+ * @param[out] out output edge list
  * @param c
  */
 template <typename value_idx = int, typename value_t = float>
 void knn_graph(const handle_t &handle, const value_t *X, size_t m, size_t n,
-               distance::DistanceType metric,
+               raft::distance::DistanceType metric,
                raft::sparse::COO<value_t, value_idx> &out, int c = 15) {
   int k = build_k(m, c);
 
diff --git a/cpp/include/raft/spatial/knn/ball_cover.hpp b/cpp/include/raft/spatial/knn/ball_cover.hpp
new file mode 100644
index 0000000000..e4b50c77e3
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/ball_cover.hpp
@@ -0,0 +1,162 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cstdint>
+
+#include <raft/linalg/distance_type.h>
+#include <thrust/transform.h>
+#include "ball_cover_common.h"
+#include "detail/ball_cover.cuh"
+#include "detail/ball_cover/common.cuh"
+
+namespace raft {
+namespace spatial {
+namespace knn {
+
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t>
+void rbc_build_index(const raft::handle_t &handle,
+                     BallCoverIndex<value_idx, value_t, value_int> &index) {
+  ASSERT(index.n == 2,
+         "Random ball cover currently only works in 2-dimensions");
+  if (index.metric == raft::distance::DistanceType::Haversine) {
+    detail::rbc_build_index(handle, index, detail::HaversineFunc());
+  } else if (index.metric == raft::distance::DistanceType::L2SqrtExpanded ||
+             index.metric == raft::distance::DistanceType::L2SqrtUnexpanded) {
+    detail::rbc_build_index(handle, index, detail::EuclideanFunc());
+  } else {
+    RAFT_FAIL("Metric not support");
+  }
+
+  index.set_index_trained();
+}
+
+/**
+ * Performs a faster exact knn in metric spaces using the triangle
+ * inequality with a number of landmark points to reduce the
+ * number of distance computations from O(n^2) to O(sqrt(n)). This
+ * performs an all neighbors knn, which can reuse memory when
+ * the index and query are the same array. This function will
+ * build the index and assumes rbc_build_index() has not already
+ * been called.
+ * @tparam value_idx knn index type
+ * @tparam value_t knn distance type
+ * @tparam value_int type for integers, such as number of rows/cols
+ * @param handle raft handle for resource management
+ * @param index ball cover index which has not yet been built
+ * @param k number of nearest neighbors to find
+ * @param perform_post_filtering if this is false, only the closest k landmarks
+ *                               are considered (which will return approximate
+ *                               results).
+ * @param[out] inds output knn indices
+ * @param[out] dists output knn distances
+ * @param weight a weight for overlap between the closest landmark and
+ *               the radius of other landmarks when pruning distances.
+ *               Setting this value below 1 can effectively turn off
+ *               computing distances against many other balls, enabling
+ *               approximate nearest neighbors. Recall can be adjusted
+ *               based on how many relevant balls are ignored. Note that
+ *               many datasets can still have great recall even by only
+ *               looking in the closest landmark.
+ */
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t>
+void rbc_all_knn_query(const raft::handle_t &handle,
+                       BallCoverIndex<value_idx, value_t, value_int> &index,
+                       value_int k, value_idx *inds, value_t *dists,
+                       bool perform_post_filtering = true, float weight = 1.0) {
+  ASSERT(index.n == 2,
+         "Random ball cover currently only works in 2-dimensions");
+  if (index.metric == raft::distance::DistanceType::Haversine) {
+    detail::rbc_all_knn_query(handle, index, k, inds, dists,
+                              detail::HaversineFunc(), perform_post_filtering,
+                              weight);
+  } else if (index.metric == raft::distance::DistanceType::L2SqrtExpanded ||
+             index.metric == raft::distance::DistanceType::L2SqrtUnexpanded) {
+    detail::rbc_all_knn_query(handle, index, k, inds, dists,
+                              detail::EuclideanFunc(), perform_post_filtering,
+                              weight);
+  } else {
+    RAFT_FAIL("Metric not supported");
+  }
+
+  index.set_index_trained();
+}
+
+/**
+ * Performs a faster exact knn in metric spaces using the triangle
+ * inequality with a number of landmark points to reduce the
+ * number of distance computations from O(n^2) to O(sqrt(n)). This
+ * function does not build the index and assumes rbc_build_index() has
+ * already been called. Use this function when the index and
+ * query arrays are different, otherwise use rbc_all_knn_query().
+ * @tparam value_idx index type
+ * @tparam value_t distances type
+ * @tparam value_int integer type for size info
+ * @param handle raft handle for resource management
+ * @param index ball cover index which has not yet been built
+ * @param k number of nearest neighbors to find
+ * @param query the
+ * @param perform_post_filtering if this is false, only the closest k landmarks
+ *                               are considered (which will return approximate
+ *                               results).
+ * @param[out] inds output knn indices
+ * @param[out] dists output knn distances
+ * @param weight a weight for overlap between the closest landmark and
+ *               the radius of other landmarks when pruning distances.
+ *               Setting this value below 1 can effectively turn off
+ *               computing distances against many other balls, enabling
+ *               approximate nearest neighbors. Recall can be adjusted
+ *               based on how many relevant balls are ignored. Note that
+ *               many datasets can still have great recall even by only
+ *               looking in the closest landmark.
+ * @param k
+ * @param inds
+ * @param dists
+ * @param n_samples
+ */
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t>
+void rbc_knn_query(const raft::handle_t &handle,
+                   BallCoverIndex<value_idx, value_t, value_int> &index,
+                   value_int k, const value_t *query, value_int n_query_pts,
+                   value_idx *inds, value_t *dists,
+                   bool perform_post_filtering = true, float weight = 1.0) {
+  ASSERT(index.n == 2,
+         "Random ball cover currently only works in 2-dimensions");
+  if (index.metric == raft::distance::DistanceType::Haversine) {
+    detail::rbc_knn_query(handle, index, k, query, n_query_pts, inds, dists,
+                          detail::HaversineFunc(), perform_post_filtering,
+                          weight);
+  } else if (index.metric == raft::distance::DistanceType::L2SqrtExpanded ||
+             index.metric == raft::distance::DistanceType::L2SqrtUnexpanded) {
+    detail::rbc_knn_query(handle, index, k, query, n_query_pts, inds, dists,
+                          detail::EuclideanFunc(), perform_post_filtering,
+                          weight);
+  } else {
+    RAFT_FAIL("Metric not supported");
+  }
+}
+
+// TODO: implement functions for:
+//  4. rbc_eps_neigh() - given a populated index, perform query against different query array
+//  5. rbc_all_eps_neigh() - populate a BallCoverIndex and query against training data
+
+}  // namespace knn
+}  // namespace spatial
+}  // namespace raft
diff --git a/cpp/include/raft/spatial/knn/ball_cover_common.h b/cpp/include/raft/spatial/knn/ball_cover_common.h
new file mode 100644
index 0000000000..ca614bb0cb
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/ball_cover_common.h
@@ -0,0 +1,99 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <raft/linalg/distance_type.h>
+#include <cstdint>
+#include <raft/handle.hpp>
+#include <rmm/device_uvector.hpp>
+
+namespace raft {
+namespace spatial {
+namespace knn {
+
+/**
+ * Stores raw index data points, sampled landmarks, the 1-nns of index points
+ * to their closest landmarks, and the ball radii of each landmark. This
+ * class is intended to be constructed once and reused across subsequent
+ * queries.
+ * @tparam value_idx
+ * @tparam value_t
+ * @tparam value_int
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t>
+class BallCoverIndex {
+ public:
+  explicit BallCoverIndex(const raft::handle_t &handle_, const value_t *X_,
+                          value_int m_, value_int n_,
+                          raft::distance::DistanceType metric_)
+    : handle(handle_),
+      X(X_),
+      m(m_),
+      n(n_),
+      metric(metric_),
+      /**
+      * the sqrt() here makes the sqrt(m)^2 a linear-time lower bound
+      *
+      * Total memory footprint of index: (2 * sqrt(m)) + (n * sqrt(m)) + (2 * m)
+      */
+      n_landmarks(sqrt(m_)),
+      R_indptr(sqrt(m_) + 1, handle.get_stream()),
+      R_1nn_cols(m_, handle.get_stream()),
+      R_1nn_dists(m_, handle.get_stream()),
+      R(sqrt(m_) * n_, handle.get_stream()),
+      R_radius(sqrt(m_), handle.get_stream()),
+      index_trained(false) {}
+
+  value_idx *get_R_indptr() { return R_indptr.data(); }
+  value_idx *get_R_1nn_cols() { return R_1nn_cols.data(); }
+  value_t *get_R_1nn_dists() { return R_1nn_dists.data(); }
+  value_t *get_R_radius() { return R_radius.data(); }
+  value_t *get_R() { return R.data(); }
+  const value_t *get_X() { return X; }
+
+  bool is_index_trained() const { return index_trained; };
+
+  // This should only be set by internal functions
+  void set_index_trained() { index_trained = true; }
+
+  const raft::handle_t &handle;
+
+  const value_int m;
+  const value_int n;
+  const value_int n_landmarks;
+
+  const value_t *X;
+
+  raft::distance::DistanceType metric;
+
+ private:
+  // CSR storing the neighborhoods for each data point
+  rmm::device_uvector<value_idx> R_indptr;
+  rmm::device_uvector<value_idx> R_1nn_cols;
+  rmm::device_uvector<value_t> R_1nn_dists;
+
+  rmm::device_uvector<value_t> R_radius;
+
+  rmm::device_uvector<value_t> R;
+
+ protected:
+  bool index_trained;
+};
+}  // namespace knn
+}  // namespace spatial
+}  // namespace raft
diff --git a/cpp/include/raft/spatial/knn/detail/ann_quantized_faiss.cuh b/cpp/include/raft/spatial/knn/detail/ann_quantized_faiss.cuh
index 77ad4afe96..0e91b5225d 100644
--- a/cpp/include/raft/spatial/knn/detail/ann_quantized_faiss.cuh
+++ b/cpp/include/raft/spatial/knn/detail/ann_quantized_faiss.cuh
@@ -17,12 +17,14 @@
 #pragma once
 
 #include "../ann_common.h"
+#include "knn_brute_force_faiss.cuh"
 
 #include "common_faiss.h"
 #include "processing.hpp"
 
 #include <raft/cudart_utils.h>
 #include <raft/cuda_utils.cuh>
+#include "processing.hpp"
 
 #include <label/classlabels.cuh>
 #include <raft/distance/distance.cuh>
diff --git a/cpp/include/raft/spatial/knn/detail/ball_cover.cuh b/cpp/include/raft/spatial/knn/detail/ball_cover.cuh
new file mode 100644
index 0000000000..46a97400e2
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/detail/ball_cover.cuh
@@ -0,0 +1,351 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <raft/handle.hpp>
+
+#include "../ball_cover_common.h"
+#include "ball_cover/common.cuh"
+#include "ball_cover/registers.cuh"
+#include "block_select_faiss.cuh"
+#include "haversine_distance.cuh"
+#include "knn_brute_force_faiss.cuh"
+#include "selection_faiss.cuh"
+
+#include <limits.h>
+#include <cstdint>
+
+#include <raft/cuda_utils.cuh>
+
+#include <raft/matrix/matrix.cuh>
+#include <raft/random/rng.cuh>
+#include <raft/sparse/convert/csr.cuh>
+#include <raft/sparse/distance/operators.cuh>
+
+#include <rmm/device_uvector.hpp>
+#include <rmm/exec_policy.hpp>
+
+#include <faiss/utils/Heap.h>
+#include <faiss/gpu/utils/Limits.cuh>
+#include <faiss/gpu/utils/Select.cuh>
+
+#include <thrust/functional.h>
+#include <thrust/reduce.h>
+#include <thrust/sequence.h>
+
+namespace raft {
+namespace spatial {
+namespace knn {
+namespace detail {
+
+/**
+ * Given a set of points in row-major order which are to be
+ * used as a set of index points, uniformly samples a subset
+ * of points to be used as landmarks.
+ * @tparam value_idx
+ * @tparam value_t
+ * @param handle
+ * @param index
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t>
+void sample_landmarks(const raft::handle_t &handle,
+                      BallCoverIndex<value_idx, value_t, value_int> &index) {
+  rmm::device_uvector<value_idx> R_1nn_cols2(index.n_landmarks,
+                                             handle.get_stream());
+  rmm::device_uvector<value_t> R_1nn_ones(index.m, handle.get_stream());
+  rmm::device_uvector<value_idx> R_indices(index.n_landmarks,
+                                           handle.get_stream());
+
+  thrust::sequence(handle.get_thrust_policy(), index.get_R_1nn_cols(),
+                   index.get_R_1nn_cols() + index.m, (value_idx)0);
+
+  thrust::fill(handle.get_thrust_policy(), R_1nn_ones.data(),
+               R_1nn_ones.data() + R_1nn_ones.size(), 1.0);
+
+  /**
+ * 1. Randomly sample sqrt(n) points from X
+ */
+  auto rng = raft::random::Rng(12345);
+  rng.sampleWithoutReplacement(handle, R_indices.data(), R_1nn_cols2.data(),
+                               index.get_R_1nn_cols(), R_1nn_ones.data(),
+                               (value_idx)index.n_landmarks, (value_idx)index.m,
+                               handle.get_stream());
+
+  raft::matrix::copyRows<value_t, value_idx, size_t>(
+    index.get_X(), index.m, index.n, index.get_R(), R_1nn_cols2.data(),
+    index.n_landmarks, handle.get_stream(), true);
+}
+
+/**
+ * Constructs a 1-nn index mapping each landmark to their closest points.
+ * @tparam value_idx
+ * @tparam value_t
+ * @param handle
+ * @param R_knn_inds_ptr
+ * @param R_knn_dists_ptr
+ * @param k
+ * @param index
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t>
+void construct_landmark_1nn(
+  const raft::handle_t &handle, const value_idx *R_knn_inds_ptr,
+  const value_t *R_knn_dists_ptr, value_int k,
+  BallCoverIndex<value_idx, value_t, value_int> &index) {
+  rmm::device_uvector<value_idx> R_1nn_inds(index.m, handle.get_stream());
+
+  value_idx *R_1nn_inds_ptr = R_1nn_inds.data();
+  value_t *R_1nn_dists_ptr = index.get_R_1nn_dists();
+
+  auto idxs = thrust::make_counting_iterator<value_idx>(0);
+  thrust::for_each(handle.get_thrust_policy(), idxs, idxs + index.m,
+                   [=] __device__(value_idx i) {
+                     R_1nn_inds_ptr[i] = R_knn_inds_ptr[i * k];
+                     R_1nn_dists_ptr[i] = R_knn_dists_ptr[i * k];
+                   });
+
+  auto keys = thrust::make_zip_iterator(
+    thrust::make_tuple(R_1nn_inds.data(), index.get_R_1nn_dists()));
+
+  // group neighborhoods for each reference landmark and sort each group by distance
+  thrust::sort_by_key(handle.get_thrust_policy(), keys, keys + index.m,
+                      index.get_R_1nn_cols(), NNComp());
+
+  // convert to CSR for fast lookup
+  raft::sparse::convert::sorted_coo_to_csr(
+    R_1nn_inds.data(), index.m, index.get_R_indptr(), index.n_landmarks + 1,
+    handle.get_stream());
+}
+
+/**
+ * Computes the k closest landmarks to a set of query points.
+ * @tparam value_idx
+ * @tparam value_t
+ * @tparam value_int
+ * @param handle
+ * @param index
+ * @param query_pts
+ * @param n_query_pts
+ * @param k
+ * @param R_knn_inds
+ * @param R_knn_dists
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t>
+void k_closest_landmarks(const raft::handle_t &handle,
+                         BallCoverIndex<value_idx, value_t, value_int> &index,
+                         const value_t *query_pts, value_int n_query_pts,
+                         value_int k, value_idx *R_knn_inds,
+                         value_t *R_knn_dists) {
+  std::vector<value_t *> input = {index.get_R()};
+  std::vector<std::uint32_t> sizes = {index.n_landmarks};
+
+  brute_force_knn_impl<std::uint32_t, std::int64_t>(
+    input, sizes, index.n, const_cast<value_t *>(query_pts), n_query_pts,
+    R_knn_inds, R_knn_dists, k, handle.get_stream(), nullptr, 0, true, true,
+    nullptr, index.metric);
+}
+
+/**
+ * Uses the sorted data points in the 1-nn landmark index to compute
+ * an array of radii for each landmark.
+ * @tparam value_idx
+ * @tparam value_t
+ * @param handle
+ * @param index
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t>
+void compute_landmark_radii(
+  const raft::handle_t &handle,
+  BallCoverIndex<value_idx, value_t, value_int> &index) {
+  auto entries = thrust::make_counting_iterator<value_idx>(0);
+
+  const value_idx *R_indptr_ptr = index.get_R_indptr();
+  const value_t *R_1nn_dists_ptr = index.get_R_1nn_dists();
+  value_t *R_radius_ptr = index.get_R_radius();
+  thrust::for_each(handle.get_thrust_policy(), entries,
+                   entries + index.n_landmarks,
+                   [=] __device__(value_idx input) {
+                     value_idx last_row_idx = R_indptr_ptr[input + 1] - 1;
+                     R_radius_ptr[input] = R_1nn_dists_ptr[last_row_idx];
+                   });
+}
+
+/**
+ * 4. Perform k-select over original KNN, using L_r to filter distances
+ *
+ * a. Map 1 row to each warp/block
+ * b. Add closest k R points to heap
+ * c. Iterate through batches of R, having each thread in the warp load a set
+ * of distances y from R (only if d(q, r) < 3 * distance to closest r) and
+ * marking the distance to be computed between x, y only
+ * if knn[k].distance >= d(x_i, R_k) + d(R_k, y)
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t, typename dist_func>
+void perform_rbc_query(const raft::handle_t &handle,
+                       BallCoverIndex<value_idx, value_t, value_int> &index,
+                       const value_t *query, value_int n_query_pts,
+                       std::uint32_t k, const value_idx *R_knn_inds,
+                       const value_t *R_knn_dists, dist_func dfunc,
+                       value_idx *inds, value_t *dists,
+                       value_int *dists_counter, value_int *post_dists_counter,
+                       float weight = 1.0, bool perform_post_filtering = true) {
+  // Compute nearest k for each neighborhood in each closest R
+  rbc_low_dim_pass_one(handle, index, query, n_query_pts, k, R_knn_inds,
+                       R_knn_dists, dfunc, inds, dists, weight, dists_counter);
+
+  if (perform_post_filtering) {
+    rbc_low_dim_pass_two(handle, index, query, n_query_pts, k, R_knn_inds,
+                         R_knn_dists, dfunc, inds, dists, weight,
+                         post_dists_counter);
+  }
+}
+
+/**
+ * Similar to a ball tree, the random ball cover algorithm
+ * uses the triangle inequality to prune distance computations
+ * in any metric space with a guarantee of sqrt(n) * c^{3/2}
+ * where `c` is an expansion constant based on the distance
+ * metric.
+ *
+ * This function variant performs an all nearest neighbors
+ * query which is useful for algorithms that need to perform
+ * A * A.T.
+ */
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t, typename distance_func>
+void rbc_build_index(const raft::handle_t &handle,
+                     BallCoverIndex<value_idx, value_t, value_int> &index,
+                     distance_func dfunc) {
+  ASSERT(index.n == 2,
+         "only 2d vectors are supported in current implementation");
+  ASSERT(!index.is_index_trained(), "index cannot be previously trained");
+
+  rmm::device_uvector<value_idx> R_knn_inds(index.m, handle.get_stream());
+  rmm::device_uvector<value_t> R_knn_dists(index.m, handle.get_stream());
+
+  /**
+   * 1. Randomly sample sqrt(n) points from X
+   */
+  sample_landmarks<value_idx, value_t>(handle, index);
+
+  /**
+   * 2. Perform knn = bfknn(X, R, k)
+   */
+  value_int k = 1;
+  k_closest_landmarks(handle, index, index.get_X(), index.m, k,
+                      R_knn_inds.data(), R_knn_dists.data());
+
+  /**
+   * 3. Create L_r = knn[:,0].T (CSR)
+   *
+   * Slice closest neighboring R
+   * Secondary sort by (R_knn_inds, R_knn_dists)
+   */
+  construct_landmark_1nn(handle, R_knn_inds.data(), R_knn_dists.data(), k,
+                         index);
+
+  /**
+   * Compute radius of each R for filtering: p(q, r) <= p(q, q_r) + radius(r)
+   * (need to take the
+   */
+  compute_landmark_radii(handle, index);
+}
+
+/**
+ * Performs an all neighbors knn query (e.g. index == query)
+ */
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t, typename distance_func>
+void rbc_all_knn_query(const raft::handle_t &handle,
+                       BallCoverIndex<value_idx, value_t, value_int> &index,
+                       value_int k, value_idx *inds, value_t *dists,
+                       distance_func dfunc,
+                       // approximate nn options
+                       bool perform_post_filtering = true, float weight = 1.0) {
+  ASSERT(index.n == 2,
+         "only 2d vectors are supported in current implementation");
+  ASSERT(index.n_landmarks >= k, "number of landmark samples must be >= k");
+  ASSERT(!index.is_index_trained(), "index cannot be previously trained");
+
+  rmm::device_uvector<value_idx> R_knn_inds(k * index.m, handle.get_stream());
+  rmm::device_uvector<value_t> R_knn_dists(k * index.m, handle.get_stream());
+
+  // For debugging / verification. Remove before releasing
+  rmm::device_uvector<value_int> dists_counter(index.m, handle.get_stream());
+  rmm::device_uvector<value_int> post_dists_counter(index.m,
+                                                    handle.get_stream());
+
+  sample_landmarks<value_idx, value_t>(handle, index);
+
+  k_closest_landmarks(handle, index, index.get_X(), index.m, k,
+                      R_knn_inds.data(), R_knn_dists.data());
+
+  construct_landmark_1nn(handle, R_knn_inds.data(), R_knn_dists.data(), k,
+                         index);
+
+  compute_landmark_radii(handle, index);
+
+  perform_rbc_query(handle, index, index.get_X(), index.m, k, R_knn_inds.data(),
+                    R_knn_dists.data(), dfunc, inds, dists,
+                    dists_counter.data(), post_dists_counter.data(), weight,
+                    perform_post_filtering);
+}
+
+/**
+ * Performs a knn query against an index. This assumes the index has
+ * already been built.
+ */
+template <typename value_idx = std::int64_t, typename value_t,
+          typename value_int = std::uint32_t, typename distance_func>
+void rbc_knn_query(const raft::handle_t &handle,
+                   BallCoverIndex<value_idx, value_t, value_int> &index,
+                   value_int k, const value_t *query, value_int n_query_pts,
+                   value_idx *inds, value_t *dists, distance_func dfunc,
+                   // approximate nn options
+                   bool perform_post_filtering = true, float weight = 1.0) {
+  ASSERT(index.n == 2,
+         "only 2d vectors are supported in current implementation");
+  ASSERT(index.n_landmarks >= k, "number of landmark samples must be >= k");
+  ASSERT(index.is_index_trained(), "index must be previously trained");
+
+  rmm::device_uvector<value_idx> R_knn_inds(k * index.m, handle.get_stream());
+  rmm::device_uvector<value_t> R_knn_dists(k * index.m, handle.get_stream());
+
+  k_closest_landmarks(handle, index, query, n_query_pts, k, R_knn_inds.data(),
+                      R_knn_dists.data());
+
+  // For debugging / verification. Remove before releasing
+  rmm::device_uvector<value_int> dists_counter(index.m, handle.get_stream());
+  rmm::device_uvector<value_int> post_dists_counter(index.m,
+                                                    handle.get_stream());
+  thrust::fill(handle.get_thrust_policy(), post_dists_counter.data(),
+               post_dists_counter.data() + index.m, 0);
+
+  perform_rbc_query(handle, index, query, n_query_pts, k, R_knn_inds.data(),
+                    R_knn_dists.data(), dfunc, inds, dists,
+                    dists_counter.data(), post_dists_counter.data(), weight,
+                    perform_post_filtering);
+}
+
+};  // namespace detail
+};  // namespace knn
+};  // namespace spatial
+};  // namespace raft
diff --git a/cpp/include/raft/spatial/knn/detail/ball_cover/common.cuh b/cpp/include/raft/spatial/knn/detail/ball_cover/common.cuh
new file mode 100644
index 0000000000..c6cb679408
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/detail/ball_cover/common.cuh
@@ -0,0 +1,91 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <thrust/functional.h>
+#include <cstdint>
+#include "../haversine_distance.cuh"
+
+namespace raft {
+namespace spatial {
+namespace knn {
+namespace detail {
+
+struct NNComp {
+  template <typename one, typename two>
+  __host__ __device__ bool operator()(const one &t1, const two &t2) {
+    // sort first by each sample's reference landmark,
+    if (thrust::get<0>(t1) < thrust::get<0>(t2)) return true;
+    if (thrust::get<0>(t1) > thrust::get<0>(t2)) return false;
+
+    // then by closest neighbor,
+    return thrust::get<1>(t1) < thrust::get<1>(t2);
+  }
+};
+
+struct HaversineFunc {
+  template <typename value_t, typename value_int = std::uint32_t>
+  __device__ __host__ __forceinline__ value_t
+  operator()(const value_t *a, const value_t *b, const value_int n_dims) {
+    return raft::spatial::knn::detail::compute_haversine(a[0], b[0], a[1],
+                                                         b[1]);
+  }
+};
+
+struct EuclideanFunc {
+  template <typename value_t, typename value_int = std::uint32_t>
+  __device__ __host__ __forceinline__ value_t
+  operator()(const value_t *a, const value_t *b, const value_int n_dims) {
+    value_t sum_sq = 0;
+    for (value_int i = 0; i < n_dims; ++i) {
+      value_t diff = a[i] - b[i];
+      sum_sq += diff * diff;
+    }
+
+    return sqrt(sum_sq);
+  }
+};
+
+/**
+ * Zeros the bit at location h in a one-hot encoded 32-bit int array
+ */
+__device__ inline void _zero_bit(std::uint32_t *arr, std::uint32_t h) {
+  int bit = h % 32;
+  int idx = h / 32;
+
+  std::uint32_t assumed;
+  std::uint32_t old = arr[idx];
+  do {
+    assumed = old;
+    old = atomicCAS(arr + idx, assumed, assumed & ~(1 << bit));
+  } while (assumed != old);
+}
+
+/**
+ * Returns whether or not bit at location h is nonzero in a one-hot
+ * encoded 32-bit in array.
+ */
+__device__ inline bool _get_val(std::uint32_t *arr, std::uint32_t h) {
+  int bit = h % 32;
+  int idx = h / 32;
+  return (arr[idx] & (1 << bit)) > 0;
+}
+
+};  // namespace detail
+};  // namespace knn
+};  // namespace spatial
+};  // namespace raft
\ No newline at end of file
diff --git a/cpp/include/raft/spatial/knn/detail/ball_cover/registers.cuh b/cpp/include/raft/spatial/knn/detail/ball_cover/registers.cuh
new file mode 100644
index 0000000000..4a476274dd
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/detail/ball_cover/registers.cuh
@@ -0,0 +1,537 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include "common.cuh"
+
+#include "../../ball_cover_common.h"
+#include "../block_select_faiss.cuh"
+#include "../haversine_distance.cuh"
+#include "../selection_faiss.cuh"
+
+#include <limits.h>
+#include <cstdint>
+
+#include <raft/cuda_utils.cuh>
+
+#include <faiss/utils/Heap.h>
+#include <faiss/gpu/utils/Limits.cuh>
+#include <faiss/gpu/utils/Select.cuh>
+
+namespace raft {
+namespace spatial {
+namespace knn {
+namespace detail {
+
+/**
+ * To find exact neighbors, we perform a post-processing stage
+ * that filters out those points which might have neighbors outside
+ * of their k closest landmarks. This is usually a very small portion
+ * of the total points.
+ * @tparam value_idx
+ * @tparam value_t
+ * @tparam value_int
+ * @tparam tpb
+ * @param X
+ * @param n_cols
+ * @param R_knn_inds
+ * @param R_knn_dists
+ * @param R_radius
+ * @param landmarks
+ * @param n_landmarks
+ * @param bitset_size
+ * @param k
+ * @param output
+ * @param weight
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t, int tpb = 32,
+          typename distance_func>
+__global__ void perform_post_filter_registers(
+  const value_t *X, value_int n_cols, const value_idx *R_knn_inds,
+  const value_t *R_knn_dists, const value_t *R_radius, const value_t *landmarks,
+  int n_landmarks, value_int bitset_size, value_int k, distance_func dfunc,
+  std::uint32_t *output, float weight = 1.0) {
+  // allocate array of size n_landmarks / 32 ints
+  extern __shared__ std::uint32_t shared_mem[];
+
+  // Start with all bits on
+  for (value_int i = threadIdx.x; i < bitset_size; i += tpb) {
+    shared_mem[i] = 0xffffffff;
+  }
+
+  __syncthreads();
+
+  // TODO: Would it be faster to use L1 for this?
+  value_t local_x_ptr[2];
+  for (value_int j = 0; j < n_cols; ++j) {
+    local_x_ptr[j] = X[n_cols * blockIdx.x + j];
+  }
+
+  value_t closest_R_dist = R_knn_dists[blockIdx.x * k + (k - 1)];
+
+  // zero out bits for closest k landmarks
+  for (value_int j = threadIdx.x; j < k; j += tpb) {
+    _zero_bit(shared_mem, (std::uint32_t)R_knn_inds[blockIdx.x * k + j]);
+  }
+
+  __syncthreads();
+
+  // Discard any landmarks where p(q, r) > p(q, r_q) + radius(r)
+  // That is, the distance between the current point and the current
+  // landmark is > the distance between the current point and
+  // its closest landmark + the radius of the current landmark.
+  for (value_int l = threadIdx.x; l < n_landmarks; l += tpb) {
+    // compute p(q, r)
+    value_t dist = dfunc(local_x_ptr, landmarks + (n_cols * l), n_cols);
+    if (dist > weight * (closest_R_dist + R_radius[l]) ||
+        dist > 3 * closest_R_dist) {
+      _zero_bit(shared_mem, l);
+    }
+  }
+
+  __syncthreads();
+
+  /**
+   * Output bitset
+   */
+  for (value_int l = threadIdx.x; l < bitset_size; l += tpb) {
+    output[blockIdx.x * bitset_size + l] = shared_mem[l];
+  }
+}
+
+/**
+ * @tparam value_idx
+ * @tparam value_t
+ * @tparam value_int
+ * @tparam bitset_type
+ * @tparam warp_q number of registers to use per warp
+ * @tparam thread_q number of registers to use within each thread
+ * @tparam tpb number of threads per block
+ * @param X
+ * @param n_cols
+ * @param bitset
+ * @param bitset_size
+ * @param R_knn_dists
+ * @param R_indptr
+ * @param R_1nn_inds
+ * @param R_1nn_dists
+ * @param knn_inds
+ * @param knn_dists
+ * @param n_landmarks
+ * @param k
+ * @param dist_counter
+ */
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t,
+          typename bitset_type = std::uint32_t, typename dist_func,
+          int warp_q = 32, int thread_q = 2, int tpb = 128, int col_q = 2>
+__global__ void compute_final_dists_registers(
+  const value_t *X_index, const value_t *X, const value_int n_cols,
+  bitset_type *bitset, value_int bitset_size, const value_t *R_knn_dists,
+  const value_idx *R_indptr, const value_idx *R_1nn_inds,
+  const value_t *R_1nn_dists, value_idx *knn_inds, value_t *knn_dists,
+  value_int n_landmarks, value_int k, dist_func dfunc,
+  value_int *dist_counter) {
+  static constexpr int kNumWarps = tpb / faiss::gpu::kWarpSize;
+
+  __shared__ value_t shared_memK[kNumWarps * warp_q];
+  __shared__ faiss::gpu::KeyValuePair<value_t, value_idx>
+    shared_memV[kNumWarps * warp_q];
+
+  const value_t *x_ptr = X + (n_cols * blockIdx.x);
+  value_t local_x_ptr[col_q];
+  for (value_int j = 0; j < n_cols; ++j) {
+    local_x_ptr[j] = x_ptr[j];
+  }
+
+  faiss::gpu::KeyValueBlockSelect<value_t, value_idx, false,
+                                  faiss::gpu::Comparator<value_t>, warp_q,
+                                  thread_q, tpb>
+    heap(faiss::gpu::Limits<value_t>::getMax(),
+         faiss::gpu::Limits<value_t>::getMax(), -1, shared_memK, shared_memV,
+         k);
+
+  const value_int n_k = faiss::gpu::utils::roundDown(k, faiss::gpu::kWarpSize);
+  value_int i = threadIdx.x;
+  for (; i < n_k; i += tpb) {
+    value_idx ind = knn_inds[blockIdx.x * k + i];
+    heap.add(knn_dists[blockIdx.x * k + i], R_knn_dists[ind * k], ind);
+  }
+
+  if (i < k) {
+    value_idx ind = knn_inds[blockIdx.x * k + i];
+    heap.addThreadQ(knn_dists[blockIdx.x * k + i], R_knn_dists[ind * k], ind);
+  }
+
+  heap.checkThreadQ();
+
+  for (value_int cur_R_ind = 0; cur_R_ind < n_landmarks; ++cur_R_ind) {
+    // if cur R overlaps cur point's closest R, it could be a
+    // candidate
+    if (_get_val(bitset + (blockIdx.x * bitset_size), cur_R_ind)) {
+      value_idx R_start_offset = R_indptr[cur_R_ind];
+      value_idx R_stop_offset = R_indptr[cur_R_ind + 1];
+      value_idx R_size = R_stop_offset - R_start_offset;
+
+      // Loop through R's neighborhood in parallel
+
+      // Round R_size to the nearest warp threads so they can
+      // all be computing in parallel.
+
+      const value_int limit =
+        faiss::gpu::utils::roundDown(R_size, faiss::gpu::kWarpSize);
+
+      i = threadIdx.x;
+      for (; i < limit; i += tpb) {
+        value_idx cur_candidate_ind = R_1nn_inds[R_start_offset + i];
+        value_t cur_candidate_dist = R_1nn_dists[R_start_offset + i];
+        value_t z = heap.warpKTopRDist == 0.00
+                      ? 0.0
+                      : (abs(heap.warpKTop - heap.warpKTopRDist) *
+                           abs(heap.warpKTopRDist - cur_candidate_dist) -
+                         heap.warpKTop * cur_candidate_dist) /
+                          heap.warpKTopRDist;
+        z = isnan(z) ? 0.0 : z;
+        // If lower bound on distance could possibly be in
+        // the closest k neighbors, compute it and add to k-select
+        value_t dist = std::numeric_limits<value_t>::max();
+        if (z <= heap.warpKTop) {
+          const value_t *y_ptr = X_index + (n_cols * cur_candidate_ind);
+          value_t local_y_ptr[col_q];
+          for (value_int j = 0; j < n_cols; ++j) {
+            local_y_ptr[j] = y_ptr[j];
+          }
+
+          dist = dfunc(local_x_ptr, local_y_ptr, n_cols);
+        }
+
+        heap.add(dist, cur_candidate_dist, cur_candidate_ind);
+      }
+
+      // second round guarantees to be only a single warp.
+      if (i < R_size) {
+        value_idx cur_candidate_ind = R_1nn_inds[R_start_offset + i];
+        value_t cur_candidate_dist = R_1nn_dists[R_start_offset + i];
+
+        value_t z = heap.warpKTopRDist == 0.00
+                      ? 0.0
+                      : (abs(heap.warpKTop - heap.warpKTopRDist) *
+                           abs(heap.warpKTopRDist - cur_candidate_dist) -
+                         heap.warpKTop * cur_candidate_dist) /
+                          heap.warpKTopRDist;
+
+        z = isnan(z) ? 0.0 : z;
+        // If lower bound on distance could possibly be in
+        // the closest k neighbors, compute it and add to k-select
+        value_t dist = std::numeric_limits<value_t>::max();
+        if (z <= heap.warpKTop) {
+          const value_t *y_ptr = X_index + (n_cols * cur_candidate_ind);
+          value_t local_y_ptr[col_q];
+          for (value_int j = 0; j < n_cols; ++j) {
+            local_y_ptr[j] = y_ptr[j];
+          }
+          dist = dfunc(local_x_ptr, local_y_ptr, n_cols);
+        }
+        heap.addThreadQ(dist, cur_candidate_dist, cur_candidate_ind);
+      }
+      heap.checkThreadQ();
+    }
+  }
+
+  heap.reduce();
+
+  for (value_int i = threadIdx.x; i < k; i += tpb) {
+    knn_dists[blockIdx.x * k + i] = shared_memK[i];
+    knn_inds[blockIdx.x * k + i] = shared_memV[i].value;
+  }
+}
+
+/**
+ * Random ball cover kernel for n_dims == 2
+ * @tparam value_idx
+ * @tparam value_t
+ * @tparam warp_q
+ * @tparam thread_q
+ * @tparam tpb
+ * @tparam value_idx
+ * @tparam value_t
+ * @param R_knn_inds
+ * @param R_knn_dists
+ * @param m
+ * @param k
+ * @param R_indptr
+ * @param R_1nn_cols
+ * @param R_1nn_dists
+ */
+template <typename value_idx = std::int64_t, typename value_t, int warp_q = 32,
+          int thread_q = 2, int tpb = 128, int col_q = 2,
+          typename value_int = std::uint32_t, typename distance_func>
+__global__ void block_rbc_kernel_registers(
+  const value_t *X_index, const value_t *X,
+  value_int n_cols,  // n_cols should be 2 or 3 dims
+  const value_idx *R_knn_inds, const value_t *R_knn_dists, value_int m,
+  value_int k, const value_idx *R_indptr, const value_idx *R_1nn_cols,
+  const value_t *R_1nn_dists, value_idx *out_inds, value_t *out_dists,
+  value_int *dist_counter, value_t *R_radius, distance_func dfunc,
+  float weight = 1.0) {
+  static constexpr value_int kNumWarps = tpb / faiss::gpu::kWarpSize;
+
+  __shared__ value_t shared_memK[kNumWarps * warp_q];
+  __shared__ faiss::gpu::KeyValuePair<value_t, value_idx>
+    shared_memV[kNumWarps * warp_q];
+
+  // TODO: Separate kernels for different widths:
+  // 1. Very small (between 3 and 32) just use registers for columns of "blockIdx.x"
+  // 2. Can fit comfortably in shared memory (32 to a few thousand?)
+  // 3. Load each time individually.
+  const value_t *x_ptr = X + (n_cols * blockIdx.x);
+
+  // Use registers only for 2d or 3d
+  value_t local_x_ptr[col_q];
+  for (value_int i = 0; i < n_cols; ++i) {
+    local_x_ptr[i] = x_ptr[i];
+  }
+
+  // Each warp works on 1 R
+  faiss::gpu::KeyValueBlockSelect<value_t, value_idx, false,
+                                  faiss::gpu::Comparator<value_t>, warp_q,
+                                  thread_q, tpb>
+    heap(faiss::gpu::Limits<value_t>::getMax(),
+         faiss::gpu::Limits<value_t>::getMax(), -1, shared_memK, shared_memV,
+         k);
+
+  value_t min_R_dist = R_knn_dists[blockIdx.x * k + (k - 1)];
+
+  value_int n_dists_computed = 0;
+
+  /**
+   * First add distances for k closest neighbors of R
+   * to the heap
+   */
+  // Start iterating through elements of each set from closest R elements,
+  // determining if the distance could even potentially be in the heap.
+  for (value_int cur_k = 0; cur_k < k; ++cur_k) {
+    // index and distance to current blockIdx.x's closest landmark
+    value_t cur_R_dist = R_knn_dists[blockIdx.x * k + cur_k];
+    value_idx cur_R_ind = R_knn_inds[blockIdx.x * k + cur_k];
+
+    // Equation (2) in Cayton's paper- prune out R's which are > 3 * p(q, r_q)
+    if (cur_R_dist > weight * (min_R_dist + R_radius[cur_R_ind])) continue;
+    if (cur_R_dist > 3 * min_R_dist) return;
+
+    // The whole warp should iterate through the elements in the current R
+    value_idx R_start_offset = R_indptr[cur_R_ind];
+    value_idx R_stop_offset = R_indptr[cur_R_ind + 1];
+
+    value_idx R_size = R_stop_offset - R_start_offset;
+
+    value_int limit =
+      faiss::gpu::utils::roundDown(R_size, faiss::gpu::kWarpSize);
+    value_int i = threadIdx.x;
+    for (; i < limit; i += tpb) {
+      // Index and distance of current candidate's nearest landmark
+      value_idx cur_candidate_ind = R_1nn_cols[R_start_offset + i];
+      value_t cur_candidate_dist = R_1nn_dists[R_start_offset + i];
+
+      // Take 2 landmarks l_1 and l_2 where l_1 is the furthest point in the heap
+      // and l_2 is the current landmark R. s is the current data point and
+      // t is the new candidate data point. We know that:
+      // d(s, t) cannot possibly be any smaller than | d(s, l_1) - d(l_1, l_2) | * | d(l_1, l_2) - d(l_2, t) | - d(s, l_1) * d(l_2, t)
+
+      // Therefore, if d(s, t) >= d(s, l_1) from the computation above, we know that the distance to the candidate point
+      // cannot possibly be in the nearest neighbors. However, if d(s, t) < d(s, l_1) then we should compute the
+      // distance because it's possible it could be smaller.
+      //
+      value_t z = heap.warpKTopRDist == 0.00
+                    ? 0.0
+                    : (abs(heap.warpKTop - heap.warpKTopRDist) *
+                         abs(heap.warpKTopRDist - cur_candidate_dist) -
+                       heap.warpKTop * cur_candidate_dist) /
+                        heap.warpKTopRDist;
+
+      z = isnan(z) ? 0.0 : z;
+      value_t dist = std::numeric_limits<value_t>::max();
+      if (i < k || z <= heap.warpKTop) {
+        const value_t *y_ptr = X_index + (n_cols * cur_candidate_ind);
+        value_t local_y_ptr[col_q];
+        for (value_int j = 0; j < n_cols; ++j) {
+          local_y_ptr[j] = y_ptr[j];
+        }
+        dist = dfunc(local_x_ptr, local_y_ptr, n_cols);
+        ++n_dists_computed;
+      }
+
+      heap.add(dist, cur_candidate_dist, cur_candidate_ind);
+    }
+
+    if (i < R_size) {
+      value_idx cur_candidate_ind = R_1nn_cols[R_start_offset + i];
+      value_t cur_candidate_dist = R_1nn_dists[R_start_offset + i];
+      value_t z = heap.warpKTopRDist == 0.0
+                    ? 0.0
+                    : (abs(heap.warpKTop - heap.warpKTopRDist) *
+                         abs(heap.warpKTopRDist - cur_candidate_dist) -
+                       heap.warpKTop * cur_candidate_dist) /
+                        heap.warpKTopRDist;
+
+      z = isnan(z) ? 0.0 : z;
+      value_t dist = std::numeric_limits<value_t>::max();
+      if (i < k || z <= heap.warpKTop) {
+        const value_t *y_ptr = X_index + (n_cols * cur_candidate_ind);
+        value_t local_y_ptr[col_q];
+        for (value_int j = 0; j < n_cols; ++j) {
+          local_y_ptr[j] = y_ptr[j];
+        }
+        dist = dfunc(local_x_ptr, local_y_ptr, n_cols);
+        ++n_dists_computed;
+      }
+
+      heap.addThreadQ(dist, cur_candidate_dist, cur_candidate_ind);
+    }
+
+    heap.checkThreadQ();
+  }
+
+  heap.reduce();
+
+  for (int i = threadIdx.x; i < k; i += tpb) {
+    out_dists[blockIdx.x * k + i] = shared_memK[i];
+    out_inds[blockIdx.x * k + i] = shared_memV[i].value;
+  }
+}
+
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t, typename dist_func>
+void rbc_low_dim_pass_one(const raft::handle_t &handle,
+                          BallCoverIndex<value_idx, value_t, value_int> &index,
+                          const value_t *query, const value_int n_query_rows,
+                          value_int k, const value_idx *R_knn_inds,
+                          const value_t *R_knn_dists, dist_func dfunc,
+                          value_idx *inds, value_t *dists, float weight,
+                          value_int *dists_counter) {
+  if (k <= 32)
+    block_rbc_kernel_registers<value_idx, value_t, 32, 2, 128, 2, value_int>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+
+  else if (k <= 64)
+    block_rbc_kernel_registers<value_idx, value_t, 64, 3, 128, 2, value_int>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+  else if (k <= 128)
+    block_rbc_kernel_registers<value_idx, value_t, 128, 3, 128, 2, value_int>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+
+  else if (k <= 256)
+    block_rbc_kernel_registers<value_idx, value_t, 256, 4, 128, 2, value_int>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+
+  else if (k <= 512)
+    block_rbc_kernel_registers<value_idx, value_t, 512, 8, 64, 2, value_int>
+      <<<n_query_rows, 64, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+
+  else if (k <= 1024)
+    block_rbc_kernel_registers<value_idx, value_t, 1024, 8, 64, 2, value_int>
+      <<<n_query_rows, 64, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, R_knn_inds, R_knn_dists, index.m, k,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, dists_counter, index.get_R_radius(), dfunc, weight);
+}
+
+template <typename value_idx, typename value_t,
+          typename value_int = std::uint32_t, typename dist_func>
+void rbc_low_dim_pass_two(const raft::handle_t &handle,
+                          BallCoverIndex<value_idx, value_t, value_int> &index,
+                          const value_t *query, const value_int n_query_rows,
+                          value_int k, const value_idx *R_knn_inds,
+                          const value_t *R_knn_dists, dist_func dfunc,
+                          value_idx *inds, value_t *dists, float weight,
+                          value_int *post_dists_counter) {
+  const value_int bitset_size = ceil(index.n_landmarks / 32.0);
+
+  rmm::device_uvector<std::uint32_t> bitset(bitset_size * index.m,
+                                            handle.get_stream());
+
+  perform_post_filter_registers<value_idx, value_t, value_int, 128, dist_func>
+    <<<n_query_rows, 128, bitset_size * sizeof(std::uint32_t),
+       handle.get_stream()>>>(index.get_X(), index.n, R_knn_inds, R_knn_dists,
+                              index.get_R_radius(), index.get_R(),
+                              index.n_landmarks, bitset_size, k, dfunc,
+                              bitset.data(), weight);
+
+  if (k <= 32)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 32, 2, 128, 2>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+  else if (k <= 64)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 64, 3, 128, 2>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+  else if (k <= 128)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 128, 3, 128, 2>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+  else if (k <= 256)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 256, 4, 128, 2>
+      <<<n_query_rows, 128, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+  else if (k <= 512)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 512, 8, 64, 2>
+      <<<n_query_rows, 64, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+  else if (k <= 1024)
+    compute_final_dists_registers<value_idx, value_t, value_int, std::uint32_t,
+                                  dist_func, 1024, 8, 64, 2>
+      <<<n_query_rows, 64, 0, handle.get_stream()>>>(
+        index.get_X(), query, index.n, bitset.data(), bitset_size, R_knn_dists,
+        index.get_R_indptr(), index.get_R_1nn_cols(), index.get_R_1nn_dists(),
+        inds, dists, index.n_landmarks, k, dfunc, post_dists_counter);
+}
+
+};  // namespace detail
+};  // namespace knn
+};  // namespace spatial
+};  // namespace raft
\ No newline at end of file
diff --git a/cpp/include/raft/spatial/knn/detail/block_select_faiss.cuh b/cpp/include/raft/spatial/knn/detail/block_select_faiss.cuh
new file mode 100644
index 0000000000..d2f7bc2210
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/detail/block_select_faiss.cuh
@@ -0,0 +1,224 @@
+/**
+ * Copyright (c) Facebook, Inc. and its affiliates.
+ *
+ * This source code is licensed under the MIT license found in the
+ * LICENSE file in the root directory of this source tree.
+ */
+
+#pragma once
+
+#include <faiss/gpu/utils/StaticUtils.h>
+#include <faiss/gpu/utils/DeviceDefs.cuh>
+#include <faiss/gpu/utils/MergeNetworkUtils.cuh>
+#include <faiss/gpu/utils/PtxUtils.cuh>
+#include <faiss/gpu/utils/Select.cuh>
+#include <faiss/gpu/utils/WarpShuffles.cuh>
+
+#include "warp_select_faiss.cuh"
+
+// TODO: Need to think further about the impact (and new boundaries created) on the registers
+// because this will change the max k that can be processed. One solution might be to break
+// up k into multiple batches for larger k.
+
+namespace faiss {
+namespace gpu {
+
+// `Dir` true, produce largest values.
+// `Dir` false, produce smallest values.
+template <typename K, typename V, bool Dir, typename Comp, int NumWarpQ,
+          int NumThreadQ, int ThreadsPerBlock>
+struct KeyValueBlockSelect {
+  static constexpr int kNumWarps = ThreadsPerBlock / kWarpSize;
+  static constexpr int kTotalWarpSortSize = NumWarpQ;
+
+  __device__ inline KeyValueBlockSelect(K initKVal, K initVKey, V initVVal,
+                                        K* smemK, KeyValuePair<K, V>* smemV,
+                                        int k)
+    : initK(initKVal),
+      initVk(initVKey),
+      initVv(initVVal),
+      numVals(0),
+      warpKTop(initKVal),
+      warpKTopRDist(initKVal),
+      sharedK(smemK),
+      sharedV(smemV),
+      kMinus1(k - 1) {
+    static_assert(utils::isPowerOf2(ThreadsPerBlock),
+                  "threads must be a power-of-2");
+    static_assert(utils::isPowerOf2(NumWarpQ), "warp queue must be power-of-2");
+
+    // Fill the per-thread queue keys with the default value
+#pragma unroll
+    for (int i = 0; i < NumThreadQ; ++i) {
+      threadK[i] = initK;
+      threadV[i].key = initVk;
+      threadV[i].value = initVv;
+    }
+
+    int laneId = getLaneId();
+    int warpId = threadIdx.x / kWarpSize;
+    warpK = sharedK + warpId * kTotalWarpSortSize;
+    warpV = sharedV + warpId * kTotalWarpSortSize;
+
+    // Fill warp queue (only the actual queue space is fine, not where
+    // we write the per-thread queues for merging)
+    for (int i = laneId; i < NumWarpQ; i += kWarpSize) {
+      warpK[i] = initK;
+      warpV[i].key = initVk;
+      warpV[i].value = initVv;
+    }
+
+    warpFence();
+  }
+
+  __device__ inline void addThreadQ(K k, K vk, V vv) {
+    if (Dir ? Comp::gt(k, warpKTop) : Comp::lt(k, warpKTop)) {
+      // Rotate right
+#pragma unroll
+      for (int i = NumThreadQ - 1; i > 0; --i) {
+        threadK[i] = threadK[i - 1];
+        threadV[i].key = threadV[i - 1].key;
+        threadV[i].value = threadV[i - 1].value;
+      }
+
+      threadK[0] = k;
+      threadV[0].key = vk;
+      threadV[0].value = vv;
+      ++numVals;
+    }
+  }
+
+  __device__ inline void checkThreadQ() {
+    bool needSort = (numVals == NumThreadQ);
+
+#if CUDA_VERSION >= 9000
+    needSort = __any_sync(0xffffffff, needSort);
+#else
+    needSort = __any(needSort);
+#endif
+
+    if (!needSort) {
+      // no lanes have triggered a sort
+      return;
+    }
+
+    // This has a trailing warpFence
+    mergeWarpQ();
+
+    // Any top-k elements have been merged into the warp queue; we're
+    // free to reset the thread queues
+    numVals = 0;
+
+#pragma unroll
+    for (int i = 0; i < NumThreadQ; ++i) {
+      threadK[i] = initK;
+      threadV[i].key = initVk;
+      threadV[i].value = initVv;
+    }
+
+    // We have to beat at least this element
+    warpKTop = warpK[kMinus1];
+    warpKTopRDist = warpV[kMinus1].key;
+
+    warpFence();
+  }
+
+  /// This function handles sorting and merging together the
+  /// per-thread queues with the warp-wide queue, creating a sorted
+  /// list across both
+  __device__ inline void mergeWarpQ() {
+    int laneId = getLaneId();
+
+    // Sort all of the per-thread queues
+    warpSortAnyRegistersKVP<K, V, NumThreadQ, !Dir, Comp>(threadK, threadV);
+
+    constexpr int kNumWarpQRegisters = NumWarpQ / kWarpSize;
+    K warpKRegisters[kNumWarpQRegisters];
+    KeyValuePair<K, V> warpVRegisters[kNumWarpQRegisters];
+
+#pragma unroll
+    for (int i = 0; i < kNumWarpQRegisters; ++i) {
+      warpKRegisters[i] = warpK[i * kWarpSize + laneId];
+      warpVRegisters[i].key = warpV[i * kWarpSize + laneId].key;
+      warpVRegisters[i].value = warpV[i * kWarpSize + laneId].value;
+    }
+
+    warpFence();
+
+    // The warp queue is already sorted, and now that we've sorted the
+    // per-thread queue, merge both sorted lists together, producing
+    // one sorted list
+    warpMergeAnyRegistersKVP<K, V, kNumWarpQRegisters, NumThreadQ, !Dir, Comp,
+                             false>(warpKRegisters, warpVRegisters, threadK,
+                                    threadV);
+
+    // Write back out the warp queue
+#pragma unroll
+    for (int i = 0; i < kNumWarpQRegisters; ++i) {
+      warpK[i * kWarpSize + laneId] = warpKRegisters[i];
+      warpV[i * kWarpSize + laneId].key = warpVRegisters[i].key;
+      warpV[i * kWarpSize + laneId].value = warpVRegisters[i].value;
+    }
+
+    warpFence();
+  }
+
+  /// WARNING: all threads in a warp must participate in this.
+  /// Otherwise, you must call the constituent parts separately.
+  __device__ inline void add(K k, K vk, V vv) {
+    addThreadQ(k, vk, vv);
+    checkThreadQ();
+  }
+
+  __device__ inline void reduce() {
+    // Have all warps dump and merge their queues; this will produce
+    // the final per-warp results
+    mergeWarpQ();
+
+    // block-wide dep; thus far, all warps have been completely
+    // independent
+    __syncthreads();
+
+    // All warp queues are contiguous in smem.
+    // Now, we have kNumWarps lists of NumWarpQ elements.
+    // This is a power of 2.
+    FinalBlockMerge<kNumWarps, ThreadsPerBlock, K, KeyValuePair<K, V>, NumWarpQ,
+                    Dir, Comp>::merge(sharedK, sharedV);
+
+    // The block-wide merge has a trailing syncthreads
+  }
+
+  // Default element key
+  const K initK;
+
+  // Default element value
+  const K initVk;
+  const V initVv;
+
+  // Number of valid elements in our thread queue
+  int numVals;
+
+  // The k-th highest (Dir) or lowest (!Dir) element
+  K warpKTop;
+
+  K warpKTopRDist;
+
+  // Thread queue values
+  K threadK[NumThreadQ];
+  KeyValuePair<K, V> threadV[NumThreadQ];
+
+  // Queues for all warps
+  K* sharedK;
+  KeyValuePair<K, V>* sharedV;
+
+  // Our warp's queue (points into sharedK/sharedV)
+  // warpK[0] is highest (Dir) or lowest (!Dir)
+  K* warpK;
+  KeyValuePair<K, V>* warpV;
+
+  // This is a cached k-1 value
+  int kMinus1;
+};
+
+}  // namespace gpu
+}  // namespace faiss
\ No newline at end of file
diff --git a/cpp/include/raft/spatial/knn/detail/knn_brute_force_faiss.cuh b/cpp/include/raft/spatial/knn/detail/knn_brute_force_faiss.cuh
index 94ace19580..43fc4a164f 100644
--- a/cpp/include/raft/spatial/knn/detail/knn_brute_force_faiss.cuh
+++ b/cpp/include/raft/spatial/knn/detail/knn_brute_force_faiss.cuh
@@ -30,6 +30,7 @@
 
 #include <raft/linalg/distance_type.h>
 #include <thrust/iterator/transform_iterator.h>
+#include <cstdint>
 #include <iostream>
 #include <raft/handle.hpp>
 #include <set>
@@ -45,8 +46,8 @@ namespace spatial {
 namespace knn {
 namespace detail {
 
-template <typename value_idx = int64_t, typename value_t = float, int warp_q,
-          int thread_q, int tpb>
+template <typename value_idx = std::int64_t, typename value_t = float,
+          int warp_q, int thread_q, int tpb>
 __global__ void knn_merge_parts_kernel(value_t *inK, value_idx *inV,
                                        value_t *outK, value_idx *outV,
                                        size_t n_samples, int n_parts,
@@ -110,8 +111,8 @@ __global__ void knn_merge_parts_kernel(value_t *inK, value_idx *inV,
   }
 }
 
-template <typename value_idx = int64_t, typename value_t = float, int warp_q,
-          int thread_q>
+template <typename value_idx = std::int64_t, typename value_t = float,
+          int warp_q, int thread_q>
 inline void knn_merge_parts_impl(value_t *inK, value_idx *inV, value_t *outK,
                                  value_idx *outV, size_t n_samples, int n_parts,
                                  int k, cudaStream_t stream,
@@ -143,7 +144,7 @@ inline void knn_merge_parts_impl(value_t *inK, value_idx *inV, value_t *outK,
  * @param stream CUDA stream to use
  * @param translations mapping of index offsets for each partition
  */
-template <typename value_idx = int64_t, typename value_t = float>
+template <typename value_idx = std::int64_t, typename value_t = float>
 inline void knn_merge_parts(value_t *inK, value_idx *inV, value_t *outK,
                             value_idx *outV, size_t n_samples, int n_parts,
                             int k, cudaStream_t stream,
@@ -195,28 +196,28 @@ inline void knn_merge_parts(value_t *inK, value_idx *inV, value_t *outK,
  * @param[in] metric corresponds to the raft::distance::DistanceType enum (default is L2Expanded)
  * @param[in] metricArg metric argument to use. Corresponds to the p arg for lp norm
  */
-template <typename IntType = int>
-void brute_force_knn_impl(std::vector<float *> &input, std::vector<int> &sizes,
-                          IntType D, float *search_items, IntType n,
-                          int64_t *res_I, float *res_D, IntType k,
-                          cudaStream_t userStream,
+template <typename IntType = int, typename IdxType = std::int64_t>
+void brute_force_knn_impl(std::vector<float *> &input,
+                          std::vector<IntType> &sizes, IntType D,
+                          float *search_items, IntType n, IdxType *res_I,
+                          float *res_D, IntType k, cudaStream_t userStream,
                           cudaStream_t *internalStreams = nullptr,
                           int n_int_streams = 0, bool rowMajorIndex = true,
                           bool rowMajorQuery = true,
-                          std::vector<int64_t> *translations = nullptr,
+                          std::vector<IdxType> *translations = nullptr,
                           raft::distance::DistanceType metric =
                             raft::distance::DistanceType::L2Expanded,
                           float metricArg = 0) {
   ASSERT(input.size() == sizes.size(),
          "input and sizes vectors should be the same size");
 
-  std::vector<int64_t> *id_ranges;
+  std::vector<IdxType> *id_ranges;
   if (translations == nullptr) {
     // If we don't have explicit translations
     // for offsets of the indices, build them
     // from the local partitions
-    id_ranges = new std::vector<int64_t>();
-    int64_t total_n = 0;
+    id_ranges = new std::vector<IdxType>();
+    IdxType total_n = 0;
     for (size_t i = 0; i < input.size(); i++) {
       id_ranges->push_back(total_n);
       total_n += sizes[i];
@@ -242,15 +243,15 @@ void brute_force_knn_impl(std::vector<float *> &input, std::vector<int> &sizes,
   int device;
   CUDA_CHECK(cudaGetDevice(&device));
 
-  rmm::device_uvector<int64_t> trans(id_ranges->size(), userStream);
+  rmm::device_uvector<std::int64_t> trans(id_ranges->size(), userStream);
   raft::update_device(trans.data(), id_ranges->data(), id_ranges->size(),
                       userStream);
 
   rmm::device_uvector<float> all_D(0, userStream);
-  rmm::device_uvector<int64_t> all_I(0, userStream);
+  rmm::device_uvector<std::int64_t> all_I(0, userStream);
 
   float *out_D = res_D;
-  int64_t *out_I = res_I;
+  IdxType *out_I = res_I;
 
   if (input.size() > 1) {
     all_D.resize(input.size() * k * n, userStream);
@@ -265,7 +266,7 @@ void brute_force_knn_impl(std::vector<float *> &input, std::vector<int> &sizes,
 
   for (size_t i = 0; i < input.size(); i++) {
     float *out_d_ptr = out_D + (i * k * n);
-    int64_t *out_i_ptr = out_I + (i * k * n);
+    IdxType *out_i_ptr = out_I + (i * k * n);
 
     cudaStream_t stream =
       raft::select_stream(userStream, internalStreams, n_int_streams, i);
diff --git a/cpp/include/raft/sparse/selection/selection.cuh b/cpp/include/raft/spatial/knn/detail/selection_faiss.cuh
similarity index 94%
rename from cpp/include/raft/sparse/selection/selection.cuh
rename to cpp/include/raft/spatial/knn/detail/selection_faiss.cuh
index 6066a36289..045edad0e6 100644
--- a/cpp/include/raft/sparse/selection/selection.cuh
+++ b/cpp/include/raft/spatial/knn/detail/selection_faiss.cuh
@@ -17,13 +17,6 @@
 #pragma once
 
 #include <raft/cudart_utils.h>
-#include <raft/sparse/cusparse_wrappers.h>
-#include <raft/cuda_utils.cuh>
-#include <raft/matrix/matrix.cuh>
-
-#include <raft/sparse/coo.cuh>
-#include <raft/sparse/csr.cuh>
-#include <raft/sparse/distance/distance.cuh>
 
 #include <faiss/gpu/GpuDistance.h>
 #include <faiss/gpu/GpuIndexFlat.h>
@@ -33,11 +26,10 @@
 #include <faiss/gpu/utils/Limits.cuh>
 #include <faiss/gpu/utils/Select.cuh>
 
-#include <cusparse_v2.h>
-
 namespace raft {
-namespace sparse {
-namespace selection {
+namespace spatial {
+namespace knn {
+namespace detail {
 
 template <typename K, typename IndexType, bool select_min, int warp_q,
           int thread_q, int tpb>
@@ -152,6 +144,7 @@ inline void select_k(value_t *inK, value_idx *inV, size_t n_rows, size_t n_cols,
                                                outV, select_min, k, stream);
 }
 
-};  // namespace selection
-};  // namespace sparse
+};  // namespace detail
+};  // namespace knn
+};  // namespace spatial
 };  // namespace raft
diff --git a/cpp/include/raft/spatial/knn/detail/warp_select_faiss.cuh b/cpp/include/raft/spatial/knn/detail/warp_select_faiss.cuh
new file mode 100644
index 0000000000..84719a0e4b
--- /dev/null
+++ b/cpp/include/raft/spatial/knn/detail/warp_select_faiss.cuh
@@ -0,0 +1,739 @@
+/**
+ * Copyright (c) Facebook, Inc. and its affiliates.
+ *
+ * This source code is licensed under the MIT license found in the
+ * LICENSE file in the root directory of this source tree.
+ */
+
+#pragma once
+
+#include <cub/cub.cuh>
+
+#include <faiss/gpu/utils/StaticUtils.h>
+#include <faiss/gpu/utils/DeviceDefs.cuh>
+#include <faiss/gpu/utils/MergeNetworkUtils.cuh>
+#include <faiss/gpu/utils/PtxUtils.cuh>
+#include <faiss/gpu/utils/WarpShuffles.cuh>
+
+namespace faiss {
+namespace gpu {
+
+template <typename _Key, typename _Value>
+struct KeyValuePair {
+  typedef _Key Key;      ///< Key data type
+  typedef _Value Value;  ///< Value data type
+
+  Key key;      ///< Item key
+  Value value;  ///< Item value
+
+  /// Constructor
+  __host__ __device__ __forceinline__ KeyValuePair() {}
+
+  /// Copy Constructors
+  __host__ __device__ __forceinline__
+  KeyValuePair(cub::KeyValuePair<_Key, _Value>& kvp)
+    : key(kvp.key), value(kvp.value) {}
+
+  __host__ __device__ __forceinline__
+  KeyValuePair(faiss::gpu::KeyValuePair<_Key, _Value>& kvp)
+    : key(kvp.key), value(kvp.value) {}
+
+  /// Constructor
+  __host__ __device__ __forceinline__ KeyValuePair(Key const& key,
+                                                   Value const& value)
+    : key(key), value(value) {}
+
+  /// Inequality operator
+  __host__ __device__ __forceinline__ bool operator!=(const KeyValuePair& b) {
+    return (value != b.value) || (key != b.key);
+  }
+};
+
+//
+// This file contains functions to:
+//
+// -perform bitonic merges on pairs of sorted lists, held in
+// registers. Each list contains N * kWarpSize (multiple of 32)
+// elements for some N.
+// The bitonic merge is implemented for arbitrary sizes;
+// sorted list A of size N1 * kWarpSize registers
+// sorted list B of size N2 * kWarpSize registers =>
+// sorted list C if size (N1 + N2) * kWarpSize registers. N1 and N2
+// are >= 1 and don't have to be powers of 2.
+//
+// -perform bitonic sorts on a set of N * kWarpSize key/value pairs
+// held in registers, by using the above bitonic merge as a
+// primitive.
+// N can be an arbitrary N >= 1; i.e., the bitonic sort here supports
+// odd sizes and doesn't require the input to be a power of 2.
+//
+// The sort or merge network is completely statically instantiated via
+// template specialization / expansion and constexpr, and it uses warp
+// shuffles to exchange values between warp lanes.
+//
+// A note about comparsions:
+//
+// For a sorting network of keys only, we only need one
+// comparison (a < b). However, what we really need to know is
+// if one lane chooses to exchange a value, then the
+// corresponding lane should also do the exchange.
+// Thus, if one just uses the negation !(x < y) in the higher
+// lane, this will also include the case where (x == y). Thus, one
+// lane in fact performs an exchange and the other doesn't, but
+// because the only value being exchanged is equivalent, nothing has
+// changed.
+// So, you can get away with just one comparison and its negation.
+//
+// If we're sorting keys and values, where equivalent keys can
+// exist, then this is a problem, since we want to treat (x, v1)
+// as not equivalent to (x, v2).
+//
+// To remedy this, you can either compare with a lexicographic
+// ordering (a.k < b.k || (a.k == b.k && a.v < b.v)), which since
+// we're predicating all of the choices results in 3 comparisons
+// being executed, or we can invert the selection so that there is no
+// middle choice of equality; the other lane will likewise
+// check that (b.k > a.k) (the higher lane has the values
+// swapped). Then, the first lane swaps if and only if the
+// second lane swaps; if both lanes have equivalent keys, no
+// swap will be performed. This results in only two comparisons
+// being executed.
+//
+// If you don't consider values as well, then this does not produce a
+// consistent ordering among (k, v) pairs with equivalent keys but
+// different values; for us, we don't really care about ordering or
+// stability here.
+//
+// I have tried both re-arranging the order in the higher lane to get
+// away with one comparison or adding the value to the check; both
+// result in greater register consumption or lower speed than just
+// perfoming both < and > comparisons with the variables, so I just
+// stick with this.
+
+// This function merges kWarpSize / 2L lists in parallel using warp
+// shuffles.
+// It works on at most size-16 lists, as we need 32 threads for this
+// shuffle merge.
+//
+// If IsBitonic is false, the first stage is reversed, so we don't
+// need to sort directionally. It's still technically a bitonic sort.
+template <typename K, typename V, int L, bool Dir, typename Comp,
+          bool IsBitonic>
+inline __device__ void warpBitonicMergeLE16KVP(K& k, KeyValuePair<K, V>& v) {
+  static_assert(utils::isPowerOf2(L), "L must be a power-of-2");
+  static_assert(L <= kWarpSize / 2, "merge list size must be <= 16");
+
+  int laneId = getLaneId();
+
+  if (!IsBitonic) {
+    // Reverse the first comparison stage.
+    // For example, merging a list of size 8 has the exchanges:
+    // 0 <-> 15, 1 <-> 14, ...
+    K otherK = shfl_xor(k, 2 * L - 1);
+    K otherVk = shfl_xor(v.key, 2 * L - 1);
+    V otherVv = shfl_xor(v.value, 2 * L - 1);
+
+    KeyValuePair<K, V> otherV = KeyValuePair(otherVk, otherVv);
+
+    // Whether we are the lesser thread in the exchange
+    bool small = !(laneId & L);
+
+    if (Dir) {
+      // See the comment above how performing both of these
+      // comparisons in the warp seems to win out over the
+      // alternatives in practice
+      bool s = small ? Comp::gt(k, otherK) : Comp::lt(k, otherK);
+      assign(s, k, otherK);
+      assign(s, v.key, otherV.key);
+      assign(s, v.value, otherV.value);
+
+    } else {
+      bool s = small ? Comp::lt(k, otherK) : Comp::gt(k, otherK);
+      assign(s, k, otherK);
+      assign(s, v.value, otherV.value);
+      assign(s, v.key, otherV.key);
+    }
+  }
+
+#pragma unroll
+  for (int stride = IsBitonic ? L : L / 2; stride > 0; stride /= 2) {
+    K otherK = shfl_xor(k, stride);
+    K otherVk = shfl_xor(v.key, stride);
+    V otherVv = shfl_xor(v.value, stride);
+
+    KeyValuePair<K, V> otherV = KeyValuePair(otherVk, otherVv);
+
+    // Whether we are the lesser thread in the exchange
+    bool small = !(laneId & stride);
+
+    if (Dir) {
+      bool s = small ? Comp::gt(k, otherK) : Comp::lt(k, otherK);
+      assign(s, k, otherK);
+      assign(s, v.key, otherV.key);
+      assign(s, v.value, otherV.value);
+
+    } else {
+      bool s = small ? Comp::lt(k, otherK) : Comp::gt(k, otherK);
+      assign(s, k, otherK);
+      assign(s, v.key, otherV.key);
+      assign(s, v.value, otherV.value);
+    }
+  }
+}
+
+// Template for performing a bitonic merge of an arbitrary set of
+// registers
+template <typename K, typename V, int N, bool Dir, typename Comp, bool Low,
+          bool Pow2>
+struct BitonicMergeStepKVP {};
+
+//
+// Power-of-2 merge specialization
+//
+
+// All merges eventually call this
+template <typename K, typename V, bool Dir, typename Comp, bool Low>
+struct BitonicMergeStepKVP<K, V, 1, Dir, Comp, Low, true> {
+  static inline __device__ void merge(K k[1], KeyValuePair<K, V> v[1]) {
+    // Use warp shuffles
+    warpBitonicMergeLE16KVP<K, V, 16, Dir, Comp, true>(k[0], v[0]);
+  }
+};
+
+template <typename K, typename V, int N, bool Dir, typename Comp, bool Low>
+struct BitonicMergeStepKVP<K, V, N, Dir, Comp, Low, true> {
+  static inline __device__ void merge(K k[N], KeyValuePair<K, V> v[N]) {
+    static_assert(utils::isPowerOf2(N), "must be power of 2");
+    static_assert(N > 1, "must be N > 1");
+
+#pragma unroll
+    for (int i = 0; i < N / 2; ++i) {
+      K& ka = k[i];
+      KeyValuePair<K, V>& va = v[i];
+
+      K& kb = k[i + N / 2];
+      KeyValuePair<K, V>& vb = v[i + N / 2];
+
+      bool s = Dir ? Comp::gt(ka, kb) : Comp::lt(ka, kb);
+      swap(s, ka, kb);
+      swap(s, va.key, vb.key);
+      swap(s, va.value, vb.value);
+    }
+
+    {
+      K newK[N / 2];
+      KeyValuePair<K, V> newV[N / 2];
+
+#pragma unroll
+      for (int i = 0; i < N / 2; ++i) {
+        newK[i] = k[i];
+        newV[i].key = v[i].key;
+        newV[i].value = v[i].value;
+      }
+
+      BitonicMergeStepKVP<K, V, N / 2, Dir, Comp, true, true>::merge(newK,
+                                                                     newV);
+
+#pragma unroll
+      for (int i = 0; i < N / 2; ++i) {
+        k[i] = newK[i];
+        v[i].key = newV[i].key;
+        v[i].value = newV[i].value;
+      }
+    }
+
+    {
+      K newK[N / 2];
+      KeyValuePair<K, V> newV[N / 2];
+
+#pragma unroll
+      for (int i = 0; i < N / 2; ++i) {
+        newK[i] = k[i + N / 2];
+        newV[i].key = v[i + N / 2].key;
+        newV[i].value = v[i + N / 2].value;
+      }
+
+      BitonicMergeStepKVP<K, V, N / 2, Dir, Comp, false, true>::merge(newK,
+                                                                      newV);
+
+#pragma unroll
+      for (int i = 0; i < N / 2; ++i) {
+        k[i + N / 2] = newK[i];
+        v[i + N / 2].key = newV[i].key;
+        v[i + N / 2].value = newV[i].value;
+      }
+    }
+  }
+};
+
+//
+// Non-power-of-2 merge specialization
+//
+
+// Low recursion
+template <typename K, typename V, int N, bool Dir, typename Comp>
+struct BitonicMergeStepKVP<K, V, N, Dir, Comp, true, false> {
+  static inline __device__ void merge(K k[N], KeyValuePair<K, V> v[N]) {
+    static_assert(!utils::isPowerOf2(N), "must be non-power-of-2");
+    static_assert(N >= 3, "must be N >= 3");
+
+    constexpr int kNextHighestPowerOf2 = utils::nextHighestPowerOf2(N);
+
+#pragma unroll
+    for (int i = 0; i < N - kNextHighestPowerOf2 / 2; ++i) {
+      K& ka = k[i];
+      KeyValuePair<K, V>& va = v[i];
+
+      K& kb = k[i + kNextHighestPowerOf2 / 2];
+      KeyValuePair<K, V>& vb = v[i + kNextHighestPowerOf2 / 2];
+
+      bool s = Dir ? Comp::gt(ka, kb) : Comp::lt(ka, kb);
+      swap(s, ka, kb);
+      swap(s, va.key, vb.key);
+      swap(s, va.value, vb.value);
+    }
+
+    constexpr int kLowSize = N - kNextHighestPowerOf2 / 2;
+    constexpr int kHighSize = kNextHighestPowerOf2 / 2;
+    {
+      K newK[kLowSize];
+      KeyValuePair<K, V> newV[kLowSize];
+
+#pragma unroll
+      for (int i = 0; i < kLowSize; ++i) {
+        newK[i] = k[i];
+        newV[i].key = v[i].key;
+        newV[i].value = v[i].value;
+      }
+
+      constexpr bool kLowIsPowerOf2 =
+        utils::isPowerOf2(N - kNextHighestPowerOf2 / 2);
+      // FIXME: compiler doesn't like this expression? compiler bug?
+      //      constexpr bool kLowIsPowerOf2 = utils::isPowerOf2(kLowSize);
+      BitonicMergeStepKVP<K, V, kLowSize, Dir, Comp,
+                          true,  // low
+                          kLowIsPowerOf2>::merge(newK, newV);
+
+#pragma unroll
+      for (int i = 0; i < kLowSize; ++i) {
+        k[i] = newK[i];
+        v[i].key = newV[i].key;
+        v[i].value = newV[i].value;
+      }
+    }
+
+    {
+      K newK[kHighSize];
+      KeyValuePair<K, V> newV[kHighSize];
+
+#pragma unroll
+      for (int i = 0; i < kHighSize; ++i) {
+        newK[i] = k[i + kLowSize];
+        newV[i].key = v[i + kLowSize].key;
+        newV[i].value = v[i + kLowSize].value;
+      }
+
+      constexpr bool kHighIsPowerOf2 =
+        utils::isPowerOf2(kNextHighestPowerOf2 / 2);
+      // FIXME: compiler doesn't like this expression? compiler bug?
+      //      constexpr bool kHighIsPowerOf2 = utils::isPowerOf2(kHighSize);
+      BitonicMergeStepKVP<K, V, kHighSize, Dir, Comp,
+                          false,  // high
+                          kHighIsPowerOf2>::merge(newK, newV);
+
+#pragma unroll
+      for (int i = 0; i < kHighSize; ++i) {
+        k[i + kLowSize] = newK[i];
+        v[i + kLowSize].key = newV[i].key;
+        v[i + kLowSize].value = newV[i].value;
+      }
+    }
+  }
+};
+
+// High recursion
+template <typename K, typename V, int N, bool Dir, typename Comp>
+struct BitonicMergeStepKVP<K, V, N, Dir, Comp, false, false> {
+  static inline __device__ void merge(K k[N], KeyValuePair<K, V> v[N]) {
+    static_assert(!utils::isPowerOf2(N), "must be non-power-of-2");
+    static_assert(N >= 3, "must be N >= 3");
+
+    constexpr int kNextHighestPowerOf2 = utils::nextHighestPowerOf2(N);
+
+#pragma unroll
+    for (int i = 0; i < N - kNextHighestPowerOf2 / 2; ++i) {
+      K& ka = k[i];
+      KeyValuePair<K, V>& va = v[i];
+
+      K& kb = k[i + kNextHighestPowerOf2 / 2];
+      KeyValuePair<K, V>& vb = v[i + kNextHighestPowerOf2 / 2];
+
+      bool s = Dir ? Comp::gt(ka, kb) : Comp::lt(ka, kb);
+      swap(s, ka, kb);
+      swap(s, va.key, vb.key);
+      swap(s, va.value, vb.value);
+    }
+
+    constexpr int kLowSize = kNextHighestPowerOf2 / 2;
+    constexpr int kHighSize = N - kNextHighestPowerOf2 / 2;
+    {
+      K newK[kLowSize];
+      KeyValuePair<K, V> newV[kLowSize];
+
+#pragma unroll
+      for (int i = 0; i < kLowSize; ++i) {
+        newK[i] = k[i];
+        newV[i].key = v[i].key;
+        newV[i].value = v[i].value;
+      }
+
+      constexpr bool kLowIsPowerOf2 =
+        utils::isPowerOf2(kNextHighestPowerOf2 / 2);
+      // FIXME: compiler doesn't like this expression? compiler bug?
+      //      constexpr bool kLowIsPowerOf2 = utils::isPowerOf2(kLowSize);
+      BitonicMergeStepKVP<K, V, kLowSize, Dir, Comp,
+                          true,  // low
+                          kLowIsPowerOf2>::merge(newK, newV);
+
+#pragma unroll
+      for (int i = 0; i < kLowSize; ++i) {
+        k[i] = newK[i];
+        v[i].key = newV[i].key;
+        v[i].value = newV[i].value;
+      }
+    }
+
+    {
+      K newK[kHighSize];
+      KeyValuePair<K, V> newV[kHighSize];
+
+#pragma unroll
+      for (int i = 0; i < kHighSize; ++i) {
+        newK[i] = k[i + kLowSize];
+        newV[i].key = v[i + kLowSize].key;
+        newV[i].value = v[i + kLowSize].value;
+      }
+
+      constexpr bool kHighIsPowerOf2 =
+        utils::isPowerOf2(N - kNextHighestPowerOf2 / 2);
+      // FIXME: compiler doesn't like this expression? compiler bug?
+      //      constexpr bool kHighIsPowerOf2 = utils::isPowerOf2(kHighSize);
+      BitonicMergeStepKVP<K, V, kHighSize, Dir, Comp,
+                          false,  // high
+                          kHighIsPowerOf2>::merge(newK, newV);
+
+#pragma unroll
+      for (int i = 0; i < kHighSize; ++i) {
+        k[i + kLowSize] = newK[i];
+        v[i + kLowSize].key = newV[i].key;
+        v[i + kLowSize].value = newV[i].value;
+      }
+    }
+  }
+};
+
+/// Merges two sets of registers across the warp of any size;
+/// i.e., merges a sorted k/v list of size kWarpSize * N1 with a
+/// sorted k/v list of size kWarpSize * N2, where N1 and N2 are any
+/// value >= 1
+template <typename K, typename V, int N1, int N2, bool Dir, typename Comp,
+          bool FullMerge = true>
+inline __device__ void warpMergeAnyRegistersKVP(K k1[N1],
+                                                KeyValuePair<K, V> v1[N1],
+                                                K k2[N2],
+                                                KeyValuePair<K, V> v2[N2]) {
+  constexpr int kSmallestN = N1 < N2 ? N1 : N2;
+
+#pragma unroll
+  for (int i = 0; i < kSmallestN; ++i) {
+    K& ka = k1[N1 - 1 - i];
+    KeyValuePair<K, V>& va = v1[N1 - 1 - i];
+
+    K& kb = k2[i];
+    KeyValuePair<K, V>& vb = v2[i];
+
+    K otherKa;
+    KeyValuePair<K, V> otherVa;
+
+    if (FullMerge) {
+      // We need the other values
+      otherKa = shfl_xor(ka, kWarpSize - 1);
+      K otherVak = shfl_xor(va.key, kWarpSize - 1);
+      V otherVav = shfl_xor(va.value, kWarpSize - 1);
+      otherVa = KeyValuePair(otherVak, otherVav);
+    }
+
+    K otherKb = shfl_xor(kb, kWarpSize - 1);
+    K otherVbk = shfl_xor(vb.key, kWarpSize - 1);
+    V otherVbv = shfl_xor(vb.value, kWarpSize - 1);
+
+    // ka is always first in the list, so we needn't use our lane
+    // in this comparison
+    bool swapa = Dir ? Comp::gt(ka, otherKb) : Comp::lt(ka, otherKb);
+    assign(swapa, ka, otherKb);
+    assign(swapa, va.key, otherVbk);
+    assign(swapa, va.value, otherVbv);
+
+    // kb is always second in the list, so we needn't use our lane
+    // in this comparison
+    if (FullMerge) {
+      bool swapb = Dir ? Comp::lt(kb, otherKa) : Comp::gt(kb, otherKa);
+      assign(swapb, kb, otherKa);
+      assign(swapb, vb.key, otherVa.key);
+      assign(swapb, vb.value, otherVa.value);
+
+    } else {
+      // We don't care about updating elements in the second list
+    }
+  }
+
+  BitonicMergeStepKVP<K, V, N1, Dir, Comp, true, utils::isPowerOf2(N1)>::merge(
+    k1, v1);
+  if (FullMerge) {
+    // Only if we care about N2 do we need to bother merging it fully
+    BitonicMergeStepKVP<K, V, N2, Dir, Comp, false,
+                        utils::isPowerOf2(N2)>::merge(k2, v2);
+  }
+}
+
+// Recursive template that uses the above bitonic merge to perform a
+// bitonic sort
+template <typename K, typename V, int N, bool Dir, typename Comp>
+struct BitonicSortStepKVP {
+  static inline __device__ void sort(K k[N], KeyValuePair<K, V> v[N]) {
+    static_assert(N > 1, "did not hit specialized case");
+
+    // Sort recursively
+    constexpr int kSizeA = N / 2;
+    constexpr int kSizeB = N - kSizeA;
+
+    K aK[kSizeA];
+    KeyValuePair<K, V> aV[kSizeA];
+
+#pragma unroll
+    for (int i = 0; i < kSizeA; ++i) {
+      aK[i] = k[i];
+      aV[i].key = v[i].key;
+      aV[i].value = v[i].value;
+    }
+
+    BitonicSortStepKVP<K, V, kSizeA, Dir, Comp>::sort(aK, aV);
+
+    K bK[kSizeB];
+    KeyValuePair<K, V> bV[kSizeB];
+
+#pragma unroll
+    for (int i = 0; i < kSizeB; ++i) {
+      bK[i] = k[i + kSizeA];
+      bV[i].key = v[i + kSizeA].key;
+      bV[i].value = v[i + kSizeA].value;
+    }
+
+    BitonicSortStepKVP<K, V, kSizeB, Dir, Comp>::sort(bK, bV);
+
+    // Merge halves
+    warpMergeAnyRegistersKVP<K, V, kSizeA, kSizeB, Dir, Comp>(aK, aV, bK, bV);
+
+#pragma unroll
+    for (int i = 0; i < kSizeA; ++i) {
+      k[i] = aK[i];
+      v[i].key = aV[i].key;
+      v[i].value = aV[i].value;
+    }
+
+#pragma unroll
+    for (int i = 0; i < kSizeB; ++i) {
+      k[i + kSizeA] = bK[i];
+      v[i + kSizeA].key = bV[i].key;
+      v[i + kSizeA].value = bV[i].value;
+    }
+  }
+};
+
+// Single warp (N == 1) sorting specialization
+template <typename K, typename V, bool Dir, typename Comp>
+struct BitonicSortStepKVP<K, V, 1, Dir, Comp> {
+  static inline __device__ void sort(K k[1], KeyValuePair<K, V> v[1]) {
+    // Update this code if this changes
+    // should go from 1 -> kWarpSize in multiples of 2
+    static_assert(kWarpSize == 32, "unexpected warp size");
+
+    warpBitonicMergeLE16KVP<K, V, 1, Dir, Comp, false>(k[0], v[0]);
+    warpBitonicMergeLE16KVP<K, V, 2, Dir, Comp, false>(k[0], v[0]);
+    warpBitonicMergeLE16KVP<K, V, 4, Dir, Comp, false>(k[0], v[0]);
+    warpBitonicMergeLE16KVP<K, V, 8, Dir, Comp, false>(k[0], v[0]);
+    warpBitonicMergeLE16KVP<K, V, 16, Dir, Comp, false>(k[0], v[0]);
+  }
+};
+
+/// Sort a list of kWarpSize * N elements in registers, where N is an
+/// arbitrary >= 1
+template <typename K, typename V, int N, bool Dir, typename Comp>
+inline __device__ void warpSortAnyRegistersKVP(K k[N],
+                                               KeyValuePair<K, V> v[N]) {
+  BitonicSortStepKVP<K, V, N, Dir, Comp>::sort(k, v);
+}
+
+// `Dir` true, produce largest values.
+// `Dir` false, produce smallest values.
+template <typename K, typename V, bool Dir, typename Comp, int NumWarpQ,
+          int NumThreadQ, int ThreadsPerBlock>
+struct KeyValueWarpSelect {
+  static constexpr int kNumWarpQRegisters = NumWarpQ / faiss::gpu::kWarpSize;
+
+  __device__ inline KeyValueWarpSelect(K initKVal,
+                                       faiss::gpu::KeyValuePair<K, V> initVVal,
+                                       int k)
+    : initK(initKVal),
+      initV(initVVal),
+      numVals(0),
+      warpKTop(initKVal),
+      warpKTopRDist(initKVal),
+      kLane((k - 1) % faiss::gpu::kWarpSize) {
+    static_assert(faiss::gpu::utils::isPowerOf2(ThreadsPerBlock),
+                  "threads must be a power-of-2");
+    static_assert(faiss::gpu::utils::isPowerOf2(NumWarpQ),
+                  "warp queue must be power-of-2");
+
+    // Fill the per-thread queue keys with the default value
+#pragma unroll
+    for (int i = 0; i < NumThreadQ; ++i) {
+      threadK[i] = initK;
+      threadV[i].key = initV.key;
+      threadV[i].value = initV.value;
+    }
+
+    // Fill the warp queue with the default value
+#pragma unroll
+    for (int i = 0; i < kNumWarpQRegisters; ++i) {
+      warpK[i] = initK;
+      warpV[i].key = initV.key;
+      warpV[i].value = initV.value;
+    }
+  }
+
+  __device__ inline void addThreadQ(K k, faiss::gpu::KeyValuePair<K, V>& v) {
+    if (Dir ? Comp::gt(k, warpKTop) : Comp::lt(k, warpKTop)) {
+      // Rotate right
+#pragma unroll
+      for (int i = NumThreadQ - 1; i > 0; --i) {
+        threadK[i] = threadK[i - 1];
+        threadV[i].key = threadV[i - 1].key;
+        threadV[i].value = threadV[i - 1].value;
+      }
+
+      threadK[0] = k;
+      threadV[0].key = v.key;
+      threadV[0].value = v.value;
+      ++numVals;
+    }
+  }
+  /// This function handles sorting and merging together the
+  /// per-thread queues with the warp-wide queue, creating a sorted
+  /// list across both
+
+  // TODO
+  __device__ inline void mergeWarpQ() {
+    // Sort all of the per-thread queues
+    faiss::gpu::warpSortAnyRegistersKVP<K, V, NumThreadQ, !Dir, Comp>(threadK,
+                                                                      threadV);
+
+    // The warp queue is already sorted, and now that we've sorted the
+    // per-thread queue, merge both sorted lists together, producing
+    // one sorted list
+    faiss::gpu::warpMergeAnyRegistersKVP<K, V, kNumWarpQRegisters, NumThreadQ,
+                                         !Dir, Comp, false>(warpK, warpV,
+                                                            threadK, threadV);
+  }
+
+  /// WARNING: all threads in a warp must participate in this.
+  /// Otherwise, you must call the constituent parts separately.
+  __device__ inline void add(K k, faiss::gpu::KeyValuePair<K, V>& v) {
+    addThreadQ(k, v);
+    checkThreadQ();
+  }
+
+  __device__ inline void reduce() {
+    // Have all warps dump and merge their queues; this will produce
+    // the final per-warp results
+    mergeWarpQ();
+  }
+
+  __device__ inline void checkThreadQ() {
+    bool needSort = (numVals == NumThreadQ);
+
+#if CUDA_VERSION >= 9000
+    needSort = __any_sync(0xffffffff, needSort);
+#else
+    needSort = __any(needSort);
+#endif
+
+    if (!needSort) {
+      // no lanes have triggered a sort
+      return;
+    }
+
+    mergeWarpQ();
+
+    // Any top-k elements have been merged into the warp queue; we're
+    // free to reset the thread queues
+    numVals = 0;
+
+#pragma unroll
+    for (int i = 0; i < NumThreadQ; ++i) {
+      threadK[i] = initK;
+      threadV[i].key = initV.key;
+      threadV[i].value = initV.value;
+    }
+
+    // We have to beat at least this element
+    warpKTopRDist = shfl(warpV[kNumWarpQRegisters - 1].key, kLane);
+    warpKTop = shfl(warpK[kNumWarpQRegisters - 1], kLane);
+  }
+
+  /// Dump final k selected values for this warp out
+  __device__ inline void writeOut(K* outK, V* outV, int k) {
+    int laneId = faiss::gpu::getLaneId();
+
+#pragma unroll
+    for (int i = 0; i < kNumWarpQRegisters; ++i) {
+      int idx = i * faiss::gpu::kWarpSize + laneId;
+
+      if (idx < k) {
+        outK[idx] = warpK[i];
+        outV[idx] = warpV[i].value;
+      }
+    }
+  }
+
+  // Default element key
+  const K initK;
+
+  // Default element value
+  const faiss::gpu::KeyValuePair<K, V> initV;
+
+  // Number of valid elements in our thread queue
+  int numVals;
+
+  // The k-th highest (Dir) or lowest (!Dir) element
+  K warpKTop;
+
+  // TopK's distance to closest landmark
+  K warpKTopRDist;
+
+  // Thread queue values
+  K threadK[NumThreadQ];
+  faiss::gpu::KeyValuePair<K, V> threadV[NumThreadQ];
+
+  // warpK[0] is highest (Dir) or lowest (!Dir)
+  K warpK[kNumWarpQRegisters];
+  faiss::gpu::KeyValuePair<K, V> warpV[kNumWarpQRegisters];
+
+  // This is what lane we should load an approximation (>=k) to the
+  // kth element from the last register in the warp queue (i.e.,
+  // warpK[kNumWarpQRegisters - 1]).
+  int kLane;
+};
+
+}  // namespace gpu
+}  // namespace faiss
diff --git a/cpp/include/raft/spatial/knn/knn.hpp b/cpp/include/raft/spatial/knn/knn.hpp
index 71c547c281..73866b35a8 100644
--- a/cpp/include/raft/spatial/knn/knn.hpp
+++ b/cpp/include/raft/spatial/knn/knn.hpp
@@ -17,6 +17,7 @@
 #pragma once
 
 #include "detail/knn_brute_force_faiss.cuh"
+#include "detail/selection_faiss.cuh"
 
 #include <raft/mr/device/buffer.hpp>
 
@@ -24,6 +25,32 @@ namespace raft {
 namespace spatial {
 namespace knn {
 
+using deviceAllocator = raft::mr::device::allocator;
+
+/**
+ * Performs a k-select across row partitioned index/distance
+ * matrices formatted like the following:
+ * row1: k0, k1, k2
+ * row2: k0, k1, k2
+ * row3: k0, k1, k2
+ * row1: k0, k1, k2
+ * row2: k0, k1, k2
+ * row3: k0, k1, k2
+ *
+ * etc...
+ *
+ * @tparam value_idx
+ * @tparam value_t
+ * @param inK
+ * @param inV
+ * @param outK
+ * @param outV
+ * @param n_samples
+ * @param n_parts
+ * @param k
+ * @param stream
+ * @param translations
+ */
 template <typename value_idx = int64_t, typename value_t = float>
 inline void knn_merge_parts(value_t *inK, value_idx *inV, value_t *outK,
                             value_idx *outV, size_t n_samples, int n_parts,
@@ -33,6 +60,34 @@ inline void knn_merge_parts(value_t *inK, value_idx *inV, value_t *outK,
                           translations);
 }
 
+/**
+ * Performs a k-select across column-partitioned index/distance
+ * matrices formatted like the following:
+ * row1: k0, k1, k2, k0, k1, k2
+ * row2: k0, k1, k2, k0, k1, k2
+ * row3: k0, k1, k2, k0, k1, k2
+ *
+ * etc...
+ *
+ * @tparam value_idx
+ * @tparam value_t
+ * @param inK
+ * @param inV
+ * @param n_rows
+ * @param n_cols
+ * @param outK
+ * @param outV
+ * @param select_min
+ * @param k
+ * @param stream
+ */
+template <typename value_idx = int, typename value_t = float>
+inline void select_k(value_t *inK, value_idx *inV, size_t n_rows, size_t n_cols,
+                     value_t *outK, value_idx *outV, bool select_min, int k,
+                     cudaStream_t stream) {
+  detail::select_k(inK, inV, n_rows, n_cols, outK, outV, select_min, k, stream);
+}
+
 /**
  * @brief Flat C++ API function to perform a brute force knn on
  * a series of input arrays and combine the results into a single
@@ -73,7 +128,6 @@ inline void brute_force_knn(
                                handle.get_num_internal_streams(), rowMajorIndex,
                                rowMajorQuery, translations, metric, metric_arg);
 }
-
 }  // namespace knn
 }  // namespace spatial
 }  // namespace raft
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index 42066061f3..43e1c65695 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -84,11 +84,12 @@ add_executable(test_raft
     test/sparse/norm.cu
     test/sparse/reduce.cu
     test/sparse/row_op.cu
-    test/sparse/selection.cu
     test/sparse/sort.cu
     test/sparse/symmetrize.cu
     test/spatial/knn.cu
     test/spatial/haversine.cu
+    test/spatial/ball_cover.cu
+    test/spatial/selection.cu
     test/spectral_matrix.cu
     test/stats/mean.cu
     test/stats/mean_center.cu
diff --git a/cpp/test/spatial/ball_cover.cu b/cpp/test/spatial/ball_cover.cu
new file mode 100644
index 0000000000..3c70ee8728
--- /dev/null
+++ b/cpp/test/spatial/ball_cover.cu
@@ -0,0 +1,271 @@
+/*
+ * Copyright (c) 2021, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <raft/cudart_utils.h>
+#include <raft/linalg/distance_type.h>
+#include <raft/spatial/knn/ball_cover.hpp>
+#include <raft/spatial/knn/detail/haversine_distance.cuh>
+#include <raft/spatial/knn/detail/knn_brute_force_faiss.cuh>
+#include <rmm/device_uvector.hpp>
+#include "../test_utils.h"
+#include "spatial_data.h"
+
+#include <thrust/transform.h>
+#include <rmm/exec_policy.hpp>
+
+#include <gtest/gtest.h>
+#include <cstdint>
+#include <iostream>
+#include <vector>
+
+namespace raft {
+namespace spatial {
+namespace knn {
+
+using namespace std;
+
+template <typename value_idx, typename value_t>
+__global__ void count_discrepancies_kernel(value_idx *actual_idx,
+                                           value_idx *expected_idx,
+                                           value_t *actual, value_t *expected,
+                                           uint32_t m, uint32_t n,
+                                           uint32_t *out, float thres = 1e-3) {
+  uint32_t row = blockDim.x * blockIdx.x + threadIdx.x;
+
+  int n_diffs = 0;
+  if (row < m) {
+    for (uint32_t i = 0; i < n; i++) {
+      value_t d = actual[row * n + i] - expected[row * n + i];
+      bool matches = fabsf(d) <= thres;
+      if (!matches) {
+        //          printf("row=%d, actual_idx=%ld, actual=%f, expected_id=%ld, expected=%f\n",
+        //                 row, actual_idx[row*n+i], actual[row*n+i], expected_idx[row*n+i], expected[row*n+i]);
+      }
+
+      n_diffs += !matches;
+      out[row] = n_diffs;
+    }
+  }
+}
+
+struct is_nonzero {
+  __host__ __device__ bool operator()(uint32_t &i) { return i > 0; }
+};
+
+template <typename value_idx, typename value_t>
+uint32_t count_discrepancies(value_idx *actual_idx, value_idx *expected_idx,
+                             value_t *actual, value_t *expected, uint32_t m,
+                             uint32_t n, uint32_t *out, cudaStream_t stream) {
+  uint32_t tpb = 256;
+  count_discrepancies_kernel<<<raft::ceildiv(m, tpb), tpb, 0, stream>>>(
+    actual_idx, expected_idx, actual, expected, m, n, out);
+
+  auto exec_policy = rmm::exec_policy(stream);
+
+  uint32_t result = thrust::count_if(exec_policy, out, out + m, is_nonzero());
+  return result;
+}
+
+struct ToRadians {
+  __device__ __host__ float operator()(float a) {
+    return a * (CUDART_PI_F / 180.0);
+  }
+};
+
+struct BallCoverInputs {
+  uint32_t k;
+  float weight;
+  raft::distance::DistanceType metric;
+};
+
+template <typename value_idx, typename value_t>
+class BallCoverKNNQueryTest : public ::testing::TestWithParam<BallCoverInputs> {
+ protected:
+  void basicTest() {
+    params = ::testing::TestWithParam<BallCoverInputs>::GetParam();
+    raft::handle_t handle;
+
+    uint32_t k = params.k;
+    float weight = params.weight;
+    auto metric = params.metric;
+
+    std::vector<value_t> h_train_inputs = spatial_data;
+
+    uint32_t n = h_train_inputs.size() / d;
+
+    rmm::device_uvector<value_idx> d_ref_I(n * k, handle.get_stream());
+    rmm::device_uvector<value_t> d_ref_D(n * k, handle.get_stream());
+
+    // Allocate input
+    rmm::device_uvector<value_t> d_train_inputs(n * d, handle.get_stream());
+    raft::update_device(d_train_inputs.data(), h_train_inputs.data(), n * d,
+                        handle.get_stream());
+
+    if (metric == raft::distance::DistanceType::Haversine) {
+      thrust::transform(handle.get_thrust_policy(), d_train_inputs.data(),
+                        d_train_inputs.data() + d_train_inputs.size(),
+                        d_train_inputs.data(), ToRadians());
+    }
+
+    cudaStream_t *int_streams = nullptr;
+    std::vector<int64_t> *translations = nullptr;
+
+    std::vector<float *> input_vec = {d_train_inputs.data()};
+    std::vector<uint32_t> sizes_vec = {n};
+
+    raft::spatial::knn::detail::brute_force_knn_impl<uint32_t, int64_t>(
+      input_vec, sizes_vec, d, d_train_inputs.data(), n, d_ref_I.data(),
+      d_ref_D.data(), k, handle.get_stream(), int_streams, 0, true, true,
+      translations, metric);
+
+    CUDA_CHECK(cudaStreamSynchronize(handle.get_stream()));
+
+    // Allocate predicted arrays
+    rmm::device_uvector<value_idx> d_pred_I(n * k, handle.get_stream());
+    rmm::device_uvector<value_t> d_pred_D(n * k, handle.get_stream());
+
+    BallCoverIndex<value_idx, value_t> index(handle, d_train_inputs.data(), n,
+                                             d, metric);
+
+    raft::spatial::knn::rbc_build_index(handle, index);
+    raft::spatial::knn::rbc_knn_query(handle, index, k, d_train_inputs.data(),
+                                      n, d_pred_I.data(), d_pred_D.data(), true,
+                                      weight);
+
+    CUDA_CHECK(cudaStreamSynchronize(handle.get_stream()));
+    // What we really want are for the distances to match exactly. The
+    // indices may or may not match exactly, depending upon the ordering which
+    // can be nondeterministic.
+
+    rmm::device_uvector<uint32_t> discrepancies(n, handle.get_stream());
+    thrust::fill(handle.get_thrust_policy(), discrepancies.data(),
+                 discrepancies.data() + discrepancies.size(), 0);
+    //
+    int res = count_discrepancies(d_ref_I.data(), d_pred_I.data(),
+                                  d_ref_D.data(), d_pred_D.data(), n, k,
+                                  discrepancies.data(), handle.get_stream());
+
+    ASSERT_TRUE(res == 0);
+  }
+
+  void SetUp() override {}
+
+  void TearDown() override {}
+
+ protected:
+  uint32_t d = 2;
+  BallCoverInputs params;
+};
+
+template <typename value_idx, typename value_t>
+class BallCoverAllKNNTest : public ::testing::TestWithParam<BallCoverInputs> {
+ protected:
+  void basicTest() {
+    params = ::testing::TestWithParam<BallCoverInputs>::GetParam();
+    raft::handle_t handle;
+
+    uint32_t k = params.k;
+    float weight = params.weight;
+    auto metric = params.metric;
+
+    std::vector<value_t> h_train_inputs = spatial_data;
+
+    uint32_t n = h_train_inputs.size() / d;
+
+    rmm::device_uvector<value_idx> d_ref_I(n * k, handle.get_stream());
+    rmm::device_uvector<value_t> d_ref_D(n * k, handle.get_stream());
+
+    // Allocate input
+    rmm::device_uvector<value_t> d_train_inputs(n * d, handle.get_stream());
+    raft::update_device(d_train_inputs.data(), h_train_inputs.data(), n * d,
+                        handle.get_stream());
+
+    if (metric == raft::distance::DistanceType::Haversine) {
+      thrust::transform(handle.get_thrust_policy(), d_train_inputs.data(),
+                        d_train_inputs.data() + d_train_inputs.size(),
+                        d_train_inputs.data(), ToRadians());
+    }
+
+    cudaStream_t *int_streams = nullptr;
+    std::vector<int64_t> *translations = nullptr;
+
+    std::vector<float *> input_vec = {d_train_inputs.data()};
+    std::vector<uint32_t> sizes_vec = {n};
+
+    raft::spatial::knn::detail::brute_force_knn_impl<uint32_t, int64_t>(
+      input_vec, sizes_vec, d, d_train_inputs.data(), n, d_ref_I.data(),
+      d_ref_D.data(), k, handle.get_stream(), int_streams, 0, true, true,
+      translations, metric);
+
+    CUDA_CHECK(cudaStreamSynchronize(handle.get_stream()));
+
+    // Allocate predicted arrays
+    rmm::device_uvector<value_idx> d_pred_I(n * k, handle.get_stream());
+    rmm::device_uvector<value_t> d_pred_D(n * k, handle.get_stream());
+
+    BallCoverIndex<value_idx, value_t> index(handle, d_train_inputs.data(), n,
+                                             d, metric);
+
+    raft::spatial::knn::rbc_all_knn_query(handle, index, k, d_pred_I.data(),
+                                          d_pred_D.data(), true, weight);
+
+    CUDA_CHECK(cudaStreamSynchronize(handle.get_stream()));
+    // What we really want are for the distances to match exactly. The
+    // indices may or may not match exactly, depending upon the ordering which
+    // can be nondeterministic.
+
+    rmm::device_uvector<uint32_t> discrepancies(n, handle.get_stream());
+    thrust::fill(handle.get_thrust_policy(), discrepancies.data(),
+                 discrepancies.data() + discrepancies.size(), 0);
+    //
+    uint32_t res = count_discrepancies(
+      d_ref_I.data(), d_pred_I.data(), d_ref_D.data(), d_pred_D.data(), n, k,
+      discrepancies.data(), handle.get_stream());
+    ASSERT_TRUE(res == 0);
+  }
+
+  void SetUp() override {}
+
+  void TearDown() override {}
+
+ protected:
+  uint32_t d = 2;
+  BallCoverInputs params;
+};
+
+typedef BallCoverAllKNNTest<int64_t, float> BallCoverAllKNNTestF;
+typedef BallCoverKNNQueryTest<int64_t, float> BallCoverKNNQueryTestF;
+
+const std::vector<BallCoverInputs> ballcover_inputs = {
+  {2, 1.0, raft::distance::DistanceType::Haversine},
+  {4, 1.0, raft::distance::DistanceType::Haversine},
+  {7, 1.0, raft::distance::DistanceType::Haversine},
+  {2, 1.0, raft::distance::DistanceType::L2SqrtUnexpanded},
+  {4, 1.0, raft::distance::DistanceType::L2SqrtUnexpanded},
+  {7, 1.0, raft::distance::DistanceType::L2SqrtUnexpanded},
+};
+
+INSTANTIATE_TEST_CASE_P(BallCoverAllKNNTest, BallCoverAllKNNTestF,
+                        ::testing::ValuesIn(ballcover_inputs));
+INSTANTIATE_TEST_CASE_P(BallCoverKNNQueryTest, BallCoverKNNQueryTestF,
+                        ::testing::ValuesIn(ballcover_inputs));
+
+TEST_P(BallCoverAllKNNTestF, Fit) { basicTest(); }
+TEST_P(BallCoverKNNQueryTestF, Fit) { basicTest(); }
+
+}  // namespace knn
+}  // namespace spatial
+}  // namespace raft
diff --git a/cpp/test/spatial/knn.cu b/cpp/test/spatial/knn.cu
index e4d05920c0..7c70f6ffac 100644
--- a/cpp/test/spatial/knn.cu
+++ b/cpp/test/spatial/knn.cu
@@ -85,15 +85,6 @@ class KNNTest : public ::testing::TestWithParam<KNNInputs> {
     build_expected_output<<<raft::ceildiv(rows_ * k_, 32), 32, 0, stream>>>(
       expected_labels_, rows_, k_, search_labels_);
 
-    raft::print_device_vector("Output indices: ", indices_, rows_ * k_,
-                              std::cout);
-    raft::print_device_vector("Output distances: ", distances_, rows_ * k_,
-                              std::cout);
-    raft::print_device_vector("Output labels: ", actual_labels_, rows_ * k_,
-                              std::cout);
-    raft::print_device_vector("Expected labels: ", expected_labels_, rows_ * k_,
-                              std::cout);
-
     ASSERT_TRUE(devArrMatch(expected_labels_, actual_labels_, rows_ * k_,
                             raft::Compare<int>()));
   }
diff --git a/cpp/test/sparse/selection.cu b/cpp/test/spatial/selection.cu
similarity index 93%
rename from cpp/test/sparse/selection.cu
rename to cpp/test/spatial/selection.cu
index 256ecfdfb7..a7f1af6034 100644
--- a/cpp/test/sparse/selection.cu
+++ b/cpp/test/spatial/selection.cu
@@ -20,10 +20,10 @@
 #include "../test_utils.h"
 
 #include <raft/sparse/utils.h>
-#include <raft/sparse/selection/selection.cuh>
+#include <raft/spatial/knn/knn.hpp>
 
 namespace raft {
-namespace sparse {
+namespace spatial {
 namespace selection {
 
 using namespace raft;
@@ -89,9 +89,8 @@ class SparseSelectionTest
 
     make_data();
 
-    raft::sparse::selection::select_k(dists, inds, n_rows, n_cols, out_dists,
-                                      out_indices, params.select_min, k,
-                                      stream);
+    raft::spatial::knn::select_k(dists, inds, n_rows, n_cols, out_dists,
+                                 out_indices, params.select_min, k, stream);
 
     CUDA_CHECK(cudaStreamSynchronize(stream));
   }
@@ -143,5 +142,5 @@ INSTANTIATE_TEST_CASE_P(SparseSelectionTest, SparseSelectionTestF,
                         ::testing::ValuesIn(inputs_i32_f));
 
 };  // end namespace selection
-};  // end namespace sparse
+};  // end namespace spatial
 };  // end namespace raft
diff --git a/cpp/test/spatial/spatial_data.h b/cpp/test/spatial/spatial_data.h
new file mode 100644
index 0000000000..87891164fc
--- /dev/null
+++ b/cpp/test/spatial/spatial_data.h
@@ -0,0 +1,27 @@
+#include <vector>
+
+namespace raft {
+namespace spatial {
+
+// Latitude and longitude coordinates of 51 US states / territories
+std::vector<float> spatial_data = {
+  63.588753, -154.493062, 32.318231, -86.902298,  35.20105,  -91.831833,
+  34.048928, -111.093731, 36.778261, -119.417932, 39.550051, -105.782067,
+  41.603221, -73.087749,  38.905985, -77.033418,  38.910832, -75.52767,
+  27.664827, -81.515754,  32.157435, -82.907123,  19.898682, -155.665857,
+  41.878003, -93.097702,  44.068202, -114.742041, 40.633125, -89.398528,
+  40.551217, -85.602364,  39.011902, -98.484246,  37.839333, -84.270018,
+  31.244823, -92.145024,  42.407211, -71.382437,  39.045755, -76.641271,
+  45.253783, -69.445469,  44.314844, -85.602364,  46.729553, -94.6859,
+  37.964253, -91.831833,  32.354668, -89.398528,  46.879682, -110.362566,
+  35.759573, -79.0193,    47.551493, -101.002012, 41.492537, -99.901813,
+  43.193852, -71.572395,  40.058324, -74.405661,  34.97273,  -105.032363,
+  38.80261,  -116.419389, 43.299428, -74.217933,  40.417287, -82.907123,
+  35.007752, -97.092877,  43.804133, -120.554201, 41.203322, -77.194525,
+  18.220833, -66.590149,  41.580095, -71.477429,  33.836081, -81.163725,
+  43.969515, -99.901813,  35.517491, -86.580447,  31.968599, -99.901813,
+  39.32098,  -111.093731, 37.431573, -78.656894,  44.558803, -72.577841,
+  47.751074, -120.740139, 43.78444,  -88.787868,  38.597626, -80.454903,
+  43.075968, -107.290284};
+};  // namespace spatial
+};  // namespace raft
\ No newline at end of file
diff --git a/cpp/test/test_utils.h b/cpp/test/test_utils.h
index b8e8fe3fa0..0f135c0121 100644
--- a/cpp/test/test_utils.h
+++ b/cpp/test/test_utils.h
@@ -21,6 +21,13 @@
 #include <memory>
 #include <raft/cuda_utils.cuh>
 
+#include <fstream>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
 namespace raft {
 
 template <typename T>
@@ -249,4 +256,39 @@ testing::AssertionResult match(const T expected, T actual, L eq_compare) {
     ms /= args.runs;                                    \
   } while (0)
 
+inline std::vector<float> read_csv(std::string filename,
+                                   bool skip_first_n_columns = 1) {
+  std::vector<float> result;
+  std::ifstream myFile(filename);
+  if (!myFile.is_open()) throw std::runtime_error("Could not open file");
+
+  std::string line, colname;
+  int val;
+
+  if (myFile.good()) {
+    std::getline(myFile, line);
+    std::stringstream ss(line);
+    while (std::getline(ss, colname, ',')) {
+    }
+  }
+
+  int n_lines = 0;
+  while (std::getline(myFile, line)) {
+    std::stringstream ss(line);
+    int colIdx = 0;
+    while (ss >> val) {
+      if (colIdx >= skip_first_n_columns) {
+        result.push_back(val);
+        if (ss.peek() == ',') ss.ignore();
+      }
+      colIdx++;
+    }
+    n_lines++;
+  }
+
+  printf("lines read: %d\n", n_lines);
+  myFile.close();
+  return result;
+}
+
 };  // end namespace raft