text: 
  /*** DO NOT EDIT THIS FILE DIRECTLY (use 'headers' to generate) ***/
  
  #ifndef hypre_UTILITIES_HPP
  #define hypre_UTILITIES_HPP
  
  #ifdef __cplusplus
  extern "C++" {
  #endif
  
  /******************************************************************************
   * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
   * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
   *
   * SPDX-License-Identifier: (Apache-2.0 OR MIT)
   ******************************************************************************/
  
  #ifndef HYPRE_UMPIRE_ALLOCATOR_H
  #define HYPRE_UMPIRE_ALLOCATOR_H
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
  #if defined(HYPRE_USING_UMPIRE_DEVICE)
  
  /*
  #include "umpire/Allocator.hpp"
  #include "umpire/ResourceManager.hpp"
  
  #include "umpire/strategy/DynamicPool.hpp"
  #include "umpire/strategy/AllocationAdvisor.hpp"
  #include "umpire/strategy/MonotonicAllocationStrategy.hpp"
  #include "umpire/util/Macros.hpp"
  */
  
  struct hypre_umpire_device_allocator
  {
     typedef char value_type;
  
     hypre_umpire_device_allocator()
     {
        // constructor
     }
  
     ~hypre_umpire_device_allocator()
     {
        // destructor
     }
  
     char *allocate(std::ptrdiff_t num_bytes)
     {
        char *ptr = NULL;
        hypre_umpire_device_pooled_allocate((void**) &ptr, num_bytes);
  
        return ptr;
     }
  
     void deallocate(char *ptr, size_t n)
     {
        hypre_umpire_device_pooled_free(ptr);
     }
  };
  
  #endif /* #ifdef HYPRE_USING_UMPIRE_DEVICE */
  #endif /* #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP) */
  
  #endif
  /******************************************************************************
   * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
   * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
   *
   * SPDX-License-Identifier: (Apache-2.0 OR MIT)
   ******************************************************************************/
  
  #ifndef HYPRE_CUDA_UTILS_H
  #define HYPRE_CUDA_UTILS_H
  
  #if defined(HYPRE_USING_GPU)
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  #include <cuda.h>
  #include <cuda_runtime.h>
  #include <cuda_profiler_api.h>
  #include <curand.h>
  #include <cublas_v2.h>
  #include <cusparse.h>
  
  #ifndef CUDART_VERSION
  #error CUDART_VERSION Undefined!
  #endif
  
  #ifndef CUDA_VERSION
  #error CUDA_VERSION Undefined!
  #endif
  
  #if CUDA_VERSION >= 11000
  #define THRUST_IGNORE_DEPRECATED_CPP11
  #define CUB_IGNORE_DEPRECATED_CPP11
  #define THRUST_IGNORE_DEPRECATED_CPP_DIALECT
  #define CUB_IGNORE_DEPRECATED_CPP_DIALECT
  #endif
  
  #define CUSPARSE_NEWAPI_VERSION 11000
  
  #elif defined(HYPRE_USING_HIP)
  
  #include <hip/hip_runtime.h>
  
  #endif // defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  
  #if defined(HYPRE_USING_ROCSPARSE)
  #include <rocsparse.h>
  #endif
  
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  #define HYPRE_CUDA_CALL(call) do {                                                           \
     cudaError_t err = call;                                                                   \
     if (cudaSuccess != err) {                                                                 \
        hypre_printf("CUDA ERROR (code = %d, %s) at %s:%d\n", err, cudaGetErrorString(err),    \
                     __FILE__, __LINE__);                                                      \
        hypre_assert(0); exit(1);                                                              \
     } } while(0)
  
  #elif defined(HYPRE_USING_HIP)
  #define HYPRE_HIP_CALL(call) do {                                                            \
     hipError_t err = call;                                                                    \
     if (hipSuccess != err) {                                                                  \
        hypre_printf("HIP ERROR (code = %d, %s) at %s:%d\n", err, hipGetErrorString(err),      \
                     __FILE__, __LINE__);                                                      \
        hypre_assert(0); exit(1);                                                              \
     } } while(0)
  
  #endif // defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  
  #define HYPRE_CUBLAS_CALL(call) do {                                                         \
     cublasStatus_t err = call;                                                                \
     if (CUBLAS_STATUS_SUCCESS != err) {                                                       \
        hypre_printf("CUBLAS ERROR (code = %d, %d) at %s:%d\n",                                \
              err, err == CUBLAS_STATUS_EXECUTION_FAILED, __FILE__, __LINE__);                 \
        hypre_assert(0); exit(1);                                                              \
     } } while(0)
  
  #define HYPRE_CUSPARSE_CALL(call) do {                                                       \
     cusparseStatus_t err = call;                                                              \
     if (CUSPARSE_STATUS_SUCCESS != err) {                                                     \
        hypre_printf("CUSPARSE ERROR (code = %d, %s) at %s:%d\n",                              \
              err, cusparseGetErrorString(err), __FILE__, __LINE__);                           \
        hypre_assert(0); exit(1);                                                              \
     } } while(0)
  
  #define HYPRE_ROCSPARSE_CALL(call) do {                                                      \
     rocsparse_status err = call;                                                              \
     if (rocsparse_status_success != err) {                                                    \
        hypre_printf("rocSPARSE ERROR (code = %d) at %s:%d\n",                                 \
              err, __FILE__, __LINE__);                                                        \
        assert(0); exit(1);                                                                    \
     } } while(0)
  
  
  #define HYPRE_CURAND_CALL(call) do {                                                         \
     curandStatus_t err = call;                                                                \
     if (CURAND_STATUS_SUCCESS != err) {                                                       \
        hypre_printf("CURAND ERROR (code = %d) at %s:%d\n", err, __FILE__, __LINE__);          \
        hypre_assert(0); exit(1);                                                              \
     } } while(0)
  
  struct hypre_cub_CachingDeviceAllocator;
  typedef struct hypre_cub_CachingDeviceAllocator hypre_cub_CachingDeviceAllocator;
  
  // HYPRE_WARP_BITSHIFT is just log2 of HYPRE_WARP_SIZE
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  #define HYPRE_WARP_SIZE       32
  #define HYPRE_WARP_BITSHIFT   5
  #elif defined(HYPRE_USING_HIP)
  #define HYPRE_WARP_SIZE       64
  #define HYPRE_WARP_BITSHIFT   6
  #endif
  
  #define HYPRE_WARP_FULL_MASK  0xFFFFFFFF
  #define HYPRE_MAX_NUM_WARPS   (64 * 64 * 32)
  #define HYPRE_FLT_LARGE       1e30
  #define HYPRE_1D_BLOCK_SIZE   512
  #define HYPRE_MAX_NUM_STREAMS 10
  
  struct hypre_CudaData
  {
  #if defined(HYPRE_USING_CURAND)
     curandGenerator_t                 curand_generator;
  #endif
  
  #if defined(HYPRE_USING_CUBLAS)
     cublasHandle_t                    cublas_handle;
  #endif
  
  #if defined(HYPRE_USING_CUSPARSE)
     cusparseHandle_t                  cusparse_handle;
  #endif
  
  #if defined(HYPRE_USING_ROCSPARSE)
     rocsparse_handle                  cusparse_handle;
  #endif
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
     cudaStream_t                      cuda_streams[HYPRE_MAX_NUM_STREAMS];
  #elif defined(HYPRE_USING_HIP)
     hipStream_t                       cuda_streams[HYPRE_MAX_NUM_STREAMS];
  #endif
  
  #ifdef HYPRE_USING_DEVICE_POOL
     hypre_uint                        cub_bin_growth;
     hypre_uint                        cub_min_bin;
     hypre_uint                        cub_max_bin;
     size_t                            cub_max_cached_bytes;
     hypre_cub_CachingDeviceAllocator *cub_dev_allocator;
     hypre_cub_CachingDeviceAllocator *cub_uvm_allocator;
  #endif
  #ifdef HYPRE_USING_UMPIRE_DEVICE
     hypre_umpire_device_allocator     umpire_device_allocator;
  #endif
     HYPRE_Int                         cuda_device;
     /* by default, hypre puts GPU computations in this stream
      * Do not be confused with the default (null) CUDA stream */
     HYPRE_Int                         cuda_compute_stream_num;
     /* work space for hypre's CUDA reducer */
     void                             *cuda_reduce_buffer;
     /* the device buffers needed to do MPI communication for struct comm */
     HYPRE_Complex*                    struct_comm_recv_buffer;
     HYPRE_Complex*                    struct_comm_send_buffer;
     HYPRE_Int                         struct_comm_recv_buffer_size;
     HYPRE_Int                         struct_comm_send_buffer_size;
     /* device spgemm options */
     HYPRE_Int                         spgemm_use_cusparse;
     HYPRE_Int                         spgemm_num_passes;
     HYPRE_Int                         spgemm_rownnz_estimate_method;
     HYPRE_Int                         spgemm_rownnz_estimate_nsamples;
     float                             spgemm_rownnz_estimate_mult_factor;
     char                              spgemm_hash_type;
     /* PMIS */
     HYPRE_Int                         use_gpu_rand;
  };
  
  #define hypre_CudaDataCubBinGrowth(data)                   ((data) -> cub_bin_growth)
  #define hypre_CudaDataCubMinBin(data)                      ((data) -> cub_min_bin)
  #define hypre_CudaDataCubMaxBin(data)                      ((data) -> cub_max_bin)
  #define hypre_CudaDataCubMaxCachedBytes(data)              ((data) -> cub_max_cached_bytes)
  #define hypre_CudaDataCubDevAllocator(data)                ((data) -> cub_dev_allocator)
  #define hypre_CudaDataCubUvmAllocator(data)                ((data) -> cub_uvm_allocator)
  #define hypre_CudaDataCudaDevice(data)                     ((data) -> cuda_device)
  #define hypre_CudaDataCudaComputeStreamNum(data)           ((data) -> cuda_compute_stream_num)
  #define hypre_CudaDataCudaReduceBuffer(data)               ((data) -> cuda_reduce_buffer)
  #define hypre_CudaDataStructCommRecvBuffer(data)           ((data) -> struct_comm_recv_buffer)
  #define hypre_CudaDataStructCommSendBuffer(data)           ((data) -> struct_comm_send_buffer)
  #define hypre_CudaDataStructCommRecvBufferSize(data)       ((data) -> struct_comm_recv_buffer_size)
  #define hypre_CudaDataStructCommSendBufferSize(data)       ((data) -> struct_comm_send_buffer_size)
  #define hypre_CudaDataSpgemmUseCusparse(data)              ((data) -> spgemm_use_cusparse)
  #define hypre_CudaDataSpgemmNumPasses(data)                ((data) -> spgemm_num_passes)
  #define hypre_CudaDataSpgemmRownnzEstimateMethod(data)     ((data) -> spgemm_rownnz_estimate_method)
  #define hypre_CudaDataSpgemmRownnzEstimateNsamples(data)   ((data) -> spgemm_rownnz_estimate_nsamples)
  #define hypre_CudaDataSpgemmRownnzEstimateMultFactor(data) ((data) -> spgemm_rownnz_estimate_mult_factor)
  #define hypre_CudaDataSpgemmHashType(data)                 ((data) -> spgemm_hash_type)
  #define hypre_CudaDataUmpireDeviceAllocator(data)          ((data) -> umpire_device_allocator)
  #define hypre_CudaDataUseGpuRand(data)                     ((data) -> use_gpu_rand)
  
  hypre_CudaData*     hypre_CudaDataCreate();
  void                hypre_CudaDataDestroy(hypre_CudaData* data);
  
  #if defined(HYPRE_USING_CURAND)
  curandGenerator_t   hypre_CudaDataCurandGenerator(hypre_CudaData *data);
  #endif
  
  #if defined(HYPRE_USING_CUBLAS)
  cublasHandle_t      hypre_CudaDataCublasHandle(hypre_CudaData *data);
  #endif
  
  #if defined(HYPRE_USING_CUSPARSE)
  cusparseHandle_t    hypre_CudaDataCusparseHandle(hypre_CudaData *data);
  #endif
  
  #if defined(HYPRE_USING_ROCSPARSE)
  rocsparse_handle    hypre_CudaDataCusparseHandle(hypre_CudaData *data);
  #endif
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_DEVICE_OPENMP)
  cudaStream_t        hypre_CudaDataCudaStream(hypre_CudaData *data, HYPRE_Int i);
  cudaStream_t        hypre_CudaDataCudaComputeStream(hypre_CudaData *data);
  #elif defined(HYPRE_USING_HIP)
  hipStream_t         hypre_CudaDataCudaStream(hypre_CudaData *data, HYPRE_Int i);
  hipStream_t         hypre_CudaDataCudaComputeStream(hypre_CudaData *data);
  #endif
  
  // Data structure and accessor routines for Cuda Sparse Triangular Matrices
  struct hypre_CsrsvData
  {
  #if defined(HYPRE_USING_CUSPARSE)
     csrsv2Info_t info_L;
     csrsv2Info_t info_U;
  #endif
     hypre_int    BufferSize;
     char        *Buffer;
  };
  
  #define hypre_CsrsvDataInfoL(data)      ((data) -> info_L)
  #define hypre_CsrsvDataInfoU(data)      ((data) -> info_U)
  #define hypre_CsrsvDataBufferSize(data) ((data) -> BufferSize)
  #define hypre_CsrsvDataBuffer(data)     ((data) -> Buffer)
  
  struct hypre_GpuMatData
  {
  #if defined(HYPRE_USING_CUSPARSE)
     cusparseMatDescr_t    mat_descr;
  #endif
  
  #if defined(HYPRE_USING_ROCSPARSE)
     rocsparse_mat_descr   mat_descr;
     rocsparse_mat_info    mat_info;
  #endif
  };
  
  #define hypre_GpuMatDataMatDecsr(data) ((data) -> mat_descr)
  #define hypre_GpuMatDataMatInfo(data)  ((data) -> mat_info)
  
  #endif //#if defined(HYPRE_USING_GPU)
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
  
  #include <thrust/execution_policy.h>
  #if defined(HYPRE_USING_CUDA)
  #include <thrust/system/cuda/execution_policy.h>
  #elif defined(HYPRE_USING_HIP)
  #include <thrust/system/hip/execution_policy.h>
  #endif
  #include <thrust/count.h>
  #include <thrust/device_ptr.h>
  #include <thrust/unique.h>
  #include <thrust/sort.h>
  #include <thrust/binary_search.h>
  #include <thrust/iterator/constant_iterator.h>
  #include <thrust/iterator/counting_iterator.h>
  #include <thrust/transform.h>
  #include <thrust/functional.h>
  #include <thrust/gather.h>
  #include <thrust/scan.h>
  #include <thrust/fill.h>
  #include <thrust/adjacent_difference.h>
  #include <thrust/inner_product.h>
  #include <thrust/logical.h>
  #include <thrust/replace.h>
  #include <thrust/sequence.h>
  
  using namespace thrust::placeholders;
  
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   * macro for launching CUDA kernels, CUDA, Thrust, Cusparse, Curand calls
   *                    NOTE: IN HYPRE'S DEFAULT STREAM
   * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   */
  
  #if defined(HYPRE_DEBUG)
  
  #if defined(HYPRE_USING_CUDA)
  #define HYPRE_CUDA_LAUNCH(kernel_name, gridsize, blocksize, ...)                                                     \
  {                                                                                                                    \
     if ( gridsize.x  == 0 || gridsize.y  == 0 || gridsize.z  == 0 ||                                                  \
          blocksize.x == 0 || blocksize.y == 0 || blocksize.z == 0 )                                                   \
     {                                                                                                                 \
        /* hypre_printf("Warning %s %d: Zero CUDA grid/block (%d %d %d) (%d %d %d)\n",                                 \
                   __FILE__, __LINE__,                                                                                 \
                   gridsize.x, gridsize.y, gridsize.z, blocksize.x, blocksize.y, blocksize.z); */                      \
     }                                                                                                                 \
     else                                                                                                              \
     {                                                                                                                 \
        (kernel_name) <<< (gridsize), (blocksize), 0, hypre_HandleCudaComputeStream(hypre_handle()) >>> (__VA_ARGS__); \
     }                                                                                                                 \
     hypre_SyncCudaComputeStream(hypre_handle());                                                                      \
     HYPRE_CUDA_CALL( cudaGetLastError() );                                                                            \
  }
  #elif defined(HYPRE_USING_HIP)
  #define HYPRE_CUDA_LAUNCH(kernel_name, gridsize, blocksize, ...)                                                     \
  {                                                                                                                    \
     if ( gridsize.x  == 0 || gridsize.y  == 0 || gridsize.z  == 0 ||                                                  \
          blocksize.x == 0 || blocksize.y == 0 || blocksize.z == 0 )                                                   \
     {                                                                                                                 \
        /* hypre_printf("Warning %s %d: Zero CUDA grid/block (%d %d %d) (%d %d %d)\n",                                 \
                   __FILE__, __LINE__,                                                                                 \
                   gridsize.x, gridsize.y, gridsize.z, blocksize.x, blocksize.y, blocksize.z); */                      \
     }                                                                                                                 \
     else                                                                                                              \
     {                                                                                                                 \
        (kernel_name) <<< (gridsize), (blocksize), 0, hypre_HandleCudaComputeStream(hypre_handle()) >>> (__VA_ARGS__); \
     }                                                                                                                 \
     hypre_SyncCudaComputeStream(hypre_handle());                                                                      \
     HYPRE_HIP_CALL( hipGetLastError() );                                                                            \
  }
  #endif //HYPRE_USING_CUDA
  
  #else // #if defined(HYPRE_DEBUG)
  
  #define HYPRE_CUDA_LAUNCH(kernel_name, gridsize, blocksize, ...)                                                     \
  {                                                                                                                    \
     if ( gridsize.x  == 0 || gridsize.y  == 0 || gridsize.z  == 0 ||                                                  \
          blocksize.x == 0 || blocksize.y == 0 || blocksize.z == 0 )                                                   \
     {                                                                                                                 \
        /* hypre_printf("Warning %s %d: Zero CUDA grid/block (%d %d %d) (%d %d %d)\n",                                 \
                   __FILE__, __LINE__,                                                                                 \
                   gridsize.x, gridsize.y, gridsize.z, blocksize.x, blocksize.y, blocksize.z); */                      \
     }                                                                                                                 \
     else                                                                                                              \
     {                                                                                                                 \
        (kernel_name) <<< (gridsize), (blocksize), 0, hypre_HandleCudaComputeStream(hypre_handle()) >>> (__VA_ARGS__); \
     }                                                                                                                 \
  }
  
  #endif // defined(HYPRE_DEBUG)
  
  /* RL: TODO Want macro HYPRE_THRUST_CALL to return value but I don't know how to do it right
   * The following one works OK for now */
  
  #ifdef HYPRE_USING_UMPIRE_DEVICE
  
  #if defined(HYPRE_USING_CUDA)
  #define HYPRE_THRUST_CALL(func_name, ...)                                                                            \
     thrust::func_name(thrust::cuda::par(hypre_HandleUmpireDeviceAllocator(hypre_handle())).on(hypre_HandleCudaComputeStream(hypre_handle())), __VA_ARGS__);
  #elif defined(HYPRE_USING_HIP)
  #define HYPRE_THRUST_CALL(func_name, ...)                                                                            \
     thrust::func_name(thrust::hip::par(hypre_HandleUmpireDeviceAllocator(hypre_handle())).on(hypre_HandleCudaComputeStream(hypre_handle())), __VA_ARGS__);
  #endif // HYPRE_USING_CUDA
  
  #elif HYPRE_USING_DEVICE_POOL
  #if defined(HYPRE_USING_CUDA)
  #define HYPRE_THRUST_CALL(func_name, ...)                                                                            \
     thrust::func_name(thrust::cuda::par(*(hypre_HandleCubDevAllocator(hypre_handle()))).on(hypre_HandleCudaComputeStream(hypre_handle())), __VA_ARGS__);
  #endif
  
  #else
  
  #if defined(HYPRE_USING_CUDA)
  #define HYPRE_THRUST_CALL(func_name, ...)                                                                            \
     thrust::func_name(thrust::cuda::par.on(hypre_HandleCudaComputeStream(hypre_handle())), __VA_ARGS__);
  #elif defined(HYPRE_USING_HIP)
  #define HYPRE_THRUST_CALL(func_name, ...)                                                                            \
     thrust::func_name(thrust::hip::par.on(hypre_HandleCudaComputeStream(hypre_handle())), __VA_ARGS__);
  #endif // HYPRE_USING_CUDA
  
  #endif // HYPRE_USING_UMPIRE_DEVICE
  
  
  /* return the number of threads in block */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_num_threads()
  {
     switch (dim)
     {
        case 1:
           return (blockDim.x);
        case 2:
           return (blockDim.x * blockDim.y);
        case 3:
           return (blockDim.x * blockDim.y * blockDim.z);
     }
  
     return -1;
  }
  
  /* return the flattened thread id in block */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_thread_id()
  {
     switch (dim)
     {
        case 1:
           return (threadIdx.x);
        case 2:
           return (threadIdx.y * blockDim.x + threadIdx.x);
        case 3:
           return (threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +
                   threadIdx.x);
     }
  
     return -1;
  }
  
  /* return the number of warps in block */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_num_warps()
  {
     return hypre_cuda_get_num_threads<dim>() >> HYPRE_WARP_BITSHIFT;
  }
  
  /* return the warp id in block */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_warp_id()
  {
     return hypre_cuda_get_thread_id<dim>() >> HYPRE_WARP_BITSHIFT;
  }
  
  /* return the thread lane id in warp */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_lane_id()
  {
     return hypre_cuda_get_thread_id<dim>() & (HYPRE_WARP_SIZE-1);
  }
  
  /* return the num of blocks in grid */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_num_blocks()
  {
     switch (dim)
     {
        case 1:
           return (gridDim.x);
        case 2:
           return (gridDim.x * gridDim.y);
        case 3:
           return (gridDim.x * gridDim.y * gridDim.z);
     }
  
     return -1;
  }
  
  /* return the flattened block id in grid */
  template <hypre_int dim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_block_id()
  {
     switch (dim)
     {
        case 1:
           return (blockIdx.x);
        case 2:
           return (blockIdx.y * gridDim.x + blockIdx.x);
        case 3:
           return (blockIdx.z * gridDim.x * gridDim.y + blockIdx.y * gridDim.x +
                   blockIdx.x);
     }
  
     return -1;
  }
  
  /* return the number of threads in grid */
  template <hypre_int bdim, hypre_int gdim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_grid_num_threads()
  {
     return hypre_cuda_get_num_blocks<gdim>() * hypre_cuda_get_num_threads<bdim>();
  }
  
  /* return the flattened thread id in grid */
  template <hypre_int bdim, hypre_int gdim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_grid_thread_id()
  {
     return hypre_cuda_get_block_id<gdim>() * hypre_cuda_get_num_threads<bdim>() +
            hypre_cuda_get_thread_id<bdim>();
  }
  
  /* return the number of warps in grid */
  template <hypre_int bdim, hypre_int gdim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_grid_num_warps()
  {
     return hypre_cuda_get_num_blocks<gdim>() * hypre_cuda_get_num_warps<bdim>();
  }
  
  /* return the flattened warp id in grid */
  template <hypre_int bdim, hypre_int gdim>
  static __device__ __forceinline__
  hypre_int hypre_cuda_get_grid_warp_id()
  {
     return hypre_cuda_get_block_id<gdim>() * hypre_cuda_get_num_warps<bdim>() +
            hypre_cuda_get_warp_id<bdim>();
  }
  
  #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
  static __device__ __forceinline__
  hypre_double atomicAdd(hypre_double* address, hypre_double val)
  {
      hypre_ulonglongint* address_as_ull = (hypre_ulonglongint*) address;
      hypre_ulonglongint old = *address_as_ull, assumed;
  
      do {
          assumed = old;
          old = atomicCAS(address_as_ull, assumed,
                          __double_as_longlong(val +
                                 __longlong_as_double(assumed)));
  
      // Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN)
      } while (assumed != old);
  
      return __longlong_as_double(old);
  }
  #endif
  
  // There are no *_sync functions in HIP
  #if defined(HYPRE_USING_HIP) || (CUDA_VERSION < 9000)
  
  template <typename T>
  static __device__ __forceinline__
  T __shfl_sync(unsigned mask, T val, hypre_int src_line, hypre_int width=HYPRE_WARP_SIZE)
  {
     return __shfl(val, src_line, width);
  }
  
  template <typename T>
  static __device__ __forceinline__
  T __shfl_down_sync(unsigned mask, T val, unsigned delta, hypre_int width=HYPRE_WARP_SIZE)
  {
     return __shfl_down(val, delta, width);
  }
  
  template <typename T>
  static __device__ __forceinline__
  T __shfl_xor_sync(unsigned mask, T val, unsigned lanemask, hypre_int width=HYPRE_WARP_SIZE)
  {
     return __shfl_xor(val, lanemask, width);
  }
  
  template <typename T>
  static __device__ __forceinline__
  T __shfl_up_sync(unsigned mask, T val, unsigned delta, hypre_int width=HYPRE_WARP_SIZE)
  {
     return __shfl_up(val, delta, width);
  }
  
  static __device__ __forceinline__
  void __syncwarp()
  {
  }
  
  #endif // #if defined(HYPRE_USING_HIP) || (CUDA_VERSION < 9000)
  
  
  // __any was technically deprecated in CUDA 7 so we don't bother
  // with this overload for CUDA, just for HIP.
  #if defined(HYPRE_USING_HIP)
  static __device__ __forceinline__
  hypre_int __any_sync(unsigned mask, hypre_int predicate)
  {
    return __any(predicate);
  }
  #endif
  
  
  template <typename T>
  static __device__ __forceinline__
  T read_only_load( const T *ptr )
  {
     return __ldg( ptr );
  }
  
  /* exclusive prefix scan */
  template <typename T>
  static __device__ __forceinline__
  T warp_prefix_sum(hypre_int lane_id, T in, T &all_sum)
  {
  #pragma unroll
     for (hypre_int d = 2; d <=HYPRE_WARP_SIZE; d <<= 1)
     {
        T t = __shfl_up_sync(HYPRE_WARP_FULL_MASK, in, d >> 1);
        if ( (lane_id & (d - 1)) == (d - 1) )
        {
           in += t;
        }
     }
  
     all_sum = __shfl_sync(HYPRE_WARP_FULL_MASK, in, HYPRE_WARP_SIZE-1);
  
     if (lane_id == HYPRE_WARP_SIZE-1)
     {
        in = 0;
     }
  
  #pragma unroll
     for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
     {
        T t = __shfl_xor_sync(HYPRE_WARP_FULL_MASK, in, d);
  
        if ( (lane_id & (d - 1)) == (d - 1))
        {
          if ( (lane_id & ((d << 1) - 1)) == ((d << 1) - 1) )
           {
              in += t;
           }
           else
           {
              in = t;
           }
        }
     }
     return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_reduce_sum(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in += __shfl_down_sync(HYPRE_WARP_FULL_MASK, in, d);
    }
    return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_allreduce_sum(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in += __shfl_xor_sync(HYPRE_WARP_FULL_MASK, in, d);
    }
    return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_reduce_max(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in = max(in, __shfl_down_sync(HYPRE_WARP_FULL_MASK, in, d));
    }
    return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_allreduce_max(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in = max(in, __shfl_xor_sync(HYPRE_WARP_FULL_MASK, in, d));
    }
    return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_reduce_min(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in = min(in, __shfl_down_sync(HYPRE_WARP_FULL_MASK, in, d));
    }
    return in;
  }
  
  template <typename T>
  static __device__ __forceinline__
  T warp_allreduce_min(T in)
  {
  #pragma unroll
    for (hypre_int d = HYPRE_WARP_SIZE/2; d > 0; d >>= 1)
    {
      in = min(in, __shfl_xor_sync(HYPRE_WARP_FULL_MASK, in, d));
    }
    return in;
  }
  
  static __device__ __forceinline__
  hypre_int next_power_of_2(hypre_int n)
  {
     if (n <= 0)
     {
        return 0;
     }
  
     /* if n is power of 2, return itself */
     if ( (n & (n - 1)) == 0 )
     {
        return n;
     }
  
     n |= (n >>  1);
     n |= (n >>  2);
     n |= (n >>  4);
     n |= (n >>  8);
     n |= (n >> 16);
     n ^= (n >>  1);
     n  = (n <<  1);
  
     return n;
  }
  
  template<typename T>
  struct absolute_value : public thrust::unary_function<T,T>
  {
    __host__ __device__ T operator()(const T &x) const
    {
      return x < T(0) ? -x : x;
    }
  };
  
  template<typename T1, typename T2>
  struct TupleComp2
  {
     typedef thrust::tuple<T1, T2> Tuple;
  
     __host__ __device__ bool operator()(const Tuple& t1, const Tuple& t2)
     {
        if (thrust::get<0>(t1) < thrust::get<0>(t2))
        {
           return true;
        }
        if (thrust::get<0>(t1) > thrust::get<0>(t2))
        {
           return false;
        }
        return hypre_abs(thrust::get<1>(t1)) > hypre_abs(thrust::get<1>(t2));
     }
  };
  
  template<typename T1, typename T2>
  struct TupleComp3
  {
     typedef thrust::tuple<T1, T2> Tuple;
  
     __host__ __device__ bool operator()(const Tuple& t1, const Tuple& t2)
     {
        if (thrust::get<0>(t1) < thrust::get<0>(t2))
        {
           return true;
        }
        if (thrust::get<0>(t1) > thrust::get<0>(t2))
        {
           return false;
        }
        if (thrust::get<0>(t2) == thrust::get<1>(t2))
        {
           return false;
        }
        return thrust::get<0>(t1) == thrust::get<1>(t1) || thrust::get<1>(t1) < thrust::get<1>(t2);
     }
  };
  
  template<typename T>
  struct is_negative : public thrust::unary_function<T,bool>
  {
     __host__ __device__ bool operator()(const T &x)
     {
        return (x < 0);
     }
  };
  
  template<typename T>
  struct is_positive : public thrust::unary_function<T,bool>
  {
     __host__ __device__ bool operator()(const T &x)
     {
        return (x > 0);
     }
  };
  
  template<typename T>
  struct is_nonnegative : public thrust::unary_function<T,bool>
  {
     __host__ __device__ bool operator()(const T &x)
     {
        return (x >= 0);
     }
  };
  
  template<typename T>
  struct in_range : public thrust::unary_function<T, bool>
  {
     T low, up;
  
     in_range(T low_, T up_) { low = low_; up = up_; }
  
     __host__ __device__ bool operator()(const T &x)
     {
        return (x >= low && x <= up);
     }
  };
  
  template<typename T>
  struct out_of_range : public thrust::unary_function<T,bool>
  {
     T low, up;
  
     out_of_range(T low_, T up_) { low = low_; up = up_; }
  
     __host__ __device__ bool operator()(const T &x)
     {
        return (x < low || x > up);
     }
  };
  
  template<typename T>
  struct less_than : public thrust::unary_function<T,bool>
  {
     T val;
  
     less_than(T val_) { val = val_; }
  
     __host__ __device__ bool operator()(const T &x)
     {
        return (x < val);
     }
  };
  
  template<typename T>
  struct equal : public thrust::unary_function<T,bool>
  {
     T val;
  
     equal(T val_) { val = val_; }
  
     __host__ __device__ bool operator()(const T &x)
     {
        return (x == val);
     }
  };
  
  /* hypre_cuda_utils.c */
  dim3 hypre_GetDefaultCUDABlockDimension();
  
  dim3 hypre_GetDefaultCUDAGridDimension( HYPRE_Int n, const char *granularity, dim3 bDim );
  
  template <typename T1, typename T2, typename T3> HYPRE_Int hypreDevice_StableSortByTupleKey(HYPRE_Int N, T1 *keys1, T2 *keys2, T3 *vals, HYPRE_Int opt);
  
  template <typename T1, typename T2, typename T3, typename T4> HYPRE_Int hypreDevice_StableSortTupleByTupleKey(HYPRE_Int N, T1 *keys1, T2 *keys2, T3 *vals1, T4 *vals2, HYPRE_Int opt);
  
  template <typename T1, typename T2, typename T3> HYPRE_Int hypreDevice_ReduceByTupleKey(HYPRE_Int N, T1 *keys1_in,  T2 *keys2_in,  T3 *vals_in, T1 *keys1_out, T2 *keys2_out, T3 *vals_out);
  
  template <typename T>
  HYPRE_Int hypreDevice_CsrRowPtrsToIndicesWithRowNum(HYPRE_Int nrows, HYPRE_Int nnz, HYPRE_Int *d_row_ptr, T *d_row_num, T *d_row_ind);
  
  template <typename T>
  HYPRE_Int hypreDevice_ScatterConstant(T *x, HYPRE_Int n, HYPRE_Int *map, T v);
  
  HYPRE_Int hypreDevice_GetRowNnz(HYPRE_Int nrows, HYPRE_Int *d_row_indices, HYPRE_Int *d_diag_ia, HYPRE_Int *d_offd_ia, HYPRE_Int *d_rownnz);
  
  HYPRE_Int hypreDevice_CopyParCSRRows(HYPRE_Int nrows, HYPRE_Int *d_row_indices, HYPRE_Int job, HYPRE_Int has_offd, HYPRE_Int first_col, HYPRE_Int *d_col_map_offd_A, HYPRE_Int *d_diag_i, HYPRE_Int *d_diag_j, HYPRE_Complex *d_diag_a, HYPRE_Int *d_offd_i, HYPRE_Int *d_offd_j, HYPRE_Complex *d_offd_a, HYPRE_Int *d_ib, HYPRE_BigInt *d_jb, HYPRE_Complex *d_ab);
  
  HYPRE_Int hypreDevice_IntegerReduceSum(HYPRE_Int m, HYPRE_Int *d_i);
  
  HYPRE_Int hypreDevice_IntegerInclusiveScan(HYPRE_Int n, HYPRE_Int *d_i);
  
  HYPRE_Int hypreDevice_IntegerExclusiveScan(HYPRE_Int n, HYPRE_Int *d_i);
  
  HYPRE_Int* hypreDevice_CsrRowPtrsToIndices(HYPRE_Int nrows, HYPRE_Int nnz, HYPRE_Int *d_row_ptr);
  
  HYPRE_Int hypreDevice_CsrRowPtrsToIndices_v2(HYPRE_Int nrows, HYPRE_Int nnz, HYPRE_Int *d_row_ptr, HYPRE_Int *d_row_ind);
  
  HYPRE_Int* hypreDevice_CsrRowIndicesToPtrs(HYPRE_Int nrows, HYPRE_Int nnz, HYPRE_Int *d_row_ind);
  
  HYPRE_Int hypreDevice_CsrRowIndicesToPtrs_v2(HYPRE_Int nrows, HYPRE_Int nnz, HYPRE_Int *d_row_ind, HYPRE_Int *d_row_ptr);
  
  HYPRE_Int hypreDevice_GenScatterAdd(HYPRE_Real *x, HYPRE_Int ny, HYPRE_Int *map, HYPRE_Real *y, char *work);
  
  HYPRE_Int hypreDevice_BigToSmallCopy(HYPRE_Int *tgt, const HYPRE_BigInt *src, HYPRE_Int size);
  
  void hypre_CudaCompileFlagCheck();
  
  #if defined(HYPRE_USING_CUDA)
  cudaError_t hypre_CachingMallocDevice(void **ptr, size_t nbytes);
  
  cudaError_t hypre_CachingMallocManaged(void **ptr, size_t nbytes);
  
  cudaError_t hypre_CachingFreeDevice(void *ptr);
  
  cudaError_t hypre_CachingFreeManaged(void *ptr);
  #endif
  
  hypre_cub_CachingDeviceAllocator * hypre_CudaDataCubCachingAllocatorCreate(hypre_uint bin_growth, hypre_uint min_bin, hypre_uint max_bin, size_t max_cached_bytes, bool skip_cleanup, bool debug, bool use_managed_memory);
  
  void hypre_CudaDataCubCachingAllocatorDestroy(hypre_CudaData *data);
  
  #endif // #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
  
  #if defined(HYPRE_USING_CUSPARSE)
  
  cudaDataType hypre_HYPREComplexToCudaDataType();
  
  cusparseIndexType_t hypre_HYPREIntToCusparseIndexType();
  
  #endif // #if defined(HYPRE_USING_CUSPARSE)
  
  #endif /* #ifndef HYPRE_CUDA_UTILS_H */
  /******************************************************************************
   * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
   * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
   *
   * SPDX-License-Identifier: (Apache-2.0 OR MIT)
   ******************************************************************************/
  
  /* CUDA reducer class */
  
  #ifndef HYPRE_CUDA_REDUCER_H
  #define HYPRE_CUDA_REDUCER_H
  
  #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
  #if !defined(HYPRE_USING_RAJA) && !defined(HYPRE_USING_KOKKOS)
  
  template<typename T> void OneBlockReduce(T *d_arr, HYPRE_Int N, T *h_out);
  
  struct HYPRE_double4
  {
     HYPRE_Real x,y,z,w;
  
     __host__ __device__
     HYPRE_double4() {}
  
     __host__ __device__
     HYPRE_double4(HYPRE_Real x1, HYPRE_Real x2, HYPRE_Real x3, HYPRE_Real x4)
     {
        x = x1;
        y = x2;
        z = x3;
        w = x4;
     }
  
     __host__ __device__
     void operator=(HYPRE_Real val)
     {
        x = y = z = w = val;
     }
  
     __host__ __device__
     void operator+=(HYPRE_double4 rhs)
     {
        x += rhs.x;
        y += rhs.y;
        z += rhs.z;
        w += rhs.w;
     }
  
  };
  
  struct HYPRE_double6
  {
     HYPRE_Real x,y,z,w,u,v;
  
     __host__ __device__
     HYPRE_double6() {}
  
     __host__ __device__
     HYPRE_double6(HYPRE_Real x1, HYPRE_Real x2, HYPRE_Real x3, HYPRE_Real x4,
                   HYPRE_Real x5, HYPRE_Real x6)
     {
        x = x1;
        y = x2;
        z = x3;
        w = x4;
        u = x5;
        v = x6;
     }
  
     __host__ __device__
     void operator=(HYPRE_Real val)
     {
        x = y = z = w = u = v = val;
     }
  
     __host__ __device__
     void operator+=(HYPRE_double6 rhs)
     {
        x += rhs.x;
        y += rhs.y;
        z += rhs.z;
        w += rhs.w;
        u += rhs.u;
        v += rhs.v;
     }
  
  };
  
  /* reduction within a warp */
  __inline__ __host__ __device__
  HYPRE_Real warpReduceSum(HYPRE_Real val)
  {
  #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
    for (HYPRE_Int offset = warpSize/2; offset > 0; offset /= 2)
    {
      val += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val, offset);
    }
  #endif
    return val;
  }
  
  __inline__ __host__ __device__
  HYPRE_double4 warpReduceSum(HYPRE_double4 val) {
  #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
    for (HYPRE_Int offset = warpSize / 2; offset > 0; offset /= 2)
    {
      val.x += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.x, offset);
      val.y += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.y, offset);
      val.z += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.z, offset);
      val.w += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.w, offset);
    }
  #endif
    return val;
  }
  
  __inline__ __host__ __device__
  HYPRE_double6 warpReduceSum(HYPRE_double6 val) {
  #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
    for (HYPRE_Int offset = warpSize / 2; offset > 0; offset /= 2)
    {
      val.x += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.x, offset);
      val.y += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.y, offset);
      val.z += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.z, offset);
      val.w += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.w, offset);
      val.u += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.u, offset);
      val.v += __shfl_down_sync(HYPRE_WARP_FULL_MASK, val.v, offset);
    }
  #endif
    return val;
  }
  
  /* reduction within a block */
  template <typename T>
  __inline__ __host__ __device__
  T blockReduceSum(T val)
  {
  #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
     //static __shared__ T shared[HYPRE_WARP_SIZE]; // Shared mem for HYPRE_WARP_SIZE partial sums
  
     __shared__ T shared[HYPRE_WARP_SIZE];        // Shared mem for HYPRE_WARP_SIZE partial sums
  
     //HYPRE_Int lane = threadIdx.x % warpSize;
     //HYPRE_Int wid  = threadIdx.x / warpSize;
     HYPRE_Int lane = threadIdx.x & (warpSize - 1);
     HYPRE_Int wid  = threadIdx.x >> HYPRE_WARP_BITSHIFT;
  
     val = warpReduceSum(val);       // Each warp performs partial reduction
  
     if (lane == 0)
     {
        shared[wid] = val;          // Write reduced value to shared memory
     }
  
     __syncthreads();               // Wait for all partial reductions
  
     //read from shared memory only if that warp existed
     if (threadIdx.x < blockDim.x / warpSize)
     {
        val = shared[lane];
     }
     else
     {
        val = 0.0;
     }
  
     if (wid == 0)
     {
        val = warpReduceSum(val); //Final reduce within first warp
     }
  
  #endif
     return val;
  }
  
  template<typename T>
  __global__ void
  OneBlockReduceKernel(T *arr, HYPRE_Int N)
  {
     T sum;
     sum = 0.0;
     if (threadIdx.x < N)
     {
        sum = arr[threadIdx.x];
     }
     sum = blockReduceSum(sum);
     if (threadIdx.x == 0)
     {
        arr[0] = sum;
     }
  }
  
  /* Reducer class */
  template <typename T>
  struct ReduceSum
  {
     using value_type = T;
  
     T init;                    /* initial value passed in */
     mutable T __thread_sum;    /* place to hold local sum of a thread,
                                   and partial sum of a block */
     T *d_buf;                  /* place to store partial sum within blocks
                                   in the 1st round, used in the 2nd round */
     HYPRE_Int nblocks;         /* number of blocks used in the 1st round */
  
     /* constructor
      * val is the initial value (added to the reduced sum) */
     __host__
     ReduceSum(T val)
     {
        init = val;
        __thread_sum = 0.0;
        nblocks = -1;
  
        if (hypre_HandleCudaReduceBuffer(hypre_handle()) == NULL)
        {
           /* allocate for the max size for reducing double6 type */
           hypre_HandleCudaReduceBuffer(hypre_handle()) = hypre_TAlloc(HYPRE_double6, 1024, HYPRE_MEMORY_DEVICE);
        }
  
        d_buf = (T*) hypre_HandleCudaReduceBuffer(hypre_handle());
     }
  
     /* copy constructor */
     __host__ __device__
     ReduceSum(const ReduceSum<T>& other)
     {
        *this = other;
     }
  
     /* reduction within blocks */
     __host__ __device__
     void BlockReduce() const
     {
  #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
        __thread_sum = blockReduceSum(__thread_sum);
        if (threadIdx.x == 0)
        {
           d_buf[blockIdx.x] = __thread_sum;
        }
  #endif
     }
  
     __host__ __device__
     void operator+=(T val) const
     {
        __thread_sum += val;
     }
  
     /* invoke the 2nd reduction at the time want the sum from the reducer */
     __host__
     operator T()
     {
        T val;
  
        HYPRE_ExecutionPolicy exec_policy = hypre_HandleStructExecPolicy(hypre_handle());
  
        if (exec_policy == HYPRE_EXEC_HOST)
        {
           val = __thread_sum;
           val += init;
        }
        else
        {
           /* 2nd reduction with only *one* block */
           hypre_assert(nblocks >= 0 && nblocks <= 1024);
           const dim3 gDim(1), bDim(1024);
           HYPRE_CUDA_LAUNCH( OneBlockReduceKernel, gDim, bDim, d_buf, nblocks );
           hypre_TMemcpy(&val, d_buf, T, 1, HYPRE_MEMORY_HOST, HYPRE_MEMORY_DEVICE);
           val += init;
        }
  
        return val;
     }
  
     /* destructor */
     __host__ __device__
     ~ReduceSum<T>()
     {
     }
  };
  
  #endif /* #if !defined(HYPRE_USING_RAJA) && !defined(HYPRE_USING_KOKKOS) */
  #endif /* #if defined(HYPRE_USING_CUDA)  || defined(HYPRE_USING_HIP) */
  #endif /* #ifndef HYPRE_CUDA_REDUCER_H */
  /******************************************************************************
   * Copyright (c) 2011, Duane Merrill.  All rights reserved.
   * Copyright (c) 2011-2018, NVIDIA CORPORATION.  All rights reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions are met:
   *     * Redistributions of source code must retain the above copyright
   *       notice, this list of conditions and the following disclaimer.
   *     * Redistributions in binary form must reproduce the above copyright
   *       notice, this list of conditions and the following disclaimer in the
   *       documentation and/or other materials provided with the distribution.
   *     * Neither the name of the NVIDIA CORPORATION nor the
   *       names of its contributors may be used to endorse or promote products
   *       derived from this software without specific prior written permission.
   *
   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
   * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED. IN NO EVENT SHALL NVIDIA CORPORATION BE LIABLE FOR ANY
   * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
   * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
   * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   *
   ******************************************************************************/
  
  /******************************************************************************
   * Simple caching allocator for device memory allocations. The allocator is
   * thread-safe and capable of managing device allocations on multiple devices.
   ******************************************************************************/
  
  #ifndef HYPRE_CUB_ALLOCATOR_HEADER
  #define HYPRE_CUB_ALLOCATOR_HEADER
  
  #if defined(HYPRE_USING_CUDA) && defined(HYPRE_USING_DEVICE_POOL)
  
  #include <set>
  #include <map>
  
  #if (__cplusplus > 199711L) || (defined(_MSC_VER) && _MSC_VER >= 1800)
      #include <mutex>
  #else
      #if defined(_WIN32) || defined(_WIN64)
          #include <intrin.h>
  
          #define WIN32_LEAN_AND_MEAN
          #define NOMINMAX
          #include <windows.h>
          #undef WIN32_LEAN_AND_MEAN
          #undef NOMINMAX
  
          /**
           * Compiler read/write barrier
           */
          #pragma intrinsic(_ReadWriteBarrier)
  
      #endif
  #endif
  
  /**
   * Simple portable mutex
   *   - Wraps std::mutex when compiled with C++11 or newer (supported on all platforms)
   *   - Uses GNU/Windows spinlock mechanisms for pre C++11 (supported on x86/x64 when compiled with cl.exe or g++)
   */
  struct hypre_cub_Mutex
  {
  #if (__cplusplus > 199711L) || (defined(_MSC_VER) && _MSC_VER >= 1800)
  
      std::mutex mtx;
  
      void Lock()
      {
          mtx.lock();
      }
  
      void Unlock()
      {
          mtx.unlock();
      }
  
      void TryLock()
      {
          mtx.try_lock();
      }
  
  #else       //__cplusplus > 199711L
  
      #if defined(_MSC_VER)
  
          // Microsoft VC++
          typedef hypre_longint Spinlock;
  
      #else
  
          // GNU g++
          typedef hypre_int Spinlock;
  
          /**
           * Compiler read/write barrier
           */
          __forceinline__ void _ReadWriteBarrier()
          {
              __sync_synchronize();
          }
  
          /**
           * Atomic exchange
           */
          __forceinline__ hypre_longint _InterlockedExchange(volatile hypre_int * const Target, const hypre_int Value)
          {
              // NOTE: __sync_lock_test_and_set would be an acquire barrier, so we force a full barrier
              _ReadWriteBarrier();
              return __sync_lock_test_and_set(Target, Value);
          }
  
          /**
           * Pause instruction to prevent excess processor bus usage
           */
          __forceinline__ void YieldProcessor()
          {
          }
  
      #endif  // defined(_MSC_VER)
  
          /// Lock member
          volatile Spinlock lock;
  
          /**
           * Constructor
           */
          hypre_cub_Mutex() : lock(0) {}
  
          /**
           * Return when the specified spinlock has been acquired
           */
          __forceinline__ void Lock()
          {
              while (1)
              {
                  if (!_InterlockedExchange(&lock, 1)) return;
                  while (lock) YieldProcessor();
              }
          }
  
  
          /**
           * Release the specified spinlock
           */
          __forceinline__ void Unlock()
          {
              _ReadWriteBarrier();
              lock = 0;
          }
  
  #endif      // __cplusplus > 199711L
  
  };
  
  #include <math.h>
  
  /******************************************************************************
   * CachingDeviceAllocator (host use)
   ******************************************************************************/
  
  /**
   * \brief A simple caching allocator for device memory allocations.
   *
   * \par Overview
   * The allocator is thread-safe and stream-safe and is capable of managing cached
   * device allocations on multiple devices.  It behaves as follows:
   *
   * \par
   * - Allocations from the allocator are associated with an \p active_stream.  Once freed,
   *   the allocation becomes available immediately for reuse within the \p active_stream
   *   with which it was associated with during allocation, and it becomes available for
   *   reuse within other streams when all prior work submitted to \p active_stream has completed.
   * - Allocations are categorized and cached by bin size.  A new allocation request of
   *   a given size will only consider cached allocations within the corresponding bin.
   * - Bin limits progress geometrically in accordance with the growth factor
   *   \p bin_growth provided during construction.  Unused device allocations within
   *   a larger bin cache are not reused for allocation requests that categorize to
   *   smaller bin sizes.
   * - Allocation requests below (\p bin_growth ^ \p min_bin) are rounded up to
   *   (\p bin_growth ^ \p min_bin).
   * - Allocations above (\p bin_growth ^ \p max_bin) are not rounded up to the nearest
   *   bin and are simply freed when they are deallocated instead of being returned
   *   to a bin-cache.
   * - %If the total storage of cached allocations on a given device will exceed
   *   \p max_cached_bytes, allocations for that device are simply freed when they are
   *   deallocated instead of being returned to their bin-cache.
   *
   * \par
   * For example, the default-constructed CachingDeviceAllocator is configured with:
   * - \p bin_growth          = 8
   * - \p min_bin             = 3
   * - \p max_bin             = 7
   * - \p max_cached_bytes    = 6MB - 1B
   *
   * \par
   * which delineates five bin-sizes: 512B, 4KB, 32KB, 256KB, and 2MB
   * and sets a maximum of 6,291,455 cached bytes per device
   *
   */
  struct hypre_cub_CachingDeviceAllocator
  {
     typedef char value_type;
  
     //---------------------------------------------------------------------
     // Constants
     //---------------------------------------------------------------------
  
     /// Out-of-bounds bin
     static const hypre_uint INVALID_BIN = (hypre_uint) -1;
  
     /// Invalid size
     static const size_t INVALID_SIZE = (size_t) -1;
  
     /// Invalid device ordinal
     static const hypre_int INVALID_DEVICE_ORDINAL = -1;
  
     //---------------------------------------------------------------------
     // Type definitions and helper types
     //---------------------------------------------------------------------
  
     /**
      * Descriptor for device memory allocations
      */
     struct BlockDescriptor
     {
        void*           d_ptr;              // Device pointer
        size_t          bytes;              // Size of allocation in bytes
        hypre_uint      bin;                // Bin enumeration
        hypre_int       device;             // device ordinal
        cudaStream_t    associated_stream;  // Associated associated_stream
        cudaEvent_t     ready_event;        // Signal when associated stream has run to the point at which this block was freed
  
        // Constructor (suitable for searching maps for a specific block, given its pointer and device)
        BlockDescriptor(void *d_ptr, hypre_int device) :
           d_ptr(d_ptr),
           bytes(0),
           bin(INVALID_BIN),
           device(device),
           associated_stream(0),
           ready_event(0)
        {}
  
        // Constructor (suitable for searching maps for a range of suitable blocks, given a device)
        BlockDescriptor(hypre_int device) :
           d_ptr(NULL),
           bytes(0),
           bin(INVALID_BIN),
           device(device),
           associated_stream(0),
           ready_event(0)
        {}
  
        // Comparison functor for comparing device pointers
        static bool PtrCompare(const BlockDescriptor &a, const BlockDescriptor &b)
        {
           if (a.device == b.device)
              return (a.d_ptr < b.d_ptr);
           else
              return (a.device < b.device);
        }
  
        // Comparison functor for comparing allocation sizes
        static bool SizeCompare(const BlockDescriptor &a, const BlockDescriptor &b)
        {
           if (a.device == b.device)
              return (a.bytes < b.bytes);
           else
              return (a.device < b.device);
        }
     };
  
     /// BlockDescriptor comparator function interface
     typedef bool (*Compare)(const BlockDescriptor &, const BlockDescriptor &);
  
     class TotalBytes {
        public:
           size_t free;
           size_t live;
           TotalBytes() { free = live = 0; }
     };
  
     /// Set type for cached blocks (ordered by size)
     typedef std::multiset<BlockDescriptor, Compare> CachedBlocks;
  
     /// Set type for live blocks (ordered by ptr)
     typedef std::multiset<BlockDescriptor, Compare> BusyBlocks;
  
     /// Map type of device ordinals to the number of cached bytes cached by each device
     typedef std::map<hypre_int, TotalBytes> GpuCachedBytes;
  
  
     //---------------------------------------------------------------------
     // Utility functions
     //---------------------------------------------------------------------
  
     /**
      * Integer pow function for unsigned base and exponent
      */
     static hypre_uint IntPow(
           hypre_uint base,
           hypre_uint exp)
     {
        hypre_uint retval = 1;
        while (exp > 0)
        {
           if (exp & 1) {
              retval = retval * base;        // multiply the result by the current base
           }
           base = base * base;                // square the base
           exp = exp >> 1;                    // divide the exponent in half
        }
        return retval;
     }
  
  
     /**
      * Round up to the nearest power-of
      */
     void NearestPowerOf(
           hypre_uint      &power,
           size_t          &rounded_bytes,
           hypre_uint       base,
           size_t           value)
     {
        power = 0;
        rounded_bytes = 1;
  
        if (value * base < value)
        {
           // Overflow
           power = sizeof(size_t) * 8;
           rounded_bytes = size_t(0) - 1;
           return;
        }
  
        while (rounded_bytes < value)
        {
           rounded_bytes *= base;
           power++;
        }
     }
  
  
     //---------------------------------------------------------------------
     // Fields
     //---------------------------------------------------------------------
  
     hypre_cub_Mutex mutex;              /// Mutex for thread-safety
  
     hypre_uint      bin_growth;         /// Geometric growth factor for bin-sizes
     hypre_uint      min_bin;            /// Minimum bin enumeration
     hypre_uint      max_bin;            /// Maximum bin enumeration
  
     size_t          min_bin_bytes;      /// Minimum bin size
     size_t          max_bin_bytes;      /// Maximum bin size
     size_t          max_cached_bytes;   /// Maximum aggregate cached bytes per device
  
     const bool      skip_cleanup;       /// Whether or not to skip a call to FreeAllCached() when destructor is called.  (The CUDA runtime may have already shut down for statically declared allocators)
     bool            debug;              /// Whether or not to print (de)allocation events to stdout
  
     GpuCachedBytes  cached_bytes;       /// Map of device ordinal to aggregate cached bytes on that device
     CachedBlocks    cached_blocks;      /// Set of cached device allocations available for reuse
     BusyBlocks      live_blocks;        /// Set of live device allocations currently in use
  
     bool            use_managed_memory; /// Whether to use managed memory or device memory
  
     //---------------------------------------------------------------------
     // Methods
     //---------------------------------------------------------------------
  
     /**
      * \brief Constructor.
      */
     hypre_cub_CachingDeviceAllocator(
           hypre_uint      bin_growth,                             ///< Geometric growth factor for bin-sizes
           hypre_uint      min_bin             = 1,                ///< Minimum bin (default is bin_growth ^ 1)
           hypre_uint      max_bin             = INVALID_BIN,      ///< Maximum bin (default is no max bin)
           size_t          max_cached_bytes    = INVALID_SIZE,     ///< Maximum aggregate cached bytes per device (default is no limit)
           bool            skip_cleanup        = false,            ///< Whether or not to skip a call to \p FreeAllCached() when the destructor is called (default is to deallocate)
           bool            debug               = false,            ///< Whether or not to print (de)allocation events to stdout (default is no stderr output)
           bool            use_managed_memory  = false)            ///< Whether to use managed memory or device memory
        :
           bin_growth(bin_growth),
           min_bin(min_bin),
           max_bin(max_bin),
           min_bin_bytes(IntPow(bin_growth, min_bin)),
           max_bin_bytes(IntPow(bin_growth, max_bin)),
           max_cached_bytes(max_cached_bytes),
           skip_cleanup(skip_cleanup),
           debug(debug),
           use_managed_memory(use_managed_memory),
           cached_blocks(BlockDescriptor::SizeCompare),
           live_blocks(BlockDescriptor::PtrCompare)
     {}
  
  
     /**
      * \brief Default constructor.
      *
      * Configured with:
      * \par
      * - \p bin_growth          = 8
      * - \p min_bin             = 3
      * - \p max_bin             = 7
      * - \p max_cached_bytes    = (\p bin_growth ^ \p max_bin) * 3) - 1 = 6,291,455 bytes
      *
      * which delineates five bin-sizes: 512B, 4KB, 32KB, 256KB, and 2MB and
      * sets a maximum of 6,291,455 cached bytes per device
      */
     hypre_cub_CachingDeviceAllocator(
           bool skip_cleanup = false,
           bool debug = false,
           bool use_managed_memory = false)
        :
           bin_growth(8),
           min_bin(3),
           max_bin(7),
           min_bin_bytes(IntPow(bin_growth, min_bin)),
           max_bin_bytes(IntPow(bin_growth, max_bin)),
           max_cached_bytes((max_bin_bytes * 3) - 1),
           skip_cleanup(skip_cleanup),
           debug(debug),
           use_managed_memory(use_managed_memory),
           cached_blocks(BlockDescriptor::SizeCompare),
           live_blocks(BlockDescriptor::PtrCompare)
     {}
  
  
     /**
      * \brief Sets the limit on the number bytes this allocator is allowed to cache per device.
      *
      * Changing the ceiling of cached bytes does not cause any allocations (in-use or
      * cached-in-reserve) to be freed.  See \p FreeAllCached().
      */
     cudaError_t SetMaxCachedBytes(
           size_t max_cached_bytes)
     {
        // Lock
        mutex.Lock();
  
        if (debug) printf("Changing max_cached_bytes (%zu -> %zu)\n", this->max_cached_bytes, max_cached_bytes);
  
        this->max_cached_bytes = max_cached_bytes;
  
        // Unlock
        mutex.Unlock();
  
        return cudaSuccess;
     }
  
  
     /**
      * \brief Provides a suitable allocation of device memory for the given size on the specified device.
      *
      * Once freed, the allocation becomes available immediately for reuse within the \p active_stream
      * with which it was associated with during allocation, and it becomes available for reuse within other
      * streams when all prior work submitted to \p active_stream has completed.
      */
     cudaError_t DeviceAllocate(
           hypre_int       device,             ///< [in] Device on which to place the allocation
           void            **d_ptr,            ///< [out] Reference to pointer to the allocation
           size_t          bytes,              ///< [in] Minimum number of bytes for the allocation
           cudaStream_t    active_stream = 0)  ///< [in] The stream to be associated with this allocation
     {
        *d_ptr                          = NULL;
        hypre_int entrypoint_device     = INVALID_DEVICE_ORDINAL;
        cudaError_t error               = cudaSuccess;
  
        if (device == INVALID_DEVICE_ORDINAL)
        {
           if ((error = cudaGetDevice(&entrypoint_device))) return error;
           device = entrypoint_device;
        }
  
        // Create a block descriptor for the requested allocation
        bool found = false;
        BlockDescriptor search_key(device);
        search_key.associated_stream = active_stream;
        NearestPowerOf(search_key.bin, search_key.bytes, bin_growth, bytes);
  
        if (search_key.bin > max_bin)
        {
           // Bin is greater than our maximum bin: allocate the request
           // exactly and give out-of-bounds bin.  It will not be cached
           // for reuse when returned.
           search_key.bin      = INVALID_BIN;
           search_key.bytes    = bytes;
        }
        else
        {
           // Search for a suitable cached allocation: lock
           mutex.Lock();
  
           if (search_key.bin < min_bin)
           {
              // Bin is less than minimum bin: round up
              search_key.bin      = min_bin;
              search_key.bytes    = min_bin_bytes;
           }
  
           // Iterate through the range of cached blocks on the same device in the same bin
           CachedBlocks::iterator block_itr = cached_blocks.lower_bound(search_key);
           while ((block_itr != cached_blocks.end())
                 && (block_itr->device == device)
                 && (block_itr->bin == search_key.bin))
           {
              // To prevent races with reusing blocks returned by the host but still
              // in use by the device, only consider cached blocks that are
              // either (from the active stream) or (from an idle stream)
              if ((active_stream == block_itr->associated_stream) ||
                    (cudaEventQuery(block_itr->ready_event) != cudaErrorNotReady))
              {
                 // Reuse existing cache block.  Insert into live blocks.
                 found = true;
                 search_key = *block_itr;
                 search_key.associated_stream = active_stream;
                 live_blocks.insert(search_key);
  
                 // Remove from free blocks
                 cached_bytes[device].free -= search_key.bytes;
                 cached_bytes[device].live += search_key.bytes;
  
                 if (debug) printf("\tDevice %d reused cached block at %p (%zu bytes) for stream %p (previously associated with stream %p).\n",
                       device, search_key.d_ptr, search_key.bytes, search_key.associated_stream, block_itr->associated_stream);
  
                 cached_blocks.erase(block_itr);
  
                 break;
              }
              block_itr++;
           }
  
           // Done searching: unlock
           mutex.Unlock();
        }
  
        // Allocate the block if necessary
        if (!found)
        {
           // Set runtime's current device to specified device (entrypoint may not be set)
           if (device != entrypoint_device)
           {
              if ((error = cudaGetDevice(&entrypoint_device))) return error;
              if ((error = cudaSetDevice(device))) return error;
           }
  
           // Attempt to allocate
  
           if (use_managed_memory)
           {
              error = cudaMallocManaged(&search_key.d_ptr, search_key.bytes);
           }
           else
           {
              error = cudaMalloc(&search_key.d_ptr, search_key.bytes);
           }
           if ((error) == cudaErrorMemoryAllocation)
           {
              // The allocation attempt failed: free all cached blocks on device and retry
              if (debug) printf("\tDevice %d failed to allocate %zu bytes for stream %p, retrying after freeing cached allocations",
                    device, search_key.bytes, search_key.associated_stream);
  
              error = cudaSuccess;    // Reset the error we will return
              cudaGetLastError();     // Reset CUDART's error
  
              // Lock
              mutex.Lock();
  
              // Iterate the range of free blocks on the same device
              BlockDescriptor free_key(device);
              CachedBlocks::iterator block_itr = cached_blocks.lower_bound(free_key);
  
              while ((block_itr != cached_blocks.end()) && (block_itr->device == device))
              {
                 // No need to worry about synchronization with the device: cudaFree is
                 // blocking and will synchronize across all kernels executing
                 // on the current device
  
                 // Free device memory and destroy stream event.
                 if ((error = cudaFree(block_itr->d_ptr))) break;
                 if ((error = cudaEventDestroy(block_itr->ready_event))) break;
  
                 // Reduce balance and erase entry
                 cached_bytes[device].free -= block_itr->bytes;
  
                 if (debug) printf("\tDevice %d freed %zu bytes.\n\t\t  %zu available blocks cached (%zu bytes), %zu live blocks (%zu bytes) outstanding.\n",
                       device, block_itr->bytes, cached_blocks.size(), cached_bytes[device].free, live_blocks.size(), cached_bytes[device].live);
  
                 cached_blocks.erase(block_itr);
  
                 block_itr++;
              }
  
              // Unlock
              mutex.Unlock();
  
              // Return under error
              if (error) return error;
  
              // Try to allocate again
  
              if (use_managed_memory)
              {
                 error = cudaMallocManaged(&search_key.d_ptr, search_key.bytes);
              }
              else
              {
                 error = cudaMalloc(&search_key.d_ptr, search_key.bytes);
              }
              if ((error)) return error;
           }
  
           // Create ready event
           if ((error = cudaEventCreateWithFlags(&search_key.ready_event, cudaEventDisableTiming)))
              return error;
  
           // Insert into live blocks
           mutex.Lock();
           live_blocks.insert(search_key);
           cached_bytes[device].live += search_key.bytes;
           mutex.Unlock();
  
           if (debug) printf("\tDevice %d allocated new device block at %p (%zu bytes associated with stream %p).\n",
                 device, search_key.d_ptr, search_key.bytes, search_key.associated_stream);
  
           // Attempt to revert back to previous device if necessary
           if ((entrypoint_device != INVALID_DEVICE_ORDINAL) && (entrypoint_device != device))
           {
              if ((error = cudaSetDevice(entrypoint_device))) return error;
           }
        }
  
        // Copy device pointer to output parameter
        *d_ptr = search_key.d_ptr;
  
        if (debug) printf("\t\t%zu available blocks cached (%zu bytes), %zu live blocks outstanding(%zu bytes).\n",
              cached_blocks.size(), cached_bytes[device].free, live_blocks.size(), cached_bytes[device].live);
  
        return error;
     }
  
  
     /**
      * \brief Provides a suitable allocation of device memory for the given size on the current device.
      *
      * Once freed, the allocation becomes available immediately for reuse within the \p active_stream
      * with which it was associated with during allocation, and it becomes available for reuse within other
      * streams when all prior work submitted to \p active_stream has completed.
      */
     cudaError_t DeviceAllocate(
           void            **d_ptr,            ///< [out] Reference to pointer to the allocation
           size_t          bytes,              ///< [in] Minimum number of bytes for the allocation
           cudaStream_t    active_stream = 0)  ///< [in] The stream to be associated with this allocation
     {
        return DeviceAllocate(INVALID_DEVICE_ORDINAL, d_ptr, bytes, active_stream);
     }
  
     char * allocate(size_t bytes)
     {
        char *ptr;
        DeviceAllocate((void **) &ptr, bytes);
  
        return ptr;
     }
  
  
     /**
      * \brief Frees a live allocation of device memory on the specified device, returning it to the allocator.
      *
      * Once freed, the allocation becomes available immediately for reuse within the \p active_stream
      * with which it was associated with during allocation, and it becomes available for reuse within other
      * streams when all prior work submitted to \p active_stream has completed.
      */
     cudaError_t DeviceFree(
           hypre_int       device,
           void*           d_ptr)
     {
        hypre_int entrypoint_device     = INVALID_DEVICE_ORDINAL;
        cudaError_t error               = cudaSuccess;
  
        if (device == INVALID_DEVICE_ORDINAL)
        {
           if ((error = cudaGetDevice(&entrypoint_device)))
              return error;
           device = entrypoint_device;
        }
  
        // Lock
        mutex.Lock();
  
        // Find corresponding block descriptor
        bool recached = false;
        BlockDescriptor search_key(d_ptr, device);
        BusyBlocks::iterator block_itr = live_blocks.find(search_key);
        if (block_itr != live_blocks.end())
        {
           // Remove from live blocks
           search_key = *block_itr;
           live_blocks.erase(block_itr);
           cached_bytes[device].live -= search_key.bytes;
  
           // Keep the returned allocation if bin is valid and we won't exceed the max cached threshold
           if ((search_key.bin != INVALID_BIN) && (cached_bytes[device].free + search_key.bytes <= max_cached_bytes))
           {
              // Insert returned allocation into free blocks
              recached = true;
              cached_blocks.insert(search_key);
              cached_bytes[device].free += search_key.bytes;
  
              if (debug) printf("\tDevice %d returned %zu bytes from associated stream %p.\n\t\t %zu available blocks cached (%zu bytes), %zu live blocks outstanding. (%zu bytes)\n",
                    device, search_key.bytes, search_key.associated_stream, cached_blocks.size(),
                    cached_bytes[device].free, live_blocks.size(), cached_bytes[device].live);
           }
        }
  
        // Unlock
        mutex.Unlock();
  
        // First set to specified device (entrypoint may not be set)
        if (device != entrypoint_device)
        {
           if ((error = cudaGetDevice(&entrypoint_device))) return error;
           if ((error = cudaSetDevice(device))) return error;
        }
  
        if (recached)
        {
           // Insert the ready event in the associated stream (must have current device set properly)
           if ((error = cudaEventRecord(search_key.ready_event, search_key.associated_stream))) return error;
        }
        else
        {
           // Free the allocation from the runtime and cleanup the event.
           if ((error = cudaFree(d_ptr))) return error;
           if ((error = cudaEventDestroy(search_key.ready_event))) return error;
  
           if (debug) printf("\tDevice %d freed %zu bytes from associated stream %p.\n\t\t  %zu available blocks cached (%zu bytes), %zu live blocks (%zu bytes) outstanding.\n",
                 device, search_key.bytes, search_key.associated_stream, cached_blocks.size(), cached_bytes[device].free, live_blocks.size(), cached_bytes[device].live);
        }
  
        // Reset device
        if ((entrypoint_device != INVALID_DEVICE_ORDINAL) && (entrypoint_device != device))
        {
           if ((error = cudaSetDevice(entrypoint_device))) return error;
        }
  
        return error;
     }
  
     /**
      * \brief Frees a live allocation of device memory on the current device, returning it to the allocator.
      *
      * Once freed, the allocation becomes available immediately for reuse within the \p active_stream
      * with which it was associated with during allocation, and it becomes available for reuse within other
      * streams when all prior work submitted to \p active_stream has completed.
      */
     cudaError_t DeviceFree(
           void*           d_ptr)
     {
        return DeviceFree(INVALID_DEVICE_ORDINAL, d_ptr);
     }
  
     void deallocate(char *ptr, size_t)
     {
        DeviceFree((void *) ptr);
     }
  
     /**
      * \brief Frees all cached device allocations on all devices
      */
     cudaError_t FreeAllCached()
     {
        cudaError_t error           = cudaSuccess;
        hypre_int entrypoint_device = INVALID_DEVICE_ORDINAL;
        hypre_int current_device    = INVALID_DEVICE_ORDINAL;
  
        mutex.Lock();
  
        while (!cached_blocks.empty())
        {
           // Get first block
           CachedBlocks::iterator begin = cached_blocks.begin();
  
           // Get entry-point device ordinal if necessary
           if (entrypoint_device == INVALID_DEVICE_ORDINAL)
           {
              if ((error = cudaGetDevice(&entrypoint_device))) break;
           }
  
           // Set current device ordinal if necessary
           if (begin->device != current_device)
           {
              if ((error = cudaSetDevice(begin->device))) break;
              current_device = begin->device;
           }
  
           // Free device memory
           if ((error = cudaFree(begin->d_ptr))) break;
           if ((error = cudaEventDestroy(begin->ready_event))) break;
  
           // Reduce balance and erase entry
           cached_bytes[current_device].free -= begin->bytes;
  
           if (debug) printf("\tDevice %d freed %zu bytes.\n\t\t  %zu available blocks cached (%zu bytes), %zu live blocks (%zu bytes) outstanding.\n",
                 current_device, begin->bytes, cached_blocks.size(), cached_bytes[current_device].free, live_blocks.size(), cached_bytes[current_device].live);
  
           cached_blocks.erase(begin);
        }
  
        mutex.Unlock();
  
        // Attempt to revert back to entry-point device if necessary
        if (entrypoint_device != INVALID_DEVICE_ORDINAL)
        {
           if ((error = cudaSetDevice(entrypoint_device))) return error;
        }
  
        return error;
     }
  
  
     /**
      * \brief Destructor
      */
     virtual ~hypre_cub_CachingDeviceAllocator()
     {
        if (!skip_cleanup)
           FreeAllCached();
     }
  };
  
  #endif // #if defined(HYPRE_USING_CUDA) && defined(HYPRE_USING_DEVICE_POOL)
  #endif // #ifndef HYPRE_CUB_ALLOCATOR_HEADER
  
  
  #ifdef __cplusplus
  }
  #endif
  
  #endif
  
  