accel.F90

Go to the documentation of this file.

!! Copyright (C) 2010-2016 X. Andrade
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
#if defined(HAVE_CUDA)
#define HAVE_ACCEL 1
#endif
 
module accel_oct_m
  use alloc_cache_oct_m
  use cuda_oct_m
  use debug_oct_m
  use global_oct_m
  use iso_c_binding, only: c_null_ptr, c_size_t
  use, intrinsic :: iso_fortran_env
  use loct_oct_m
  use math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use types_oct_m
  use parser_oct_m
  use profiling_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use string_oct_m
 
  implicit none
 
  private
 
  public ::                              &
    accel_context_t,                     &
    accel_device_t,                      &
    accel_mem_t,                         &
    accel_kernel_t,                      &
    accel_t,                             &
    accel_is_enabled,                    &
    accel_allow_cpu_only,                &
    accel_init,                          &
    accel_end,                           &
    accel_padded_size,                   &
    accel_kernel_start_call,             &
    accel_kernel_build,                  &
    accel_create_buffer,                 &
    accel_write_buffer,                  &
    accel_read_buffer,                   &
    accel_free_buffer,                   &
    accel_move_buffer,                   &
    accel_detach_buffer,                 &
    accel_ensure_buffer_size,            &
    accel_buffer_is_allocated,           &
    accel_grid_size,                     &
    accel_grid_size_extend_dim,          &
    accel_finish,                        &
    accel_set_kernel_arg,                &
    accel_max_block_size,                &
    accel_kernel_block_size,             &
    accel_kernel_run,                    &
    accel_set_buffer_to_zero,            &
    accel_copy_buffer,                   &
    accel_shared_memory_size,            &
    accel_global_memory_size,            &
    accel_max_size_per_dim,              &
    accel_get_device_pointer,            &
    daccel_get_pointer_with_offset,      &
    zaccel_get_pointer_with_offset,      &
    accel_clean_pointer,                 &
    accel_set_stream,                    &
    accel_get_stream,                    &
    accel_synchronize_all_streams,       &
    accel_get_unfolded_size
 
  integer, public, parameter ::                 &
    ACCEL_MEM_READ_ONLY  = 0,                   &
    accel_mem_read_write = 1,                   &
    accel_mem_write_only = 2
 
  type accel_context_t
    ! Components are public by default
#if defined(HAVE_CUDA)
    type(c_ptr)      :: cuda_context
#else
    integer          :: dummy
#endif
  end type accel_context_t
 
  type accel_device_t
    ! Components are public by default
#if defined(HAVE_CUDA)
    type(c_ptr)      :: cuda_device
#else
    integer         :: dummy
#endif
  end type accel_device_t
 
  type accel_t
    ! Components are public by default
    type(accel_context_t)  :: context
    type(accel_device_t)   :: device
    type(c_ptr)            :: cuda_stream
    type(c_ptr)            :: program_map
    integer                :: max_block_size
    integer(int64)         :: shared_memory_size
    integer(int64)         :: global_memory_size
    logical                :: enabled
    logical                :: allow_CPU_only
    logical                :: cuda_mpi
    integer                :: warp_size
    integer(int64)         :: initialize_buffers
    character(len=32)      :: debug_flag
    integer(int64)         :: max_block_dim(3)
    integer(int64)         :: max_grid_dim(3)
  end type accel_t
 
  type accel_mem_t
    ! Components are public by default
    type(c_ptr)            :: mem
    integer(c_size_t)      :: size = 0
    type(type_t)           :: type
    integer                :: flags = 0
    logical                :: allocated = .false.
  end type accel_mem_t
 
  type accel_kernel_t
    ! Components are public by default
#ifdef HAVE_CUDA
    type(c_ptr)                   :: cuda_kernel
    type(c_ptr)                   :: arguments
