[INFRA] Import NVIDIA/CCCL upstream as optimization reference library

CCCL (CUDA C++ Core Libraries) provides:
- CUB: device/block/warp-level GPU primitives (reduce, scan, sort, topk)
- Thrust: high-level parallel algorithms (transform_reduce, sort, scan)
- libcudacxx: CUDA C++ standard library (atomics, barriers, memory)
- cudax: experimental features (memory resources, allocators)
- Tuning policies: per-SM hardware-specific algorithm parameters

Competition optimization vectors mapped to CCCL:
- Output TPS (83% weight): warp_reduce, block_reduce, device_topk
- Input TPS (14% weight): device_scan, block_load, prefetch
- Cache TPS (3% weight): prefix caching strategy patterns
- Memory (0.9 util): pooled/cached/buddy allocators

Source: https://github.com/NVIDIA/cccl (shallow clone, HEAD only)
License: Apache-2.0

cccl_upstream/thrust/examples/CMakeLists.txt

@@ -0,0 +1,179 @@
+# Setup FileCheck if requested and available:
+option(
+  THRUST_ENABLE_EXAMPLE_FILECHECK
+  "Check example output with the LLVM FileCheck utility."
+  OFF
+)
+set(filecheck_data_path "${Thrust_SOURCE_DIR}/internal/test")
+
+if (THRUST_ENABLE_EXAMPLE_FILECHECK)
+  # TODO this should go into a find module
+  find_program(
+    THRUST_FILECHECK_EXECUTABLE
+    DOC "Path to the LLVM FileCheck utility."
+    NAMES
+      FileCheck
+      FileCheck-3.9
+      FileCheck-4.0
+      FileCheck-5.0
+      FileCheck-6.0
+      FileCheck-7
+      FileCheck-8
+      FileCheck-9
+  )
+
+  if (NOT THRUST_FILECHECK_EXECUTABLE)
+    message(
+      FATAL_ERROR
+      "Could not find the LLVM FileCheck utility. Set THRUST_FILECHECK_EXECUTABLE manually, "
+      "or disable THRUST_ENABLE_EXAMPLE_FILECHECK."
+    )
+  endif()
+
+  execute_process(
+    COMMAND
+      "${THRUST_FILECHECK_EXECUTABLE}"
+      "${filecheck_data_path}/thrust.smoke.filecheck"
+    INPUT_FILE "${Thrust_SOURCE_DIR}/cmake/filecheck_smoke_test"
+    RESULT_VARIABLE exit_code
+  )
+
+  if (0 EQUAL exit_code)
+    message(STATUS "FileCheck enabled: ${THRUST_FILECHECK_EXECUTABLE}")
+  else()
+    message(
+      FATAL_ERROR
+      "The current THRUST_FILECHECK_EXECUTABLE ('${THRUST_FILECHECK_EXECUTABLE}') "
+      "does not seem to be a valid FileCheck executable."
+    )
+  endif()
+endif()
+
+## thrust_add_example
+#
+# Add an example executable and register it with ctest.
+#
+# target_name_var: Variable name to overwrite with the name of the example
+#   target. Useful for post-processing target information per-backend.
+# example_name: The name of the example minus "<config_prefix>.example." For
+#   instance, examples/vector.cu will be "vector", and examples/cuda/copy.cu
+#   would be "cuda.copy".
+# example_src: The source file that implements the example.
+# thrust_target: The reference thrust target with configuration information.
+#
+function(
+  thrust_add_example
+  target_name_var
+  example_name
+  example_src
+  thrust_target
+)
+  thrust_get_target_property(config_host ${thrust_target} HOST)
+  thrust_get_target_property(config_device ${thrust_target} DEVICE)
+  thrust_get_target_property(config_prefix ${thrust_target} PREFIX)
+
+  # Wrap the .cu file in .cpp for non-CUDA backends
+  if ("CUDA" STREQUAL "${config_device}")
+    set(real_example_src "${example_src}")
+  else()
+    thrust_wrap_cu_in_cpp(real_example_src "${example_src}" ${thrust_target})
+  endif()
+
+  # The actual name of the test's target:
+  set(example_target ${config_prefix}.example.${example_name})
+  set(${target_name_var} ${example_target} PARENT_SCOPE)
+
+  cccl_add_executable(${example_target} SOURCES "${real_example_src}")
+  target_link_libraries(${example_target} PRIVATE ${thrust_target})
+  target_include_directories(
+    ${example_target}
+    PRIVATE "${Thrust_SOURCE_DIR}/examples"
+  )
+
+  if ("CUDA" STREQUAL "${config_device}")
+    thrust_configure_cuda_target(${example_target} RDC ${THRUST_FORCE_RDC})
+  endif()
+
+  # We do not want to explicitly include `host_device.h` if not needed,
+  # so force include the file for non CUDA targets.
+  set(gx_msvc "$<CXX_COMPILER_ID:MSVC>")
+  set(gx_cxx "$<COMPILE_LANGUAGE:CXX>")
+  set(gx_cxx_msvc "$<AND:${gx_cxx},${gx_msvc}>")
+  set(gx_cxx_not_msvc "$<AND:${gx_cxx},$<NOT:${gx_msvc}>>")
+  target_compile_options(
+    ${example_target}
+    PRIVATE
+      "$<${gx_cxx_msvc}:SHELL:/FI include/host_device.h>"
+      "$<${gx_cxx_not_msvc}:SHELL:-include include/host_device.h>"
+  )
+  target_compile_definitions(
+    ${example_target}
+    PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>:THRUST_EXAMPLE_DEVICE_SIDE>
+  )
+
+  # Get the name of FileCheck input by stripping out the config name.
+  # (e.g. "thrust.cpp.cuda.cpp14.example.xxx" -> "thrust.example.xxx.filecheck")
+  string(
+    REPLACE
+    "${config_prefix}"
+    "thrust"
+    filecheck_reference_file
+    "${example_target}.filecheck"
+  )
+
+  add_test(
+    NAME ${example_target}
+    # gersemi: off
+    COMMAND
+      "${CMAKE_COMMAND}"
+        "-DEXAMPLE_EXECUTABLE=$<TARGET_FILE:${example_target}>"
+        "-DFILECHECK_ENABLED=${THRUST_ENABLE_EXAMPLE_FILECHECK}"
+        "-DFILECHECK_EXECUTABLE=${THRUST_FILECHECK_EXECUTABLE}"
+        "-DREFERENCE_FILE=${filecheck_data_path}/${filecheck_reference_file}"
+        -P "${Thrust_SOURCE_DIR}/cmake/ThrustRunExample.cmake"
+    # gersemi: on
+  )
+
+  # Run OMP/TBB tests in serial. Multiple OMP processes will massively
+  # oversubscribe the machine with GCC's OMP, and we want to test these with
+  # the full CPU available to each unit test.
+  set(config_systems ${config_host} ${config_device})
+  if (("OMP" IN_LIST config_systems) OR ("TBB" IN_LIST config_systems))
+    set_tests_properties(${example_target} PROPERTIES RUN_SERIAL ON)
+  endif()
+
+  # Check for per-example script. Script will be included in the current scope
+  # to allow custom property modifications.
+  get_filename_component(example_cmake_script "${example_src}" NAME_WLE)
+  set(
+    example_cmake_script
+    "${CMAKE_CURRENT_LIST_DIR}/${example_cmake_script}.cmake"
+  )
+  # Use a glob so we can detect if this changes:
+  file(
+    GLOB example_cmake_script
+    RELATIVE "${CMAKE_CURRENT_LIST_DIR}"
+    CONFIGURE_DEPENDS
+    "${example_cmake_script}"
+  )
+  if (example_cmake_script) # Will be non-empty only if the script exists
+    include("${example_cmake_script}")
+  endif()
+endfunction()
+
+file(
+  GLOB example_srcs
+  RELATIVE "${CMAKE_CURRENT_LIST_DIR}"
+  CONFIGURE_DEPENDS
+  *.cu
+  *.cpp
+)
+
+foreach (thrust_target IN LISTS THRUST_TARGETS)
+  foreach (example_src IN LISTS example_srcs)
+    get_filename_component(example_name "${example_src}" NAME_WLE)
+    thrust_add_example(example_target ${example_name} "${example_src}" ${thrust_target})
+  endforeach()
+endforeach()
+
+add_subdirectory(cuda)

cccl_upstream/thrust/examples/README.md