#endif
    logical                       :: initialized = .false.
    type(accel_kernel_t), pointer :: next
    integer                       :: arg_count
    character(len=128)            :: kernel_name
  end type accel_kernel_t
 
  type(accel_t), public :: accel
 
  ! the kernels
  type(accel_kernel_t), public, target, save :: kernel_vpsi
  type(accel_kernel_t), public, target, save :: kernel_vpsi_complex
  type(accel_kernel_t), public, target, save :: kernel_vpsi_spinors
  type(accel_kernel_t), public, target, save :: kernel_vpsi_spinors_complex
  type(accel_kernel_t), public, target, save :: kernel_daxpy
  type(accel_kernel_t), public, target, save :: kernel_zaxpy
  type(accel_kernel_t), public, target, save :: kernel_saxpy
  type(accel_kernel_t), public, target, save :: kernel_caxpy
  type(accel_kernel_t), public, target, save :: kernel_copy
  type(accel_kernel_t), public, target, save :: kernel_scopy
  type(accel_kernel_t), public, target, save :: kernel_copy_complex_to_real
  type(accel_kernel_t), public, target, save :: kernel_copy_real_to_complex
  type(accel_kernel_t), public, target, save :: dpack
  type(accel_kernel_t), public, target, save :: zpack
  type(accel_kernel_t), public, target, save :: spack
  type(accel_kernel_t), public, target, save :: cpack
  type(accel_kernel_t), public, target, save :: dunpack
  type(accel_kernel_t), public, target, save :: zunpack
  type(accel_kernel_t), public, target, save :: sunpack
  type(accel_kernel_t), public, target, save :: cunpack
  type(accel_kernel_t), public, target, save :: kernel_ghost_reorder
  type(accel_kernel_t), public, target, save :: kernel_density_real
  type(accel_kernel_t), public, target, save :: kernel_density_complex
  type(accel_kernel_t), public, target, save :: kernel_density_spinors
  type(accel_kernel_t), public, target, save :: kernel_phase
  type(accel_kernel_t), public, target, save :: kernel_phase_spiral
  type(accel_kernel_t), public, target, save :: zkernel_dot_matrix_spinors
  type(accel_kernel_t), public, target, save :: dzmul_batch
  type(accel_kernel_t), public, target, save :: zzmul_batch
 
  interface accel_grid_size
    module procedure accel_grid_size_i8, accel_grid_size_i4, accel_grid_size_array_i8, accel_grid_size_array_i4
  end interface accel_grid_size
 
  interface accel_grid_size_extend_dim
    module procedure accel_grid_size_extend_dim_i8, accel_grid_size_extend_dim_i4
  end interface accel_grid_size_extend_dim
 
  interface accel_padded_size
    module procedure accel_padded_size_i8, accel_padded_size_i4
  end interface accel_padded_size
 
  interface accel_create_buffer
    module procedure accel_create_buffer_4, accel_create_buffer_8
  end interface accel_create_buffer
 
  interface accel_ensure_buffer_size
    module procedure accel_ensure_buffer_size_4, accel_ensure_buffer_size_8
  end interface accel_ensure_buffer_size
 
  interface accel_kernel_run
    module procedure accel_kernel_run_4, accel_kernel_run_8
  end interface accel_kernel_run
 
  interface accel_set_buffer_to_zero
    module procedure accel_set_buffer_to_zero_i8, accel_set_buffer_to_zero_i4
  end interface accel_set_buffer_to_zero
 
  interface accel_write_buffer
    module procedure iaccel_write_buffer_single, laccel_write_buffer_single, daccel_write_buffer_single, zaccel_write_buffer_single
    module procedure iaccel_write_buffer_0, laccel_write_buffer_0, daccel_write_buffer_0, zaccel_write_buffer_0
    module procedure iaccel_write_buffer_1, laccel_write_buffer_1, daccel_write_buffer_1, zaccel_write_buffer_1
    module procedure iaccel_write_buffer_2, laccel_write_buffer_2, daccel_write_buffer_2, zaccel_write_buffer_2
    module procedure iaccel_write_buffer_3, laccel_write_buffer_3, daccel_write_buffer_3, zaccel_write_buffer_3
    module procedure iaccel_write_buffer_4, laccel_write_buffer_4, daccel_write_buffer_4, zaccel_write_buffer_4
    module procedure iaccel_write_buffer_5, laccel_write_buffer_5, daccel_write_buffer_5, zaccel_write_buffer_5
    module procedure iaccel_write_buffer_6, laccel_write_buffer_6, daccel_write_buffer_6, zaccel_write_buffer_6
    module procedure iaccel_write_buffer_0_int32, laccel_write_buffer_0_int32, daccel_write_buffer_0_int32, &
      zaccel_write_buffer_0_int32
    module procedure iaccel_write_buffer_1_int32, laccel_write_buffer_1_int32, daccel_write_buffer_1_int32, &
      zaccel_write_buffer_1_int32
    module procedure iaccel_write_buffer_2_int32, laccel_write_buffer_2_int32, daccel_write_buffer_2_int32, &
      zaccel_write_buffer_2_int32
    module procedure iaccel_write_buffer_3_int32, laccel_write_buffer_3_int32, daccel_write_buffer_3_int32, &
      zaccel_write_buffer_3_int32
    module procedure iaccel_write_buffer_4_int32, laccel_write_buffer_4_int32, daccel_write_buffer_4_int32, &
      zaccel_write_buffer_4_int32
    module procedure iaccel_write_buffer_5_int32, laccel_write_buffer_5_int32, daccel_write_buffer_5_int32, &
      zaccel_write_buffer_5_int32
    module procedure iaccel_write_buffer_6_int32, laccel_write_buffer_6_int32, daccel_write_buffer_6_int32, &
      zaccel_write_buffer_6_int32
    module procedure saccel_write_buffer_single, caccel_write_buffer_single
    module procedure saccel_write_buffer_0, caccel_write_buffer_0, saccel_write_buffer_1, caccel_write_buffer_1
    module procedure saccel_write_buffer_2, caccel_write_buffer_2, saccel_write_buffer_3, caccel_write_buffer_3
    module procedure saccel_write_buffer_4, caccel_write_buffer_4, saccel_write_buffer_5, caccel_write_buffer_5
    module procedure saccel_write_buffer_6, caccel_write_buffer_6
    module procedure saccel_write_buffer_0_int32, caccel_write_buffer_0_int32
    module procedure saccel_write_buffer_1_int32, caccel_write_buffer_1_int32
    module procedure saccel_write_buffer_2_int32, caccel_write_buffer_2_int32
    module procedure saccel_write_buffer_3_int32, caccel_write_buffer_3_int32
    module procedure saccel_write_buffer_4_int32, caccel_write_buffer_4_int32
    module procedure saccel_write_buffer_5_int32, caccel_write_buffer_5_int32
    module procedure saccel_write_buffer_6_int32, caccel_write_buffer_6_int32
  end interface accel_write_buffer
 
  interface accel_read_buffer
    module procedure iaccel_read_buffer_0, laccel_read_buffer_0, daccel_read_buffer_0, zaccel_read_buffer_0
    module procedure iaccel_read_buffer_1, laccel_read_buffer_1, daccel_read_buffer_1, zaccel_read_buffer_1
    module procedure iaccel_read_buffer_2, laccel_read_buffer_2, daccel_read_buffer_2, zaccel_read_buffer_2
    module procedure iaccel_read_buffer_3, laccel_read_buffer_3, daccel_read_buffer_3, zaccel_read_buffer_3
    module procedure iaccel_read_buffer_4, laccel_read_buffer_4, daccel_read_buffer_4, zaccel_read_buffer_4
    module procedure iaccel_read_buffer_5, laccel_read_buffer_5, daccel_read_buffer_5, zaccel_read_buffer_5
    module procedure iaccel_read_buffer_6, laccel_read_buffer_6, daccel_read_buffer_6, zaccel_read_buffer_6
    module procedure iaccel_read_buffer_0_int32, laccel_read_buffer_0_int32, daccel_read_buffer_0_int32, zaccel_read_buffer_0_int32
    module procedure iaccel_read_buffer_1_int32, laccel_read_buffer_1_int32, daccel_read_buffer_1_int32, zaccel_read_buffer_1_int32
    module procedure iaccel_read_buffer_2_int32, laccel_read_buffer_2_int32, daccel_read_buffer_2_int32, zaccel_read_buffer_2_int32
    module procedure iaccel_read_buffer_3_int32, laccel_read_buffer_3_int32, daccel_read_buffer_3_int32, zaccel_read_buffer_3_int32
    module procedure iaccel_read_buffer_4_int32, laccel_read_buffer_4_int32, daccel_read_buffer_4_int32, zaccel_read_buffer_4_int32
    module procedure iaccel_read_buffer_5_int32, laccel_read_buffer_5_int32, daccel_read_buffer_5_int32, zaccel_read_buffer_5_int32
    module procedure iaccel_read_buffer_6_int32, laccel_read_buffer_6_int32, daccel_read_buffer_6_int32, zaccel_read_buffer_6_int32
    module procedure saccel_read_buffer_0, caccel_read_buffer_0, saccel_read_buffer_1, caccel_read_buffer_1
    module procedure saccel_read_buffer_2, caccel_read_buffer_2, saccel_read_buffer_3, caccel_read_buffer_3
    module procedure saccel_read_buffer_4, caccel_read_buffer_4, saccel_read_buffer_5, caccel_read_buffer_5
    module procedure saccel_read_buffer_6, caccel_read_buffer_6
    module procedure saccel_read_buffer_0_int32, caccel_read_buffer_0_int32
    module procedure saccel_read_buffer_1_int32, caccel_read_buffer_1_int32
    module procedure saccel_read_buffer_2_int32, caccel_read_buffer_2_int32
    module procedure saccel_read_buffer_3_int32, caccel_read_buffer_3_int32
    module procedure saccel_read_buffer_4_int32, caccel_read_buffer_4_int32
    module procedure saccel_read_buffer_5_int32, caccel_read_buffer_5_int32
    module procedure saccel_read_buffer_6_int32, caccel_read_buffer_6_int32
  end interface accel_read_buffer
 
  interface accel_set_kernel_arg
    module procedure                       &
      accel_set_kernel_arg_buffer,  &
      iaccel_set_kernel_arg_data,   &
      laccel_set_kernel_arg_data,   &
      daccel_set_kernel_arg_data,   &
      zaccel_set_kernel_arg_data,   &
      saccel_set_kernel_arg_data,   &
      caccel_set_kernel_arg_data
  end interface accel_set_kernel_arg
 
  interface accel_get_device_pointer
    module procedure iaccel_get_device_pointer_1, laccel_get_device_pointer_1
    module procedure iaccel_get_device_pointer_2, laccel_get_device_pointer_2
    module procedure iaccel_get_device_pointer_3, laccel_get_device_pointer_3
    module procedure daccel_get_device_pointer_1, zaccel_get_device_pointer_1
    module procedure daccel_get_device_pointer_2, zaccel_get_device_pointer_2
    module procedure daccel_get_device_pointer_3, zaccel_get_device_pointer_3
    module procedure iaccel_get_device_pointer_1l, laccel_get_device_pointer_1l
    module procedure iaccel_get_device_pointer_2l, laccel_get_device_pointer_2l
    module procedure iaccel_get_device_pointer_3l, laccel_get_device_pointer_3l
    module procedure daccel_get_device_pointer_1l, zaccel_get_device_pointer_1l
    module procedure daccel_get_device_pointer_2l, zaccel_get_device_pointer_2l
    module procedure daccel_get_device_pointer_3l, zaccel_get_device_pointer_3l
    module procedure saccel_get_device_pointer_1, caccel_get_device_pointer_1
    module procedure saccel_get_device_pointer_2, caccel_get_device_pointer_2
    module procedure saccel_get_device_pointer_3, caccel_get_device_pointer_3
    module procedure saccel_get_device_pointer_1l, caccel_get_device_pointer_1l
    module procedure saccel_get_device_pointer_2l, caccel_get_device_pointer_2l
    module procedure saccel_get_device_pointer_3l, caccel_get_device_pointer_3l
  end interface accel_get_device_pointer
 
  integer :: buffer_alloc_count
  integer(int64) :: allocated_mem
  type(accel_kernel_t), pointer :: head
  type(alloc_cache_t) :: memcache
 
contains
 
  pure logical function accel_is_enabled() result(enabled)
#ifdef HAVE_ACCEL
    enabled = accel%enabled
#else
    enabled = .false.
#endif
  end function accel_is_enabled
 
  ! ------------------------------------------
 
  pure logical function accel_allow_cpu_only() result(allow)
#ifdef HAVE_ACCEL
    allow = accel%allow_CPU_only
#else
    allow = .true.
#endif
  end function accel_allow_cpu_only
 
  ! ------------------------------------------
 
  subroutine accel_init(base_grp, namespace)
    type(mpi_grp_t),     intent(inout) :: base_grp
    type(namespace_t),   intent(in)    :: namespace
 
    logical  :: disable, default, run_benchmark
    integer  :: idevice
#ifdef HAVE_CUDA
    integer :: dim
#ifdef HAVE_MPI
    character(len=256) :: sys_name
#endif
#endif
 
    push_sub(accel_init)
 
    buffer_alloc_count = 0
 
    !%Variable DisableAccel
    !%Type logical
    !%Default yes
    !%Section Execution::Accel
    !%Description
    !% If Octopus was compiled with CUDA support, it will
    !% try to initialize and use an accelerator device. By setting this
    !% variable to <tt>yes</tt> you force Octopus not to use an accelerator even it is available.
    !%End
    call messages_obsolete_variable(namespace, 'DisableOpenCL', 'DisableAccel')
#ifdef HAVE_ACCEL
    default = .false.
#else
    default = .true.
#endif
    call parse_variable(namespace, 'DisableAccel', default, disable)
    accel%enabled = .not. disable
 
#ifndef HAVE_ACCEL
    if (accel%enabled) then
      message(1) = 'Octopus was compiled without Cuda support.'
      call messages_fatal(1)
    end if
#endif
 
    if (.not. accel_is_enabled()) then
      pop_sub(accel_init)
      return
    end if
 
    call messages_obsolete_variable(namespace, 'AccelPlatform')
    call messages_obsolete_variable(namespace, 'OpenCLPlatform', 'AccelPlatform')
 
    !%Variable AccelDevice
    !%Type integer
    !%Default 0
    !%Section Execution::Accel
    !%Description
    !% This variable selects the GPU that Octopus will use. You can specify a
    !% numerical id to select a specific device.
    !%
    !% In case of MPI enabled runs devices are distributed in a round robin fashion,
    !% starting at this value.
    !%End
    call parse_variable(namespace, 'AccelDevice', 0, idevice)
 
    call messages_obsolete_variable(namespace, 'OpenCLDevice', 'AccelDevice')
 
    if (idevice < 0) then
      call messages_write('Invalid AccelDevice')
      call messages_fatal()
    end if
 
    call messages_print_with_emphasis(msg="GPU acceleration", namespace=namespace)
 
#ifdef HAVE_CUDA
    if (idevice<0) idevice = 0
    call cuda_init(accel%context%cuda_context, accel%device%cuda_device, accel%cuda_stream, &
      idevice, base_grp%rank)
#ifdef HAVE_MPI
    call loct_sysname(sys_name)
    write(message(1), '(A,I5,A,I5,2A)') "Rank ", base_grp%rank, " uses device number ", idevice, &
      " on ", trim(sys_name)
    call messages_info(1, all_nodes = .true.)
#endif
 
    call cublas_init(accel%cuda_stream)
#endif
 
 
    ! Get some device information that we will need later
#ifdef HAVE_CUDA
    call cuda_device_total_memory(accel%device%cuda_device, accel%global_memory_size)
    call cuda_device_shared_memory(accel%device%cuda_device, accel%shared_memory_size)
    call cuda_device_max_threads_per_block(accel%device%cuda_device, accel%max_block_size)
    call cuda_device_get_warpsize(accel%device%cuda_device, accel%warp_size)
    call cuda_device_max_block_dim_x(accel%device%cuda_device, dim)
    accel%max_block_dim(1) = int(dim, int64)
    call cuda_device_max_block_dim_y(accel%device%cuda_device, dim)
    accel%max_block_dim(2) = int(dim, int64)
    call cuda_device_max_block_dim_z(accel%device%cuda_device, dim)
    accel%max_block_dim(3) = int(dim, int64)
    call cuda_device_max_grid_dim_x(accel%device%cuda_device, dim)
    accel%max_grid_dim(1) = int(dim, int64)
    call cuda_device_max_grid_dim_y(accel%device%cuda_device, dim)
    accel%max_grid_dim(2) = int(dim, int64)
    call cuda_device_max_grid_dim_z(accel%device%cuda_device, dim)
    accel%max_grid_dim(3) = int(dim, int64)
#endif
 
    if (base_grp%is_root()) call device_info()
 
    ! initialize the cache used to speed up allocations
    call alloc_cache_init(memcache, nint(0.25_real64*accel%global_memory_size, int64))
 
    ! now initialize the kernels
    call accel_kernel_global_init()
 
#if defined(HAVE_HIP)
    accel%debug_flag = "-g"
#elif defined(HAVE_CUDA)
    accel%debug_flag = "-lineinfo"
#endif
 
    call accel_kernel_start_call(kernel_vpsi, 'vpsi.cu', "vpsi<double>")
    call accel_kernel_start_call(kernel_vpsi_complex, 'vpsi.cu', "vpsi_complex<double>")
    call accel_kernel_start_call(kernel_vpsi_spinors, 'vpsi.cu', "vpsi_spinors<double>")
    call accel_kernel_start_call(kernel_vpsi_spinors_complex, 'vpsi.cu', "vpsi_spinors_complex<double>")
    call accel_kernel_start_call(kernel_daxpy, 'axpy.cu', "axpy<double>")
    call accel_kernel_start_call(kernel_zaxpy, 'axpy.cu', "axpy<complex<double>>")
    call accel_kernel_start_call(kernel_saxpy, 'axpy.cu', "axpy<float>")
    call accel_kernel_start_call(kernel_caxpy, 'axpy.cu', "axpy<complex<float>>")
    call accel_kernel_start_call(dpack, 'pack.cu', "pack<double>")
    call accel_kernel_start_call(zpack, 'pack.cu', "pack<complex<double>>")
    call accel_kernel_start_call(spack, 'pack.cu', "pack<float>")
    call accel_kernel_start_call(cpack, 'pack.cu', "pack<complex<float>>")
    call accel_kernel_start_call(dunpack, 'pack.cu', "unpack<double>")
    call accel_kernel_start_call(zunpack, 'pack.cu', "unpack<complex<double>>")
    call accel_kernel_start_call(sunpack, 'pack.cu', "unpack<float>")
    call accel_kernel_start_call(cunpack, 'pack.cu', "unpack<complex<float>>")
    call accel_kernel_start_call(kernel_copy, 'copy.cu', "copy<double>")
    call accel_kernel_start_call(kernel_scopy, 'copy.cu', "copy<float>")
    call accel_kernel_start_call(kernel_copy_complex_to_real, 'copy.cu', "copy_complex_to_real<double>")
    call accel_kernel_start_call(kernel_copy_real_to_complex, 'copy.cu', "copy_real_to_complex<double>")
    call accel_kernel_start_call(kernel_ghost_reorder, 'ghost.cu', "ghost_reorder<double>")
    call accel_kernel_start_call(kernel_density_real, 'density.cu', "density<double, double>")
    call accel_kernel_start_call(kernel_density_complex, 'density.cu', "density<double, complex<double>>")
    call accel_kernel_start_call(kernel_density_spinors, 'density.cu', "density_spinors<double>")
    call accel_kernel_start_call(kernel_phase, 'phase.cu', "phase<double>")
    call accel_kernel_start_call(zkernel_dot_matrix_spinors, 'mesh_batch.cu', "zdot_matrix_spinors<double>")
 
 
    call accel_kernel_start_call(dzmul_batch, 'mul.cu', 'zmul_batch<double, double>')
    call accel_kernel_start_call(zzmul_batch, 'mul.cu', 'zmul_batch<double, complex<double>>')
 
 
    !%Variable AccelBenchmark
    !%Type logical
    !%Default no
    !%Section Execution::Accel
    !%Description
    !% If this variable is set to yes, Octopus will run some
    !% routines to benchmark the performance of the accelerator device.
    !%End
    call parse_variable(namespace, 'AccelBenchmark', .false., run_benchmark)
 
    call messages_obsolete_variable(namespace, 'OpenCLBenchmark', 'AccelBenchmark')
 
    if (run_benchmark) then
      call accel_check_bandwidth()
    end if
 
    !%Variable GPUAwareMPI
    !%Type logical
    !%Section Execution::Accel
    !%Description
    !% If Octopus was compiled with GPU support and MPI support and if the MPI
    !% implementation is GPU-aware (i.e., it supports communication using device pointers),
    !% this switch can be set to true to use the GPU-aware MPI features. The advantage
    !% of this approach is that it can do, e.g., peer-to-peer copies between devices without
    !% going through the host memory.
    !% The default is false, except when the configure switch --enable-cudampi is set, in which
    !% case this variable is set to true.
    !%End
#ifdef HAVE_CUDA_MPI
    default = .true.
#else
    default = .false.
#endif
    call parse_variable(namespace, 'GPUAwareMPI', default, accel%cuda_mpi)
    if (accel%cuda_mpi) then
#ifndef HAVE_CUDA_MPI
      call messages_write("Warning: trying to use GPU-aware MPI, but we have not detected support in the linked MPI library.")
      call messages_warning()
#endif
      call messages_write("Using GPU-aware MPI.")
      call messages_info()
    end if
 
 
    !%Variable AllowCPUonly
    !%Type logical
    !%Section Execution::Accel
    !%Description
    !% In order to prevent waste of resources, the code will normally stop when the GPU is disabled due to
    !% incomplete implementations or incompatibilities. AllowCPUonly = yes overrides this and allows the
    !% code execution also in these cases.
    !%End
#if defined (HAVE_ACCEL)
    default = .false.
#else
    default = .true.
#endif
    call parse_variable(namespace, 'AllowCPUonly', default, accel%allow_CPU_only)
 
 
    !%Variable InitializeGPUBuffers
    !%Type integer
    !%Default no
    !%Section Execution::Accel
    !%Description
    !% Initialize new GPU buffers to zero on creation (use only for debugging, as it has a performance impact!).
    !%Option no 0
    !% Do not initialize GPU buffers.
    !%Option yes 1
    !% Initialize GPU buffers to zero.
    !%Option nan 2
    !% Initialize GPU buffers to nan.
    !%End
    call parse_variable(namespace, 'InitializeGPUBuffers', option__initializegpubuffers__no, accel%initialize_buffers)
    if (.not. varinfo_valid_option('InitializeGPUBuffers', accel%initialize_buffers)) then
      call messages_input_error(namespace, 'InitializeGPUBuffers')
    end if
 
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(accel_init)
 
  contains
 
    subroutine device_info()
#ifdef HAVE_CUDA
      integer :: version
      character(kind=c_char) :: cval_str(257)
#endif
      integer :: major, minor
      character(len=256) :: val_str
 
      push_sub(accel_init.device_info)
 
      call messages_new_line()
      call messages_write('Selected device:')
      call messages_new_line()
 
#ifdef HAVE_CUDA
#ifdef __HIP_PLATFORM_AMD__
      call messages_write('      Framework              : ROCm')
#else
      call messages_write('      Framework              : CUDA')
#endif
#endif
      call messages_info()
 
#ifdef HAVE_CUDA
      call messages_write('      Device type            : GPU', new_line = .true.)
#ifdef __HIP_PLATFORM_AMD__
      call messages_write('      Device vendor          : AMD Corporation', new_line = .true.)
#else
      call messages_write('      Device vendor          : NVIDIA Corporation', new_line = .true.)
#endif
#endif
 
#ifdef HAVE_CUDA
      cval_str = c_null_char
      call cuda_device_name(accel%device%cuda_device, cval_str)
      call string_c_to_f(cval_str, val_str)