@@ -0,0 +1,13 @@
+Once Thrust has been installed, these example programs can be compiled
+directly with nvcc.  For example, the following command will compile the
+`norm` example.
+
+```bash
+$ nvcc norm.cu -o norm
+```
+
+These examples are also available online:
+https://github.com/NVIDIA/cccl/tree/main/thrust/examples
+
+For any serious experimentation, we recommend using CMake and [CCCL from GitHub](https://github.com/NVIDIA/cccl).
+We also provide consistent and convenient development environments as [devcontainers](../../.devcontainers/README.md).

cccl_upstream/thrust/examples/arbitrary_transformation.cu

@@ -0,0 +1,93 @@
+#include <thrust/detail/config.h>
+
+#include <thrust/device_vector.h>
+#include <thrust/for_each.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/zip_function.h>
+
+#include <iostream>
+
+// This example shows how to implement an arbitrary transformation of
+// the form output[i] = F(first[i], second[i], third[i], ... ).
+// In this example, we use a function with 3 inputs and 1 output.
+//
+// Iterators for all four vectors (3 inputs + 1 output) are "zipped"
+// into a single sequence of tuples with the zip_iterator.
+//
+// The arbitrary_functor receives a tuple that contains four elements,
+// which are references to values in each of the four sequences. When we
+// access the tuple 't' with the get() function,
+//      get<0>(t) returns a reference to A[i],
+//      get<1>(t) returns a reference to B[i],
+//      get<2>(t) returns a reference to C[i],
+//      get<3>(t) returns a reference to D[i].
+//
+// In this example, we can implement the transformation,
+//      D[i] = A[i] + B[i] * C[i];
+// by invoking arbitrary_functor() on each of the tuples using for_each.
+//
+// If we are using a functor that is not designed for zip iterators by taking a
+// tuple instead of individual arguments we can adapt this function using the
+// zip_function adaptor (C++11 only).
+//
+// Note that we could extend this example to implement functions with an
+// arbitrary number of input arguments by zipping more sequence together.
+// With the same approach we can have multiple *output* sequences, if we
+// wanted to implement something like
+//      D[i] = A[i] + B[i] * C[i];
+//      E[i] = A[i] + B[i] + C[i];
+//
+// The possibilities are endless! :)
+
+struct arbitrary_functor1
+{
+  template <typename Tuple>
+  __host__ __device__ void operator()(Tuple t)
+  {
+    // D[i] = A[i] + B[i] * C[i];
+    cuda::std::get<3>(t) = cuda::std::get<0>(t) + cuda::std::get<1>(t) * cuda::std::get<2>(t);
+  }
+};
+
+struct arbitrary_functor2
+{
+  __host__ __device__ void operator()(const float& a, const float& b, const float& c, float& d)
+  {
+    // D[i] = A[i] + B[i] * C[i];
+    d = a + b * c;
+  }
+};
+
+int main()
+{
+  // allocate and initialize
+  thrust::device_vector<float> A{3, 4, 0, 8, 2};
+  thrust::device_vector<float> B{6, 7, 2, 1, 8};
+  thrust::device_vector<float> C{2, 5, 7, 4, 3};
+  thrust::device_vector<float> D1(5);
+
+  // apply the transformation
+  thrust::for_each(thrust::make_zip_iterator(A.begin(), B.begin(), C.begin(), D1.begin()),
+                   thrust::make_zip_iterator(A.end(), B.end(), C.end(), D1.end()),
+                   arbitrary_functor1());
+
+  // print the output
+  std::cout << "Tuple functor" << '\n';
+  for (size_t i = 0; i < A.size(); i++)
+  {
+    std::cout << A[i] << " + " << B[i] << " * " << C[i] << " = " << D1[i] << '\n';
+  }
+
+  // apply the transformation using zip_function
+  thrust::device_vector<float> D2(5);
+  thrust::for_each(thrust::make_zip_iterator(A.begin(), B.begin(), C.begin(), D2.begin()),
+                   thrust::make_zip_iterator(A.end(), B.end(), C.end(), D2.end()),
+                   thrust::make_zip_function(arbitrary_functor2()));
+
+  // print the output
+  std::cout << "N-ary functor" << '\n';
+  for (size_t i = 0; i < A.size(); i++)
+  {
+    std::cout << A[i] << " + " << B[i] << " * " << C[i] << " = " << D2[i] << '\n';
+  }
+}

cccl_upstream/thrust/examples/basic_vector.cu

@@ -0,0 +1,40 @@
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+
+#include <iostream>
+
+int main()
+{
+  // H holds 4 integers
+  thrust::host_vector<int> H{14, 20, 38, 46};
+
+  // H.size() returns the size of vector H
+  std::cout << "H has size " << H.size() << '\n';
+
+  // print contents of H
+  for (size_t i = 0; i < H.size(); i++)
+  {
+    std::cout << "H[" << i << "] = " << H[i] << '\n';
+  }
+
+  // resize H
+  H.resize(2);
+
+  std::cout << "H now has size " << H.size() << '\n';
+
+  // Copy host_vector H to device_vector D
+  thrust::device_vector<int> D = H;
+
+  // elements of D can be modified
+  D[0] = 99;
+  D[1] = 88;
+
+  // print contents of D
+  for (size_t i = 0; i < D.size(); i++)
+  {
+    std::cout << "D[" << i << "] = " << D[i] << '\n';
+  }
+
+  // H and D are automatically deleted when the function returns
+  return 0;
+}

cccl_upstream/thrust/examples/bounding_box.cu

@@ -0,0 +1,101 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/random.h>
+#include <thrust/transform_reduce.h>
+
+#include <cuda/std/utility>
+
+#include <iostream>
+
+// This example shows how to compute a bounding box
+// for a set of points in two dimensions.
+
+struct point2d
+{
+  float x, y;
+
+  __host__ __device__ point2d()
+      : x(0)
+      , y(0)
+  {}
+
+  __host__ __device__ point2d(float _x, float _y)
+      : x(_x)
+      , y(_y)
+  {}
+};
+
+// bounding box type
+struct bbox
+{
+  // construct an empty box
+  bbox() = default;
+
+  // construct a box from a single point
+  __host__ __device__ bbox(const point2d& point)
+      : lower_left(point)
+      , upper_right(point)
+  {}
+
+  // construct a box from a single point
+  __host__ __device__ bbox& operator=(const point2d& point)
+  {
+    lower_left  = point;
+    upper_right = point;
+    return *this;
+  }
+
+  // construct a box from a pair of points
+  __host__ __device__ bbox(const point2d& ll, const point2d& ur)
+      : lower_left(ll)
+      , upper_right(ur)
+  {}
+
+  point2d lower_left, upper_right;
+};
+
+// reduce a pair of bounding boxes (a,b) to a bounding box containing a and b
+struct bbox_union
+{
+  __host__ __device__ bbox operator()(bbox a, bbox b)
+  {
+    // lower left corner
+    point2d ll(thrust::min(a.lower_left.x, b.lower_left.x), thrust::min(a.lower_left.y, b.lower_left.y));
+
+    // upper right corner
+    point2d ur(thrust::max(a.upper_right.x, b.upper_right.x), thrust::max(a.upper_right.y, b.upper_right.y));
+
+    return bbox(ll, ur);
+  }
+};
+
+int main()
+{
+  const size_t N = 40;
+
+  // allocate storage for points
+  thrust::device_vector<point2d> points(N);
+
+  // generate some random points in the unit square
+  thrust::default_random_engine rng;
+  thrust::uniform_real_distribution<float> u01(0.0f, 1.0f);
+  for (size_t i = 0; i < N; i++)
+  {
+    float x   = u01(rng);
+    float y   = u01(rng);
+    points[i] = point2d(x, y);
+  }
+
+  // initial bounding box contains first point
+  bbox init(points[0], points[0]);
+
+  // compute the bounding box for the point set
+  bbox result = thrust::reduce(points.begin(), points.end(), init, bbox_union{});
+
+  // print output
+  std::cout << "bounding box " << std::fixed;
+  std::cout << "(" << result.lower_left.x << "," << result.lower_left.y << ") ";
+  std::cout << "(" << result.upper_right.x << "," << result.upper_right.y << ")" << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/bucket_sort2d.cu

@@ -0,0 +1,106 @@
+#include <thrust/binary_search.h>
+#include <thrust/device_vector.h>
+#include <thrust/generate.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random.h>
+#include <thrust/sort.h>
+
+#include <iomanip>
+#include <iostream>
+
+// define a 2d float vector
+using vec2 = cuda::std::tuple<float, float>;
+
+// return a random vec2 in [0,1)^2
+vec2 make_random_vec2()
+{
+  static thrust::default_random_engine rng;
+  static thrust::uniform_real_distribution<float> u01(0.0f, 1.0f);
+  float x = u01(rng);
+  float y = u01(rng);
+  return vec2(x, y);
+}
+
+// hash a point in the unit square to the index of
+// the grid bucket that contains it
+struct point_to_bucket_index
+{
+  unsigned int width; // buckets in the x dimension (grid spacing = 1/width)
+  unsigned int height; // buckets in the y dimension (grid spacing = 1/height)
+
+  __host__ __device__ point_to_bucket_index(unsigned int width, unsigned int height)
+      : width(width)
+      , height(height)
+  {}
+
+  __host__ __device__ unsigned int operator()(const vec2& v) const
+  {
+    // find the raster indices of p's bucket
+    unsigned int x = static_cast<unsigned int>(cuda::std::get<0>(v) * static_cast<float>(width));
+    unsigned int y = static_cast<unsigned int>(cuda::std::get<1>(v) * static_cast<float>(height));
+
+    // return the bucket's linear index
+    return y * width + x;
+  }
+};
+
+int main()
+{
+  const size_t N = 1000000;
+
+  // allocate some random points in the unit square on the host
+  thrust::host_vector<vec2> h_points(N);
+  thrust::generate(h_points.begin(), h_points.end(), make_random_vec2);
+
+  // transfer to device
+  thrust::device_vector<vec2> points = h_points;
+
+  // allocate storage for a 2D grid
+  // of dimensions w x h
+  unsigned int w = 200, h = 100;
+
+  // the grid data structure keeps a range per grid bucket:
+  // each bucket_begin[i] indexes the first element of bucket i's list of points
+  // each bucket_end[i] indexes one past the last element of bucket i's list of points
+  thrust::device_vector<unsigned int> bucket_begin(w * h);
+  thrust::device_vector<unsigned int> bucket_end(w * h);
+
+  // allocate storage for each point's bucket index
+  thrust::device_vector<unsigned int> bucket_indices(N);
+
+  // transform the points to their bucket indices
+  thrust::transform(points.begin(), points.end(), bucket_indices.begin(), point_to_bucket_index(w, h));
+
+  // sort the points by their bucket index
+  thrust::sort_by_key(bucket_indices.begin(), bucket_indices.end(), points.begin());
+
+  // find the beginning of each bucket's list of points
+  thrust::counting_iterator<unsigned int> search_begin(0);
+  thrust::lower_bound(
+    bucket_indices.begin(),
+    bucket_indices.end(),
+    search_begin,
+    search_begin + static_cast<decltype(search_begin)::difference_type>(w) * h,
+    bucket_begin.begin());
+
+  // find the end of each bucket's list of points
+  thrust::upper_bound(
+    bucket_indices.begin(),
+    bucket_indices.end(),
+    search_begin,
+    search_begin + static_cast<decltype(search_begin)::difference_type>(w) * h,
+    bucket_end.begin());
+
+  // write out bucket (150, 50)'s list of points
+  unsigned int bucket_idx = 50 * w + 150;
+  std::cout << "bucket (150, 50)'s list of points:" << '\n';
+  std::cout << std::fixed << std::setprecision(6);
+  for (unsigned int point_idx = bucket_begin[bucket_idx]; point_idx != bucket_end[bucket_idx]; ++point_idx)
+  {
+    vec2 p = points[point_idx];
+    std::cout << "(" << cuda::std::get<0>(p) << "," << cuda::std::get<1>(p) << ")" << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/constant_iterator.cu

@@ -0,0 +1,26 @@
+#define CCCL_IGNORE_DEPRECATED_API
+
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform.h>
+
+#include <cuda/iterator>
+
+#include <iostream>
+#include <iterator>
+
+int main()
+{
+  thrust::device_vector<int> data{3, 7, 2, 5};
+
+  // add 10 to all values in data
+  thrust::transform(data.begin(), data.end(), cuda::constant_iterator<int>(10), data.begin(), cuda::std::plus<int>());
+
+  // data is now [13, 17, 12, 15]
+
+  // print result
+  thrust::copy(data.begin(), data.end(), std::ostream_iterator<int>(std::cout, "\n"));
+
+  return 0;
+}

cccl_upstream/thrust/examples/counting_iterator.cu

@@ -0,0 +1,34 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+
+#include <iostream>
+#include <iterator>
+
+int main()
+{
+  // this example computes indices for all the nonzero values in a sequence
+
+  // sequence of zero and nonzero values
+  thrust::device_vector<int> stencil{0, 1, 1, 0, 0, 1, 0, 1};
+
+  // storage for the nonzero indices
+  thrust::device_vector<int> indices(8);
+
+  // counting iterators define a sequence [0, 8)
+  thrust::counting_iterator<int> first(0);
+  thrust::counting_iterator<int> last = first + 8;
+
+  // compute indices of nonzero elements
+  using IndexIterator = thrust::device_vector<int>::iterator;
+
+  IndexIterator indices_end = thrust::copy_if(first, last, stencil.begin(), indices.begin(), cuda::std::identity{});
+  // indices now contains [1,2,5,7]
+
+  // print result
+  std::cout << "found " << cuda::std::distance(indices.begin(), indices_end) << " nonzero values at indices:\n";
+  thrust::copy(indices.begin(), indices_end, std::ostream_iterator<int>(std::cout, "\n"));
+
+  return 0;
+}

cccl_upstream/thrust/examples/cpp_integration/README

@@ -0,0 +1,20 @@
+This example shows how to link a Thrust program contained in
+a .cu file with a C++ program contained in a .cpp file.  Note
+that device_vector only appears in the .cu file while host_vector
+appears in both.  This relects the fact that algorithms on device
+vectors are only available when the contents of the program are
+located in a .cu file and compiled with the nvcc compiler.
+
+On a Linux system where Thrust is installed in the default location
+we can use the following procedure to compile the two parts of the
+program and link them together.
+
+  $ nvcc -O2 -c device.cu
+  $ g++  -O2 -c host.cpp   -I/usr/local/cuda/include/
+  $ nvcc -o tester device.o host.o
+
+Alternatively, we can use g++ to perform final linking step.
+
+  $ nvcc -O2 -c device.cu
+  $ g++  -O2 -c host.cpp   -I/usr/local/cuda/include/
+  $ g++ -o tester device.o host.o -L/usr/local/cuda/lib64 -lcudart

cccl_upstream/thrust/examples/cpp_integration/device.cu

@@ -0,0 +1,17 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/sort.h>
+
+#include "device.h"
+
+void sort_on_device(thrust::host_vector<int>& h_vec)
+{
+  // transfer data to the device
+  thrust::device_vector<int> d_vec = h_vec;
+
+  // sort data on the device
+  thrust::sort(d_vec.begin(), d_vec.end());
+
+  // transfer data back to host
+  thrust::copy(d_vec.begin(), d_vec.end(), h_vec.begin());
+}

cccl_upstream/thrust/examples/cpp_integration/device.h

@@ -0,0 +1,6 @@
+#pragma once
+
+#include <thrust/host_vector.h>
+
+// function prototype
+void sort_on_device(thrust::host_vector<int>& V);

cccl_upstream/thrust/examples/cpp_integration/host.cpp

@@ -0,0 +1,27 @@
+#include <thrust/generate.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/sort.h>
+
+#include <cstdlib>
+#include <iostream>
+#include <iterator>
+
+// defines the function prototype
+#include "device.h"
+
+int main()
+{
+  // generate 20 random numbers on the host
+  thrust::host_vector<int> h_vec(20);
+  thrust::default_random_engine rng;
+  thrust::generate(h_vec.begin(), h_vec.end(), rng);
+
+  // interface to CUDA code
+  sort_on_device(h_vec);
+
+  // print sorted array
+  thrust::copy(h_vec.begin(), h_vec.end(), std::ostream_iterator<int>(std::cout, "\n"));
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/CMakeLists.txt

@@ -0,0 +1,20 @@
+file(
+  GLOB example_srcs
+  RELATIVE "${CMAKE_CURRENT_LIST_DIR}"
+  CONFIGURE_DEPENDS
+  *.cu
+  *.cpp
+)
+
+foreach (thrust_target IN LISTS THRUST_TARGETS)
+  thrust_get_target_property(config_device ${thrust_target} DEVICE)
+  if (NOT config_device STREQUAL "CUDA")
+    continue()
+  endif()
+
+  foreach (example_src IN LISTS example_srcs)
+    get_filename_component(example_name "${example_src}" NAME_WLE)
+    string(PREPEND example_name "cuda.")
+    thrust_add_example(example_target ${example_name} "${example_src}" ${thrust_target})
+  endforeach()
+endforeach()

cccl_upstream/thrust/examples/cuda/async_reduce.cu

@@ -0,0 +1,78 @@
+#include <thrust/detail/config.h>
+
+#include <thrust/device_vector.h>
+#include <thrust/reduce.h>
+#include <thrust/system/cuda/execution_policy.h>
+
+#include <cassert>
+#include <future>
+
+// This example demonstrates two ways to achieve algorithm invocations that are asynchronous with
+// the calling thread.
+//
+// The first method wraps a call to thrust::reduce inside a __global__ function. Since __global__ function
+// launches are asynchronous with the launching thread, this achieves asynchrony. The result of the reduction
+// is stored to a pointer to CUDA global memory. The calling thread waits for the result of the reduction to
+// be ready by synchronizing with the CUDA stream on which the __global__ function is launched.
+//
+// The second method uses the C++11 library function, std::async, to create concurrency. The lambda function
+// given to std::async returns the result of thrust::reduce to a std::future. The calling thread can use the
+// std::future to wait for the result of the reduction. This method requires a compiler which supports
+// C++11-capable language and library constructs.
+
+#ifdef THRUST_EXAMPLE_DEVICE_SIDE
+template <typename Iterator, typename T, typename BinaryOperation, typename Pointer>
+__global__ void reduce_kernel(Iterator first, Iterator last, T init, BinaryOperation binary_op, Pointer result)
+{
+  *result = thrust::reduce(thrust::cuda::par, first, last, init, binary_op);
+}
+#endif
+
+int main()
+{
+  size_t n = 1 << 20;
+  thrust::device_vector<unsigned int> data(n, 1);
+  thrust::device_vector<unsigned int> result(1, 0);
+
+  // method 1: call thrust::reduce from an asynchronous CUDA kernel launch
+
+  // create a CUDA stream
+  cudaStream_t s;
+  cudaStreamCreate(&s);
+
+  // launch a CUDA kernel with only 1 thread on our stream
+#ifdef THRUST_EXAMPLE_DEVICE_SIDE
+  reduce_kernel<<<1, 1, 0, s>>>(data.begin(), data.end(), 0, cuda::std::plus<int>(), result.data());
+#else
+  result[0] = thrust::reduce(thrust::cuda::par, data.begin(), data.end(), 0, cuda::std::plus<int>());
+#endif
+
+  // wait for the stream to finish
+  cudaStreamSynchronize(s);
+
+  // our result should be ready
+  assert(result[0] == n);
+
+  cudaStreamDestroy(s);
+
+  // reset the result
+  result[0] = 0;
+
+  // method 2: use std::async to create asynchrony
+  // copy all the algorithm parameters
+  auto begin        = data.begin();
+  auto end          = data.end();
+  unsigned int init = 0;
+  auto binary_op    = cuda::std::plus<unsigned int>();
+
+  // std::async captures the algorithm parameters by value
+  // use std::launch::async to ensure the creation of a new thread
+  std::future<unsigned int> future_result = std::async(std::launch::async, [=] {
+    return thrust::reduce(begin, end, init, binary_op);
+  });
+
+  // wait on the result and check that it is correct
+  assert(future_result.get() == n);
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/custom_temporary_allocation.cu

@@ -0,0 +1,154 @@
+#include <thrust/generate.h>
+#include <thrust/host_vector.h>
+#include <thrust/sort.h>
+#include <thrust/system/cuda/execution_policy.h>
+#include <thrust/system/cuda/vector.h>
+
+#include <cuda/std/utility>
+
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <stdexcept>
+
+// This example demonstrates how to control how Thrust allocates temporary
+// storage during algorithms such as thrust::sort. The idea will be to create a
+// simple cache of allocations to search when temporary storage is requested.
+// If a hit is found in the cache, we quickly return the cached allocation
+// instead of resorting to the more expensive thrust::cuda::malloc.
+
+// Note: Thrust now has its own caching allocator layer; if you just need a
+// caching allocator, you ought to use that. This example is still useful
+// as a demonstration of how to use a Thrust custom allocator.
+
+// Note: this implementation cached_allocator is not thread-safe. If multiple
+// (host) threads use the same cached_allocator then they should gain exclusive
+// access to the allocator before accessing its methods.
+
+struct not_my_pointer_exception : std::exception
+{
+  explicit not_my_pointer_exception(void* p)
+  {
+    std::stringstream s;
+    s << "Pointer `" << p << "` was not allocated by this allocator.";
+    message = s.str();
+  }
+
+  const char* what() const noexcept override
+  {
+    return message.c_str();
+  }
+
+private:
+  std::string message;
+};
+
+// A simple allocator for caching cudaMalloc allocations.
+// A minimum allocator needs to provide at least `value_type`, `allocate` and `deallocate`.
+struct cached_allocator
+{
+  using value_type = char;
+
+  ~cached_allocator()
+  {
+    free_all();
+  }
+
+  char* allocate(std::ptrdiff_t num_bytes)
+  {
+    std::cout << "cached_allocator::allocate(): num_bytes == " << num_bytes << '\n';
+
+    char* result = nullptr;
+
+    // Search the cache for a free block.
+    auto free_block_it = free_blocks.find(num_bytes);
+    if (free_block_it != free_blocks.end())
+    {
+      std::cout << "cached_allocator::allocate(): found a free block" << '\n';
+      result = free_block_it->second;
+      free_blocks.erase(free_block_it);
+    }
+    else
+    {
+      // No allocation of the right size exists, so create a new one with `thrust::cuda::malloc`.
+      std::cout << "cached_allocator::allocate(): allocating new block" << '\n';
+      // Allocate memory and convert the resulting `thrust::cuda::pointer` to a raw pointer.
+      result = thrust::cuda::malloc<char>(num_bytes).get();
+    }
+
+    // Insert the allocated pointer into the `allocated_blocks` map.
+    allocated_blocks.insert(std::pair{result, num_bytes});
+
+    return result;
+  }
+
+  void deallocate(char* ptr, size_t)
+  {
+    std::cout << "cached_allocator::deallocate(): ptr == " << reinterpret_cast<void*>(ptr) << '\n';
+
+    // Erase the allocated block from the allocated blocks map.
+    auto it = allocated_blocks.find(ptr);
+    if (it == allocated_blocks.end())
+    {
+      throw not_my_pointer_exception(ptr);
+    }
+
+    const std::ptrdiff_t num_bytes = it->second;
+    allocated_blocks.erase(it);
+
+    // Insert the block into the free blocks map.
+    free_blocks.insert(std::make_pair(num_bytes, ptr));
+  }
+
+private:
+  std::multimap<std::ptrdiff_t, char*> free_blocks;
+  std::map<char*, std::ptrdiff_t> allocated_blocks;
+
+  void free_all()
+  {
+    std::cout << "cached_allocator::free_all()" << '\n';
+
+    // Deallocate all outstanding blocks in both lists.
+    for (auto [bytes, ptr] : free_blocks)
+    {
+      thrust::cuda::free(thrust::cuda::pointer<char>(ptr));
+    }
+
+    for (auto [ptr, bytes] : allocated_blocks)
+    {
+      thrust::cuda::free(thrust::cuda::pointer<char>(ptr));
+    }
+  }
+};
+
+int main()
+{
+  std::size_t num_elements = 32768;
+
+  thrust::host_vector<int> h_input(num_elements);
+
+  // Generate random input.
+  thrust::generate(h_input.begin(), h_input.end(), rand);
+
+  thrust::cuda::vector<int> d_input = h_input;
+  thrust::cuda::vector<int> d_result(num_elements);
+
+  std::size_t num_trials = 5;
+
+  cached_allocator alloc;
+
+  for (std::size_t i = 0; i < num_trials; ++i)
+  {
+    d_result = d_input;
+
+    // Pass alloc to execution policy cuda::par. It will handle allocations needed inside sort.
+    thrust::sort(thrust::cuda::par(alloc), d_result.begin(), d_result.end());
+
+    // Ensure the result is sorted.
+    assert(thrust::is_sorted(d_result.begin(), d_result.end()));
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/explicit_cuda_stream.cu

@@ -0,0 +1,77 @@
+#include <thrust/device_vector.h>
+#include <thrust/execution_policy.h> // For thrust::device
+#include <thrust/reduce.h>
+#include <thrust/sequence.h>
+
+#include <iostream>
+
+#include <cuda_runtime.h>
+
+// This example shows how to execute a Thrust device algorithm on an explicit
+// CUDA stream. The simple program below fills a vector with the numbers
+// [0, 1000) (thrust::sequence) and then performs a scan operation
+// (thrust::inclusive_scan) on them. Both algorithms are executed on the same
+// custom CUDA stream using the CUDA execution policies.
+//
+// Thrust provides two execution policies that accept CUDA streams that differ
+// in when/if they synchronize the stream:
+// 1. thrust::cuda::par.on(stream)
+//      - `stream` will *always* be synchronized before an algorithm returns.
+//      - This is the default `thrust::device` policy when compiling with the
+//        CUDA device backend.
+// 2. thrust::cuda::par_nosync.on(stream)
+//      - `stream` will only be synchronized when necessary for correctness
+//        (e.g., returning a result from `thrust::reduce`). This is a hint that
+//        may be ignored by an algorithm's implementation.
+
+int main()
+{
+  thrust::device_vector<int> d_vec(1000);
+
+  // Create the stream:
+  cudaStream_t custom_stream;
+  cudaError_t err = cudaStreamCreate(&custom_stream);
+  if (err != cudaSuccess)
+  {
+    std::cerr << "Error creating stream: " << cudaGetErrorString(err) << "\n";
+    return 1;
+  }
+
+  // Construct a new `nosync` execution policy with the custom stream
+  auto nosync_exec_policy = thrust::cuda::par_nosync.on(custom_stream);
+
+  // Fill the vector with sequential data.
+  // This will execute using the custom stream and the stream will *not* be
+  // synchronized before the function returns, meaning asynchronous work may
+  // still be executing after returning and the contents of `d_vec` are
+  // undefined. Synchronization is not needed here because the following
+  // `inclusive_scan` is executed on the same stream and is therefore guaranteed
+  // to be ordered after the `sequence`
+  thrust::sequence(nosync_exec_policy, d_vec.begin(), d_vec.end());
+
+  // Construct a new *synchronous* execution policy with the same custom stream
+  auto sync_exec_policy = thrust::cuda::par.on(custom_stream);
+
+  // Compute in-place inclusive sum scan of data in the vector.
+  // This also executes in the custom stream, but the execution policy ensures
+  // the stream is synchronized before the algorithm returns. This guarantees
+  // there is no pending asynchronous work and the contents of `d_vec` are
+  // immediately accessible.
+  thrust::inclusive_scan(sync_exec_policy, d_vec.cbegin(), d_vec.cend(), d_vec.begin());
+
+  // This access is only valid because the stream has been synchronized
+  int sum = d_vec.back();
+
+  // Free the stream:
+  err = cudaStreamDestroy(custom_stream);
+  if (err != cudaSuccess)
+  {
+    std::cerr << "Error destroying stream: " << cudaGetErrorString(err) << "\n";
+    return 1;
+  }
+
+  // Print the sum:
+  std::cout << "sum is " << sum << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/global_device_vector.cu

@@ -0,0 +1,34 @@
+#include <thrust/detail/config.h>
+
+#include <thrust/device_vector.h>
+
+// If you create a global `thrust::device_vector` with the default allocator,
+// you'll get an error during program termination when the memory of the vector
+// is freed, as the CUDA runtime cannot be used during program termination.
+//
+// To get around this, you can create your own allocator which ignores
+// deallocation failures that occur because the CUDA runtime is shut down.
+
+extern "C" cudaError_t cudaFreeIgnoreShutdown(void* ptr)
+{
+  cudaError_t const err = cudaFree(ptr);
+  if (cudaSuccess == err || cudaErrorCudartUnloading == err)
+  {
+    return cudaSuccess;
+  }
+  return err;
+}
+
+using device_ignore_shutdown_memory_resource =
+  thrust::system::cuda::detail::cuda_memory_resource<cudaMalloc, cudaFreeIgnoreShutdown, thrust::cuda::pointer<void>>;
+
+template <typename T>
+using device_ignore_shutdown_allocator =
+  thrust::mr::stateless_resource_allocator<T, thrust::device_ptr_memory_resource<device_ignore_shutdown_memory_resource>>;
+
+thrust::device_vector<double, device_ignore_shutdown_allocator<double>> d_vec;
+
+int main()
+{
+  d_vec.resize(25);
+}

cccl_upstream/thrust/examples/cuda/range_view.cu

@@ -0,0 +1,209 @@
+#include <thrust/device_vector.h>
+#include <thrust/execution_policy.h>
+#include <thrust/for_each.h>
+#include <thrust/iterator/counting_iterator.h>
+
+#include <cuda/std/iterator>
+
+#include <iostream>
+
+// This example demonstrates the use of a view: a non-owning wrapper for an
+// iterator range which presents a container-like interface to the user.
+//
+// For example, a view of a device_vector's data can be helpful when we wish to
+// access that data from a device function. Even though device_vectors are not
+// accessible from device functions, the range_view class allows us to access
+// and manipulate its data as if we were manipulating a real container.
+
+template <class Iterator>
+class range_view
+{
+public:
+  using iterator        = Iterator;
+  using value_type      = typename cuda::std::iterator_traits<iterator>::value_type;
+  using pointer         = typename cuda::std::iterator_traits<iterator>::pointer;
+  using difference_type = typename cuda::std::iterator_traits<iterator>::difference_type;
+  using reference       = typename cuda::std::iterator_traits<iterator>::reference;
+
+private:
+  const iterator first;
+  const iterator last;
+
+public:
+  __host__ __device__ range_view(Iterator first, Iterator last)
+      : first(first)
+      , last(last)
+  {}
+  ~range_view() = default;
+
+  __host__ __device__ difference_type size() const
+  {
+    return cuda::std::distance(first, last);
+  }
+
+  __host__ __device__ reference operator[](difference_type n)
+  {
+    return *(first + n);
+  }
+  __host__ __device__ const reference operator[](difference_type n) const
+  {
+    return *(first + n);
+  }
+
+  __host__ __device__ iterator begin()
+  {
+    return first;
+  }
+  __host__ __device__ const iterator cbegin() const
+  {
+    return first;
+  }
+  __host__ __device__ iterator end()
+  {
+    return last;
+  }
+  __host__ __device__ const iterator cend() const
+  {
+    return last;
+  }
+
+  __host__ __device__ cuda::std::reverse_iterator<iterator> rbegin()
+  {
+    return cuda::std::reverse_iterator<iterator>(end());
+  }
+  __host__ __device__ const cuda::std::reverse_iterator<const iterator> crbegin() const
+  {
+    return cuda::std::reverse_iterator<const iterator>(cend());
+  }
+  __host__ __device__ cuda::std::reverse_iterator<iterator> rend()
+  {
+    return cuda::std::reverse_iterator<iterator>(begin());
+  }
+  __host__ __device__ const cuda::std::reverse_iterator<const iterator> crend() const
+  {
+    return cuda::std::reverse_iterator<const iterator>(cbegin());
+  }
+  __host__ __device__ reference front()
+  {
+    return *begin();
+  }
+  __host__ __device__ const reference front() const
+  {
+    return *cbegin();
+  }
+
+  __host__ __device__ reference back()
+  {
+    return *end();
+  }
+  __host__ __device__ const reference back() const
+  {
+    return *cend();
+  }
+
+  __host__ __device__ bool empty() const
+  {
+    return size() == 0;
+  }
+};
+
+// This helper function creates a range_view from iterator and the number of
+// elements
+template <class Iterator, class Size>
+range_view<Iterator> __host__ __device__ make_range_view(Iterator first, Size n)
+{
+  return range_view<Iterator>(first, first + n);
+}
+
+// This helper function creates a range_view from a pair of iterators
+template <class Iterator>
+range_view<Iterator> __host__ __device__ make_range_view(Iterator first, Iterator last)
+{
+  return range_view<Iterator>(first, last);
+}
+
+// This helper function creates a range_view from a Vector
+template <class Vector>
+range_view<typename Vector::iterator> __host__ make_range_view(Vector& v)
+{
+  return range_view<typename Vector::iterator>(v.begin(), v.end());
+}
+
+// This saxpy functor stores view of X, Y, Z array, and accesses them in
+// vector-like way
+template <class View1, class View2, class View3>
+struct saxpy_functor
+{
+  const float a;
+  View1 x;
+  View2 y;
+  View3 z;
+
+  __host__ __device__ saxpy_functor(float _a, View1 _x, View2 _y, View3 _z)
+      : a(_a)
+      , x(_x)
+      , y(_y)
+      , z(_z)
+  {}
+
+  __host__ __device__ void operator()(int i)
+  {
+    z[i] = a * x[i] + y[i];
+  }
+};
+
+// saxpy function, which can either be called form host or device
+// The views are passed by value
+template <class View1, class View2, class View3>
+__host__ __device__ void saxpy(float A, View1 X, View2 Y, View3 Z)
+{
+  // Z = A * X + Y
+  const int size = X.size();
+  thrust::for_each(thrust::device,
+                   thrust::make_counting_iterator(0),
+                   thrust::make_counting_iterator(size),
+                   saxpy_functor<View1, View2, View3>(A, X, Y, Z));
+}
+
+struct f1
+{
+  __host__ __device__ float operator()(float x) const
+  {
+    return x * 3;
+  }
+};
+
+int main()
+{
+  using std::cout;
+
+  // initialize host arrays
+  float x[4] = {1.0, 1.0, 1.0, 1.0};
+  float y[4] = {1.0, 2.0, 3.0, 4.0};
+  float z[4] = {0.0};
+
+  thrust::device_vector<float> X(x, x + 4);
+  thrust::device_vector<float> Y(y, y + 4);
+  thrust::device_vector<float> Z(z, z + 4);
+
+  saxpy(
+    2.0,
+
+    // make a range view of a pair of transform_iterators
+    make_range_view(thrust::make_transform_iterator(X.cbegin(), f1()), thrust::make_transform_iterator(X.cend(), f1())),
+
+    // range view of normal_iterators
+    make_range_view(Y.begin(), cuda::std::distance(Y.begin(), Y.end())),
+
+    // range view of naked pointers
+    make_range_view(Z.data().get(), 4));
+
+  // print values from original device_vector<float> Z
+  // to ensure that range view was mapped to this vector
+  for (std::size_t i = 0, n = Z.size(); i < n; ++i)
+  {
+    cout << "z[" << i << "]= " << Z[i] << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/unwrap_pointer.cu

@@ -0,0 +1,31 @@
+#include <thrust/device_free.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_ptr.h>
+#include <thrust/device_vector.h>
+
+#include <cuda.h>
+
+int main()
+{
+  size_t N = 10;
+
+  // create a device_ptr
+  thrust::device_ptr<int> dev_ptr = thrust::device_malloc<int>(N);
+
+  // extract raw pointer from device_ptr
+  int* raw_ptr = thrust::raw_pointer_cast(dev_ptr);
+
+  // use raw_ptr in CUDA API functions
+  cudaMemset(raw_ptr, 0, N * sizeof(int));
+
+  // free memory
+  thrust::device_free(dev_ptr);
+
+  // we can use the same approach for device_vector
+  thrust::device_vector<int> d_vec(N);
+
+  // note: d_vec.data() returns a device_ptr
+  raw_ptr = thrust::raw_pointer_cast(d_vec.data());
+
+  return 0;
+}

cccl_upstream/thrust/examples/cuda/wrap_pointer.cu

@@ -0,0 +1,27 @@
+#include <thrust/device_ptr.h>
+#include <thrust/fill.h>
+
+#include <cuda.h>
+
+int main()
+{
+  size_t N = 10;
+
+  // obtain raw pointer to device memory
+  int* raw_ptr;
+  cudaMalloc((void**) &raw_ptr, N * sizeof(int));
+
+  // wrap raw pointer with a device_ptr
+  thrust::device_ptr<int> dev_ptr = thrust::device_pointer_cast(raw_ptr);
+
+  // use device_ptr in Thrust algorithms
+  thrust::fill(dev_ptr, dev_ptr + static_cast<std::ptrdiff_t>(N), (int) 0);
+
+  // access device memory transparently through device_ptr
+  dev_ptr[0] = 1;
+
+  // free memory
+  cudaFree(raw_ptr);
+
+  return 0;
+}

cccl_upstream/thrust/examples/device_ptr.cu

@@ -0,0 +1,44 @@
+#include <thrust/device_free.h>
+#include <thrust/device_malloc.h>
+#include <thrust/device_ptr.h>
+#include <thrust/reduce.h>
+#include <thrust/sequence.h>
+
+#include <cassert>
+#include <iostream>
+
+int main()
+{
+  // allocate memory buffer to store 10 integers on the device
+  thrust::device_ptr<int> d_ptr = thrust::device_malloc<int>(10);
+
+  // device_ptr supports pointer arithmetic
+  thrust::device_ptr<int> first = d_ptr;
+  thrust::device_ptr<int> last  = d_ptr + 10;
+  std::cout << "device array contains " << cuda::std::distance(first, last) << " values\n";
+
+  // algorithms work as expected
+  thrust::sequence(first, last);
+  std::cout << "sum of values is " << thrust::reduce(first, last) << "\n";
+
+  // device memory can be read and written transparently
+  d_ptr[0] = 10;
+  d_ptr[1] = 11;
+  d_ptr[2] = d_ptr[0] + d_ptr[1];
+
+  // device_ptr can be converted to a "raw" pointer for use in other APIs and kernels, etc.
+  int* raw_ptr = thrust::raw_pointer_cast(d_ptr);
+
+  // note: raw_ptr cannot necessarily be accessed by the host!
+
+  // conversely, raw pointers can be wrapped
+  [[maybe_unused]] thrust::device_ptr<int> wrapped_ptr = thrust::device_pointer_cast(raw_ptr);
+
+  // back to where we started
+  assert(wrapped_ptr == d_ptr);
+
+  // deallocate device memory
+  thrust::device_free(d_ptr);
+
+  return 0;
+}

cccl_upstream/thrust/examples/discrete_voronoi.cu

@@ -0,0 +1,259 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random.h>
+
+#include <cuda/std/tuple>
+
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "include/timer.h"
+
+// Compute an approximate Voronoi Diagram with a Jump Flooding Algorithm (JFA)
+//
+// References
+//   http://en.wikipedia.org/wiki/Voronoi_diagram
+//   http://www.comp.nus.edu.sg/~tants/jfa.html
+//   http://www.utdallas.edu/~guodongrong/Papers/Dissertation.pdf
+//
+// Thanks to David Coeurjolly for contributing this example
+
+// minFunctor
+// Tuple  = <seeds,seeds + k,seeds + m*k, seeds - k,
+//           seeds - m*k, seeds+ k+m*k,seeds + k-m*k,
+//           seeds- k+m*k,seeds - k+m*k, i>
+struct voronoi_site_selector
+{
+  int m, n, k;
+
+  __host__ __device__ voronoi_site_selector(int m, int n, int k)
+      : m(m)
+      , n(n)
+      , k(k)
+  {}
+
+  // To decide I have to change my current Voronoi site
+  __host__ __device__ int minVoro(int x_i, int y_i, int p, int q)
+  {
+    if (q == m * n)
+    {
+      return p;
+    }
+
+    // coordinates of points p and q
+    int y_q = q / m;
+    int x_q = q - y_q * m;
+    int y_p = p / m;
+    int x_p = p - y_p * m;
+
+    // squared distances
+    int d_iq = (x_i - x_q) * (x_i - x_q) + (y_i - y_q) * (y_i - y_q);
+    int d_ip = (x_i - x_p) * (x_i - x_p) + (y_i - y_p) * (y_i - y_p);
+
+    if (d_iq < d_ip)
+    {
+      return q; // q is closer
+    }
+    else
+    {
+      return p;
+    }
+  }
+
+  // For each point p+{-k,0,k}, we keep the Site with minimum distance
+  template <typename Tuple>
+  __host__ __device__ int operator()(const Tuple& t)
+  {
+    // Current point and site
+    int i = cuda::std::get<9>(t);
+    int v = cuda::std::get<0>(t);
+
+    // Current point coordinates
+    int y = i / m;
+    int x = i - y * m;
+
+    if (x >= k)
+    {
+      v = minVoro(x, y, v, cuda::std::get<3>(t));
+
+      if (y >= k)
+      {
+        v = minVoro(x, y, v, cuda::std::get<8>(t));
+      }
+
+      if (y + k < n)
+      {
+        v = minVoro(x, y, v, cuda::std::get<7>(t));
+      }
+    }
+
+    if (x + k < m)
+    {
+      v = minVoro(x, y, v, cuda::std::get<1>(t));
+
+      if (y >= k)
+      {
+        v = minVoro(x, y, v, cuda::std::get<6>(t));
+      }
+      if (y + k < n)
+      {
+        v = minVoro(x, y, v, cuda::std::get<5>(t));
+      }
+    }
+
+    if (y >= k)
+    {
+      v = minVoro(x, y, v, cuda::std::get<4>(t));
+    }
+    if (y + k < n)
+    {
+      v = minVoro(x, y, v, cuda::std::get<2>(t));
+    }
+
+    // global return
+    return v;
+  }
+};
+
+// print an M-by-N array
+template <typename T>
+void print(int m, int n, const thrust::device_vector<T>& d_data)
+{
+  thrust::host_vector<T> h_data = d_data;
+
+  for (int i = 0; i < m; i++)
+  {
+    for (int j = 0; j < n; j++)
+    {
+      std::cout << std::setw(4) << h_data[i * n + j] << " ";
+    }
+    std::cout << "\n";
+  }
+}
+
+void generate_random_sites(thrust::host_vector<int>& t, int Nb, int m, int n)
+{
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, m * n - 1);
+
+  for (int k = 0; k < Nb; k++)
+  {
+    int index = dist(rng);
+    t[index]  = index + 1;
+  }
+}
+
+// Export the tab to PGM image format
+void vector_to_pgm(thrust::host_vector<int>& t, int m, int n, const char* out)
+{
+  assert(static_cast<int>(t.size()) == m * n && "Vector size does not match image dims.");
+
+  std::fstream f(out, std::fstream::out);
+  f << "P2\n";
+  f << m << " " << n << "\n";
+  f << "253\n";
+
+  // Hash function to map values to [0,255]
+  auto to_grey_level = [](int in_value) -> int {
+    return (71 * in_value) % 253;
+  };
+
+  for (int value : t)
+  {
+    f << to_grey_level(value) << " ";
+  }
+  f << "\n";
+  f.close();
+}
+
+/************Main Jfa loop********************/
+// Perform a jump with step k
+void jfa(thrust::device_vector<int>& in, thrust::device_vector<int>& out, unsigned int k, int m, int n)
+{
+  thrust::transform(
+    thrust::make_zip_iterator(
+      in.begin(),
+      in.begin() + k,
+      in.begin() + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k,
+      in.begin() - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() + k + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() + k - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      thrust::counting_iterator<int>(0)),
+    thrust::make_zip_iterator(
+      in.begin(),
+      in.begin() + k,
+      in.begin() + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k,
+      in.begin() - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() + k + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() + k - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k + m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      in.begin() - k - m * k, // NOLINT(bugprone-misplaced-widening-cast)
+      thrust::counting_iterator<int>(0))
+      + n * m, // NOLINT(bugprone-misplaced-widening-cast)
+    out.begin(),
+    voronoi_site_selector(m, n, static_cast<int>(k)));
+}
+/********************************************/
+
+void display_time(timer& t)
+{
+  std::cout << "  ( " << 1e3 * t.elapsed() << "ms )" << '\n';
+}
+
+int main()
+{
+  int m = 2048; // number of rows
+  int n = 2048; // number of columns
+  int s = 1000; // number of sites
+
+  timer t;
+
+  // Host vector to encode a 2D image
+  std::cout << "[Initialize " << m << "x" << n << " Image]" << '\n';
+  t.restart();
+  thrust::host_vector<int> seeds_host(m * n, m * n);
+  generate_random_sites(seeds_host, s, m, n);
+  display_time(t);
+
+  std::cout << "[Copy to Device]" << '\n';
+  t.restart();
+  thrust::device_vector<int> seeds = seeds_host;
+  thrust::device_vector<int> temp(seeds);
+  display_time(t);
+
+  // JFA+1  : before entering the log(n) loop, we perform a jump with k=1
+  std::cout << "[JFA stepping]" << '\n';
+  t.restart();
+  jfa(seeds, temp, 1, m, n);
+  seeds.swap(temp);
+
+  // JFA : main loop with k=n/2, n/4, ..., 1
+  for (int k = thrust::max(m, n) / 2; k > 0; k /= 2)
+  {
+    jfa(seeds, temp, k, m, n);
+    seeds.swap(temp);
+  }
+
+  display_time(t);
+  std::cout << "  ( " << static_cast<double>(seeds.size()) / (1e6 * t.elapsed()) << " MPixel/s ) " << '\n';
+
+  std::cout << "[Device to Host Copy]" << '\n';
+  t.restart();
+  seeds_host = seeds;
+  display_time(t);
+
+  std::cout << "[PGM Export]" << '\n';
+  t.restart();
+  vector_to_pgm(seeds_host, m, n, "discrete_voronoi.pgm");
+  display_time(t);
+
+  return 0;
+}

cccl_upstream/thrust/examples/dot_products_with_zip.cu

@@ -0,0 +1,122 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/random.h>
+#include <thrust/transform.h>
+
+#include <iostream>
+
+// This example shows how thrust::zip_iterator can be used to create a
+// 'virtual' array of structures.  In this case the structure is a 3d
+// vector type (Float3) whose (x,y,z) components will be stored in
+// three separate float arrays.  The zip_iterator "zips" these arrays
+// into a single virtual Float3 array.
+
+// We'll use a 3-tuple to store our 3d vector type
+using Float3 = cuda::std::tuple<float, float, float>;
+
+// This functor implements the dot product between 3d vectors
+struct DotProduct
+{
+  __host__ __device__ float operator()(const Float3& a, const Float3& b) const
+  {
+    return cuda::std::get<0>(a) * cuda::std::get<0>(b) + // x components
+           cuda::std::get<1>(a) * cuda::std::get<1>(b) + // y components
+           cuda::std::get<2>(a) * cuda::std::get<2>(b); // z components
+  }
+};
+
+// Return a host vector with random values in the range [0,1)
+thrust::host_vector<float> random_vector(const size_t N, unsigned int seed = thrust::default_random_engine::default_seed)
+{
+  thrust::default_random_engine rng(seed);
+  thrust::uniform_real_distribution<float> u01(0.0f, 1.0f);
+  thrust::host_vector<float> temp(N);
+  for (size_t i = 0; i < N; ++i)
+  {
+    temp[i] = u01(rng);
+  }
+  return temp;
+}
+
+int main()
+{
+  // number of vectors
+  const size_t N = 1000;
+
+  // We'll store the components of the 3d vectors in separate arrays. One set of
+  // arrays will store the 'A' vectors and another set will store the 'B' vectors.
+
+  // This 'structure of arrays' (SoA) approach is usually more efficient than the
+  // 'array of structures' (AoS) approach.  The primary reason is that structures,
+  // like Float3, don't always obey the memory coalescing rules, so they are not
+  // efficiently transferred to and from memory.  Another reason to prefer SoA to
+  // AoS is that we don't always want to process all members of the structure.  For
+  // example, if we only need to look at first element of the structure then it
+  // is wasteful to load the entire structure from memory.  With the SoA approach,
+  // we can chose which elements of the structure we wish to read.
+
+  thrust::device_vector<float> A0 = random_vector(N); // x components of the 'A' vectors
+  thrust::device_vector<float> A1 = random_vector(N); // y components of the 'A' vectors
+  thrust::device_vector<float> A2 = random_vector(N); // z components of the 'A' vectors
+
+  thrust::device_vector<float> B0 = random_vector(N); // x components of the 'B' vectors
+  thrust::device_vector<float> B1 = random_vector(N); // y components of the 'B' vectors
+  thrust::device_vector<float> B2 = random_vector(N); // z components of the 'B' vectors
+
+  // Storage for result of each dot product
+  thrust::device_vector<float> result(N);
+
+  // We'll now illustrate two ways to use zip_iterator to compute the dot
+  // products.  The first method is verbose but shows how the parts fit together.
+  // The second method hides these details and is more concise.
+
+  // METHOD #1
+  // Defining a zip_iterator type can be a little cumbersome ...
+  using FloatIterator      = thrust::device_vector<float>::iterator;
+  using FloatIteratorTuple = cuda::std::tuple<FloatIterator, FloatIterator, FloatIterator>;
+  using Float3Iterator     = thrust::zip_iterator<FloatIteratorTuple>;
+
+  // Now we'll create some zip_iterators for A and B
+  Float3Iterator A_first = thrust::make_zip_iterator(A0.begin(), A1.begin(), A2.begin());
+  Float3Iterator A_last  = thrust::make_zip_iterator(A0.end(), A1.end(), A2.end());
+  Float3Iterator B_first = thrust::make_zip_iterator(B0.begin(), B1.begin(), B2.begin());
+
+  // Finally, we pass the zip_iterators into transform() as if they
+  // were 'normal' iterators for a device_vector<Float3>.
+  thrust::transform(A_first, A_last, B_first, result.begin(), DotProduct());
+
+  // METHOD #2
+  // Alternatively, we can avoid creating variables for X_first, X_last,
+  // and Y_first and invoke transform() directly.
+  thrust::transform(
+    thrust::make_zip_iterator(A0.begin(), A1.begin(), A2.begin()),
+    thrust::make_zip_iterator(A0.end(), A1.end(), A2.end()),
+    thrust::make_zip_iterator(B0.begin(), B1.begin(), B2.begin()),
+    result.begin(),
+    DotProduct());
+
+  // Finally, we'll print a few results
+
+  // Example output
+  // (0.840188,0.45724,0.0860517) * (0.0587587,0.456151,0.322409) = 0.285683
+  // (0.394383,0.640368,0.180886) * (0.0138811,0.24875,0.0221609) = 0.168775
+  // (0.783099,0.717092,0.426423) * (0.622212,0.0699601,0.234811) = 0.63755
+  // (0.79844,0.460067,0.0470658) * (0.0391351,0.742097,0.354747) = 0.389358
+  std::cout << std::fixed;
+  for (size_t i = 0; i < 4; i++)
+  {
+    Float3 a  = A_first[static_cast<std::ptrdiff_t>(i)];
+    Float3 b  = B_first[static_cast<std::ptrdiff_t>(i)];
+    float dot = result[i];
+
+    std::cout << "(" << cuda::std::get<0>(a) << "," << cuda::std::get<1>(a) << "," << cuda::std::get<2>(a) << ")";
+    std::cout << " * ";
+    std::cout << "(" << cuda::std::get<0>(b) << "," << cuda::std::get<1>(b) << "," << cuda::std::get<2>(b) << ")";
+    std::cout << " = ";
+    std::cout << dot << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/expand.cu

@@ -0,0 +1,80 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/gather.h>
+#include <thrust/reduce.h>
+#include <thrust/scan.h>
+
+#include <iostream>
+#include <iterator>
+
+// This example demonstrates how to expand an input sequence by
+// replicating each element a variable number of times. For example,
+//
+//   expand([2,2,2],[A,B,C]) -> [A,A,B,B,C,C]
+//   expand([3,0,1],[A,B,C]) -> [A,A,A,C]
+//   expand([1,3,2],[A,B,C]) -> [A,B,B,B,C,C]
+//
+// The element counts are assumed to be non-negative integers
+
+template <typename InputIterator1, typename InputIterator2, typename OutputIterator>
+OutputIterator expand(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, OutputIterator output)
+{
+  using difference_type = typename cuda::std::iterator_traits<InputIterator1>::difference_type;
+
+  difference_type input_size  = cuda::std::distance(first1, last1);
+  difference_type output_size = thrust::reduce(first1, last1);
+
+  // scan the counts to obtain output offsets for each input element
+  thrust::device_vector<difference_type> output_offsets(input_size, 0);
+  thrust::exclusive_scan(first1, last1, output_offsets.begin());
+
+  // scatter the nonzero counts into their corresponding output positions
+  thrust::device_vector<difference_type> output_indices(output_size, 0);
+  thrust::scatter_if(
+    thrust::counting_iterator<difference_type>(0),
+    thrust::counting_iterator<difference_type>(input_size),
+    output_offsets.begin(),
+    first1,
+    output_indices.begin());
+
+  // compute max-scan over the output indices, filling in the holes
+  thrust::inclusive_scan(
+    output_indices.begin(), output_indices.end(), output_indices.begin(), cuda::maximum<difference_type>{});
+
+  // gather input values according to index array (output = first2[output_indices])
+  thrust::gather(output_indices.begin(), output_indices.end(), first2, output);
+
+  // return output + output_size
+  cuda::std::advance(output, output_size);
+  return output;
+}
+
+template <typename Vector>
+void print(const std::string& s, const Vector& v)
+{
+  using T = typename Vector::value_type;
+
+  std::cout << s;
+  thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
+  std::cout << '\n';
+}
+
+int main()
+{
+  thrust::device_vector<int> d_counts = {3, 5, 2, 0, 1, 3, 4, 2, 4};
+  thrust::device_vector<int> d_values = {1, 2, 3, 4, 5, 6, 7, 8, 9};
+
+  const size_t output_size = thrust::reduce(d_counts.begin(), d_counts.end());
+  thrust::device_vector<int> d_output(output_size);
+
+  // expand values according to counts
+  expand(d_counts.begin(), d_counts.end(), d_values.begin(), d_output.begin());
+
+  std::cout << "Expanding values according to counts" << '\n';
+  print(" counts ", d_counts);
+  print(" values ", d_values);
+  print(" output ", d_output);
+
+  return 0;
+}

cccl_upstream/thrust/examples/fill_copy_sequence.cu

@@ -0,0 +1,33 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/host_vector.h>
+#include <thrust/sequence.h>
+
+#include <iostream>
+
+int main()
+{
+  // initialize all ten integers of a device_vector to 1
+  thrust::device_vector<int> D(10, 1);
+
+  // set the first seven elements of a vector to 9
+  thrust::fill(D.begin(), D.begin() + 7, 9);
+
+  // initialize a host_vector with the first five elements of D
+  thrust::host_vector<int> H(D.begin(), D.begin() + 5);
+
+  // set the elements of H to 0, 1, 2, 3, ...
+  thrust::sequence(H.begin(), H.end());
+
+  // copy all of H back to the beginning of D
+  thrust::copy(H.begin(), H.end(), D.begin());
+
+  // print D
+  for (size_t i = 0; i < D.size(); i++)
+  {
+    std::cout << "D[" << i << "] = " << D[i] << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/histogram.cu

@@ -0,0 +1,182 @@
+#include <thrust/adjacent_difference.h>
+#include <thrust/binary_search.h>
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/inner_product.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random.h>
+#include <thrust/sort.h>
+
+#include <cuda/iterator>
+
+#include <iomanip>
+#include <iostream>
+#include <iterator>
+
+// This example illustrates several methods for computing a
+// histogram [1] with Thrust.  We consider standard "dense"
+// histograms, where some bins may have zero entries, as well
+// as "sparse" histograms, where only the nonzero bins are
+// stored.  For example, histograms for the data set
+//    [2 1 0 0 2 2 1 1 1 1 4]
+// which contains 2 zeros, 5 ones, and 3 twos and 1 four, is
+//    [2 5 3 0 1]
+// using the dense method and
+//    [(0,2), (1,5), (2,3), (4,1)]
+// using the sparse method. Since there are no threes, the
+// sparse histogram representation does not contain a bin
+// for that value.
+//
+// Note that we choose to store the sparse histogram in two
+// separate arrays, one array of keys and one array of bin counts,
+//    [0 1 2 4] - keys
+//    [2 5 3 1] - bin counts
+// This "structure of arrays" format is generally faster and
+// more convenient to process than the alternative "array
+// of structures" layout.
+//
+// The best histogramming methods depends on the application.
+// If the number of bins is relatively small compared to the
+// input size, then the binary search-based dense histogram
+// method is probably best.  If the number of bins is comparable
+// to the input size, then the reduce_by_key-based sparse method
+// ought to be faster.  When in doubt, try both and see which
+// is fastest.
+//
+// [1] http://en.wikipedia.org/wiki/Histogram
+
+// simple routine to print contents of a vector
+template <typename Vector>
+void print_vector(const std::string& name, const Vector& v)
+{
+  using T = typename Vector::value_type;
+  std::cout << "  " << std::setw(20) << name << "  ";
+  thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
+  std::cout << '\n';
+}
+
+// dense histogram using binary search
+template <typename Vector1, typename Vector2>
+void dense_histogram(const Vector1& input, Vector2& histogram)
+{
+  using ValueType = typename Vector1::value_type; // input value type
+  using IndexType = typename Vector2::value_type; // histogram index type
+
+  // copy input data (could be skipped if input is allowed to be modified)
+  thrust::device_vector<ValueType> data(input);
+
+  // print the initial data
+  print_vector("initial data", data);
+
+  // sort data to bring equal elements together
+  thrust::sort(data.begin(), data.end());
+
+  // print the sorted data
+  print_vector("sorted data", data);
+
+  // number of histogram bins is equal to the maximum value plus one
+  IndexType num_bins = data.back() + 1;
+
+  // resize histogram storage
+  histogram.resize(num_bins);
+
+  // find the end of each bin of values
+  thrust::counting_iterator<IndexType> search_begin(0);
+  thrust::upper_bound(data.begin(), data.end(), search_begin, search_begin + num_bins, histogram.begin());
+
+  // print the cumulative histogram
+  print_vector("cumulative histogram", histogram);
+
+  // compute the histogram by taking differences of the cumulative histogram
+  thrust::adjacent_difference(histogram.begin(), histogram.end(), histogram.begin());
+
+  // print the histogram
+  print_vector("histogram", histogram);
+}
+
+// sparse histogram using reduce_by_key
+template <typename Vector1, typename Vector2, typename Vector3>
+void sparse_histogram(const Vector1& input, Vector2& histogram_values, Vector3& histogram_counts)
+{
+  using ValueType = typename Vector1::value_type; // input value type
+  using IndexType = typename Vector3::value_type; // histogram index type
+
+  // copy input data (could be skipped if input is allowed to be modified)
+  thrust::device_vector<ValueType> data(input);
+
+  // print the initial data
+  print_vector("initial data", data);
+
+  // sort data to bring equal elements together
+  thrust::sort(data.begin(), data.end());
+
+  // print the sorted data
+  print_vector("sorted data", data);
+
+  // number of histogram bins is equal to number of unique values (assumes data.size() > 0)
+  IndexType num_bins = thrust::inner_product(
+    data.begin(),
+    data.end() - 1,
+    data.begin() + 1,
+    IndexType(1),
+    cuda::std::plus<IndexType>(),
+    cuda::std::not_equal_to<ValueType>());
+
+  // resize histogram storage
+  histogram_values.resize(num_bins);
+  histogram_counts.resize(num_bins);
+
+  // compact find the end of each bin of values
+  thrust::reduce_by_key(
+    data.begin(), data.end(), cuda::constant_iterator<IndexType>(1), histogram_values.begin(), histogram_counts.begin());
+
+  // print the sparse histogram
+  print_vector("histogram values", histogram_values);
+  print_vector("histogram counts", histogram_counts);
+}
+
+int main()
+{
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, 9);
+
+  const int N = 40;
+  const int S = 4;
+
+  // generate random data on the host
+  thrust::host_vector<int> input(N);
+  for (int i = 0; i < N; i++)
+  {
+    int sum = 0;
+    for (int j = 0; j < S; j++)
+    {
+      sum += dist(rng);
+    }
+    input[i] = sum / S;
+  }
+
+  // demonstrate dense histogram method
+  {
+    std::cout << "Dense Histogram" << '\n';
+    thrust::device_vector<int> histogram;
+    dense_histogram(input, histogram);
+  }
+
+  // demonstrate sparse histogram method
+  {
+    std::cout << "Sparse Histogram" << '\n';
+    thrust::device_vector<int> histogram_values;
+    thrust::device_vector<int> histogram_counts;
+    sparse_histogram(input, histogram_values, histogram_counts);
+  }
+
+  // Note:
+  // A dense histogram can be converted to a sparse histogram
+  // using stream compaction (i.e. thrust::copy_if).
+  // A sparse histogram can be expanded into a dense histogram
+  // by initializing the dense histogram to zero (with thrust::fill)
+  // and then scattering the histogram counts (with thrust::scatter).
+
+  return 0;
+}

cccl_upstream/thrust/examples/include/host_device.h

@@ -0,0 +1,18 @@
+// SPDX-FileCopyrightText: Copyright (c) 2008-2009, NVIDIA Corporation. All rights reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#pragma once
+
+#include <cuda/__cccl_config>
+
+#if !_CCCL_CUDA_COMPILATION()
+
+#  ifndef __host__
+#    define __host__
+#  endif
+
+#  ifndef __device__
+#    define __device__
+#  endif
+
+#endif

cccl_upstream/thrust/examples/include/timer.h

@@ -0,0 +1,100 @@
+// SPDX-FileCopyrightText: Copyright (c) 2008-2009, NVIDIA Corporation. All rights reserved.
+// SPDX-License-Identifier: Apache-2.0
+
+#pragma once
+
+// A simple timer class
+
+#ifdef __CUDACC__
+
+// use CUDA's high-resolution timers when possible
+#  include <thrust/system/cuda/error.h>
+#  include <thrust/system_error.h>
+
+#  include <string>
+
+#  include <cuda_runtime_api.h>
+
+void cuda_safe_call(cudaError_t error, const std::string& message = "")
+{
+  if (error)
+  {
+    throw thrust::system_error(error, thrust::cuda_category(), message);
+  }
+}
+
+struct timer
+{
+  cudaEvent_t start;
+  cudaEvent_t end;
+
+  timer()
+  {
+    cuda_safe_call(cudaEventCreate(&start));
+    cuda_safe_call(cudaEventCreate(&end));
+    restart();
+  }
+
+  ~timer()
+  {
+    static_cast<void>(cudaEventDestroy(start));
+    static_cast<void>(cudaEventDestroy(end));
+  }
+
+  void restart()
+  {
+    cuda_safe_call(cudaEventRecord(start, nullptr));
+  }
+
+  double elapsed()
+  {
+    cuda_safe_call(cudaEventRecord(end, nullptr));
+    cuda_safe_call(cudaEventSynchronize(end));
+
+    float ms_elapsed;
+    cuda_safe_call(cudaEventElapsedTime(&ms_elapsed, start, end));
+    return ms_elapsed / 1e3;
+  }
+
+  double epsilon()
+  {
+    return 0.5e-6;
+  }
+};
+
+#else
+
+// fallback to clock()
+#  include <ctime>
+
+struct timer
+{
+  clock_t start;
+  clock_t end;
+
+  timer()
+  {
+    restart();
+  }
+
+  ~timer() = default;
+
+  void restart()
+  {
+    start = clock();
+  }
+
+  double elapsed()
+  {
+    end = clock();
+
+    return static_cast<double>(end - start) / static_cast<double>(CLOCKS_PER_SEC);
+  }
+
+  double epsilon()
+  {
+    return 1.0 / static_cast<double>(CLOCKS_PER_SEC);
+  }
+};
+
+#endif

cccl_upstream/thrust/examples/lambda.cu

@@ -0,0 +1,85 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/transform.h>
+
+#include <iostream>
+
+// This example demonstrates the use of placeholders to implement
+// the SAXPY operation (i.e. Y[i] = a * X[i] + Y[i]).
+//
+// Placeholders enable developers to write concise inline expressions
+// instead of full functors for many simple operations.  For example,
+// the placeholder expression "_1 + _2" means to add the first argument,
+// represented by _1, to the second argument, represented by _2.
+// The names _1, _2, _3, _4 ... _10 represent the first ten arguments
+// to the function.
+//
+// In this example, the placeholder expression "a * _1 + _2" is used
+// to implement the SAXPY operation.  Note that the placeholder
+// implementation is considerably shorter and written inline.
+
+// allows us to use "_1" instead of "thrust::placeholders::_1"
+using namespace thrust::placeholders;
+
+// implementing SAXPY with a functor is cumbersome and verbose
+struct saxpy_functor
+{
+  float a;
+
+  saxpy_functor(float a)
+      : a(a)
+  {}
+
+  __host__ __device__ float operator()(float x, float y)
+  {
+    return a * x + y;
+  }
+};
+
+int main()
+{
+  // input data
+  float a                             = 2.0f;
+  thrust::device_vector<float> x_data = {1, 2, 3, 4};
+  thrust::device_vector<float> y_data = {1, 1, 1, 1};
+
+  // SAXPY implemented with a functor (function object)
+  {
+    thrust::device_vector<float> X = x_data;
+    thrust::device_vector<float> Y = y_data;
+
+    thrust::transform(
+      X.begin(),
+      X.end(), // input range #1
+      Y.begin(), // input range #2
+      Y.begin(), // output range
+      saxpy_functor(a)); // functor
+
+    std::cout << "SAXPY (functor method)" << '\n';
+    for (size_t i = 0; i < Y.size(); i++)
+    {
+      std::cout << a << " * " << x_data[i] << " + " << y_data[i] << " = " << Y[i] << '\n';
+    }
+  }
+
+  // SAXPY implemented with a placeholders
+  {
+    thrust::device_vector<float> X = x_data;
+    thrust::device_vector<float> Y = y_data;
+
+    thrust::transform(
+      X.begin(),
+      X.end(), // input range #1
+      Y.begin(), // input range #2
+      Y.begin(), // output range
+      a * _1 + _2); // placeholder expression
+
+    std::cout << "SAXPY (placeholder method)" << '\n';
+    for (size_t i = 0; i < Y.size(); i++)
+    {
+      std::cout << a << " * " << x_data[i] << " + " << y_data[i] << " = " << Y[i] << '\n';
+    }
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/lexicographical_sort.cu

@@ -0,0 +1,91 @@
+#include <thrust/device_vector.h>
+#include <thrust/gather.h>
+#include <thrust/generate.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/sequence.h>
+#include <thrust/sort.h>
+
+#include <iostream>
+
+// This example shows how to perform a lexicographical sort on multiple keys.
+//
+// http://en.wikipedia.org/wiki/Lexicographical_order
+
+template <typename KeyVector, typename PermutationVector>
+void update_permutation(KeyVector& keys, PermutationVector& permutation)
+{
+  // temporary storage for keys
+  KeyVector temp(keys.size());
+
+  // permute the keys with the current reordering
+  thrust::gather(permutation.begin(), permutation.end(), keys.begin(), temp.begin());
+
+  // stable_sort the permuted keys and update the permutation
+  thrust::stable_sort_by_key(temp.begin(), temp.end(), permutation.begin());
+}
+
+template <typename KeyVector, typename PermutationVector>
+void apply_permutation(KeyVector& keys, PermutationVector& permutation)
+{
+  // copy keys to temporary vector
+  KeyVector temp(keys.begin(), keys.end());
+
+  // permute the keys
+  thrust::gather(permutation.begin(), permutation.end(), temp.begin(), keys.begin());
+}
+
+thrust::host_vector<int> random_vector(size_t N)
+{
+  thrust::host_vector<int> vec(N);
+  static thrust::default_random_engine rng;
+  static thrust::uniform_int_distribution<int> dist(0, 9);
+
+  for (size_t i = 0; i < N; i++)
+  {
+    vec[i] = dist(rng);
+  }
+
+  return vec;
+}
+
+int main()
+{
+  size_t N = 20;
+
+  // generate three arrays of random values
+  thrust::device_vector<int> upper  = random_vector(N);
+  thrust::device_vector<int> middle = random_vector(N);
+  thrust::device_vector<int> lower  = random_vector(N);
+
+  std::cout << "Unsorted Keys" << '\n';
+  for (size_t i = 0; i < upper.size(); i++)
+  {
+    std::cout << "(" << upper[i] << "," << middle[i] << "," << lower[i] << ")" << '\n';
+  }
+
+  // initialize permutation to [0, 1, 2, ... ,N-1]
+  thrust::device_vector<int> permutation(N);
+  thrust::sequence(permutation.begin(), permutation.end());
+
+  // sort from least significant key to most significant keys
+  update_permutation(lower, permutation);
+  update_permutation(middle, permutation);
+  update_permutation(upper, permutation);
+
+  // Note: keys have not been modified
+  // Note: permutation now maps unsorted keys to sorted order
+
+  // permute the key arrays by the final permutation
+  apply_permutation(lower, permutation);
+  apply_permutation(middle, permutation);
+  apply_permutation(upper, permutation);
+
+  std::cout << "Sorted Keys" << '\n';
+  for (size_t i = 0; i < upper.size(); i++)
+  {
+    std::cout << "(" << upper[i] << "," << middle[i] << "," << lower[i] << ")" << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/max_abs_diff.cu

@@ -0,0 +1,37 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/inner_product.h>
+
+#include <cuda/functional>
+
+#include <iostream>
+
+// this example computes the maximum absolute difference
+// between the elements of two vectors
+
+template <typename T>
+struct abs_diff
+{
+  __host__ __device__ T operator()(const T& a, const T& b)
+  {
+    return fabsf(b - a);
+  }
+};
+
+int main()
+{
+  thrust::device_vector<float> d_a = {1.0, 2.0, 3.0, 4.0};
+  thrust::device_vector<float> d_b = {2.0, 4.0, 3.0, 0.0};
+
+  // initial value of the reduction
+  float init = 0;
+
+  // binary operations
+  cuda::maximum<float> binary_op1{};
+  abs_diff<float> binary_op2;
+
+  float max_abs_diff = thrust::inner_product(d_a.begin(), d_a.end(), d_b.begin(), init, binary_op1, binary_op2);
+
+  std::cout << "maximum absolute difference: " << max_abs_diff << '\n';
+  return 0;
+}

cccl_upstream/thrust/examples/minimal_custom_backend.cu

@@ -0,0 +1,59 @@
+#include <thrust/device_vector.h>
+#include <thrust/execution_policy.h>
+#include <thrust/for_each.h>
+#include <thrust/functional.h>
+#include <thrust/transform.h>
+
+#include <iostream>
+
+// This example demonstrates how to build a minimal custom
+// Thrust backend by intercepting for_each's dispatch.
+
+// We begin by defining a "system", which distinguishes our novel
+// backend from other Thrust backends.
+// We'll derive my_system from thrust::device_execution_policy to inherit
+// the functionality of the default device backend.
+// Note that we pass the name of our system as a template parameter
+// to thrust::device_execution_policy.
+struct my_system : thrust::device_execution_policy<my_system>
+{};
+
+// Next, we'll create a novel version of for_each which only
+// applies to algorithm invocations executed with my_system.
+// Our version of for_each will print a message and then call
+// the regular device version of for_each.
+
+// The first parameter to our version for_each is my_system. This allows
+// Thrust to locate it when dispatching thrust::for_each.
+// The following parameters are as normal.
+template <typename Iterator, typename Function>
+Iterator for_each(my_system, Iterator first, Iterator last, Function f)
+{
+  // output a message
+  std::cout << "Hello, world from for_each(my_system)!" << '\n';
+
+  // to call the normal device version of for_each, pass thrust::device as the first parameter.
+  return thrust::for_each(thrust::device, first, last, f);
+}
+
+int main()
+{
+  thrust::device_vector<int> vec(1);
+
+  // create an instance of our system
+  my_system sys;
+
+  // To invoke our version of for_each, pass sys as the first parameter
+  thrust::for_each(sys, vec.begin(), vec.end(), cuda::std::negate{});
+
+  // Other algorithms that Thrust implements with thrust::for_each will also
+  // cause our version of for_each to be invoked when we pass an instance of my_system as the first parameter.
+  // Even though we did not define a special version of transform, Thrust dispatches the version it knows
+  // for thrust::device_execution_policy, which my_system inherits.
+  thrust::transform(sys, vec.begin(), vec.end(), vec.begin(), cuda::std::identity{});
+
+  // Invocations without my_system are handled normally.
+  thrust::for_each(vec.begin(), vec.end(), cuda::std::negate{});
+
+  return 0;
+}

cccl_upstream/thrust/examples/minmax.cu

@@ -0,0 +1,91 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/transform_reduce.h>
+
+#include <iostream>
+
+// compute minimum and maximum values in a single reduction
+
+// minmax_pair stores the minimum and maximum
+// values that have been encountered so far
+template <typename T>
+struct minmax_pair
+{
+  T min_val;
+  T max_val;
+};
+
+// minmax_unary_op is a functor that takes in a value x and
+// returns a minmax_pair whose minimum and maximum values
+// are initialized to x.
+template <typename T>
+struct minmax_unary_op
+{
+  __host__ __device__ minmax_pair<T> operator()(const T& x) const
+  {
+    minmax_pair<T> result;
+    result.min_val = x;
+    result.max_val = x;
+    return result;
+  }
+};
+
+// minmax_binary_op is a functor that accepts two minmax_pair
+// structs and returns a new minmax_pair whose minimum and
+// maximum values are the min() and max() respectively of
+// the minimums and maximums of the input pairs
+template <typename T>
+struct minmax_binary_op
+{
+  __host__ __device__ minmax_pair<T> operator()(const minmax_pair<T>& x, const minmax_pair<T>& y) const
+  {
+    minmax_pair<T> result;
+    result.min_val = thrust::min(x.min_val, y.min_val);
+    result.max_val = thrust::max(x.max_val, y.max_val);
+    return result;
+  }
+};
+
+int main()
+{
+  // input size
+  size_t N = 10;
+
+  // initialize random number generator
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(10, 99);
+
+  // initialize data on host
+  thrust::host_vector<int> host_data(N);
+  for (auto& e : host_data)
+  {
+    e = dist(rng);
+  }
+  thrust::device_vector<int> data = host_data;
+
+  // setup arguments
+  minmax_unary_op<int> unary_op;
+  minmax_binary_op<int> binary_op;
+
+  // initialize reduction with the first value
+  minmax_pair<int> init = unary_op(data[0]);
+
+  // compute minimum and maximum values
+  minmax_pair<int> result = thrust::transform_reduce(data.begin(), data.end(), unary_op, init, binary_op);
+
+  // print results
+  std::cout << "[ ";
+  for (auto& e : host_data)
+  {
+    std::cout << e << " ";
+  }
+  std::cout << "]" << '\n';
+
+  std::cout << "minimum = " << result.min_val << '\n';
+  std::cout << "maximum = " << result.max_val << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/mode.cu

@@ -0,0 +1,80 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/random.h>
+#include <thrust/reduce.h>
+#include <thrust/sort.h>
+#include <thrust/unique.h>
+
+#include <cuda/iterator>
+
+#include <iostream>
+#include <iterator>
+
+// This example compute the mode [1] of a set of numbers.  If there
+// are multiple modes, one with the smallest value it returned.
+//
+// [1] http://en.wikipedia.org/wiki/Mode_(statistics)
+
+int main()
+{
+  const size_t N = 30;
+  const size_t M = 10;
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, M - 1);
+
+  // generate random data on the host
+  thrust::host_vector<int> h_data(N);
+  for (auto& e : h_data)
+  {
+    e = dist(rng);
+  }
+
+  // transfer data to device
+  thrust::device_vector<int> d_data(h_data);
+
+  // print the initial data
+  std::cout << "initial data" << '\n';
+  thrust::copy(d_data.begin(), d_data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // sort data to bring equal elements together
+  thrust::sort(d_data.begin(), d_data.end());
+
+  // print the sorted data
+  std::cout << "sorted data" << '\n';
+  thrust::copy(d_data.begin(), d_data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // count number of unique keys
+  size_t num_unique = thrust::unique_count(d_data.begin(), d_data.end());
+
+  // count multiplicity of each key
+  thrust::device_vector<int> d_output_keys(num_unique);
+  thrust::device_vector<int> d_output_counts(num_unique);
+  thrust::reduce_by_key(
+    d_data.begin(), d_data.end(), cuda::constant_iterator<int>(1), d_output_keys.begin(), d_output_counts.begin());
+
+  // print the counts
+  std::cout << "values" << '\n';
+  thrust::copy(d_output_keys.begin(), d_output_keys.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // print the counts
+  std::cout << "counts" << '\n';
+  thrust::copy(d_output_counts.begin(), d_output_counts.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // find the index of the maximum count
+  thrust::device_vector<int>::iterator mode_iter;
+  mode_iter = thrust::max_element(d_output_counts.begin(), d_output_counts.end());
+
+  int mode        = d_output_keys[cuda::std::distance(d_output_counts.begin(), mode_iter)];
+  int occurrences = *mode_iter;
+
+  std::cout << "Modal value " << mode << " occurs " << occurrences << " times " << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/monte_carlo.cu

@@ -0,0 +1,76 @@
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random.h>
+#include <thrust/transform_reduce.h>
+
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+
+// we could vary M & N to find the perf sweet spot
+
+__host__ __device__ unsigned int hash(unsigned int a)
+{
+  a = (a + 0x7ed55d16) + (a << 12);
+  a = (a ^ 0xc761c23c) ^ (a >> 19);
+  a = (a + 0x165667b1) + (a << 5);
+  a = (a + 0xd3a2646c) ^ (a << 9);
+  a = (a + 0xfd7046c5) + (a << 3);
+  a = (a ^ 0xb55a4f09) ^ (a >> 16);
+  return a;
+}
+
+struct estimate_pi
+{
+  __host__ __device__ float operator()(unsigned int thread_id)
+  {
+    float sum      = 0;
+    unsigned int N = 10000; // samples per thread
+
+    unsigned int seed = hash(thread_id);
+
+    // seed a random number generator
+    thrust::default_random_engine rng(seed);
+
+    // create a mapping from random numbers to [0,1)
+    thrust::uniform_real_distribution<float> u01(0, 1);
+
+    // take N samples in a quarter circle
+    for (unsigned int i = 0; i < N; ++i)
+    {
+      // draw a sample from the unit square
+      float x = u01(rng);
+      float y = u01(rng);
+
+      // measure distance from the origin
+      float dist = sqrtf(x * x + y * y);
+
+      // add 1.0f if (u0,u1) is inside the quarter circle
+      if (dist <= 1.0f)
+      {
+        sum += 1.0f;
+      }
+    }
+
+    // multiply by 4 to get the area of the whole circle
+    sum *= 4.0f;
+
+    // divide by N
+    return sum / static_cast<float>(N);
+  }
+};
+
+int main()
+{
+  // use 30K independent seeds
+  int M = 30000;
+
+  float estimate = thrust::transform_reduce(
+    thrust::counting_iterator<int>(0), thrust::counting_iterator<int>(M), estimate_pi(), 0.0f, cuda::std::plus<float>());
+  estimate /= static_cast<float>(M);
+
+  std::cout << std::setprecision(3);
+  std::cout << "pi is approximately " << estimate << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/monte_carlo_disjoint_sequences.cu

@@ -0,0 +1,82 @@
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/random.h>
+#include <thrust/transform_reduce.h>
+
+#include <cmath>
+#include <iostream>
+
+// The technique demonstrated in the example monte_carlo.cu
+// assigns an independently seeded random number generator to each
+// of 30K threads, and uses a hashing scheme based on thread index to
+// seed each RNG. This technique, while simple, may be succeptible
+// to correlation among the streams of numbers generated by each RNG
+// because there is no guarantee that the streams are not disjoint.
+// This example demonstrates a slightly more sophisticated technique
+// which ensures that the subsequences generated in each thread are
+// disjoint. To achieve this, we use a single common stream
+// of random numbers, but partition it among threads to ensure no overlap
+// of substreams. The substreams are generated procedurally using
+// default_random_engine's discard(n) member function, which skips
+// past n states of the RNG. This function is accelerated and executes
+// in O(lg n) time.
+
+struct estimate_pi
+{
+  __host__ __device__ float operator()(unsigned int thread_id)
+  {
+    float sum      = 0;
+    unsigned int N = 5000; // samples per stream
+
+    // note that M * N <= default_random_engine::max,
+    // which is also the period of this particular RNG
+    // this ensures the substreams are disjoint
+
+    // create a random number generator
+    // note that each thread uses an RNG with the same seed
+    thrust::default_random_engine rng;
+
+    // jump past the numbers used by the subsequences before me
+    rng.discard(N * thread_id); // NOLINT(bugprone-misplaced-widening-cast)
+
+    // create a mapping from random numbers to [0,1)
+    thrust::uniform_real_distribution<float> u01(0, 1);
+
+    // take N samples in a quarter circle
+    for (unsigned int i = 0; i < N; ++i)
+    {
+      // draw a sample from the unit square
+      float x = u01(rng);
+      float y = u01(rng);
+
+      // measure distance from the origin
+      float dist = sqrtf(x * x + y * y);
+
+      // add 1.0f if (u0,u1) is inside the quarter circle
+      if (dist <= 1.0f)
+      {
+        sum += 1.0f;
+      }
+    }
+
+    // multiply by 4 to get the area of the whole circle
+    sum *= 4.0f;
+
+    // divide by N
+    return sum / static_cast<float>(N);
+  }
+};
+
+int main()
+{
+  // use 30K subsequences of random numbers
+  int M = 30000;
+
+  float estimate = thrust::transform_reduce(
+    thrust::counting_iterator<int>(0), thrust::counting_iterator<int>(M), estimate_pi(), 0.0f, cuda::std::plus<float>());
+  estimate /= static_cast<float>(M);
+
+  std::cout << "pi is around " << estimate << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/mr_basic.cu

@@ -0,0 +1,77 @@
+#include <thrust/device_ptr.h>
+#include <thrust/device_vector.h>
+#include <thrust/mr/allocator.h>
+#include <thrust/mr/disjoint_pool.h>
+#include <thrust/mr/new.h>
+#include <thrust/mr/pool.h>
+
+#include <cassert>
+
+template <typename Vec>
+void do_stuff_with_vector(typename Vec::allocator_type alloc)
+{
+  Vec v1(alloc);
+  v1.push_back(1);
+  assert(v1.back() == 1);
+
+  Vec v2(alloc);
+  v2 = v1;
+
+  v1.swap(v2);
+
+  v1.clear();
+  v1.resize(2);
+  assert(v1.size() == 2);
+}
+
+int main()
+{
+  thrust::mr::new_delete_resource memres;
+
+  {
+    // no virtual calls will be issued
+    using Alloc = thrust::mr::allocator<int, thrust::mr::new_delete_resource>;
+    Alloc alloc(&memres);
+
+    do_stuff_with_vector<thrust::host_vector<int, Alloc>>(alloc);
+  }
+
+  {
+    // virtual calls will be issued - wrapping in a polymorphic wrapper
+    thrust::mr::polymorphic_adaptor_resource<void*> adaptor(&memres);
+    using Alloc = thrust::mr::polymorphic_allocator<int, void*>;
+    Alloc alloc(&adaptor);
+
+    do_stuff_with_vector<thrust::host_vector<int, Alloc>>(alloc);
+  }
+
+  {
+    // use the global device_ptr-flavored device memory resource
+    using Resource = thrust::device_ptr_memory_resource<thrust::device_memory_resource>;
+    thrust::mr::polymorphic_adaptor_resource<thrust::device_ptr<void>> adaptor(
+      thrust::mr::get_global_resource<Resource>());
+    using Alloc = thrust::mr::polymorphic_allocator<int, thrust::device_ptr<void>>;
+    Alloc alloc(&adaptor);
+
+    do_stuff_with_vector<thrust::device_vector<int, Alloc>>(alloc);
+  }
+
+  using Pool = thrust::mr::unsynchronized_pool_resource<thrust::mr::new_delete_resource>;
+  Pool pool(&memres);
+  {
+    using Alloc = thrust::mr::allocator<int, Pool>;
+    Alloc alloc(&pool);
+
+    do_stuff_with_vector<thrust::host_vector<int, Alloc>>(alloc);
+  }
+
+  using DisjointPool =
+    thrust::mr::disjoint_unsynchronized_pool_resource<thrust::mr::new_delete_resource, thrust::mr::new_delete_resource>;
+  DisjointPool disjoint_pool(&memres, &memres);
+  {
+    using Alloc = thrust::mr::allocator<int, DisjointPool>;
+    Alloc alloc(&disjoint_pool);
+
+    do_stuff_with_vector<thrust::host_vector<int, Alloc>>(alloc);
+  }
+}

cccl_upstream/thrust/examples/norm.cu

@@ -0,0 +1,46 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/transform_reduce.h>
+
+#include <cmath>
+#include <iostream>
+
+//   This example computes the norm [1] of a vector.  The norm is
+// computed by squaring all numbers in the vector, summing the
+// squares, and taking the square root of the sum of squares.  In
+// Thrust this operation is efficiently implemented with the
+// transform_reduce() algorithm.  Specifically, we first transform
+// x -> x^2 and the compute a standard plus reduction.  Since there
+// is no built-in functor for squaring numbers, we define our own
+// square functor.
+//
+// [1] http://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm
+
+// square<T> computes the square of a number f(x) -> x*x
+template <typename T>
+struct square
+{
+  __host__ __device__ T operator()(const T& x) const
+  {
+    return x * x;
+  }
+};
+
+int main()
+{
+  // initialize device vector directly
+  thrust::device_vector<float> d_x = {1.0, 2.0, 3.0, 4.0};
+
+  // setup arguments
+  square<float> unary_op;
+  cuda::std::plus<float> binary_op;
+  float init = 0;
+
+  // compute norm
+  float norm = std::sqrt(thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op));
+
+  std::cout << "norm is " << norm << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/padded_grid_reduction.cu

@@ -0,0 +1,120 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/zip_iterator.h>
+#include <thrust/random.h>
+#include <thrust/transform_reduce.h>
+
+#include <cfloat>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+
+// This example computes the minimum and maximum values
+// over a padded grid.  The padded values are not considered
+// during the reduction operation.
+
+// transform a tuple (int,value) into a tuple (bool,value,value)
+// where the bool is true for valid grid values and false for
+// values in the padded region of the grid
+template <typename IndexType, typename ValueType>
+struct transform_tuple
+{
+  using InputTuple  = typename cuda::std::tuple<IndexType, ValueType>;
+  using OutputTuple = typename cuda::std::tuple<bool, ValueType, ValueType>;
+
+  IndexType n, N;
+
+  transform_tuple(IndexType n, IndexType N)
+      : n(n)
+      , N(N)
+  {}
+
+  __host__ __device__ OutputTuple operator()(const InputTuple& t) const
+  {
+    bool is_valid = (cuda::std::get<0>(t) % N) < n;
+    return OutputTuple(is_valid, cuda::std::get<1>(t), cuda::std::get<1>(t));
+  }
+};
+
+// reduce two tuples (bool,value,value) into a single tuple such that output
+// contains the smallest and largest *valid* values.
+template <typename IndexType, typename ValueType>
+struct reduce_tuple
+{
+  using Tuple = typename cuda::std::tuple<bool, ValueType, ValueType>;
+
+  __host__ __device__ Tuple operator()(const Tuple& t0, const Tuple& t1) const
+  {
+    if (cuda::std::get<0>(t0) && cuda::std::get<0>(t1)) // both valid
+    {
+      return Tuple(true,
+                   thrust::min(cuda::std::get<1>(t0), cuda::std::get<1>(t1)),
+                   thrust::max(cuda::std::get<2>(t0), cuda::std::get<2>(t1)));
+    }
+    else if (cuda::std::get<0>(t0))
+    {
+      return t0;
+    }
+    else
+    {
+      return t1; // if t0 is not valid, return t1 whether it is valid or not
+    }
+  }
+};
+
+int main()
+{
+  int M = 10; // number of rows
+  int n = 11; // number of columns excluding padding
+  int N = 16; // number of columns including padding
+
+  thrust::default_random_engine rng(12345);
+  thrust::uniform_real_distribution<float> dist(0.0f, 1.0f);
+
+  thrust::device_vector<float> data(M * N, -1);
+
+  // initialize valid values in grid
+  for (int i = 0; i < M; i++)
+  {
+    for (int j = 0; j < n; j++)
+    {
+      data[static_cast<std::size_t>(i) * N + j] = dist(rng);
+    }
+  }
+
+  // print full grid
+  std::cout << "padded grid" << '\n';
+  std::cout << std::fixed << std::setprecision(4);
+  for (int i = 0; i < M; i++)
+  {
+    std::cout << " ";
+    for (int j = 0; j < N; j++)
+    {
+      std::cout << data[(static_cast<std::size_t>(i) * N) + j] << " ";
+    }
+    std::cout << "\n";
+  }
+  std::cout << "\n";
+
+  // compute min & max over valid region of the 2d grid
+  using result_type = cuda::std::tuple<bool, float, float>;
+
+  result_type init(true, FLT_MAX, -FLT_MAX); // initial value
+  transform_tuple<int, float> unary_op(n, N); // transformation operator
+  reduce_tuple<int, float> binary_op; // reduction operator
+
+  result_type result = thrust::transform_reduce(
+    thrust::make_zip_iterator(thrust::counting_iterator<int>(0), data.begin()),
+    thrust::make_zip_iterator(cuda::std::tuple(thrust::counting_iterator<int>(0), data.begin()))
+      + static_cast<std::ptrdiff_t>(data.size()),
+    unary_op,
+    init,
+    binary_op);
+
+  std::cout << "minimum value: " << cuda::std::get<1>(result) << '\n';
+  std::cout << "maximum value: " << cuda::std::get<2>(result) << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/permutation_iterator.cu

@@ -0,0 +1,27 @@
+#include <thrust/device_vector.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/reduce.h>
+
+#include <iostream>
+
+// this example fuses a gather operation with a reduction for
+// greater efficiency than separate gather() and reduce() calls
+
+int main()
+{
+  // gather locations
+  thrust::device_vector<int> map = {3, 1, 0, 5};
+
+  // array to gather from
+  thrust::device_vector<int> source = {10, 20, 30, 40, 50, 60};
+
+  // fuse gather with reduction:
+  //   sum = source[map[0]] + source[map[1]] + ...
+  int sum = thrust::reduce(thrust::make_permutation_iterator(source.begin(), map.begin()),
+                           thrust::make_permutation_iterator(source.begin(), map.end()));
+
+  // print sum
+  std::cout << "sum is " << sum << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/print_version.cu

@@ -0,0 +1,15 @@
+#include <thrust/version.h>
+
+#include <iostream>
+
+int main()
+{
+  int major    = THRUST_MAJOR_VERSION;
+  int minor    = THRUST_MINOR_VERSION;
+  int subminor = THRUST_SUBMINOR_VERSION;
+  int patch    = THRUST_PATCH_NUMBER;
+
+  std::cout << "Thrust v" << major << "." << minor << "." << subminor << "-" << patch << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/raw_reference_cast.cu

@@ -0,0 +1,106 @@
+#include <thrust/detail/raw_reference_cast.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/sequence.h>
+
+#include <iostream>
+
+// This example illustrates how to use the raw_reference_cast to convert
+// system-specific reference wrappers into native references.
+//
+// Using iterators in the manner described here is generally discouraged.
+// Users should only resort to this technique if there is no viable
+// implementation of a given operation in terms of Thrust algorithms.
+// For example this particular example is better solved with thrust::copy,
+// which is safer and potentially faster.  Only use this approach after all
+// safer alternatives have been exhausted.
+//
+// When a Thrust iterator is referenced (e.g. *iter) the result is not
+// a native or "raw" reference like int& or float&.  Instead,
+// the result is a type such as thrust::system::cuda::reference<int>
+// or thrust::system::tbb::reference<float>, depending on the system
+// to which the data belongs.  These reference wrappers are necessary
+// to make expressions like *iter1 = *iter2; work correctly when
+// iter1 and iter2 refer to data in different memory spaces on
+// heterogenous systems.
+//
+// The raw_reference_cast function essentially strips away the system-specific
+// meta-data so it should only be used when the code is guaranteed to be
+// executed within an appropriate context.
+
+__host__ __device__ void assign_reference_to_reference(int& x, int& y)
+{
+  y = x;
+}
+
+__host__ __device__ void assign_value_to_reference(int x, int& y)
+{
+  y = x;
+}
+
+template <typename InputIterator, typename OutputIterator>
+struct copy_iterators
+{
+  InputIterator input;
+  OutputIterator output;
+
+  copy_iterators(InputIterator input, OutputIterator output)
+      : input(input)
+      , output(output)
+  {}
+
+  __host__ __device__ void operator()(int i)
+  {
+    InputIterator in   = input + i;
+    OutputIterator out = output + i;
+
+    // invalid - reference<int> is not convertible to int&
+    // assign_reference_to_reference(*in, *out);
+
+    // valid - reference<int> explicitly converted to int&
+    assign_reference_to_reference(thrust::raw_reference_cast(*in), thrust::raw_reference_cast(*out));
+
+    // valid - since reference<int> is convertible to int
+    assign_value_to_reference(*in, thrust::raw_reference_cast(*out));
+  }
+};
+
+template <typename Vector>
+void print(const std::string& name, const Vector& v)
+{
+  using T = typename Vector::value_type;
+
+  std::cout << name << ": ";
+  thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
+  std::cout << "\n";
+}
+
+int main()
+{
+  using Vector   = thrust::device_vector<int>;
+  using Iterator = Vector::iterator;
+  using System   = thrust::device_system_tag;
+
+  // allocate device memory
+  Vector A(5);
+  Vector B(5);
+
+  // initialize A and B
+  thrust::sequence(A.begin(), A.end());
+  thrust::fill(B.begin(), B.end(), 0);
+
+  std::cout << "Before A->B Copy" << '\n';
+  print("A", A);
+  print("B", B);
+
+  // note: we must specify the System to ensure correct execution
+  thrust::for_each(thrust::counting_iterator<int, System>(0),
+                   thrust::counting_iterator<int, System>(5),
+                   copy_iterators<Iterator, Iterator>(A.begin(), B.begin()));
+
+  std::cout << "After A->B Copy" << '\n';
+  print("A", A);
+  print("B", B);
+
+  return 0;
+}

cccl_upstream/thrust/examples/remove_points2d.cu

@@ -0,0 +1,78 @@
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/remove.h>
+
+#include <iostream>
+
+// This example generates random points in the
+// unit square [0,1)x[0,1) and then removes all
+// points where x^2 + y^2 > 1
+//
+// The x and y coordinates are stored in separate arrays
+// and a zip_iterator is used to combine them together
+
+template <typename T>
+struct is_outside_circle
+{
+  template <typename Tuple>
+  inline __host__ __device__ bool operator()(const Tuple& tuple) const
+  {
+    // unpack the tuple into x and y coordinates
+    const T x = cuda::std::get<0>(tuple);
+    const T y = cuda::std::get<1>(tuple);
+
+    if (x * x + y * y > 1)
+    {
+      return true;
+    }
+    else
+    {
+      return false;
+    }
+  }
+};
+
+int main()
+{
+  const size_t N = 20;
+
+  // generate random points in the unit square on the host
+  thrust::default_random_engine rng;
+  thrust::uniform_real_distribution<float> u01(0.0f, 1.0f);
+  thrust::host_vector<float> x(N);
+  thrust::host_vector<float> y(N);
+  for (size_t i = 0; i < N; i++)
+  {
+    x[i] = u01(rng);
+    y[i] = u01(rng);
+  }
+
+  // print the initial points
+  std::cout << std::fixed;
+  std::cout << "Generated " << N << " points" << '\n';
+  for (size_t i = 0; i < x.size(); i++)
+  {
+    std::cout << "(" << x[i] << "," << y[i] << ")" << '\n';
+  }
+  std::cout << '\n';
+
+  // remove points where x^2 + y^2 > 1 and determine new array sizes
+  size_t new_size =
+    thrust::remove_if(thrust::make_zip_iterator(x.begin(), y.begin()),
+                      thrust::make_zip_iterator(x.end(), y.end()),
+                      is_outside_circle<float>())
+    - thrust::make_zip_iterator(x.begin(), y.begin());
+
+  // resize the vectors (note: this does not free any memory)
+  x.resize(new_size);
+  y.resize(new_size);
+
+  // print the filtered points
+  std::cout << "After stream compaction, " << new_size << " points remain" << '\n';
+  for (size_t i = 0; i < x.size(); i++)
+  {
+    std::cout << "(" << x[i] << "," << y[i] << ")" << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/repeated_range.cu

@@ -0,0 +1,92 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+
+#include <iostream>
+
+// this example illustrates how to make repeated access to a range of values
+// examples:
+//   repeated_range([0, 1, 2, 3], 1) -> [0, 1, 2, 3]
+//   repeated_range([0, 1, 2, 3], 2) -> [0, 0, 1, 1, 2, 2, 3, 3]
+//   repeated_range([0, 1, 2, 3], 3) -> [0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3]
+//   ...
+
+template <typename Iterator>
+class repeated_range
+{
+public:
+  using difference_type = typename cuda::std::iterator_traits<Iterator>::difference_type;
+
+  struct repeat_functor
+  {
+    difference_type repeats;
+
+    repeat_functor(difference_type repeats)
+        : repeats(repeats)
+    {}
+
+    __host__ __device__ difference_type operator()(const difference_type& i) const
+    {
+      return i / repeats;
+    }
+  };
+
+  using CountingIterator    = typename thrust::counting_iterator<difference_type>;
+  using TransformIterator   = typename thrust::transform_iterator<repeat_functor, CountingIterator>;
+  using PermutationIterator = typename thrust::permutation_iterator<Iterator, TransformIterator>;
+
+  // type of the repeated_range iterator
+  using iterator = PermutationIterator;
+
+  // construct repeated_range for the range [first,last)
+  repeated_range(Iterator first, Iterator last, difference_type repeats)
+      : first(first)
+      , last(last)
+      , repeats(repeats)
+  {}
+
+  iterator begin() const
+  {
+    return PermutationIterator(first, TransformIterator(CountingIterator(0), repeat_functor(repeats)));
+  }
+
+  iterator end() const
+  {
+    return begin() + repeats * (last - first);
+  }
+
+protected:
+  Iterator first;
+  Iterator last;
+  difference_type repeats;
+};
+
+int main()
+{
+  thrust::device_vector<int> data{10, 20, 30, 40};
+
+  // print the initial data
+  std::cout << "range        ";
+  thrust::copy(data.begin(), data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  using Iterator = thrust::device_vector<int>::iterator;
+
+  // create repeated_range with elements repeated twice
+  repeated_range<Iterator> twice(data.begin(), data.end(), 2);
+  std::cout << "repeated x2: ";
+  thrust::copy(twice.begin(), twice.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // create repeated_range with elements repeated x3
+  repeated_range<Iterator> thrice(data.begin(), data.end(), 3);
+  std::cout << "repeated x3: ";
+  thrust::copy(thrice.begin(), thrice.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/run_length_decoding.cu

@@ -0,0 +1,63 @@
+#include <thrust/binary_search.h>
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/gather.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/scan.h>
+
+#include <iostream>
+#include <iterator>
+
+// This example decodes a run-length code [1] for an array of characters.
+//
+// [1] http://en.wikipedia.org/wiki/Run-length_encoding
+
+int main()
+{
+  // allocate storage for compressed input and run lengths
+  thrust::device_vector<char> input(6);
+  thrust::device_vector<int> lengths(6);
+
+  // clang-format off
+  input[0] = 'a';  lengths[0] = 3;
+  input[1] = 'b';  lengths[1] = 5;
+  input[2] = 'c';  lengths[2] = 1;
+  input[3] = 'd';  lengths[3] = 2;
+  input[4] = 'e';  lengths[4] = 9;
+  input[5] = 'f';  lengths[5] = 2;
+  // clang-format on
+
+  // print the initial data
+  std::cout << "run-length encoded input:" << '\n';
+  for (size_t i = 0; i < 6; i++)
+  {
+    std::cout << "(" << input[i] << "," << lengths[i] << ")";
+  }
+  std::cout << '\n' << '\n';
+
+  // scan the lengths
+  thrust::inclusive_scan(lengths.begin(), lengths.end(), lengths.begin());
+
+  // output size is sum of the run lengths
+  int N = lengths.back();
+
+  // compute input index for each output element
+  thrust::device_vector<int> indices(N);
+  thrust::lower_bound(
+    lengths.begin(),
+    lengths.end(),
+    thrust::counting_iterator<int>(1),
+    thrust::counting_iterator<int>(N + 1),
+    indices.begin());
+
+  // gather input elements
+  thrust::device_vector<char> output(N);
+  thrust::gather(indices.begin(), indices.end(), input.begin(), output.begin());
+
+  // print the initial data
+  std::cout << "decoded output:" << '\n';
+  thrust::copy(output.begin(), output.end(), std::ostream_iterator<char>(std::cout, ""));
+  std::cout << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/run_length_encoding.cu

@@ -0,0 +1,54 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/reduce.h>
+
+#include <cuda/iterator>
+
+#include <iostream>
+#include <iterator>
+
+// This example computes a run-length code [1] for an array of characters.
+//
+// [1] http://en.wikipedia.org/wiki/Run-length_encoding
+
+int main()
+{
+  // input data on the host
+  const char data[] = "aaabbbbbcddeeeeeeeeeff";
+
+  const size_t N = (sizeof(data) / sizeof(char)) - 1;
+
+  // copy input data to the device
+  thrust::device_vector<char> input(data, data + N);
+
+  // allocate storage for output data and run lengths
+  thrust::device_vector<char> output(N);
+  thrust::device_vector<int> lengths(N);
+
+  // print the initial data
+  std::cout << "input data:" << '\n';
+  thrust::copy(input.begin(), input.end(), std::ostream_iterator<char>(std::cout, ""));
+  std::cout << '\n' << '\n';
+
+  // compute run lengths
+  size_t num_runs =
+    thrust::reduce_by_key(
+      input.begin(),
+      input.end(), // input key sequence
+      cuda::constant_iterator<int>(1), // input value sequence
+      output.begin(), // output key sequence
+      lengths.begin() // output value sequence
+      )
+      .first
+    - output.begin(); // compute the output size
+
+  // print the output
+  std::cout << "run-length encoded output:" << '\n';
+  for (size_t i = 0; i < num_runs; i++)
+  {
+    std::cout << "(" << output[i] << "," << lengths[i] << ")";
+  }
+  std::cout << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/saxpy.cu

@@ -0,0 +1,78 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/transform.h>
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+
+// This example illustrates how to implement the SAXPY
+// operation (Y[i] = a * X[i] + Y[i]) using Thrust.
+// The saxpy_slow function demonstrates the most
+// straightforward implementation using a temporary
+// array and two separate transformations, one with
+// multiplies and one with plus.  The saxpy_fast function
+// implements the operation with a single transformation
+// and represents "best practice".
+
+struct saxpy_functor
+{
+  const float a;
+
+  saxpy_functor(float _a)
+      : a(_a)
+  {}
+
+  __host__ __device__ float operator()(const float& x, const float& y) const
+  {
+    return a * x + y;
+  }
+};
+
+void saxpy_fast(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
+{
+  // Y <- A * X + Y
+  thrust::transform(X.begin(), X.end(), Y.begin(), Y.begin(), saxpy_functor(A));
+}
+
+void saxpy_slow(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
+{
+  thrust::device_vector<float> temp(X.size());
+
+  // temp <- A
+  thrust::fill(temp.begin(), temp.end(), A);
+
+  // temp <- A * X
+  thrust::transform(X.begin(), X.end(), temp.begin(), temp.begin(), cuda::std::multiplies<float>());
+
+  // Y <- A * X + Y
+  thrust::transform(temp.begin(), temp.end(), Y.begin(), Y.begin(), cuda::std::plus<float>());
+}
+
+int main()
+{
+  // initialize host arrays
+  thrust::host_vector<float> x{1.0, 1.0, 1.0, 1.0};
+  thrust::host_vector<float> y{1.0, 2.0, 3.0, 4.0};
+
+  {
+    // transfer to device
+    thrust::device_vector<float> X(x);
+    thrust::device_vector<float> Y(y);
+
+    // slow method
+    saxpy_slow(2.0, X, Y);
+  }
+
+  {
+    // transfer to device
+    thrust::device_vector<float> X(x);
+    thrust::device_vector<float> Y(y);
+
+    // fast method
+    saxpy_fast(2.0, X, Y);
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/scan_by_key.cu

@@ -0,0 +1,91 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/scan.h>
+
+#include <iostream>
+
+// BinaryPredicate for the head flag segment representation
+// equivalent to cuda::std::not_fn(thrust::project2nd<int,int>()));
+template <typename HeadFlagType>
+struct head_flag_predicate
+{
+  __host__ __device__ bool operator()(HeadFlagType, HeadFlagType right) const
+  {
+    return !right;
+  }
+};
+
+template <typename Vector>
+void print(const Vector& v)
+{
+  for (const auto& e : v)
+  {
+    std::cout << e << " ";
+  }
+  std::cout << '\n';
+}
+
+int main()
+{
+  int keys[]   = {0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 4, 4, 5, 5, 5}; // segments represented with keys
+  int flags[]  = {1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0}; // segments represented with head flags
+  int values[] = {2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}; // values corresponding to each key
+
+  int N = sizeof(keys) / sizeof(int); // number of elements
+
+  // copy input data to device
+  thrust::device_vector<int> d_keys(keys, keys + N);
+  thrust::device_vector<int> d_flags(flags, flags + N);
+  thrust::device_vector<int> d_values(values, values + N);
+
+  // allocate storage for output
+  thrust::device_vector<int> d_output(N);
+
+  // inclusive scan using keys
+  thrust::inclusive_scan_by_key(d_keys.begin(), d_keys.end(), d_values.begin(), d_output.begin());
+
+  std::cout << "Inclusive Segmented Scan w/ Key Sequence\n";
+  std::cout << " keys          : ";
+  print(d_keys);
+  std::cout << " input values  : ";
+  print(d_values);
+  std::cout << " output values : ";
+  print(d_output);
+
+  // inclusive scan using head flags
+  thrust::inclusive_scan_by_key(
+    d_flags.begin(), d_flags.end(), d_values.begin(), d_output.begin(), head_flag_predicate<int>());
+
+  std::cout << "\nInclusive Segmented Scan w/ Head Flag Sequence\n";
+  std::cout << " head flags    : ";
+  print(d_flags);
+  std::cout << " input values  : ";
+  print(d_values);
+  std::cout << " output values : ";
+  print(d_output);
+
+  // exclusive scan using keys
+  thrust::exclusive_scan_by_key(d_keys.begin(), d_keys.end(), d_values.begin(), d_output.begin());
+
+  std::cout << "\nExclusive Segmented Scan w/ Key Sequence\n";
+  std::cout << " keys          : ";
+  print(d_keys);
+  std::cout << " input values  : ";
+  print(d_values);
+  std::cout << " output values : ";
+  print(d_output);
+
+  // exclusive scan using head flags
+  thrust::exclusive_scan_by_key(
+    d_flags.begin(), d_flags.end(), d_values.begin(), d_output.begin(), 0, head_flag_predicate<int>());
+
+  std::cout << "\nExclusive Segmented Scan w/ Head Flag Sequence\n";
+  std::cout << " head flags    : ";
+  print(d_flags);
+  std::cout << " input values  : ";
+  print(d_values);
+  std::cout << " output values : ";
+  print(d_output);
+
+  return 0;
+}

cccl_upstream/thrust/examples/scan_matrix_by_rows.cu

@@ -0,0 +1,81 @@
+#include <thrust/device_vector.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/scan.h>
+#include <thrust/sequence.h>
+
+#include <cuda/std/cassert>
+
+// We have a matrix stored in a `thrust::device_vector`. We want to perform a
+// scan on each row of a matrix.
+
+__host__ void scan_matrix_by_rows0(thrust::device_vector<int>& u, int n, int m)
+{
+  // Here, we launch a separate scan for each row in the matrix. This works,
+  // but each kernel only does a small amount of work. It would be better if we
+  // could launch one big kernel for the entire matrix.
+  for (int i = 0; i < n; ++i)
+  {
+    thrust::inclusive_scan(
+      u.begin() + m * i, u.begin() + m * (i + 1), u.begin() + m * i); // NOLINT(bugprone-misplaced-widening-cast)
+  }
+}
+
+// We can batch the operation using `thrust::inclusive_scan_by_key`, which
+// scans each group of consecutive equal keys. All we need to do is generate
+// the right key sequence. We want the keys for elements on the same row to
+// be identical.
+
+// So first, we define an unary function object which takes the index of an
+// element and returns the row that it belongs to.
+
+struct which_row
+{
+  int row_length;
+
+  __host__ __device__ which_row(int row_length_)
+      : row_length(row_length_)
+  {}
+
+  __host__ __device__ int operator()(int idx) const
+  {
+    return idx / row_length;
+  }
+};
+
+__host__ void scan_matrix_by_rows1(thrust::device_vector<int>& u, int n, int m)
+{
+  // This `thrust::counting_iterator` represents the index of the element.
+  thrust::counting_iterator<int> c_first(0);
+
+  // We construct a `thrust::transform_iterator` which applies the `which_row`
+  // function object to the index of each element.
+  thrust::transform_iterator<which_row, thrust::counting_iterator<int>> t_first(c_first, which_row(m));
+
+  // Finally, we use our `thrust::transform_iterator` as the key sequence to
+  // `thrust::inclusive_scan_by_key`.
+  thrust::inclusive_scan_by_key(
+    t_first, t_first + n * m, u.begin(), u.begin()); // NOLINT(bugprone-misplaced-widening-cast)
+}
+
+int main()
+{
+  int const n = 4;
+  int const m = 5;
+
+  thrust::device_vector<int> u0(n * m);
+  thrust::sequence(u0.begin(), u0.end());
+  scan_matrix_by_rows0(u0, n, m);
+
+  thrust::device_vector<int> u1(n * m);
+  thrust::sequence(u1.begin(), u1.end());
+  scan_matrix_by_rows1(u1, n, m);
+
+  for (int i = 0; i < n; ++i)
+  {
+    for (int j = 0; j < m; ++j)
+    {
+      assert(u0[j + m * i] == u1[j + m * i]);
+    }
+  }
+}

cccl_upstream/thrust/examples/set_operations.cu

@@ -0,0 +1,157 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/iterator/discard_iterator.h>
+#include <thrust/merge.h>
+#include <thrust/set_operations.h>
+
+#include <iostream>
+
+// This example illustrates use of the set operation algorithms
+//  - merge
+//  - set_union
+//  - set_intersection
+//  - set_difference
+//  - set_symmetric_difference
+//
+// In this context a "set" is simply a sequence of sorted values,
+// allowing the standard set operations to be performed more efficiently
+// than on unsorted data.  Since the output of a set operation is a valid
+// set (i.e. a sorted sequence) it is possible to apply the set operations
+// in a nested fashion to compute arbitrary set expressions.
+//
+// Set operation usage notes:
+//   - The output set size is variable (except for thrust::merge),
+//     so the return value is important.
+//   - Generally one would conservatively allocate storage for the output
+//     and then resize or shrink an output container as necessary.
+//     Alternatively, one can compute the exact output size by
+//     outputting to a discard_iterator.  This approach is more computationally
+//     expensive (approximately 2x), but conserves memory capacity.
+//     Refer to the SetIntersectionSize function for implementation details.
+//   - Sets are allowed to have duplicate elements, which are carried
+//     through to the output in a algorithm-specific manner.  Refer
+//     to the full documentation for precise semantics.
+
+// helper routine
+template <typename String, typename Vector>
+void print(const String& s, const Vector& v)
+{
+  std::cout << s << " [";
+  for (const auto& e : v)
+  {
+    std::cout << " " << e;
+  }
+  std::cout << " ]" << '\n';
+}
+
+template <typename Vector>
+void Merge(const Vector& A, const Vector& B)
+{
+  // merged output is always exactly A.size() + B.size()
+  Vector C(A.size() + B.size());
+
+  thrust::merge(A.begin(), A.end(), B.begin(), B.end(), C.begin());
+
+  print("Merge(A,B)", C);
+}
+
+template <typename Vector>
+void SetUnion(const Vector& A, const Vector& B)
+{
+  // union output is at most A.size() + B.size()
+  Vector C(A.size() + B.size());
+
+  // set_union returns an iterator C_end denoting the end of input
+  typename Vector::iterator C_end;
+
+  C_end = thrust::set_union(A.begin(), A.end(), B.begin(), B.end(), C.begin());
+
+  // shrink C to exactly fit output
+  C.erase(C_end, C.end());
+
+  print("Union(A,B)", C);
+}
+
+template <typename Vector>
+void SetIntersection(const Vector& A, const Vector& B)
+{
+  // intersection output is at most min(A.size(), B.size())
+  Vector C(thrust::min(A.size(), B.size()));
+
+  // set_union returns an iterator C_end denoting the end of input
+  typename Vector::iterator C_end;
+
+  C_end = thrust::set_intersection(A.begin(), A.end(), B.begin(), B.end(), C.begin());
+
+  // shrink C to exactly fit output
+  C.erase(C_end, C.end());
+
+  print("Intersection(A,B)", C);
+}
+
+template <typename Vector>
+void SetDifference(const Vector& A, const Vector& B)
+{
+  // difference output is at most A.size()
+  Vector C(A.size());
+
+  // set_union returns an iterator C_end denoting the end of input
+  typename Vector::iterator C_end;
+
+  C_end = thrust::set_difference(A.begin(), A.end(), B.begin(), B.end(), C.begin());
+
+  // shrink C to exactly fit output
+  C.erase(C_end, C.end());
+
+  print("Difference(A,B)", C);
+}
+
+template <typename Vector>
+void SetSymmetricDifference(const Vector& A, const Vector& B)
+{
+  // symmetric difference output is at most A.size() + B.size()
+  Vector C(A.size() + B.size());
+
+  // set_union returns an iterator C_end denoting the end of input
+  typename Vector::iterator C_end;
+
+  C_end = thrust::set_symmetric_difference(A.begin(), A.end(), B.begin(), B.end(), C.begin());
+
+  // shrink C to exactly fit output
+  C.erase(C_end, C.end());
+
+  print("SymmetricDifference(A,B)", C);
+}
+
+template <typename Vector>
+void SetIntersectionSize(const Vector& A, const Vector& B)
+{
+  // computes the exact size of the intersection without allocating output
+  thrust::discard_iterator<> C_begin, C_end;
+
+  C_end = thrust::set_intersection(A.begin(), A.end(), B.begin(), B.end(), C_begin);
+
+  std::cout << "SetIntersectionSize(A,B) " << (C_end - C_begin) << '\n';
+}
+
+int main()
+{
+  int a[] = {0, 2, 4, 5, 6, 8, 9};
+  int b[] = {0, 1, 2, 3, 5, 7, 8};
+
+  thrust::device_vector<int> A(a, a + sizeof(a) / sizeof(int));
+  thrust::device_vector<int> B(b, b + sizeof(b) / sizeof(int));
+
+  print("Set A", A);
+  print("Set B", B);
+
+  Merge(A, B);
+  SetUnion(A, B);
+  SetIntersection(A, B);
+  SetDifference(A, B);
+  SetSymmetricDifference(A, B);
+
+  SetIntersectionSize(A, B);
+
+  return 0;
+}

cccl_upstream/thrust/examples/simple_moving_average.cu

@@ -0,0 +1,101 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/scan.h>
+#include <thrust/sequence.h>
+#include <thrust/transform.h>
+
+#include <iomanip>
+#include <iostream>
+
+// Efficiently computes the simple moving average (SMA) [1] of a data series
+// using a parallel prefix-sum or "scan" operation.
+//
+// Note: additional numerical precision should be used in the cumulative summing
+// stage when computing the SMA of large data series.  The most straightforward
+// remedy is to replace 'float' with 'double'.   Alternatively a Kahan or
+// "compensated" summation algorithm could be applied [2].
+//
+// [1] http://en.wikipedia.org/wiki/Moving_average#Simple_moving_average
+// [2] http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+
+// compute the difference of two positions in the cumumulative sum and
+// divide by the SMA window size w.
+template <typename T>
+struct minus_and_divide
+{
+  T w;
+
+  minus_and_divide(T w)
+      : w(w)
+  {}
+
+  __host__ __device__ T operator()(const T& a, const T& b) const
+  {
+    return (a - b) / w;
+  }
+};
+
+template <typename InputVector, typename OutputVector>
+void simple_moving_average(const InputVector& data, size_t w, OutputVector& output)
+{
+  using T = typename InputVector::value_type;
+
+  if (data.size() < w)
+  {
+    return;
+  }
+
+  // allocate storage for cumulative sum
+  thrust::device_vector<T> temp(data.size() + 1);
+
+  // compute cumulative sum
+  thrust::exclusive_scan(data.begin(), data.end(), temp.begin());
+  temp[data.size()] = data.back() + temp[data.size() - 1];
+
+  // compute moving averages from cumulative sum
+  thrust::transform(temp.begin() + w, temp.end(), temp.begin(), output.begin(), minus_and_divide<T>(T(w)));
+}
+
+int main()
+{
+  // length of data series
+  size_t n = 30;
+
+  // window size of the moving average
+  size_t w = 4;
+
+  // generate random data series
+  thrust::host_vector<float> host_data(n);
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, 10);
+  for (auto& e : host_data)
+  {
+    e = static_cast<float>(dist(rng));
+  }
+  thrust::device_vector<float> data = host_data;
+
+  // allocate storage for averages
+  thrust::device_vector<float> averages(data.size() - (w - 1));
+
+  // compute SMA using standard summation
+  simple_moving_average(data, w, averages);
+
+  // print data series
+  std::cout << "data series: [ ";
+  for (const auto& value : data)
+  {
+    std::cout << value << " ";
+  }
+  std::cout << "]" << '\n';
+
+  // print moving averages
+  std::cout << "simple moving averages (window = " << w << ")" << '\n';
+  for (size_t i = 0; i < averages.size(); i++)
+  {
+    std::cout << "  [" << std::setw(2) << i << "," << std::setw(2) << (i + w) << ") = " << averages[i] << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/sort.cu

@@ -0,0 +1,194 @@
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/sequence.h>
+#include <thrust/sort.h>
+
+#include <cuda/std/utility>
+
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+
+// Helper routines
+
+void initialize(thrust::device_vector<int>& v)
+{
+  thrust::default_random_engine rng(123456);
+  thrust::uniform_int_distribution<int> dist(10, 99);
+  thrust::host_vector<int> host_data(v.size());
+  for (auto& e : host_data)
+  {
+    e = dist(rng);
+  }
+  v = host_data;
+}
+
+void initialize(thrust::device_vector<float>& v)
+{
+  thrust::default_random_engine rng(123456);
+  thrust::uniform_int_distribution<int> dist(2, 19);
+  thrust::host_vector<float> host_data(v.size());
+  for (auto& e : host_data)
+  {
+    e = static_cast<float>(dist(rng)) / 2.0f;
+  }
+  v = host_data;
+}
+
+void initialize(thrust::device_vector<cuda::std::pair<int, int>>& v)
+{
+  thrust::default_random_engine rng(123456);
+  thrust::uniform_int_distribution<int> dist(0, 9);
+  thrust::host_vector<cuda::std::pair<int, int>> host_data(v.size());
+  for (auto& e : host_data)
+  {
+    int a = dist(rng);
+    int b = dist(rng);
+    e     = cuda::std::make_pair(a, b);
+  }
+  v = host_data;
+}
+
+void initialize(thrust::device_vector<int>& v1, thrust::device_vector<int>& v2)
+{
+  thrust::default_random_engine rng(123456);
+  thrust::uniform_int_distribution<int> dist(10, 99);
+  thrust::host_vector<int> host_data(v1.size());
+  for (auto& e : host_data)
+  {
+    e = dist(rng);
+  }
+  v1 = host_data;
+  thrust::sequence(v2.begin(), v2.end(), 0);
+}
+
+void print(const thrust::device_vector<int>& v)
+{
+  for (const auto& value : v)
+  {
+    std::cout << " " << value;
+  }
+  std::cout << "\n";
+}
+
+void print(const thrust::device_vector<float>& v)
+{
+  for (const auto& value : v)
+  {
+    std::cout << " " << std::fixed << std::setprecision(1) << value;
+  }
+  std::cout << "\n";
+}
+
+void print(const thrust::device_vector<cuda::std::pair<int, int>>& v)
+{
+  for (const auto& p : v)
+  {
+    cuda::std::pair<int, int> local_p = p;
+    std::cout << " (" << local_p.first << "," << local_p.second << ")";
+  }
+  std::cout << "\n";
+}
+
+void print(thrust::device_vector<int>& v1, thrust::device_vector<int> v2)
+{
+  for (size_t i = 0; i < v1.size(); i++)
+  {
+    std::cout << " (" << v1[i] << "," << std::setw(2) << v2[i] << ")";
+  }
+  std::cout << "\n";
+}
+
+// user-defined comparison operator that acts like less<int>,
+// except even numbers are considered to be smaller than odd numbers
+struct evens_before_odds
+{
+  __host__ __device__ bool operator()(int x, int y)
+  {
+    if (x % 2 == y % 2)
+    {
+      return x < y;
+    }
+    else if (x % 2)
+    {
+      return false;
+    }
+    else
+    {
+      return true;
+    }
+  }
+};
+
+int main()
+{
+  size_t N = 16;
+
+  std::cout << "sorting integers\n";
+  {
+    thrust::device_vector<int> keys(N);
+    initialize(keys);
+    print(keys);
+    thrust::sort(keys.begin(), keys.end());
+    print(keys);
+  }
+
+  std::cout << "\nsorting integers (descending)\n";
+  {
+    thrust::device_vector<int> keys(N);
+    initialize(keys);
+    print(keys);
+    thrust::sort(keys.begin(), keys.end(), cuda::std::greater<int>());
+    print(keys);
+  }
+
+  std::cout << "\nsorting integers (user-defined comparison)\n";
+  {
+    thrust::device_vector<int> keys(N);
+    initialize(keys);
+    print(keys);
+    thrust::sort(keys.begin(), keys.end(), evens_before_odds());
+    print(keys);
+  }
+
+  std::cout << "\nsorting floats\n";
+  {
+    thrust::device_vector<float> keys(N);
+    initialize(keys);
+    print(keys);
+    thrust::sort(keys.begin(), keys.end());
+    print(keys);
+  }
+
+  std::cout << "\nsorting pairs\n";
+  {
+    thrust::device_vector<cuda::std::pair<int, int>> keys(N);
+    initialize(keys);
+    print(keys);
+    thrust::sort(keys.begin(), keys.end());
+    print(keys);
+  }
+
+  std::cout << "\nkey-value sorting\n";
+  {
+    thrust::device_vector<int> keys(N);
+    thrust::device_vector<int> values(N);
+    initialize(keys, values);
+    print(keys, values);
+    thrust::sort_by_key(keys.begin(), keys.end(), values.begin());
+    print(keys, values);
+  }
+
+  std::cout << "\nkey-value sorting (descending)\n";
+  {
+    thrust::device_vector<int> keys(N);
+    thrust::device_vector<int> values(N);
+    initialize(keys, values);
+    print(keys, values);
+    thrust::sort_by_key(keys.begin(), keys.end(), values.begin(), cuda::std::greater<int>());
+    print(keys, values);
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/sorting_aos_vs_soa.cu

@@ -0,0 +1,96 @@
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/sort.h>
+
+#include <cuda/std/cassert>
+
+#include <iostream>
+
+#include "include/timer.h"
+
+// This examples compares sorting performance using Array of Structures (AoS)
+// and Structure of Arrays (SoA) data layout.  Legacy applications will often
+// store data in C/C++ structs, such as MyStruct defined below.  Although
+// Thrust can process array of structs, it is typically less efficient than
+// the equivalent structure of arrays layout.  In this particular example,
+// the optimized SoA approach is approximately *five times faster* than the
+// traditional AoS method.  Therefore, it is almost always worthwhile to
+// convert AoS data structures to SoA.
+
+struct MyStruct
+{
+  int key;
+  float value;
+
+  __host__ __device__ bool operator<(const MyStruct other) const
+  {
+    return key < other.key;
+  }
+};
+
+void initialize_keys(thrust::device_vector<int>& keys)
+{
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, 2147483647);
+
+  thrust::host_vector<int> h_keys(keys.size());
+
+  for (auto& e : h_keys)
+  {
+    e = dist(rng);
+  }
+
+  keys = h_keys;
+}
+
+void initialize_keys(thrust::device_vector<MyStruct>& structures)
+{
+  thrust::default_random_engine rng;
+  thrust::uniform_int_distribution<int> dist(0, 2147483647);
+
+  thrust::host_vector<MyStruct> h_structures(structures.size());
+
+  for (auto& s : h_structures)
+  {
+    s.key = dist(rng);
+  }
+
+  structures = h_structures;
+}
+
+int main()
+{
+  size_t N = 2 * 1024 * 1024;
+
+  // Sort Key-Value pairs using Array of Structures (AoS) storage
+  {
+    thrust::device_vector<MyStruct> structures(N);
+
+    initialize_keys(structures);
+
+    timer t;
+
+    thrust::sort(structures.begin(), structures.end());
+    assert(thrust::is_sorted(structures.begin(), structures.end()));
+
+    std::cout << "AoS sort took " << 1e3 * t.elapsed() << " milliseconds" << '\n';
+  }
+
+  // Sort Key-Value pairs using Structure of Arrays (SoA) storage
+  {
+    thrust::device_vector<int> keys(N);
+    thrust::device_vector<float> values(N);
+
+    initialize_keys(keys);
+
+    timer t;
+
+    thrust::sort_by_key(keys.begin(), keys.end(), values.begin());
+    assert(thrust::is_sorted(keys.begin(), keys.end()));
+
+    std::cout << "SoA sort took " << 1e3 * t.elapsed() << " milliseconds" << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/sparse_vector.cu

@@ -0,0 +1,125 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/inner_product.h>
+#include <thrust/merge.h>
+#include <thrust/reduce.h>
+
+#include <cassert>
+#include <iostream>
+
+template <typename IndexVector, typename ValueVector>
+void print_sparse_vector(const IndexVector& A_index, const ValueVector& A_value)
+{
+  assert(A_index.size() == A_value.size());
+
+  for (size_t i = 0; i < A_index.size(); i++)
+  {
+    std::cout << "(" << A_index[i] << "," << A_value[i] << ") ";
+  }
+  std::cout << '\n';
+}
+
+template <typename IndexVector1,
+          typename ValueVector1,
+          typename IndexVector2,
+          typename ValueVector2,
+          typename IndexVector3,
+          typename ValueVector3>
+void sum_sparse_vectors(
+  const IndexVector1& A_index,
+  const ValueVector1& A_value,
+  const IndexVector2& B_index,
+  const ValueVector2& B_value,
+  IndexVector3& C_index,
+  ValueVector3& C_value)
+{
+  using IndexType = typename IndexVector3::value_type;
+  using ValueType = typename ValueVector3::value_type;
+
+  assert(A_index.size() == A_value.size());
+  assert(B_index.size() == B_value.size());
+
+  size_t A_size = A_index.size();
+  size_t B_size = B_index.size();
+
+  // allocate storage for the combined contents of sparse vectors A and B
+  IndexVector3 temp_index(A_size + B_size);
+  ValueVector3 temp_value(A_size + B_size);
+
+  // merge A and B by index
+  thrust::merge_by_key(
+    A_index.begin(),
+    A_index.end(),
+    B_index.begin(),
+    B_index.end(),
+    A_value.begin(),
+    B_value.begin(),
+    temp_index.begin(),
+    temp_value.begin());
+
+  // compute number of unique indices
+  size_t C_size =
+    thrust::inner_product(
+      temp_index.begin(),
+      temp_index.end() - 1,
+      temp_index.begin() + 1,
+      size_t(0),
+      cuda::std::plus<size_t>(),
+      cuda::std::not_equal_to<IndexType>())
+    + 1;
+
+  // allocate space for output
+  C_index.resize(C_size);
+  C_value.resize(C_size);
+
+  // sum values with the same index
+  thrust::reduce_by_key(
+    temp_index.begin(),
+    temp_index.end(),
+    temp_value.begin(),
+    C_index.begin(),
+    C_value.begin(),
+    cuda::std::equal_to<IndexType>(),
+    cuda::std::plus<ValueType>());
+}
+
+int main()
+{
+  // initialize sparse vector A with 4 elements
+  thrust::device_vector<int> A_index(4);
+  thrust::device_vector<float> A_value(4);
+
+  // clang-format off
+  A_index[0] = 2;  A_value[0] = 10;
+  A_index[1] = 3;  A_value[1] = 60;
+  A_index[2] = 5;  A_value[2] = 20;
+  A_index[3] = 8;  A_value[3] = 40;
+  // clang-format on
+
+  // initialize sparse vector B with 6 elements
+  thrust::device_vector<int> B_index(6);
+  thrust::device_vector<float> B_value(6);
+
+  // clang-format off
+  B_index[0] = 1;  B_value[0] = 50;
+  B_index[1] = 2;  B_value[1] = 30;
+  B_index[2] = 4;  B_value[2] = 80;
+  B_index[3] = 5;  B_value[3] = 30;
+  B_index[4] = 7;  B_value[4] = 90;
+  B_index[5] = 8;  B_value[5] = 10;
+  // clang-format on
+
+  // compute sparse vector C = A + B
+  thrust::device_vector<int> C_index;
+  thrust::device_vector<float> C_value;
+
+  sum_sparse_vectors(A_index, A_value, B_index, B_value, C_index, C_value);
+
+  std::cout << "Computing C = A + B for sparse vectors A and B" << '\n';
+  std::cout << "A ";
+  print_sparse_vector(A_index, A_value);
+  std::cout << "B ";
+  print_sparse_vector(B_index, B_value);
+  std::cout << "C ";
+  print_sparse_vector(C_index, C_value);
+}

cccl_upstream/thrust/examples/stream_compaction.cu

@@ -0,0 +1,76 @@
+#include <thrust/copy.h>
+#include <thrust/count.h>
+#include <thrust/device_vector.h>
+#include <thrust/remove.h>
+#include <thrust/sequence.h>
+
+#include <cuda/std/iterator>
+
+#include <iostream>
+#include <string>
+
+// this functor returns true if the argument is odd, and false otherwise
+template <typename T>
+struct is_odd
+{
+  __host__ __device__ bool operator()(T x)
+  {
+    return x % 2;
+  }
+};
+
+template <typename Iterator>
+void print_range(const std::string& name, Iterator first, Iterator last)
+{
+  using T = cuda::std::iter_value_t<Iterator>;
+
+  std::cout << name << ": ";
+  thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));
+  std::cout << "\n";
+}
+
+int main()
+{
+  // input size
+  size_t N = 10;
+
+  // define some types
+  using Vector   = thrust::device_vector<int>;
+  using Iterator = Vector::iterator;
+
+  // allocate storage for array
+  Vector values(N);
+
+  // initialize array to [0, 1, 2, ... ]
+  thrust::sequence(values.begin(), values.end());
+
+  print_range("values", values.begin(), values.end());
+
+  // allocate output storage, here we conservatively assume all values will be copied
+  Vector output(values.size());
+
+  // copy odd numbers to separate array
+  Iterator output_end = thrust::copy_if(values.begin(), values.end(), output.begin(), is_odd<int>());
+
+  print_range("output", output.begin(), output_end);
+
+  // another approach is to count the number of values that will
+  // be copied, and allocate an array of the right size
+  size_t N_odd = thrust::count_if(values.begin(), values.end(), is_odd<int>());
+
+  Vector small_output(N_odd);
+
+  thrust::copy_if(values.begin(), values.end(), small_output.begin(), is_odd<int>());
+
+  print_range("small_output", small_output.begin(), small_output.end());
+
+  // we can also compact sequences with the remove functions, which do the opposite of copy
+  Iterator values_end = thrust::remove_if(values.begin(), values.end(), is_odd<int>());
+
+  // since the values after values_end are garbage, we'll resize the vector
+  values.resize(values_end - values.begin());
+
+  print_range("values", values.begin(), values.end());
+
+  return 0;
+}

cccl_upstream/thrust/examples/strided_range.cu

@@ -0,0 +1,102 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+
+#include <iostream>
+
+// this example illustrates how to make strided access to a range of values
+// examples:
+//   strided_range([0, 1, 2, 3, 4, 5, 6], 1) -> [0, 1, 2, 3, 4, 5, 6]
+//   strided_range([0, 1, 2, 3, 4, 5, 6], 2) -> [0, 2, 4, 6]
+//   strided_range([0, 1, 2, 3, 4, 5, 6], 3) -> [0, 3, 6]
+//   ...
+
+template <typename Iterator>
+class strided_range
+{
+public:
+  using difference_type = typename cuda::std::iterator_traits<Iterator>::difference_type;
+
+  struct stride_functor
+  {
+    difference_type stride;
+
+    stride_functor(difference_type stride)
+        : stride(stride)
+    {}
+
+    __host__ __device__ difference_type operator()(const difference_type& i) const
+    {
+      return stride * i;
+    }
+  };
+
+  using CountingIterator    = typename thrust::counting_iterator<difference_type>;
+  using TransformIterator   = typename thrust::transform_iterator<stride_functor, CountingIterator>;
+  using PermutationIterator = typename thrust::permutation_iterator<Iterator, TransformIterator>;
+
+  // type of the strided_range iterator
+  using iterator = PermutationIterator;
+
+  // construct strided_range for the range [first,last)
+  strided_range(Iterator first, Iterator last, difference_type stride)
+      : first(first)
+      , last(last)
+      , stride(stride)
+  {}
+
+  iterator begin() const
+  {
+    return PermutationIterator(first, TransformIterator(CountingIterator(0), stride_functor(stride)));
+  }
+
+  iterator end() const
+  {
+    return begin() + ((last - first) + (stride - 1)) / stride;
+  }
+
+protected:
+  Iterator first;
+  Iterator last;
+  difference_type stride;
+};
+
+int main()
+{
+  thrust::device_vector<int> data(8);
+  data[0] = 10;
+  data[1] = 20;
+  data[2] = 30;
+  data[3] = 40;
+  data[4] = 50;
+  data[5] = 60;
+  data[6] = 70;
+  data[7] = 80;
+
+  // print the initial data
+  std::cout << "data: ";
+  thrust::copy(data.begin(), data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  using Iterator = thrust::device_vector<int>::iterator;
+
+  // create strided_range with indices [0,2,4,6]
+  strided_range<Iterator> evens(data.begin(), data.end(), 2);
+  std::cout << "sum of even indices: " << thrust::reduce(evens.begin(), evens.end()) << '\n';
+
+  // create strided_range with indices [1,3,5,7]
+  strided_range<Iterator> odds(data.begin() + 1, data.end(), 2);
+  std::cout << "sum of odd indices:  " << thrust::reduce(odds.begin(), odds.end()) << '\n';
+
+  // set odd elements to 0 with fill()
+  std::cout << "setting odd indices to zero: ";
+  thrust::fill(odds.begin(), odds.end(), 0);
+  thrust::copy(data.begin(), data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/sum.cu

@@ -0,0 +1,39 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/generate.h>
+#include <thrust/host_vector.h>
+#include <thrust/random.h>
+#include <thrust/reduce.h>
+
+#include <iostream>
+
+int my_rand()
+{
+  static thrust::default_random_engine rng;
+  static thrust::uniform_int_distribution<int> dist(0, 9999);
+  return dist(rng);
+}
+
+int main()
+{
+  // generate random data on the host
+  thrust::host_vector<int> h_vec(100);
+  thrust::generate(h_vec.begin(), h_vec.end(), my_rand);
+
+  // transfer to device and compute sum
+  thrust::device_vector<int> d_vec = h_vec;
+
+  // initial value of the reduction
+  int init = 0;
+
+  // binary operation used to reduce values
+  cuda::std::plus<int> binary_op;
+
+  // compute sum on the device
+  int sum = thrust::reduce(d_vec.begin(), d_vec.end(), init, binary_op);
+
+  // print the sum
+  std::cout << "sum is " << sum << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/sum_columns.cmake

@@ -0,0 +1,31 @@
+# GCC 7 has issues with capturing lambdas in non-CUDA backends triggering unused parameter warnings
+# in libcudacxx's __destroy_at. Disable this example for GCC 7.
+if (
+  "GNU" STREQUAL "${CMAKE_CXX_COMPILER_ID}"
+  AND CMAKE_CXX_COMPILER_VERSION VERSION_LESS 8
+)
+  set_target_properties(${example_target} PROPERTIES EXCLUDE_FROM_ALL TRUE)
+  set_tests_properties(${example_target} PROPERTIES DISABLED TRUE)
+  return()
+endif()
+
+target_compile_options(
+  ${example_target}
+  PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: --extended-lambda>
+)
+# This check is actually not correct, because we must check the host compiler, not the CXX compiler.
+# We rely on these usually being the same ;)
+if (
+  "Clang" STREQUAL "${CMAKE_CXX_COMPILER_ID}"
+  AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 13
+)
+  # When clang >= 13 is used as host compiler, we get the following warning:
+  #   nvcc_internal_extended_lambda_implementation:312:22: error: definition of implicit copy constructor for '__nv_hdl_wrapper_t<false, true, false, __nv_dl_tag<void (*)(), &TestAddressStabilityLambda, 2>, int (const int &)>' is deprecated because it has a user-declared copy assignment operator [-Werror,-Wdeprecated-copy]
+  #   312 | __nv_hdl_wrapper_t & operator=(const __nv_hdl_wrapper_t &in) = delete;
+  #       |                      ^
+  # Let's suppress it until NVBug 4980157 is resolved.
+  target_compile_options(
+    ${example_target}
+    PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: -Wno-deprecated-copy>
+  )
+endif()

cccl_upstream/thrust/examples/sum_columns.cu

@@ -0,0 +1,78 @@
+#include <thrust/execution_policy.h>
+#include <thrust/for_each.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/discard_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/random.h>
+#include <thrust/reduce.h>
+#include <thrust/tabulate.h>
+#include <thrust/universal_vector.h>
+
+#include <cuda/iterator>
+#include <cuda/std/mdspan>
+
+#include <iomanip>
+#include <iostream>
+
+#include "include/host_device.h"
+
+int main()
+{
+  const int rows = 32;
+  const int cols = 16;
+
+  // Create a 2D multidimensional array of ints.
+  thrust::universal_vector<int> data(rows * cols, 42);
+  cuda::std::mdspan M(thrust::raw_pointer_cast(data.data()), rows, cols);
+
+  // Create an iterator to the flat linear index space for the multidimensional array.
+  auto flat_idx = cuda::counting_iterator(0);
+
+  // Fill the array with pseudorandom inputs in parallel on the device.
+  thrust::tabulate(thrust::device, M.data_handle(), M.data_handle() + M.size(), [] __host__ __device__(int flat) {
+    thrust::default_random_engine rng;
+    thrust::uniform_int_distribution<int> dist(0, 3);
+    rng.discard(flat); // Advance to the current element's position.
+    return dist(rng);
+  });
+
+  // Create a range to the column index of each element.
+  auto col_idx_begin = thrust::make_transform_iterator(flat_idx, [=] __host__ __device__(int flat) {
+    return flat / rows;
+  });
+  auto col_idx_end   = col_idx_begin + static_cast<std::ptrdiff_t>(M.size());
+
+  // Create a transposed view of the multidimensional array.
+  auto M_transposed = thrust::make_permutation_iterator(
+    M.data_handle(), thrust::make_transform_iterator(cuda::counting_iterator(0), [=] __host__ __device__(int flat) {
+      int i = flat / cols;
+      int j = flat % cols;
+      return i + j * rows;
+    }));
+
+  // Sum each column, storing the result in a new vector.
+  thrust::universal_vector<int> sums(cols);
+  thrust::reduce_by_key(
+    thrust::device, col_idx_begin, col_idx_end, M_transposed, thrust::make_discard_iterator(), sums.begin());
+
+  // Output the result.
+  thrust::for_each_n(thrust::seq, flat_idx, rows, [&](int i) {
+    std::cout << "[ ";
+    thrust::for_each_n(thrust::seq, flat_idx, cols, [&](int j) {
+      std::cout << std::setw(2) << M(i, j) << " ";
+    });
+    std::cout << "]\n";
+  });
+
+  std::cout << "  ";
+  thrust::for_each_n(thrust::seq, flat_idx, cols, [&](int) {
+    std::cout << " = ";
+  });
+  std::cout << "\n";
+
+  std::cout << "  ";
+  thrust::for_each_n(thrust::seq, flat_idx, cols, [&](int j) {
+    std::cout << std::setw(2) << sums[j] << " ";
+  });
+  std::cout << "\n";
+}

cccl_upstream/thrust/examples/sum_rows.cmake

@@ -0,0 +1,31 @@
+# GCC 7 has issues with capturing lambdas in non-CUDA backends triggering unused parameter warnings
+# in libcudacxx's __destroy_at. Disable this example for GCC 7.
+if (
+  "GNU" STREQUAL "${CMAKE_CXX_COMPILER_ID}"
+  AND CMAKE_CXX_COMPILER_VERSION VERSION_LESS 8
+)
+  set_target_properties(${example_target} PROPERTIES EXCLUDE_FROM_ALL TRUE)
+  set_tests_properties(${example_target} PROPERTIES DISABLED TRUE)
+  return()
+endif()
+
+target_compile_options(
+  ${example_target}
+  PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: --extended-lambda>
+)
+# This check is actually not correct, because we must check the host compiler, not the CXX compiler.
+# We rely on these usually being the same ;)
+if (
+  "Clang" STREQUAL "${CMAKE_CXX_COMPILER_ID}"
+  AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 13
+)
+  # When clang >= 13 is used as host compiler, we get the following warning:
+  #   nvcc_internal_extended_lambda_implementation:312:22: error: definition of implicit copy constructor for '__nv_hdl_wrapper_t<false, true, false, __nv_dl_tag<void (*)(), &TestAddressStabilityLambda, 2>, int (const int &)>' is deprecated because it has a user-declared copy assignment operator [-Werror,-Wdeprecated-copy]
+  #   312 | __nv_hdl_wrapper_t & operator=(const __nv_hdl_wrapper_t &in) = delete;
+  #       |                      ^
+  # Let's suppress it until NVBug 4980157 is resolved.
+  target_compile_options(
+    ${example_target}
+    PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: -Wno-deprecated-copy>
+  )
+endif()

cccl_upstream/thrust/examples/sum_rows.cu

@@ -0,0 +1,58 @@
+#include <thrust/execution_policy.h>
+#include <thrust/for_each.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/discard_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/random.h>
+#include <thrust/reduce.h>
+#include <thrust/tabulate.h>
+#include <thrust/universal_vector.h>
+
+#include <cuda/iterator>
+#include <cuda/std/mdspan>
+
+#include <iomanip>
+#include <iostream>
+
+#include "include/host_device.h"
+
+int main()
+{
+  const int rows = 32;
+  const int cols = 16;
+
+  // Create a 2D multidimensional array of ints.
+  thrust::universal_vector<int> data(rows * cols, 42);
+  cuda::std::mdspan M(thrust::raw_pointer_cast(data.data()), rows, cols);
+
+  // Create an iterator to the flat linear index space for the multidimensional array.
+  auto flat_idx = cuda::counting_iterator(0);
+
+  // Fill the array with pseudorandom inputs in parallel on the device.
+  thrust::tabulate(thrust::device, M.data_handle(), M.data_handle() + M.size(), [] __host__ __device__(int flat) {
+    thrust::default_random_engine rng;
+    thrust::uniform_int_distribution<int> dist(0, 3);
+    rng.discard(flat); // Advance to the current element's position.
+    return dist(rng);
+  });
+
+  // Create a range to the row index of each element.
+  auto row_idx_begin = thrust::make_transform_iterator(flat_idx, [=] __host__ __device__(int flat) {
+    return flat / cols;
+  });
+  auto row_idx_end   = row_idx_begin + static_cast<std::ptrdiff_t>(M.size());
+
+  // Sum each row, storing the result in a new vector.
+  thrust::universal_vector<int> sums(rows);
+  thrust::reduce_by_key(
+    thrust::device, row_idx_begin, row_idx_end, M.data_handle(), thrust::make_discard_iterator(), sums.begin());
+
+  // Output the result.
+  thrust::for_each_n(thrust::seq, flat_idx, rows, [&](int i) {
+    std::cout << "[ ";
+    thrust::for_each_n(thrust::seq, flat_idx, cols, [&](int j) {
+      std::cout << std::setw(2) << M(i, j) << " ";
+    });
+    std::cout << "] = " << sums[i] << "\n";
+  });
+}

cccl_upstream/thrust/examples/summary_statistics.cu

@@ -0,0 +1,166 @@
+#include <thrust/device_vector.h>
+#include <thrust/extrema.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/transform_reduce.h>
+
+#include <cuda/std/iterator>
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+// This example computes several statistical properties of a data
+// series in a single reduction.  The algorithm is described in detail here:
+// http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
+//
+// Thanks to Joseph Rhoads for contributing this example
+
+// structure used to accumulate the moments and other
+// statistical properties encountered so far.
+template <typename T>
+struct summary_stats_data
+{
+  T n;
+  T min;
+  T max;
+  T mean;
+  T M2;
+  T M3;
+  T M4;
+
+  // initialize to the identity element
+  void initialize()
+  {
+    n = mean = M2 = M3 = M4 = 0;
+    min                     = std::numeric_limits<T>::max();
+    max                     = std::numeric_limits<T>::min();
+  }
+
+  T variance()
+  {
+    return M2 / (n - 1);
+  }
+  T variance_n()
+  {
+    return M2 / n;
+  }
+  T skewness()
+  {
+    return std::sqrt(n) * M3 / std::pow(M2, (T) 1.5);
+  }
+  T kurtosis()
+  {
+    return n * M4 / (M2 * M2);
+  }
+};
+
+// stats_unary_op is a functor that takes in a value x and
+// returns a variace_data whose mean value is initialized to x.
+template <typename T>
+struct summary_stats_unary_op
+{
+  __host__ __device__ summary_stats_data<T> operator()(const T& x) const
+  {
+    summary_stats_data<T> result;
+    result.n    = 1;
+    result.min  = x;
+    result.max  = x;
+    result.mean = x;
+    result.M2   = 0;
+    result.M3   = 0;
+    result.M4   = 0;
+
+    return result;
+  }
+};
+
+// summary_stats_binary_op is a functor that accepts two summary_stats_data
+// structs and returns a new summary_stats_data which are an
+// approximation to the summary_stats for
+// all values that have been aggregated so far
+template <typename T>
+struct summary_stats_binary_op
+{
+  __host__ __device__ summary_stats_data<T>
+  operator()(const summary_stats_data<T>& x, const summary_stats_data<T>& y) const
+  {
+    summary_stats_data<T> result;
+
+    // precompute some common subexpressions
+    T n  = x.n + y.n;
+    T n2 = n * n;
+    T n3 = n2 * n;
+
+    T delta  = y.mean - x.mean;
+    T delta2 = delta * delta;
+    T delta3 = delta2 * delta;
+    T delta4 = delta3 * delta;
+
+    // Basic number of samples (n), min, and max
+    result.n   = n;
+    result.min = thrust::min(x.min, y.min);
+    result.max = thrust::max(x.max, y.max);
+
+    result.mean = x.mean + delta * y.n / n;
+
+    result.M2 = x.M2 + y.M2;
+    result.M2 += delta2 * x.n * y.n / n;
+
+    result.M3 = x.M3 + y.M3;
+    result.M3 += delta3 * x.n * y.n * (x.n - y.n) / n2;
+    result.M3 += (T) 3.0 * delta * (x.n * y.M2 - y.n * x.M2) / n;
+
+    result.M4 = x.M4 + y.M4;
+    result.M4 += delta4 * x.n * y.n * (x.n * x.n - x.n * y.n + y.n * y.n) / n3;
+    result.M4 += (T) 6.0 * delta2 * (x.n * x.n * y.M2 + y.n * y.n * x.M2) / n2;
+    result.M4 += (T) 4.0 * delta * (x.n * y.M3 - y.n * x.M3) / n;
+
+    return result;
+  }
+};
+
+template <typename Iterator>
+void print_range(const std::string& name, Iterator first, Iterator last)
+{
+  using T = cuda::std::iter_value_t<Iterator>;
+
+  std::cout << name << ": ";
+  thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));
+  std::cout << "\n";
+}
+
+int main()
+{
+  using T = float;
+
+  // initialize host array
+  thrust::host_vector<T> h_x{4, 7, 13, 16};
+
+  // transfer to device
+  thrust::device_vector<T> d_x(h_x);
+
+  // setup arguments
+  summary_stats_unary_op<T> unary_op;
+  summary_stats_binary_op<T> binary_op;
+  summary_stats_data<T> init;
+
+  init.initialize();
+
+  // compute summary statistics
+  summary_stats_data<T> result = thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op);
+
+  std::cout << "******Summary Statistics Example*****" << '\n';
+  print_range("The data", d_x.begin(), d_x.end());
+
+  std::cout << "Count              : " << result.n << '\n';
+  std::cout << "Minimum            : " << result.min << '\n';
+  std::cout << "Maximum            : " << result.max << '\n';
+  std::cout << "Mean               : " << result.mean << '\n';
+  std::cout << "Variance           : " << result.variance() << '\n';
+  std::cout << "Standard Deviation : " << std::sqrt(result.variance_n()) << '\n';
+  std::cout << "Skewness           : " << result.skewness() << '\n';
+  std::cout << "Kurtosis           : " << result.kurtosis() << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/summed_area_table.cu

@@ -0,0 +1,119 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/gather.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/scan.h>
+
+#include <iomanip>
+#include <iostream>
+
+// This example computes a summed area table using segmented scan
+// http://en.wikipedia.org/wiki/Summed_area_table
+
+// convert a linear index to a linear index in the transpose
+struct transpose_index
+{
+  size_t m, n;
+
+  __host__ __device__ transpose_index(size_t _m, size_t _n)
+      : m(_m)
+      , n(_n)
+  {}
+
+  __host__ __device__ size_t operator()(size_t linear_index)
+  {
+    size_t i = linear_index / n;
+    size_t j = linear_index % n;
+
+    return m * j + i;
+  }
+};
+
+// convert a linear index to a row index
+struct row_index
+{
+  size_t n;
+
+  __host__ __device__ row_index(size_t _n)
+      : n(_n)
+  {}
+
+  __host__ __device__ size_t operator()(size_t i)
+  {
+    return i / n;
+  }
+};
+
+// transpose an M-by-N array
+template <typename T>
+void transpose(size_t m, size_t n, thrust::device_vector<T>& src, thrust::device_vector<T>& dst)
+{
+  thrust::counting_iterator<size_t> indices(0);
+
+  thrust::gather(thrust::make_transform_iterator(indices, transpose_index(n, m)),
+                 thrust::make_transform_iterator(indices, transpose_index(n, m)) + dst.size(),
+                 src.begin(),
+                 dst.begin());
+}
+
+// scan the rows of an M-by-N array
+template <typename T>
+void scan_horizontally(size_t n, thrust::device_vector<T>& d_data)
+{
+  thrust::counting_iterator<size_t> indices(0);
+
+  thrust::inclusive_scan_by_key(
+    thrust::make_transform_iterator(indices, row_index(n)),
+    thrust::make_transform_iterator(indices, row_index(n)) + d_data.size(),
+    d_data.begin(),
+    d_data.begin());
+}
+
+// print an M-by-N array
+template <typename T>
+void print(size_t m, size_t n, thrust::device_vector<T>& d_data)
+{
+  thrust::host_vector<T> h_data = d_data;
+
+  for (size_t i = 0; i < m; i++)
+  {
+    for (size_t j = 0; j < n; j++)
+    {
+      std::cout << std::setw(8) << h_data[i * n + j] << " ";
+    }
+    std::cout << "\n";
+  }
+}
+
+int main()
+{
+  size_t m = 3; // number of rows
+  size_t n = 4; // number of columns
+
+  // 2d array stored in row-major order [(0,0), (0,1), (0,2) ... ]
+  thrust::device_vector<int> data(m * n, 1);
+
+  std::cout << "[step 0] initial array" << '\n';
+  print(m, n, data);
+
+  std::cout << "[step 1] scan horizontally" << '\n';
+  scan_horizontally(n, data);
+  print(m, n, data);
+
+  std::cout << "[step 2] transpose array" << '\n';
+  thrust::device_vector<int> temp(m * n);
+  transpose(m, n, data, temp);
+  print(n, m, temp);
+
+  std::cout << "[step 3] scan transpose horizontally" << '\n';
+  scan_horizontally(m, temp);
+  print(n, m, temp);
+
+  std::cout << "[step 4] transpose the transpose" << '\n';
+  transpose(n, m, temp, data);
+  print(m, n, data);
+
+  return 0;
+}

cccl_upstream/thrust/examples/tiled_range.cu

@@ -0,0 +1,92 @@
+#include <thrust/copy.h>
+#include <thrust/device_vector.h>
+#include <thrust/fill.h>
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/permutation_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+
+#include <iostream>
+
+// this example illustrates how to tile a range multiple times
+// examples:
+//   tiled_range([0, 1, 2, 3], 1) -> [0, 1, 2, 3]
+//   tiled_range([0, 1, 2, 3], 2) -> [0, 1, 2, 3, 0, 1, 2, 3]
+//   tiled_range([0, 1, 2, 3], 3) -> [0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]
+//   ...
+
+template <typename Iterator>
+class tiled_range
+{
+public:
+  using difference_type = typename cuda::std::iterator_traits<Iterator>::difference_type;
+
+  struct tile_functor
+  {
+    difference_type tile_size;
+
+    tile_functor(difference_type tile_size)
+        : tile_size(tile_size)
+    {}
+
+    __host__ __device__ difference_type operator()(const difference_type& i) const
+    {
+      return i % tile_size;
+    }
+  };
+
+  using CountingIterator    = typename thrust::counting_iterator<difference_type>;
+  using TransformIterator   = typename thrust::transform_iterator<tile_functor, CountingIterator>;
+  using PermutationIterator = typename thrust::permutation_iterator<Iterator, TransformIterator>;
+
+  // type of the tiled_range iterator
+  using iterator = PermutationIterator;
+
+  // construct repeated_range for the range [first,last)
+  tiled_range(Iterator first, Iterator last, difference_type tiles)
+      : first(first)
+      , last(last)
+      , tiles(tiles)
+  {}
+
+  iterator begin() const
+  {
+    return PermutationIterator(first, TransformIterator(CountingIterator(0), tile_functor(last - first)));
+  }
+
+  iterator end() const
+  {
+    return begin() + tiles * (last - first);
+  }
+
+protected:
+  Iterator first;
+  Iterator last;
+  difference_type tiles;
+};
+
+int main()
+{
+  thrust::device_vector<int> data{10, 20, 30, 40};
+
+  // print the initial data
+  std::cout << "range        ";
+  thrust::copy(data.begin(), data.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  using Iterator = thrust::device_vector<int>::iterator;
+
+  // create tiled_range with two tiles
+  tiled_range<Iterator> two(data.begin(), data.end(), 2);
+  std::cout << "two tiles:   ";
+  thrust::copy(two.begin(), two.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  // create tiled_range with three tiles
+  tiled_range<Iterator> three(data.begin(), data.end(), 3);
+  std::cout << "three tiles: ";
+  thrust::copy(three.begin(), three.end(), std::ostream_iterator<int>(std::cout, " "));
+  std::cout << '\n';
+
+  return 0;
+}

cccl_upstream/thrust/examples/transform_input_output_iterator.cu

@@ -0,0 +1,117 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/gather.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/transform_input_output_iterator.h>
+#include <thrust/sequence.h>
+
+#include <iostream>
+
+// Base 2 fixed point
+class ScaledInteger
+{
+  int value_;
+  int scale_;
+
+public:
+  __host__ __device__ ScaledInteger(int value, int scale)
+      : value_{value}
+      , scale_{scale}
+  {}
+
+  __host__ __device__ int value() const
+  {
+    return value_;
+  }
+
+  __host__ __device__ ScaledInteger rescale(int scale) const
+  {
+    int shift  = scale - scale_;
+    int result = shift < 0 ? value_ << (-shift) : value_ >> shift;
+    return ScaledInteger{result, scale};
+  }
+
+  __host__ __device__ friend ScaledInteger operator+(ScaledInteger a, ScaledInteger b)
+  {
+    // Rescale inputs to the lesser of the two scales
+    if (b.scale_ < a.scale_)
+    {
+      a = a.rescale(b.scale_);
+    }
+    else if (a.scale_ < b.scale_)
+    {
+      b = b.rescale(a.scale_);
+    }
+    return ScaledInteger{a.value_ + b.value_, a.scale_};
+  }
+};
+
+struct ValueToScaledInteger
+{
+  int scale;
+
+  __host__ __device__ ScaledInteger operator()(const int& value) const
+  {
+    return ScaledInteger{value, scale};
+  }
+};
+
+struct ScaledIntegerToValue
+{
+  int scale;
+
+  __host__ __device__ int operator()(const ScaledInteger& scaled) const
+  {
+    return scaled.rescale(scale).value();
+  }
+};
+
+int main()
+{
+  const size_t size = 4;
+  thrust::device_vector<int> A(size);
+  thrust::device_vector<int> B(size);
+  thrust::device_vector<int> C(size);
+
+  thrust::sequence(A.begin(), A.end(), 1);
+  thrust::sequence(B.begin(), B.end(), 5);
+
+  const int A_scale = 16; // Values in A are left shifted by 16
+  const int B_scale = 8; // Values in B are left shifted by 8
+  const int C_scale = 4; // Values in C are left shifted by 4
+
+  auto A_begin = thrust::make_transform_input_output_iterator(
+    A.begin(), ValueToScaledInteger{A_scale}, ScaledIntegerToValue{A_scale});
+  auto A_end =
+    thrust::make_transform_input_output_iterator(A.end(), ValueToScaledInteger{A_scale}, ScaledIntegerToValue{A_scale});
+  auto B_begin = thrust::make_transform_input_output_iterator(
+    B.begin(), ValueToScaledInteger{B_scale}, ScaledIntegerToValue{B_scale});
+  auto C_begin = thrust::make_transform_input_output_iterator(
+    C.begin(), ValueToScaledInteger{C_scale}, ScaledIntegerToValue{C_scale});
+
+  // Sum A and B as ScaledIntegers, storing the scaled result in C
+  thrust::transform(A_begin, A_end, B_begin, C_begin, cuda::std::plus<ScaledInteger>{});
+
+  thrust::host_vector<int> A_h(A);
+  thrust::host_vector<int> B_h(B);
+  thrust::host_vector<int> C_h(C);
+
+  std::cout << std::hex;
+
+  std::cout << "Expected [ ";
+  for (size_t i = 0; i < size; i++)
+  {
+    const int expected = ((A_h[i] << A_scale) + (B_h[i] << B_scale)) >> C_scale;
+    std::cout << expected << " ";
+  }
+  std::cout << "] \n";
+
+  std::cout << "Result   [ ";
+  for (const auto& value : C_h)
+  {
+    std::cout << value << " ";
+  }
+  std::cout << "] \n";
+
+  return 0;
+}

cccl_upstream/thrust/examples/transform_iterator.cu

@@ -0,0 +1,131 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/iterator/counting_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/reduce.h>
+
+#include <cuda/std/iterator>
+
+#include <iostream>
+#include <string>
+
+// this functor clamps a value to the range [lo, hi]
+template <typename T>
+struct clamp
+{
+  T lo, hi;
+
+  __host__ __device__ clamp(T _lo, T _hi)
+      : lo(_lo)
+      , hi(_hi)
+  {}
+
+  __host__ __device__ T operator()(T x)
+  {
+    if (x < lo)
+    {
+      return lo;
+    }
+    else if (x < hi)
+    {
+      return x;
+    }
+    else
+    {
+      return hi;
+    }
+  }
+};
+
+template <typename T>
+struct simple_negate
+{
+  __host__ __device__ T operator()(T x)
+  {
+    return -x;
+  }
+};
+
+template <typename Iterator>
+void print_range(const std::string& name, Iterator first, Iterator last)
+{
+  using T = cuda::std::iter_value_t<Iterator>;
+
+  std::cout << name << ": ";
+  thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));
+  std::cout << "\n";
+}
+
+int main()
+{
+  // clamp values to the range [1, 5]
+  int lo = 1;
+  int hi = 5;
+
+  // define some types
+  using Vector         = thrust::device_vector<int>;
+  using VectorIterator = Vector::iterator;
+
+  // initialize values
+  Vector values(8);
+
+  values[0] = 2;
+  values[1] = 5;
+  values[2] = 7;
+  values[3] = 1;
+  values[4] = 6;
+  values[5] = 0;
+  values[6] = 3;
+  values[7] = 8;
+
+  print_range("values         ", values.begin(), values.end());
+
+  // define some more types
+  using ClampedVectorIterator = thrust::transform_iterator<clamp<int>, VectorIterator>;
+
+  // create a transform_iterator that applies clamp() to the values array
+  ClampedVectorIterator cv_begin = thrust::make_transform_iterator(values.begin(), clamp<int>(lo, hi));
+  ClampedVectorIterator cv_end   = cv_begin + static_cast<std::ptrdiff_t>(values.size());
+
+  // now [clamped_begin, clamped_end) defines a sequence of clamped values
+  print_range("clamped values ", cv_begin, cv_end);
+
+  ////
+  // compute the sum of the clamped sequence with reduce()
+  std::cout << "sum of clamped values : " << thrust::reduce(cv_begin, cv_end) << "\n";
+
+  ////
+  // combine transform_iterator with other fancy iterators like counting_iterator
+  using CountingIterator        = thrust::counting_iterator<int>;
+  using ClampedCountingIterator = thrust::transform_iterator<clamp<int>, CountingIterator>;
+
+  CountingIterator count_begin(0);
+  CountingIterator count_end(10);
+
+  print_range("sequence         ", count_begin, count_end);
+
+  ClampedCountingIterator cs_begin = thrust::make_transform_iterator(count_begin, clamp<int>(lo, hi));
+  ClampedCountingIterator cs_end   = thrust::make_transform_iterator(count_end, clamp<int>(lo, hi));
+
+  print_range("clamped sequence ", cs_begin, cs_end);
+
+  ////
+  // combine transform_iterator with another transform_iterator
+  using NegatedClampedCountingIterator = thrust::transform_iterator<cuda::std::negate<int>, ClampedCountingIterator>;
+
+  NegatedClampedCountingIterator ncs_begin = thrust::make_transform_iterator(cs_begin, cuda::std::negate<int>());
+  NegatedClampedCountingIterator ncs_end   = thrust::make_transform_iterator(cs_end, cuda::std::negate<int>());
+
+  print_range("negated sequence ", ncs_begin, ncs_end);
+
+  ////
+  // when a functor does not define result_type, a third template argument must be provided
+  using NegatedVectorIterator = thrust::transform_iterator<simple_negate<int>, VectorIterator, int>;
+
+  NegatedVectorIterator nv_begin(values.begin(), simple_negate<int>());
+  NegatedVectorIterator nv_end(values.end(), simple_negate<int>());
+
+  print_range("negated values ", nv_begin, nv_end);
+
+  return 0;
+}

cccl_upstream/thrust/examples/transform_output_iterator.cu

@@ -0,0 +1,46 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/gather.h>
+#include <thrust/host_vector.h>
+#include <thrust/iterator/transform_output_iterator.h>
+
+#include <iostream>
+
+struct Functor
+{
+  template <class Tuple>
+  __host__ __device__ float operator()(const Tuple& tuple) const
+  {
+    const float x = cuda::std::get<0>(tuple);
+    const float y = cuda::std::get<1>(tuple);
+    return x * y * 2.0f / 3.0f;
+  }
+};
+
+int main()
+{
+  thrust::host_vector<float> u{4, 3, 2, 1};
+  thrust::host_vector<float> v{-1, 1, 1, -1};
+  thrust::host_vector<int> idx{3, 0, 1};
+  thrust::host_vector<float> w{0, 0, 0};
+
+  thrust::device_vector<float> U(u);
+  thrust::device_vector<float> V(v);
+  thrust::device_vector<int> IDX(idx);
+  thrust::device_vector<float> W(w);
+
+  // gather multiple elements and apply a function before writing result in memory
+  thrust::gather(IDX.begin(),
+                 IDX.end(),
+                 thrust::make_zip_iterator(U.begin(), V.begin()),
+                 thrust::make_transform_output_iterator(W.begin(), Functor()));
+
+  std::cout << "result= [ ";
+  for (const auto& value : W)
+  {
+    std::cout << value << " ";
+  }
+  std::cout << "] \n";
+
+  return 0;
+}

cccl_upstream/thrust/examples/uninitialized_vector.cu

@@ -0,0 +1,73 @@
+// Occasionally, it is advantageous to avoid initializing the individual
+// elements of a device_vector. For example, the default behavior of
+// zero-initializing numeric data may introduce undesirable overhead.
+// This example demonstrates how to avoid default construction of a
+// device_vector's data by using a custom allocator.
+
+#include <thrust/device_allocator.h>
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/logical.h>
+
+#include <cassert>
+
+// uninitialized_allocator is an allocator which
+// derives from device_allocator and which has a
+// no-op construct member function
+template <typename T>
+struct uninitialized_allocator : thrust::device_allocator<T>
+{
+  // the default generated constructors and destructors are implicitly
+  // marked __host__ __device__, but the current Thrust device_allocator
+  // can only be constructed and destroyed on the host; therefore, we
+  // define these as host only
+  __host__ uninitialized_allocator() {} // NOLINT(modernize-use-equals-default)
+  __host__ uninitialized_allocator(const uninitialized_allocator& other)
+      : thrust::device_allocator<T>(other)
+  {}
+  __host__ ~uninitialized_allocator() {} // NOLINT(modernize-use-equals-default)
+
+  uninitialized_allocator& operator=(const uninitialized_allocator&) = default;
+
+  // for correctness, you should also redefine rebind when you inherit
+  // from an allocator type; this way, if the allocator is rebound somewhere,
+  // it's going to be rebound to the correct type - and not to its base
+  // type for U
+  template <typename U>
+  struct rebind
+  {
+    using other = uninitialized_allocator<U>;
+  };
+
+  // note that construct is annotated as
+  // a __host__ __device__ function
+  __host__ __device__ void construct(T*)
+  {
+    // no-op
+  }
+};
+
+// to make a device_vector which does not initialize its elements,
+// use uninitialized_allocator as the 2nd template parameter
+using uninitialized_vector = thrust::device_vector<float, uninitialized_allocator<float>>;
+
+int main()
+{
+  uninitialized_vector vec(10);
+
+  // the initial value of vec's 10 elements is undefined
+
+  // resize without default value does not initialize elements
+  vec.resize(20);
+
+  // resize with default value does initialize elements
+  vec.resize(30, 13);
+
+  // the value of elements [0,20) is still undefined
+  // but the value of elements [20,30) is 13:
+
+  using namespace thrust::placeholders;
+  assert(thrust::all_of(vec.begin() + 20, vec.end(), _1 == 13));
+
+  return 0;
+}

cccl_upstream/thrust/examples/weld_vertices.cu

@@ -0,0 +1,84 @@
+#include <thrust/binary_search.h>
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/remove.h>
+#include <thrust/sort.h>
+#include <thrust/unique.h>
+
+#include <iostream>
+
+/*
+ * This example "welds" triangle vertices together by taking as
+ * input "triangle soup" and eliminating redundant vertex positions
+ * and shared edges.  A connected mesh is the result.
+ *
+ *
+ * Input: 9 vertices representing a mesh with 3 triangles
+ *
+ *  Mesh              Vertices
+ *    ------           (2)      (5)--(4)    (8)
+ *    | \ 2| \          | \       \   |      | \
+ *    |  \ |  \   <->   |  \       \  |      |  \
+ *    | 0 \| 1 \        |   \       \ |      |   \
+ *    -----------      (0)--(1)      (3)    (6)--(7)
+ *
+ *   (vertex 1 equals vertex 3, vertex 2 equals vertex 5, ...)
+ *
+ * Output: mesh representation with 5 vertices and 9 indices
+ *
+ *  Vertices            Indices
+ *   (1)--(3)            [(0,2,1),
+ *    | \  | \            (2,3,1),
+ *    |  \ |  \           (2,4,3)]
+ *    |   \|   \
+ *   (0)--(2)--(4)
+ */
+
+// define a 2d float vector
+using vec2 = cuda::std::tuple<float, float>;
+
+int main()
+{
+  // allocate memory for input mesh representation
+  thrust::device_vector<vec2> input(9);
+
+  thrust::host_vector<vec2> h_input{
+    vec2(0, 0),
+    vec2(1, 0),
+    vec2(0, 1), // First Triangle
+    vec2(1, 0),
+    vec2(1, 1),
+    vec2(0, 1), // Second Triangle
+    vec2(1, 0),
+    vec2(2, 0),
+    vec2(1, 1) // Third Triangle
+  };
+  input = h_input;
+
+  // allocate space for output mesh representation
+  thrust::device_vector<vec2> vertices = input;
+  thrust::device_vector<unsigned int> indices(input.size());
+
+  // sort vertices to bring duplicates together
+  thrust::sort(vertices.begin(), vertices.end());
+
+  // find unique vertices and erase redundancies
+  vertices.erase(thrust::unique(vertices.begin(), vertices.end()), vertices.end());
+
+  // find index of each input vertex in the list of unique vertices
+  thrust::lower_bound(vertices.begin(), vertices.end(), input.begin(), input.end(), indices.begin());
+
+  // print output mesh representation
+  std::cout << "Output Representation" << '\n';
+  for (size_t i = 0; i < vertices.size(); i++)
+  {
+    vec2 v = vertices[i];
+    std::cout << " vertices[" << i << "] = (" << cuda::std::get<0>(v) << "," << cuda::std::get<1>(v) << ")" << '\n';
+  }
+  for (size_t i = 0; i < indices.size(); i++)
+  {
+    std::cout << " indices[" << i << "] = " << indices[i] << '\n';
+  }
+
+  return 0;
+}

cccl_upstream/thrust/examples/word_count.cmake

@@ -0,0 +1,21 @@
+target_compile_options(
+  ${example_target}
+  PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: --extended-lambda>
+)
+
+# This check is actually not correct, because we must check the host compiler, not the CXX compiler.
+# We rely on these usually being the same ;)
+if (
+  "Clang" STREQUAL "${CMAKE_CXX_COMPILER_ID}"
+  AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 13
+)
+  # When clang >= 13 is used as host compiler, we get the following warning:
+  #   nvcc_internal_extended_lambda_implementation:312:22: error: definition of implicit copy constructor for '__nv_hdl_wrapper_t<false, true, false, __nv_dl_tag<void (*)(), &TestAddressStabilityLambda, 2>, int (const int &)>' is deprecated because it has a user-declared copy assignment operator [-Werror,-Wdeprecated-copy]
+  #   312 | __nv_hdl_wrapper_t & operator=(const __nv_hdl_wrapper_t &in) = delete;
+  #       |                      ^
+  # Let's suppress it until NVBug 4980157 is resolved.
+  target_compile_options(
+    ${example_target}
+    PRIVATE $<$<COMPILE_LANG_AND_ID:CUDA,NVIDIA>: -Wno-deprecated-copy>
+  )
+endif()

cccl_upstream/thrust/examples/word_count.cu

@@ -0,0 +1,79 @@
+#include <thrust/device_vector.h>
+#include <thrust/functional.h>
+#include <thrust/host_vector.h>
+#include <thrust/inner_product.h>
+#include <thrust/reduce.h>
+
+#include <cuda/std/iterator> // Required for std::begin/std::end
+
+#include <iostream>
+
+// This example computes the number of words in a text sample
+// with a single call to thrust::inner_product.  The algorithm
+// counts the number of characters which start a new word, i.e.
+// the number of characters where input[i] is an alphabetical
+// character and input[i-1] is not an alphabetical character.
+
+// determines whether the character is alphabetical
+__host__ __device__ bool is_alpha(const char c)
+{
+  return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
+}
+
+int word_count(const thrust::device_vector<char>& input)
+{
+  // check for empty string
+  if (input.empty())
+  {
+    return 0;
+  }
+
+  // determines whether the right character begins a new word
+  auto is_word_start = [] __host__ __device__(const char& left, const char& right) {
+    return is_alpha(right) && !is_alpha(left);
+  };
+
+  // compute the number characters that start a new word
+  int wc = thrust::inner_product(
+    input.begin(),
+    input.end() - 1, // sequence of left characters
+    input.begin() + 1, // sequence of right characters
+    0, // initialize sum to 0
+    cuda::std::plus<int>{}, // sum values together
+    is_word_start // how to compare the left and right characters
+  );
+
+  // if the first character is alphabetical, then it also begins a word
+  if (is_alpha(input.front()))
+  {
+    wc++;
+  }
+
+  return wc;
+}
+
+int main()
+{
+  // Paragraph from 'The Raven' by Edgar Allan Poe
+  // http://en.wikipedia.org/wiki/The_Raven
+  const char raw_input[] =
+    "  But the raven, sitting lonely on the placid bust, spoke only,\n"
+    "  That one word, as if his soul in that one word he did outpour.\n"
+    "  Nothing further then he uttered - not a feather then he fluttered -\n"
+    "  Till I scarcely more than muttered `Other friends have flown before -\n"
+    "  On the morrow he will leave me, as my hopes have flown before.'\n"
+    "  Then the bird said, `Nevermore.'\n";
+
+  std::cout << "Text sample:\n";
+  std::cout << raw_input << "\n";
+
+  // transfer to device
+  thrust::device_vector<char> input(cuda::std::begin(raw_input), cuda::std::end(raw_input));
+
+  // count words
+  int wc = word_count(input);
+
+  std::cout << "Text sample contains " << wc << " words\n";
+
+  return 0;
+}