#endif
      call messages_write('      Device name            : '//trim(val_str))
      call messages_new_line()
 
#ifdef HAVE_CUDA
      call cuda_device_capability(accel%device%cuda_device, major, minor)
#endif
      call messages_write('      Cuda capabilities      :')
      call messages_write(major, fmt = '(i2)')
      call messages_write('.')
      call messages_write(minor, fmt = '(i1)')
      call messages_new_line()
 
      ! VERSION
#ifdef HAVE_CUDA
      call cuda_driver_version(version)
      call messages_write('      Driver version         : ')
      call messages_write(version)
#endif
      call messages_new_line()
 
 
      call messages_write('      Device memory          :')
      call messages_write(accel%global_memory_size, units=unit_megabytes)
      call messages_new_line()
 
      call messages_write('      Shared memory          :')
      call messages_write(accel%shared_memory_size, units=unit_kilobytes)
      call messages_new_line()
 
      call messages_write('      Max. block size        :')
      call messages_write(accel%max_block_size)
      call messages_new_line()
 
      call messages_info()
 
      pop_sub(accel_init.device_info)
    end subroutine device_info
 
  end subroutine accel_init
 
  ! ------------------------------------------
  subroutine accel_end(namespace)
    type(namespace_t), intent(in) :: namespace
 
    integer(int64) :: hits, misses
    real(real64) :: volume_hits, volume_misses
    logical :: found
    type(accel_mem_t) :: tmp
 
    push_sub(accel_end)
 
    if (accel_is_enabled()) then
 
      ! Release global buffers
 
      do
        call alloc_cache_get(memcache, alloc_cache_any_size, found, tmp%mem)
        if (.not. found) exit
 
#ifdef HAVE_CUDA
        call cuda_mem_free(tmp%mem)
#endif
      end do
 
      call alloc_cache_end(memcache, hits, misses, volume_hits, volume_misses)
 
      call messages_print_with_emphasis(msg="Acceleration-device allocation cache", namespace=namespace)
 
      call messages_new_line()
      call messages_write('    Number of allocations    =')
      call messages_write(hits + misses, new_line = .true.)
      call messages_write('    Volume of allocations    =')
      call messages_write(volume_hits + volume_misses, fmt = 'f18.1', units = unit_gigabytes, align_left = .true., &
        new_line = .true.)
      call messages_write('    Hit ratio                =')
      if (hits + misses > 0) then
        call messages_write(hits/real(hits + misses, real64)*100, fmt='(f6.1)', align_left = .true.)
      else
        call messages_write(m_zero, fmt='(f6.1)', align_left = .true.)
      end if
      call messages_write('%', new_line = .true.)
      call messages_write('    Volume hit ratio         =')
      if (volume_hits + volume_misses > 0) then
        call messages_write(volume_hits/(volume_hits + volume_misses)*100, fmt='(f6.1)', align_left = .true.)
      else
        call messages_write(m_zero, fmt='(f6.1)', align_left = .true.)
      end if
      call messages_write('%')
      call messages_new_line()
      call messages_info()
 
      call messages_print_with_emphasis(namespace=namespace)
    end if
 
    call accel_kernel_global_end()
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cublas_end()
      if (.not. accel%cuda_mpi) then ! CUDA aware MPI finalize will do the cleanup
        call cuda_end(accel%context%cuda_context, accel%device%cuda_device)
      end if
#endif
 
      if (buffer_alloc_count /= 0) then
        call messages_write('Accel:')
        call messages_write(real(allocated_mem, real64) , fmt = 'f12.1', units = unit_megabytes, align_left = .true.)
        call messages_write(' in ')
        call messages_write(buffer_alloc_count)
        call messages_write(' buffers were not deallocated.')
        call messages_fatal()
      end if
 
    end if
 
    pop_sub(accel_end)
  end subroutine accel_end
 
  ! ------------------------------------------
 
  subroutine accel_grid_size_array_i8(n, blocksizes, gridsizes)
    integer(int64), intent(in)  :: n(:)
    integer(int64), intent(in)  :: blocksizes(:)
    integer(int64), intent(out) :: gridsizes(:)
 
    integer :: dim, i
 
    dim = ubound(n, dim=1)
    assert(dim == ubound(blocksizes, dim=1))
    assert(dim == ubound(gridsizes, dim=1))
 
    do i = 1, dim
      gridsizes(i) = (n(i) + blocksizes(i) - 1_int64) / blocksizes(i)
      gridsizes(i) = min(gridsizes(i), accel%max_grid_dim(i))
    end do
  end subroutine accel_grid_size_array_i8
 
  ! ------------------------------------------
 
  subroutine accel_grid_size_array_i4(n, blocksizes, gridsizes)
    integer,       intent(in)  :: n(:)
    integer,       intent(in)  :: blocksizes(:)
    integer,       intent(out) :: gridsizes(:)
 
    integer(int64) :: gridsizes64(size(gridsizes))
 
    call accel_grid_size_array_i8(int(n(:), int64), int(blocksizes(:), int64), gridsizes64)
 
    gridsizes = int(gridsizes64, int32)
  end subroutine accel_grid_size_array_i4
 
  ! ------------------------------------------
 
  subroutine accel_grid_size_i8(n, blocksizes, gridsizes)
    integer(int64), intent(in)  :: n
    integer(int64), intent(in)  :: blocksizes
    integer(int64), intent(out) :: gridsizes
 
    integer(int64) :: temp(1)
 
    call accel_grid_size_array_i8( (/n/), (/blocksizes/), temp)
 
    gridsizes = temp(1)
  end subroutine accel_grid_size_i8
 
  ! ------------------------------------------
 
  subroutine accel_grid_size_i4(n, blocksizes, gridsizes)
    integer,       intent(in)  :: n
    integer,       intent(in)  :: blocksizes
    integer,       intent(out) :: gridsizes
 
    integer(int64) :: temp(1)
 
    call accel_grid_size_array_i8(int((/n/), int64), int((/blocksizes/), int64), temp)
 
    gridsizes = int(temp(1), int32)
  end subroutine accel_grid_size_i4
 
! ------------------------------------------
 
  subroutine accel_grid_size_extend_dim_i8(n, pack_size, gridsizes, blocksizes, kernel)
    integer(int64),                   intent(in)  :: n
    integer(int64),                   intent(in)  :: pack_size
    integer(int64), dimension(3),     intent(out) :: gridsizes
    integer(int64), dimension(3),     intent(out) :: blocksizes
    type(accel_kernel_t), optional, intent(inout) :: kernel
 
    integer(int64) :: bsize, dim2, dim3
    integer(int64), dimension(3) :: nn
 
    if(present(kernel)) then
      bsize = accel_kernel_block_size(kernel)/pack_size
    else
      bsize = accel_max_block_size()/pack_size
    end if
 
    dim3 = n/(accel_max_size_per_dim(2)*bsize) + 1
    dim2 = min(accel_max_size_per_dim(2)*bsize, pad(n, bsize))
 
    nn = (/pack_size, dim2, dim3/)
    blocksizes = (/pack_size, bsize, 1_int64/)
 
    call accel_grid_size(nn, blocksizes, gridsizes)
  end subroutine accel_grid_size_extend_dim_i8
 
  ! ------------------------------------------
 
  subroutine accel_grid_size_extend_dim_i4(n, pack_size, gridsizes, blocksizes, kernel)
    integer,                        intent(in)  :: n
    integer,                        intent(in)  :: pack_size
    integer, dimension(3),          intent(out) :: gridsizes
    integer, dimension(3),          intent(out) :: blocksizes
    type(accel_kernel_t), optional, intent(inout) :: kernel
 
    integer(int64) :: gridsizes64(3), blocksizes64(3)
 
    call accel_grid_size_extend_dim_i8(int(n, int64), int(pack_size, int64), &
      gridsizes64, blocksizes64, kernel=kernel)
 
    gridsizes  = int(gridsizes64, int32)
    blocksizes = int(blocksizes64, int32)
  end subroutine accel_grid_size_extend_dim_i4
 
  ! ------------------------------------------
 
  integer(int64) function accel_padded_size_i8(nn) result(psize)
    integer(int64), intent(in) :: nn
 
    integer(int64) :: modnn, bsize
 
    psize = nn
 
    if (accel_is_enabled()) then
 
      bsize = accel_max_block_size()
 
      psize = nn
      modnn = mod(nn, bsize)
      if (modnn /= 0) psize = psize + bsize - modnn
 
    end if
 
  end function accel_padded_size_i8
 
  ! ------------------------------------------
 
  integer(int32) function accel_padded_size_i4(nn) result(psize)
    integer(int32), intent(in) :: nn
 
    psize = int(accel_padded_size_i8(int(nn, int64)), int32)
 
  end function accel_padded_size_i4
 
  ! ------------------------------------------
 
  subroutine accel_create_buffer_4(this, flags, type, size, set_zero, async)
    type(accel_mem_t),  intent(inout) :: this
    integer,            intent(in)    :: flags
    type(type_t),       intent(in)    :: type
    integer,            intent(in)    :: size
    logical,  optional, intent(in)    :: set_zero
    logical,  optional, intent(in)    :: async
 
    call accel_create_buffer_8(this, flags, type, int(size, int64), set_zero, async)
  end subroutine accel_create_buffer_4
 
  ! ------------------------------------------
 
  subroutine accel_create_buffer_8(this, flags, type, size, set_zero, async)
    type(accel_mem_t),  intent(inout) :: this
    integer,            intent(in)    :: flags
    type(type_t),       intent(in)    :: type
    integer(int64),     intent(in)    :: size
    logical,  optional, intent(in)    :: set_zero
    logical,  optional, intent(in)    :: async
 
    integer(int64) :: fsize
    logical    :: found
    integer(int64) :: initialize_buffers
 
    push_sub(accel_create_buffer_8)
 
    this%type = type
    this%size = size
    this%flags = flags
    fsize = int(size, int64)*types_get_size(type)
    this%allocated = .true.
 
    if (fsize > 0) then
 
      call alloc_cache_get(memcache, fsize, found, this%mem)
 
      if (.not. found) then
#ifdef HAVE_CUDA
        if(optional_default(async, .false.)) then
          call cuda_mem_alloc_async(this%mem, fsize)
        else
          call cuda_mem_alloc(this%mem, fsize)
        end if
#endif
      end if
 
      buffer_alloc_count = buffer_alloc_count + 1
      allocated_mem = allocated_mem + fsize
 
    end if
 
    if (present(set_zero)) then
      initialize_buffers = merge(option__initializegpubuffers__yes, option__initializegpubuffers__no, set_zero)
    else
      initialize_buffers = accel%initialize_buffers
    end if
    select case (initialize_buffers)
    case (option__initializegpubuffers__yes)
      call accel_set_buffer_to(this, type, int(z'00', int8), size, async=async)
    case (option__initializegpubuffers__nan)
      call accel_set_buffer_to(this, type, int(z'FF', int8), size, async=async)
    end select
 
    pop_sub(accel_create_buffer_8)
  end subroutine accel_create_buffer_8
 
  ! ------------------------------------------
 
  subroutine accel_free_buffer(this, async)
    type(accel_mem_t), intent(inout) :: this
    logical, optional, intent(in) :: async
 
    logical :: put
    integer(int64) :: fsize
 
    push_sub(accel_free_buffer)
 
    if (this%size > 0) then
 
      fsize = int(this%size, int64)*types_get_size(this%type)
 
      call alloc_cache_put(memcache, fsize, this%mem, put)
 
      if (.not. put) then
#ifdef HAVE_CUDA
        if (optional_default(async, .false.)) then
          call cuda_mem_free_async(this%mem)
        else
          call cuda_mem_free(this%mem)
        end if
#endif
      end if
 
      buffer_alloc_count = buffer_alloc_count - 1
      allocated_mem = allocated_mem + fsize
 
    end if
 
    this%size = 0
    this%flags = 0
 
    this%allocated = .false.
 
    pop_sub(accel_free_buffer)
  end subroutine accel_free_buffer
 
  ! ------------------------------------------
 
  subroutine accel_move_buffer(buffer_from, buffer_to)
    type(accel_mem_t), intent(inout) :: buffer_from
    type(accel_mem_t), intent(inout) :: buffer_to
 
    push_sub(accel_move_buffer)
 
    ! Can not move an allocated buffer, the destination buffer should be empty
    assert(.not. accel_buffer_is_allocated(buffer_to))
 
    buffer_to%mem = buffer_from%mem
    buffer_to%size = buffer_from%size
    buffer_to%type = buffer_from%type
    buffer_to%flags = buffer_from%flags
    buffer_to%allocated = buffer_from%allocated
 
    call accel_detach_buffer(buffer_from)
 
    pop_sub(accel_move_buffer)
  end subroutine accel_move_buffer
 
  ! ------------------------------------------
 
  subroutine accel_detach_buffer(this)
    type(accel_mem_t), intent(inout) :: this
 
    push_sub(accel_detach_buffer)
 
    this%mem = c_null_ptr
    this%size = 0
    this%type = type_none
    this%flags = 0
    this%allocated = .false.
 
    pop_sub(accel_detach_buffer)
  end subroutine accel_detach_buffer
 
  ! ------------------------------------------------------
 
  ! Check if the temporary buffers are the right size, if not reallocate them
  subroutine accel_ensure_buffer_size_4(buffer, flags, type, required_size, set_zero, async)
    type(accel_mem_t), intent(inout) :: buffer
    integer,           intent(in)    :: flags
    type(type_t),      intent(in)    :: type
    integer,           intent(in)    :: required_size
    logical,           intent(in)    :: set_zero
    logical, optional, intent(in)    :: async
 
    call accel_ensure_buffer_size_8(buffer, flags, type, int(required_size, int64), set_zero, async)
  end subroutine accel_ensure_buffer_size_4
 
  ! ------------------------------------------
 
  subroutine accel_ensure_buffer_size_8(buffer, flags, type, required_size, set_zero, async)
    type(accel_mem_t), intent(inout) :: buffer
    integer,           intent(in)    :: flags
    type(type_t),      intent(in)    :: type
    integer(int64),    intent(in)    :: required_size
    logical,           intent(in)    :: set_zero
    logical, optional, intent(in)    :: async
 
    push_sub(accel_ensure_buffer_size_8)
 
 
    if (accel_buffer_is_allocated(buffer) .and. buffer%size < required_size) then
      call accel_free_buffer(buffer, async=optional_default(async, .false.))
    end if
 
    if (.not. accel_buffer_is_allocated(buffer)) then
      call accel_create_buffer(buffer, flags, type, required_size, set_zero=set_zero, async=optional_default(async, .false.))
    end if
 
    pop_sub(accel_ensure_buffer_size_8)
  end subroutine accel_ensure_buffer_size_8
 
  ! ------------------------------------------
 
  logical pure function accel_buffer_is_allocated(this) result(allocated)
    type(accel_mem_t), intent(in) :: this
 
    allocated = this%allocated
  end function accel_buffer_is_allocated
 
  ! -----------------------------------------
 
  subroutine accel_finish()
    ! no push_sub, called too frequently
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cuda_context_synchronize()
#endif
    end if
  end subroutine accel_finish
 
  ! ------------------------------------------
 
  subroutine accel_set_kernel_arg_buffer(kernel, narg, buffer)
    type(accel_kernel_t), intent(inout) :: kernel
    integer,              intent(in)    :: narg
    type(accel_mem_t),    intent(in)    :: buffer
 
    assert(accel_buffer_is_allocated(buffer))
 
    ! no push_sub, called too frequently
#ifdef HAVE_CUDA
    call cuda_kernel_set_arg_buffer(kernel%arguments, buffer%mem, narg)
#endif
 
  end subroutine accel_set_kernel_arg_buffer
 
  ! ------------------------------------------
 
 
  subroutine accel_kernel_run_8(kernel, gridsizes, blocksizes, shared_memory_size)
    type(accel_kernel_t),  intent(inout) :: kernel
    integer(int64),           intent(in) :: gridsizes(:)
    integer(int64),           intent(in) :: blocksizes(:)
    integer(int64), optional, intent(in) :: shared_memory_size
 
    integer :: dim
    integer(int64) :: gsizes(1:3)
    integer(int64) :: bsizes(1:3)
 
    ! no push_sub, called too frequently
 
    ! CUDA needs all dimensions
    gsizes = 1
    bsizes = 1
 
    dim = ubound(gridsizes, dim=1)
 
    assert(dim == ubound(blocksizes, dim=1))
 
    ! if one size is zero, there is nothing to do
    if (any(gridsizes == 0)) return
 
    assert(all(blocksizes > 0))
 
    gsizes(1:dim) = gridsizes(1:dim)
    bsizes(1:dim) = blocksizes(1:dim)
 
#ifdef HAVE_CUDA
    ! Maximum dimension of a block
    if (any(bsizes(1:3) > accel%max_block_dim(1:3))) then
      message(1) = "Maximum dimension of a block too large in kernel "//trim(kernel%kernel_name)
      message(2) = "The following conditions should be fulfilled:"
      write(message(3), "(A, I8, A, I8)") "Dim 1: ", bsizes(1), " <= ", accel%max_block_dim(1)
      write(message(4), "(A, I8, A, I8)") "Dim 2: ", bsizes(2), " <= ", accel%max_block_dim(2)
      write(message(5), "(A, I8, A, I8)") "Dim 3: ", bsizes(3), " <= ", accel%max_block_dim(3)
      message(6) = "This is an internal error, please contact the developers."
      call messages_fatal(6)
    end if
 
 
    ! Maximum number of threads per block
    if (product(bsizes) > accel_max_block_size()) then
      message(1) = "Maximum number of threads per block too large in kernel "//trim(kernel%kernel_name)
      message(2) = "The following condition should be fulfilled:"
      write(message(3), "(I8, A, I8)") product(bsizes), " <= ", accel_max_block_size()
      message(4) = "This is an internal error, please contact the developers."
      call messages_fatal(4)
    end if
 
    ! Maximum dimensions of the grid of thread block
    if (any(gsizes(1:3) > accel%max_grid_dim(1:3))) then
      message(1) = "Maximum dimension of grid too large in kernel "//trim(kernel%kernel_name)
      message(2) = "The following conditions should be fulfilled:"
      write(message(3), "(A, I8, A, I10)") "Dim 1: ", gsizes(1), " <= ", accel%max_grid_dim(1)
      write(message(4), "(A, I8, A, I10)") "Dim 2: ", gsizes(2), " <= ", accel%max_grid_dim(2)
      write(message(5), "(A, I8, A, I10)") "Dim 3: ", gsizes(3), " <= ", accel%max_grid_dim(3)
      message(6) = "This is an internal error, please contact the developers."
      call messages_fatal(6)
    end if
 
    if(present(shared_memory_size)) then
 
      if (shared_memory_size > accel%shared_memory_size) then
        message(1) = "Shared memory too large in kernel "//trim(kernel%kernel_name)
        message(2) = "The following condition should be fulfilled:"
        message(3) = "Requested shared memory <= Available shared memory"
        write(message(4), '(a,f12.6,a)') "Requested shared memory: ", real(shared_memory_size, real64) /1024.0, " Kb"
        write(message(5), '(a,f12.6,a)') "Available shared memory: ", real(accel%shared_memory_size, real64) /1024.0, " Kb"
        message(6) = "This is an internal error, please contact the developers."
        call messages_fatal(6)
      else if (shared_memory_size <= 0) then
        message(1) = "Invalid shared memory size in kernel "//trim(kernel%kernel_name)
        write(message(2), '(a,f12.6,a)') "Shared memory size requested: ", real(shared_memory_size, real64) /1024.0, " Kb"
        message(3) = "This is an internal error, please contact the developers."
        call messages_fatal(3)
      end if
    end if
 
    call cuda_launch_kernel(kernel%cuda_kernel, gsizes(1), bsizes(1), &
      optional_default(shared_memory_size, 0_int64), kernel%arguments)
#endif
 
  end subroutine accel_kernel_run_8
 
  ! -----------------------------------------------
 
 
  subroutine accel_kernel_run_4(kernel, gridsizes, blocksizes, shared_memory_size)
    type(accel_kernel_t),  intent(inout) :: kernel
    integer,                  intent(in) :: gridsizes(:)
    integer,                  intent(in) :: blocksizes(:)
    integer(int64), optional, intent(in) :: shared_memory_size
 
    call accel_kernel_run_8(kernel, int(gridsizes, int64), int(blocksizes, int64), shared_memory_size)
 
  end subroutine accel_kernel_run_4
 
  ! -----------------------------------------------
 
  integer pure function accel_max_block_size() result(max_block_size)
    max_block_size = accel%max_block_size
  end function accel_max_block_size
 
  ! -----------------------------------------------
 
  integer function accel_kernel_block_size(kernel) result(block_size)
    type(accel_kernel_t), intent(inout) :: kernel
 
#ifdef HAVE_CUDA
    integer :: max_block_size
#endif
 
    block_size = 0
 
#ifdef HAVE_CUDA
    call cuda_kernel_max_threads_per_block(kernel%cuda_kernel, max_block_size)
    if (debug%info .and. max_block_size /= accel%max_block_size) then
      write(message(1), "(A, I5, A)") "A kernel can use only less threads per block (", max_block_size, ")", &
        "than available on the device (", accel%max_block_size, ")"
      call messages_info(1)
    end if
 
    ! recommended number of threads per block is 256 according to the CUDA best practice guide
    ! see https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/index.html#thread-and-block-heuristics
    block_size = 256
 
    ! make sure we do not use more threads per block than available for this kernel
    block_size = min(block_size, max_block_size)
#endif
 
  end function accel_kernel_block_size
 
  ! ----------------------------------------------------
 
  subroutine accel_copy_buffer(dest, src, type, nval, offset_dest, offset_src, async)
    type(accel_mem_t),        intent(inout) :: dest
    type(accel_mem_t),        intent(in)    :: src
    type(type_t),             intent(in)    :: type
    integer(int64),           intent(in)    :: nval
    integer(int64), optional, intent(in)    :: offset_dest
    integer(int64), optional, intent(in)    :: offset_src
    logical,        optional, intent(in)    :: async
 
    integer(int64) :: type_size, offset_dest_, offset_src_
 
    push_sub(accel_copy_buffer)
    call profiling_in("ACCEL_COPY_BUFFER")
 
    if (nval == 0) then
      call profiling_out("ACCEL_COPY_BUFFER")
      pop_sub(accel_copy_buffer)
      return
    end if
    assert(nval > 0)
 
    type_size = types_get_size(type)
 
    offset_dest_ = 0_int64
    if (present(offset_dest)) offset_dest_ = offset_dest*type_size
    offset_src_ = 0_int64
    if (present(offset_src)) offset_src_ = offset_src*type_size
 
    assert(offset_dest_ + nval*type_size <= int(dest%size, int64)*types_get_size(dest%type))
    assert(offset_src_ + nval*type_size <= int(src%size, int64)*types_get_size(src%type))
 
    call cuda_memcpy_dtod(dest%mem, src%mem, nval*type_size, offset_dest_, offset_src_)
 
    if (.not. optional_default(async, .false.)) call accel_finish()
 
    call profiling_out("ACCEL_COPY_BUFFER")
    pop_sub(accel_copy_buffer)
  end subroutine accel_copy_buffer
 
  ! ----------------------------------------------------
 
  subroutine accel_set_buffer_to(buffer, type, val, nval, offset, async)
    type(accel_mem_t),        intent(inout) :: buffer
    type(type_t),             intent(in)    :: type
    integer(int8),            intent(in)    :: val
    integer(int64),           intent(in)    :: nval
    integer(int64), optional, intent(in)    :: offset
    logical,        optional, intent(in)    :: async
 
    integer(int64) :: nval_, offset_, type_size
 
    push_sub(accel_set_buffer_to)
 
    if (nval == 0) then
      pop_sub(accel_set_buffer_to)
      return
    end if
    assert(nval > 0)
 
    if (present(offset)) then
      assert(offset >= 0)
      if(offset > buffer%size) then
        pop_sub(accel_set_buffer_to)
        return
      end if
    end if
 
    type_size = types_get_size(type)
 
    nval_ = nval*type_size
 
    offset_ = 0_int64
    if (present(offset)) offset_ = offset*type_size
 
    call cuda_mem_set_async(buffer%mem, val, nval_, offset_)
    if(.not. optional_default(async, .false.)) call accel_finish()
 
    pop_sub(accel_set_buffer_to)
  end subroutine accel_set_buffer_to
 
  ! ----------------------------------------------------
 
  subroutine accel_set_buffer_to_zero_i8(buffer, type, nval, offset, async)
    type(accel_mem_t),        intent(inout) :: buffer
    type(type_t),             intent(in)    :: type
    integer(int64),           intent(in)    :: nval
    integer(int64), optional, intent(in)    :: offset
    logical,        optional, intent(in)    :: async
 
    push_sub(accel_set_buffer_to_zero_i8)
 
    call accel_set_buffer_to(buffer, type, int(z'00', int8), nval, offset, async)
 
    pop_sub(accel_set_buffer_to_zero_i8)
  end subroutine accel_set_buffer_to_zero_i8
 
  ! ----------------------------------------------------
 
  subroutine accel_set_buffer_to_zero_i4(buffer, type, nval, offset, async)
    type(accel_mem_t),        intent(inout) :: buffer
    type(type_t),             intent(in)    :: type
    integer(int32),           intent(in)    :: nval
    integer(int32), optional, intent(in)    :: offset
    logical,        optional, intent(in)    :: async
 
    push_sub(accel_set_buffer_to_zero_i4)
 
    if (present(offset)) then
      call accel_set_buffer_to_zero_i8(buffer, type, int(nval, int64), int(offset, int64), async=async)
    else
      call accel_set_buffer_to_zero_i8(buffer, type, int(nval, int64), async=async)
    end if
 
    pop_sub(accel_set_buffer_to_zero_i4)
  end subroutine accel_set_buffer_to_zero_i4
 
  ! ----------------------------------------------------
 
  subroutine accel_check_bandwidth()
    integer :: itime
    integer, parameter :: times = 10
    integer :: size
    real(real64)   :: time, stime
    real(real64)   :: read_bw, write_bw
    type(accel_mem_t) :: buff
    real(real64), allocatable :: data(:)
 
    call messages_new_line()
    call messages_write('Info: Benchmarking the bandwidth between main memory and device memory')
    call messages_new_line()
    call messages_info()
 
    call messages_write(' Buffer size   Read bw  Write bw')
    call messages_new_line()
    call messages_write('       [MiB]   [MiB/s]   [MiB/s]')
    call messages_info()
 
    size = 15000
    do
      safe_allocate(data(1:size))
      call accel_create_buffer(buff, accel_mem_read_write, type_float, size)
 
      stime = loct_clock()
      do itime = 1, times
        call accel_write_buffer(buff, size, data)
        call accel_finish()
      end do
      time = (loct_clock() - stime)/real(times, real64)
 
      write_bw = real(size, real64) *8.0_real64/time
 
      stime = loct_clock()
      do itime = 1, times
        call accel_read_buffer(buff, size, data)
      end do
      call accel_finish()
 
      time = (loct_clock() - stime)/real(times, real64)
      read_bw = real(size, real64) *8.0_real64/time
 
      call messages_write(size*8.0_real64/1024.0_real64**2)
      call messages_write(write_bw/1024.0_real64**2, fmt = '(f10.1)')
      call messages_write(read_bw/1024.0_real64**2, fmt = '(f10.1)')
      call messages_info()
 
      call accel_free_buffer(buff)
 
      safe_deallocate_a(data)
 
      size = int(size*2.0)
 
      if (size > 50000000) exit
    end do
  end subroutine accel_check_bandwidth
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_global_init()
 
    push_sub(accel_kernel_global_init)
 
    nullify(head)
 
    call cuda_program_map_init(accel%program_map)
 
    call accel_kernel_list()
 
    pop_sub(accel_kernel_global_init)
  end subroutine accel_kernel_global_init
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_list()
 
    push_sub(accel_kernel_list)
 
    call accel_program_add('axpy.cu', 'axpy<double>')
    call accel_program_add('axpy.cu', 'axpy<complex<double>>')
    call accel_program_add('axpy.cu', 'axpy<float>')
    call accel_program_add('axpy.cu', 'axpy<complex<float>>')
    call accel_program_add('axpy.cu', 'axpby<double>')
    call accel_program_add('axpy.cu', 'axpby<complex<double>>')
    call accel_program_add('axpy.cu', 'scal2v<double>')
    call accel_program_add('axpy.cu', 'scal2v<complex<double>, Conjugation::No>')
    call accel_program_add('axpy.cu', 'scal2v<complex<double>, Conjugation::Yes>')
    call accel_program_add('axpy.cu', 'axpy_vec<double>')
    call accel_program_add('axpy.cu', 'axpy_vec<complex<double>>')
    call accel_program_add('axpy.cu', 'axpy_vec<float>')
    call accel_program_add('axpy.cu', 'axpy_vec<complex<float>>')
    call accel_program_add('axpy.cu', 'scal_vec<double>')
    call accel_program_add('axpy.cu', 'scal_vec<complex<double>>')
    call accel_program_add('axpy.cu', 'scal_vec<float>')
    call accel_program_add('axpy.cu', 'scal_vec<complex<float>>')
    call accel_program_add('axpy.cu', 'xpay_vec<double>')
    call accel_program_add('axpy.cu', 'xpay_vec<complex<double>>')
    call accel_program_add('axpy.cu', 'xpay_vec<float>')
    call accel_program_add('axpy.cu', 'xpay_vec<complex<float>>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<double, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<double, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<complex<double>, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<complex<double>, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<float, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<float, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<complex<float>, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_axpy_function<complex<float>, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<double, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<double, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<complex<double>, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<complex<double>, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<float, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<float, DimSpin::D2>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<complex<float>, DimSpin::D1>')
    call accel_program_add('axpy.cu', &
      'batch_ax_function_py<complex<float>, DimSpin::D2>')
 
    call accel_program_add('batch_mul.cu', 'mul_batch<double>')
    call accel_program_add('batch_mul.cu', 'mul_mf<double, double>')
    call accel_program_add('batch_mul.cu', 'mul_batch<complex<double>, Conjugation::No>')
    call accel_program_add('batch_mul.cu', 'mul_batch<complex<double>, Conjugation::Yes>')
    call accel_program_add('batch_mul.cu', 'mul_mf<complex<double>, double, Conjugation::No>')
    call accel_program_add('batch_mul.cu', 'mul_mf<complex<double>, double, Conjugation::Yes>')
    call accel_program_add('batch_mul.cu', 'mul_mf<complex<double>, complex<double>, Conjugation::No>')
    call accel_program_add('batch_mul.cu', 'mul_mf<complex<double>, complex<double>, Conjugation::Yes>')
 
    call accel_program_add('boundaries.cu', 'boundaries_periodic<double>')
    call accel_program_add('boundaries.cu', 'boundaries_periodic_corr<double>')
    call accel_program_add('boundaries.cu', 'boundaries_periodic_recv<double>')
    call accel_program_add('boundaries.cu', 'boundaries_periodic_recv_corr<double>')
    call accel_program_add('boundaries.cu', 'boundaries_periodic_send<double>')
 
    call accel_program_add('copy.cu', 'add_with_map<double>')
    call accel_program_add('copy.cu', 'copy<double>')
    call accel_program_add('copy.cu', 'copy<float>')
    call accel_program_add('copy.cu', 'copy_complex_to_real<double>')
    call accel_program_add('copy.cu', 'copy_real_to_complex<double>')
    call accel_program_add('copy.cu', 'copy_with_map<double>')
 
    call accel_program_add('curl.cu', 'curl<double>')
    call accel_program_add('curl.cu', 'curl<complex<double>>')
 
    call accel_program_add('density.cu', 'density<double, double>')
    call accel_program_add('density.cu', 'density<double, complex<double>>')
    call accel_program_add('density.cu', 'current_accumulate<double>')
    call accel_program_add('density.cu', 'density_spinors<double>')
 
    call accel_program_add('dftu_projector.cu', &
      'dftu_pos_mat_elem<double, double, Indexing::Indirect>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_pos_mat_elem<double, complex<double>, Indexing::Indirect>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_pos_mat_elem<double, complex<double>, Indexing::Direct>')
    call accel_program_add('dftu_projector.cu', 'dftu_pos_mat_elem_phase<double>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_pos_mat_elem<double, double, Indexing::Direct>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_bra<double>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_bra<complex<double>>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_bra_submesh<complex<double>>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_bra_submesh<double>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_ket<double>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_ket<complex<double>>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_ket_submesh<complex<double>>')
    call accel_program_add('dftu_projector.cu', &
      'dftu_projector_ket_submesh<double>')
 
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<double, double, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<double, double, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<double, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<double, double, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<double, double, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<double, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<double, complex<double>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<double, complex<double>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<complex<double>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<double, complex<double>, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<double, complex<double>, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<complex<double>, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<float, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<float, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<complex<float>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_codensity_batch<complex<float>, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<float, float, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<float, float, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<float, float, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<float, float, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<float, complex<float>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<float, complex<float>, DimSpin::D1>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_batch<float, complex<float>, DimSpin::D2>')
    call accel_program_add('exchange_local.cu', &
      'exchange_accumulate_return_batch<float, complex<float>, DimSpin::D2>')
 
    call accel_program_add('forces.cu', 'density_gradient<double, double>')
    call accel_program_add('forces.cu', 'density_gradient<double, complex<double>>')
 
    call accel_program_add('get_points.cu', 'get_selected_points<double>')
 
    call accel_program_add('ghost.cu', 'ghost_reorder<double>')
 
    call accel_program_add('mesh_batch.cu', 'zdot_matrix_spinors<double>')
 
    call accel_program_add('mesh_batch_single.cu', 'batch_dotpv_partial<double>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_dotpv_reduce<double, double, DimSpin::D1>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_dotpv_reduce<double, double, DimSpin::D2>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_mf_dotp<double, DimSpin::D1>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_mf_dotp<double, DimSpin::D2>')
    call accel_program_add('mesh_batch_single.cu', 'batch_dotpv_partial<complex<double>>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_dotpv_reduce<double, complex<double>, DimSpin::D1>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_dotpv_reduce<double, complex<double>, DimSpin::D2>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_mf_dotp<complex<double>, DimSpin::D1>')
    call accel_program_add('mesh_batch_single.cu', &
      'batch_mf_dotp<complex<double>, DimSpin::D2>')
 
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, double, Dim::D1>')
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, double, Dim::D2>')
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, double, Dim::D3>')
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, complex<double>, Dim::D1>')
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, complex<double>, Dim::D2>')
    call accel_program_add('mesh_interpolate.cu', &
      'mesh_interpolate<double, complex<double>, Dim::D3>')
 
    call accel_program_add('mesh_to_cube.cu', 'cube_to_mesh_batch<double>')
    call accel_program_add('mesh_to_cube.cu', 'mesh_to_cube_batch<double>')
    call accel_program_add('mesh_to_cube.cu', 'cube_to_mesh_batch<complex<double>>')
    call accel_program_add('mesh_to_cube.cu', 'mesh_to_cube_batch<complex<double>>')
    call accel_program_add('mesh_to_cube.cu', 'submesh_to_cube<double>')
    call accel_program_add('mesh_to_cube.cu', 'cube_to_submesh<double>')
    call accel_program_add('mesh_to_cube.cu', 'submesh_to_cube<complex<double>>')
    call accel_program_add('mesh_to_cube.cu', 'cube_to_submesh<complex<double>>')
    call accel_program_add('mesh_to_cube.cu', 'cube_to_mesh_batch<float>')
    call accel_program_add('mesh_to_cube.cu', 'mesh_to_cube_batch<float>')
    call accel_program_add('mesh_to_cube.cu', 'cube_to_mesh_batch<complex<float>>')
    call accel_program_add('mesh_to_cube.cu', 'mesh_to_cube_batch<complex<float>>')
 
    call accel_program_add('mul.cu', 'zmul_batch<double, double>')
    call accel_program_add('mul.cu', 'zmul_batch<double, complex<double>>')
    call accel_program_add('mul.cu', 'zmul_batch<float, float>')
    call accel_program_add('mul.cu', 'zmul_batch<float, complex<float>>')
 
    call accel_program_add('operate.cu', 'operate<double>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Direct>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Indirect>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Direct, Symmetry::Antisymmetric>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Direct, Symmetry::Symmetric>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Indirect, Symmetry::Antisymmetric>')
    call accel_program_add('operate.cu', 'operate_map<double, Indexing::Indirect, Symmetry::Symmetric>')
 
    call accel_program_add('pack.cu', 'pack<double>')
    call accel_program_add('pack.cu', 'unpack<double>')
    call accel_program_add('pack.cu', 'pack<complex<double>>')
    call accel_program_add('pack.cu', 'unpack<complex<double>>')
    call accel_program_add('pack.cu', 'pack<float>')
    call accel_program_add('pack.cu', 'unpack<float>')
    call accel_program_add('pack.cu', 'pack<complex<float>>')
    call accel_program_add('pack.cu', 'unpack<complex<float>>')
 
    call accel_program_add('pes.cu', 'zpes_flux_sph_integrate<double>')
    call accel_program_add('pes.cu', 'zpes_flux_sph_integxx<double>')
    call accel_program_add('pes.cu', 'zpes_flux_sph_integxx_2<double>')
 
    call accel_program_add('phase.cu', 'phase<double>')
    call accel_program_add('phase.cu', 'phase_hamiltonian<double>')
    call accel_program_add('phase.cu', 'update_phases<double, Dim::D1>')
    call accel_program_add('phase.cu', 'update_phases<double, Dim::D2>')
    call accel_program_add('phase.cu', 'update_phases<double, Dim::D3>')
 
    call accel_program_add('phase_spiral.cu', 'phase_spiral_apply<double>')
 
    call accel_program_add('pml.cu', 'pml_apply<double>')
    call accel_program_add('pml.cu', 'pml_apply_new<double>')
    call accel_program_add('pml.cu', 'pml_copy<double>')
    call accel_program_add('pml.cu', 'pml_update_conv<double>')
    call accel_program_add('pml.cu', 'pml_update_new<double>')
 
    call accel_program_add('points.cu', 'get_points<double>')
    call accel_program_add('points.cu', 'set_points<double>')
 
    call accel_program_add('precision_copy.cu', 'precision_copy<double, float>')
    call accel_program_add('precision_copy.cu', 'precision_copy<float, double>')
    call accel_program_add('precision_copy.cu', &
      'precision_copy<complex<double>, complex<float>>')
    call accel_program_add('precision_copy.cu', &
      'precision_copy<complex<float>, complex<double>>')
 
    call accel_program_add('projector.cu', 'projector_bra_force<double, double>')
    call accel_program_add('projector.cu', 'projector_bra_force_phase<double, double>')
    call accel_program_add('projector.cu', 'projector_ket<double>')
    call accel_program_add('projector.cu', 'projector_ket_phase<double, double>')
    call accel_program_add('projector.cu', 'projector_ket_phase_spiral<double, double>')
    call accel_program_add('projector.cu', 'dprojector_mix<double>')
    call accel_program_add('projector.cu', 'projector_r_vnl_bra<double, double>')
    call accel_program_add('projector.cu', 'projector_r_vnl_bra_phase<double, double>')
    call accel_program_add('projector.cu', 'projector_r_vnl_ket<double, double>')
    call accel_program_add('projector.cu', 'projector_r_vnl_ket_phase<double, double>')
    call accel_program_add('projector.cu', 'projector_bra_phase_spiral_gather_all<double>')
    call accel_program_add('projector.cu', 'projector_bra_force<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_bra_force_phase<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_ket<complex<double>>')
    call accel_program_add('projector.cu', 'projector_ket_phase<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_ket_phase_spiral<double, complex<double>>')
    call accel_program_add('projector.cu', 'zprojector_mix<double>')
    call accel_program_add('projector.cu', 'projector_r_vnl_bra<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_r_vnl_bra_phase<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_r_vnl_ket<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_r_vnl_ket_phase<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_bra_gather_all<double>')
    call accel_program_add('projector.cu', 'projector_bra_gather_all<complex<double>>')
    call accel_program_add('projector.cu', 'projector_bra_phase_gather_all<double>')
    call accel_program_add('projector.cu', 'projector_bra_scatter_all<double, double>')
    call accel_program_add('projector.cu', 'projector_bra_scatter_all<double, complex<double>>')
    call accel_program_add('projector.cu', 'projector_bra_scatter_complex_all<double, complex<double>>')
 
    call accel_program_add('projector_commutator.cu', 'projector_commutator_bra<double, double>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_bra_phase<double, double>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_ket<double, double>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_ket_phase<double, double>')
    call accel_program_add('projector_commutator.cu', 'dprojector_mix_commutator<double>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_bra<double, complex<double>>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_bra_phase<double, complex<double>>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_ket<double, complex<double>>')
    call accel_program_add('projector_commutator.cu', 'projector_commutator_ket_phase<double, complex<double>>')
    call accel_program_add('projector_commutator.cu', 'zprojector_mix_commutator<double>')
 
    call accel_program_add('split.cu', 'merge_complex<double>')
    call accel_program_add('split.cu', 'merge_complex<float>')
    call accel_program_add('split.cu', 'split_complex<double>')
    call accel_program_add('split.cu', 'split_complex<float>')
 
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<double, Conjugation::No>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<double, Conjugation::Yes>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<complex<double>, Conjugation::No>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<complex<double>, Conjugation::Yes>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<float, Conjugation::No>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<float, Conjugation::Yes>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<complex<float>, Conjugation::No>')
    call accel_program_add('symmetrize_batch.cu', &
      'symmetrizer_apply_single_batch<complex<float>, Conjugation::Yes>')
 
    call accel_program_add('uvw_to_xyz.cu', &
      'uvw_to_xyz<double, Dim::D1, double*, int, double*, int>')
    call accel_program_add('uvw_to_xyz.cu', &
      'uvw_to_xyz<double, Dim::D2, double*, int, double*, int, double*, int, double*, int>')
    call accel_program_add('uvw_to_xyz.cu', &
      'uvw_to_xyz<double, Dim::D3, double*, int, double*, int, double*, int, double*, int, double*, int, double*, int>')
 
    call accel_program_add('vpsi.cu', 'vpsi<double>')
    call accel_program_add('vpsi.cu', 'vpsi_complex<double>')
    call accel_program_add('vpsi.cu', 'vpsi_spinors<double>')
    call accel_program_add('vpsi.cu', 'vpsi_spinors_complex<double>')
 
    call accel_program_add('xc_dens_block.cu', 'xc_dens_apply_corrections<double>')
    call accel_program_add('xc_dens_block.cu', 'xc_dens_extract_block<double>')
 
    pop_sub(accel_kernel_list)
  end subroutine accel_kernel_list
 
  !------------------------------------------------------------
 
  subroutine accel_program_add(file_name, name)
    character(len=*), intent(in) :: file_name
    character(len=*), intent(in) :: name
 
    push_sub(accel_program_add)
 
    call cuda_program_register(accel%program_map, string_f_to_c(trim(file_name)), &
      string_f_to_c(trim(name)))
 
    pop_sub(accel_program_add)
  end subroutine accel_program_add
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_global_end()
    type(accel_kernel_t), pointer :: next_head
 
    push_sub(accel_kernel_global_end)
 
    do
      if (.not. associated(head)) exit
      next_head => head%next
      call accel_kernel_end(head)
      head => next_head
    end do
 
    if (accel_is_enabled()) then
      call cuda_program_map_end(accel%program_map)
    end if
 
    pop_sub(accel_kernel_global_end)
  end subroutine accel_kernel_global_end
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_build(this, file_name, kernel_name, flags)
    type(accel_kernel_t),        intent(inout) :: this
    character(len=*),            intent(in)    :: file_name
    character(len=*),            intent(in)    :: kernel_name
    character(len=*), optional,  intent(in)    :: flags
 
#ifdef HAVE_CUDA
    character(len=1000) :: all_flags
#endif
 
    push_sub(accel_kernel_build)
 
    call profiling_in("ACCEL_COMPILE", exclude = .true.)
 
#ifdef HAVE_CUDA
    all_flags = '-I'//trim(conf%share)//'/kernels/'//" "//trim(accel%debug_flag)
 
    if (present(flags)) then
      all_flags = trim(all_flags)//' '//trim(flags)
    end if
 
    call cuda_program_build(accel%program_map, accel%device%cuda_device, &
      string_f_to_c(trim(file_name)), string_f_to_c(trim(all_flags)))
 
    call cuda_program_get_kernel(this%cuda_kernel, accel%program_map, &
      string_f_to_c(trim(file_name)), string_f_to_c(trim(all_flags)), &
      string_f_to_c(trim(kernel_name)))
    call cuda_alloc_arg_array(this%arguments)
#endif
 
    this%initialized = .true.
    this%kernel_name = trim(kernel_name)
 
    call profiling_out("ACCEL_COMPILE")
 
    pop_sub(accel_kernel_build)
  end subroutine accel_kernel_build
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_end(this)
    type(accel_kernel_t), intent(inout) :: this
 
    push_sub(accel_kernel_end)
 
#ifdef HAVE_CUDA
    call cuda_free_arg_array(this%arguments)
    call cuda_release_kernel(this%cuda_kernel)
    ! modules are not released here, since they are not associated to a kernel
#endif
 
    this%initialized = .false.
 
    pop_sub(accel_kernel_end)
  end subroutine accel_kernel_end
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_start_call(this, file_name, kernel_name, flags)
    type(accel_kernel_t), target, intent(inout) :: this
    character(len=*),             intent(in)    :: file_name
    character(len=*),             intent(in)    :: kernel_name
    character(len=*), optional,   intent(in)    :: flags
 
    push_sub(accel_kernel_start_call)
 
    if (.not. this%initialized) then
      call accel_kernel_build(this, file_name, kernel_name, flags)
      this%next => head
      head => this
    end if
 
    pop_sub(accel_kernel_start_call)
  end subroutine accel_kernel_start_call
 
  !--------------------------------------------------------------
 
  integer(int64) pure function accel_global_memory_size() result(size)
 
    size = accel%global_memory_size
 
  end function accel_global_memory_size
 
  !--------------------------------------------------------------
 
  integer(int64) pure function accel_shared_memory_size() result(size)
 
    size = accel%shared_memory_size
 
  end function accel_shared_memory_size
  !--------------------------------------------------------------
 
  integer pure function accel_max_size_per_dim(dim) result(size)
    integer, intent(in) :: dim
 
    size = 0
#ifdef HAVE_CUDA
    size = 32768
    if (dim == 1) size = 2**30
#endif
  end function accel_max_size_per_dim
 
  ! ------------------------------------------------------
 
  subroutine accel_set_stream(stream_number)
    integer, intent(in) :: stream_number
 
    push_sub(accel_set_stream)
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cuda_set_stream(accel%cuda_stream, stream_number)
      ! The stream number selects which of the per-stream BLAS handles the
      ! wrappers use; see cublas_init() in src/basic/cublas.cc.
      call cublas_set_stream(accel%cuda_stream, stream_number)
#endif
    end if
 
    pop_sub(accel_set_stream)
  end subroutine accel_set_stream
 
  ! ------------------------------------------------------
 
  subroutine accel_get_stream(stream_number)
    integer, intent(inout) :: stream_number
 
    push_sub(accel_get_stream)
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cuda_get_stream(stream_number)
#endif
    end if
 
    pop_sub(accel_get_stream)
  end subroutine accel_get_stream
 
  ! ------------------------------------------------------
 
  subroutine accel_synchronize_all_streams()
    push_sub(accel_synchronize_all_streams)
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cuda_synchronize_all_streams()
#endif
    end if
 
    pop_sub(accel_synchronize_all_streams)
  end subroutine accel_synchronize_all_streams
 
  function daccel_get_pointer_with_offset(buffer, offset) result(buffer_offset)
    type(c_ptr), intent(in) :: buffer
    integer(int64), intent(in) :: offset
    type(c_ptr) :: buffer_offset
 
    push_sub(daccel_get_pointer_with_offset)
#ifdef HAVE_CUDA
    call cuda_get_pointer_with_offset(buffer, offset, buffer_offset)
#else
    ! this is needed to make the compiler happy for non-GPU compilations
    buffer_offset = buffer
#endif
    pop_sub(daccel_get_pointer_with_offset)
  end function daccel_get_pointer_with_offset
 
  function zaccel_get_pointer_with_offset(buffer, offset) result(buffer_offset)
    type(c_ptr), intent(in) :: buffer
    integer(int64), intent(in) :: offset
    type(c_ptr) :: buffer_offset
 
    push_sub(zaccel_get_pointer_with_offset)
#ifdef HAVE_CUDA
    call cuda_get_pointer_with_offset(buffer, 2_int64*offset, buffer_offset)
#else
    ! this is needed to make the compiler happy for non-GPU compilations
    buffer_offset = buffer
#endif
    pop_sub(zaccel_get_pointer_with_offset)
  end function zaccel_get_pointer_with_offset
 
  subroutine accel_clean_pointer(buffer)
    type(c_ptr), intent(in) :: buffer
 
    push_sub(accel_clean_pointer)
#ifdef HAVE_CUDA
    call cuda_clean_pointer(buffer)
#endif
    pop_sub(accel_clean_pointer)
  end subroutine accel_clean_pointer
 
  subroutine accel_get_unfolded_size(size, grid_size, thread_block_size)
    integer(int64), intent(in)  :: size
    integer(int64), intent(out) :: grid_size
    integer(int64), intent(out) :: thread_block_size
 
    push_sub(accel_get_unfolded_size)
#ifdef __HIP_PLATFORM_AMD__
    ! not benefitial for AMD chips
    grid_size = 1_int64
    thread_block_size = size
#else
    grid_size = size
    thread_block_size = accel%warp_size
#endif
    pop_sub(accel_get_unfolded_size)
  end subroutine accel_get_unfolded_size
 
#include "undef.F90"
#include "real.F90"
#include "accel_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "accel_inc.F90"
 
#include "undef.F90"
#include "real32.F90"
#include "accel_inc.F90"
 
#include "undef.F90"
#include "complex32.F90"
#include "accel_inc.F90"
 
#include "undef.F90"
#include "integer.F90"
#include "accel_inc.F90"
 
#include "undef.F90"
#include "integer8.F90"
#include "accel_inc.F90"
 
end module accel_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: