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_OPENCL) && defined(HAVE_CUDA)
#error "Cannot compile with OpenCL and Cuda support at the same time"
#endif
 
#if defined(HAVE_OPENCL) || defined(HAVE_CUDA)
#define HAVE_ACCEL 1
#endif
 
module accel_oct_m
  use alloc_cache_oct_m
#ifdef HAVE_OPENCL
  use cl
#endif
#if defined(HAVE_CLBLAS) || defined(HAVE_CLBLAST)
  use clblas_oct_m
#endif
  use cuda_oct_m
#ifdef HAVE_CLFFT
  use clfft
#endif
  use debug_oct_m
  use global_oct_m
  use iso_c_binding, only: c_size_t
  use, intrinsic :: iso_fortran_env
  use loct_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
 
  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_release_buffer,         &
    accel_buffer_is_allocated,    &
    accel_finish,                 &
    accel_set_kernel_arg,         &
    accel_max_workgroup_size,     &
    accel_kernel_workgroup_size,  &
    accel_kernel_run,             &
    accel_set_buffer_to_zero,     &
    accel_use_shared_mem,         &
    clblas_print_error,           &
    clfft_print_error,            &
    accel_local_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
 
#ifdef HAVE_OPENCL
  integer, public, parameter ::                 &
    ACCEL_MEM_READ_ONLY  = cl_mem_read_only,    &
    accel_mem_read_write = cl_mem_read_write,   &
    accel_mem_write_only = cl_mem_write_only
#else
  integer, public, parameter ::                 &
    ACCEL_MEM_READ_ONLY  = 0,                   &
    accel_mem_read_write = 1,                   &
    accel_mem_write_only = 2
#endif
 
  type accel_context_t
    ! Components are public by default
#ifdef HAVE_OPENCL
    type(cl_context) :: cl_context
#elif 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
#ifdef HAVE_OPENCL
    type(cl_device_id) :: cl_device
#elif 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
#ifdef HAVE_OPENCL
    type(cl_command_queue) :: command_queue
#endif
    type(c_ptr)            :: cublas_handle
    type(c_ptr)            :: cuda_stream
    type(c_ptr)            :: module_map
    integer                :: max_workgroup_size
    integer(int64)            :: local_memory_size
    integer(int64)            :: global_memory_size
    logical                :: enabled
    logical                :: allow_CPU_only
    logical                :: shared_mem
    logical                :: cuda_mpi
    integer                :: warp_size
    logical                :: 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
#ifdef HAVE_OPENCL
    type(cl_mem)           :: mem
#else
    type(c_ptr)            :: mem
#endif
    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_OPENCL
    type(cl_kernel)               :: kernel
#endif
#ifdef HAVE_CUDA
    type(c_ptr)                   :: cuda_kernel
    type(c_ptr)                   :: cuda_module
    type(c_ptr)                   :: arguments
#endif
    integer(int64)                   :: cuda_shared_mem
    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_copy
  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 :: dunpack
  type(accel_kernel_t), public, target, save :: zunpack
  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 :: dkernel_dot_matrix
  type(accel_kernel_t), public, target, save :: zkernel_dot_matrix
  type(accel_kernel_t), public, target, save :: zkernel_dot_matrix_spinors
  type(accel_kernel_t), public, target, save :: dkernel_batch_axpy
  type(accel_kernel_t), public, target, save :: zkernel_batch_axpy
  type(accel_kernel_t), public, target, save :: dkernel_ax_function_py
  type(accel_kernel_t), public, target, save :: zkernel_ax_function_py
  type(accel_kernel_t), public, target, save :: dkernel_batch_dotp
  type(accel_kernel_t), public, target, save :: zkernel_batch_dotp
  type(accel_kernel_t), public, target, save :: dzmul
  type(accel_kernel_t), public, target, save :: zzmul
  type(accel_kernel_t), public, target, save :: set_one
 
  ! kernels used locally
  type(accel_kernel_t), save :: set_zero
 
  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_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_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
  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_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
  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,   &
      accel_set_kernel_arg_local
  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
  end interface accel_get_device_pointer
 
 
  integer, parameter  ::      &
    opencl_gpu         = -1,  &
    opencl_cpu         = -2,  &
    opencl_accelerator = -3,  &
    opencl_default     = -4
 
 
  integer, parameter  ::      &
    cl_plat_invalid   = -1,   &
    cl_plat_amd       = -2,   &
    cl_plat_nvidia    = -3,   &
    cl_plat_ati       = -4,   &
    cl_plat_intel     = -5
 
  ! a "convenience" public variable
  integer, public :: cl_status
 
  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, iplatform
#ifdef HAVE_OPENCL
    integer  :: device_type
    integer :: cl_status, idev
    integer  :: ndevices, ret_devices, nplatforms, iplat
    character(len=256) :: device_name
    type(cl_platform_id) :: platform_id
    type(cl_program) :: prog
    type(cl_platform_id), allocatable :: allplatforms(:)
    type(cl_device_id), allocatable :: alldevices(:)
#endif
#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 OpenCL or 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 OpenCL or Cuda support.'
      call messages_fatal(1)
    end if
#endif
 
    if (.not. accel_is_enabled()) then
      pop_sub(accel_init)
      return
    end if
 
    !%Variable AccelPlatform
    !%Type integer
    !%Default 0
    !%Section Execution::Accel
    !%Description
    !% This variable selects the OpenCL platform that Octopus will
    !% use. You can give an explicit platform number or use one of
    !% the options that select a particular vendor
    !% implementation. Platform 0 is used by default.
    !%
    !% This variable has no effect for CUDA.
    !%Option amd -2
    !% Use the AMD OpenCL platform.
    !%Option nvidia -3
    !% Use the Nvidia OpenCL platform.
    !%Option ati -4
    !% Use the ATI (old AMD) OpenCL platform.
    !%Option intel -5
    !% Use the Intel OpenCL platform.
    !%End
    call parse_variable(namespace, 'AccelPlatform', 0, iplatform)
 
    call messages_obsolete_variable(namespace, 'OpenCLPlatform', 'AccelPlatform')
 
    !%Variable AccelDevice
    !%Type integer
    !%Default gpu
    !%Section Execution::Accel
    !%Description
    !% This variable selects the OpenCL or CUDA accelerator device
    !% that Octopus will use. You can specify one of the options below
    !% or a numerical id to select a specific device.
    !%
    !% Values >= 0 select the device to be used. In case of MPI enabled runs
    !% devices are distributed in a round robin fashion, starting at this value.
    !%Option gpu -1
    !% If available, Octopus will use a GPU.
    !%Option cpu -2
    !% If available, Octopus will use a CPU (only for OpenCL).
    !%Option accelerator -3
    !% If available, Octopus will use an accelerator (only for OpenCL).
    !%Option accel_default -4
    !% Octopus will use the default device specified by the implementation.
    !% implementation.
    !%End
    call parse_variable(namespace, 'AccelDevice', opencl_gpu, idevice)
 
    call messages_obsolete_variable(namespace, 'OpenCLDevice', 'AccelDevice')
 
    if (idevice < opencl_default) 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
 
    ! no shared mem support in our cuda interface (for the moment)
    accel%shared_mem = .true.
 
    call cublas_init(accel%cublas_handle, accel%cuda_stream)
#endif
 
#ifdef HAVE_OPENCL
    call profiling_in('CL_INIT')
 
    call clgetplatformids(nplatforms, cl_status)
    if (cl_status /= cl_success) call opencl_print_error(cl_status, "GetPlatformIDs")
 
    safe_allocate(allplatforms(1:nplatforms))
 
    call clgetplatformids(allplatforms, iplat, cl_status)
    if (cl_status /= cl_success) call opencl_print_error(cl_status, "GetPlatformIDs")
 
    call messages_write('Info: Available CL platforms: ')
    call messages_write(nplatforms)
    call messages_info()
 
    do iplat = 1, nplatforms
 
      call clgetplatforminfo(allplatforms(iplat), cl_platform_name, device_name, cl_status)
 
      if (iplatform < 0) then
        if (iplatform == get_platform_id(device_name)) iplatform = iplat - 1
      end if
 
      if (iplatform == iplat - 1) then
        call messages_write('    * Platform ')
      else
        call messages_write('      Platform ')
      end if
 
      call messages_write(iplat - 1)
      call messages_write(' : '//device_name)
      call clgetplatforminfo(allplatforms(iplat), cl_platform_version, device_name, cl_status)
      call messages_write(' ('//trim(device_name)//')')
      call messages_info()
    end do
 
    call messages_info()
 
    if (iplatform >= nplatforms .or. iplatform < 0) then
      call messages_write('Requested CL platform does not exist')
      if (iplatform > 0) then
        call messages_write('(platform = ')
        call messages_write(iplatform)
        call messages_write(').')
      end if
      call messages_fatal()
    end if
 
    platform_id = allplatforms(iplatform + 1)
 
    safe_deallocate_a(allplatforms)
 
    call clgetdeviceids(platform_id, cl_device_type_all, ndevices, cl_status)
 
    call messages_write('Info: Available CL devices: ')
    call messages_write(ndevices)
    call messages_info()
 
    safe_allocate(alldevices(1:ndevices))
 
    ! list all devices
 
    call clgetdeviceids(platform_id, cl_device_type_all, alldevices, ret_devices, cl_status)
 
    do idev = 1, ndevices
      call messages_write('      Device ')
      call messages_write(idev - 1)
      call clgetdeviceinfo(alldevices(idev), cl_device_name, device_name, cl_status)
      call messages_write(' : '//device_name)
      call messages_info()
    end do
 
    select case (idevice)
    case (opencl_gpu)
      device_type = cl_device_type_gpu
    case (opencl_cpu)
      device_type = cl_device_type_cpu
    case (opencl_accelerator)
      device_type = cl_device_type_accelerator
    case (opencl_default)
      device_type = cl_device_type_default
    case default
      device_type = cl_device_type_all
    end select
 
    ! now get a list of the selected type
    call clgetdeviceids(platform_id, device_type, alldevices, ret_devices, cl_status)
 
    if (ret_devices < 1) then
      ! we didnt find a device of the selected type, we ask for the default device
      call clgetdeviceids(platform_id, cl_device_type_default, alldevices, ret_devices, cl_status)
 
      if (ret_devices < 1) then
        ! if this does not work, we ask for all devices
        call clgetdeviceids(platform_id, cl_device_type_all, alldevices, ret_devices, cl_status)
      end if
 
      if (ret_devices < 1) then
        call messages_write('Cannot find an OpenCL device')
        call messages_fatal()
      end if
    end if
 
    ! the number of devices can be smaller
    ndevices = ret_devices
 
    if (idevice < 0) then
      if (base_grp%size > 1) then
        ! with MPI we have to select the device so multiple GPUs in one
        ! node are correctly distributed
        call select_device(idevice)
      else
        idevice = 0
      end if
    end if
 
    if (idevice >= ndevices) then
      call messages_write('Requested CL device does not exist (device = ')
      call messages_write(idevice)
      call messages_write(', platform = ')
      call messages_write(iplatform)
      call messages_write(').')
      call messages_fatal()
    end if
 
    accel%device%cl_device = alldevices(idevice + 1)
 
    ! create the context
    accel%context%cl_context = clcreatecontext(platform_id, accel%device%cl_device, cl_status)
    if (cl_status /= cl_success) call opencl_print_error(cl_status, "CreateContext")
 
    safe_deallocate_a(alldevices)
 
    accel%command_queue = clcreatecommandqueue(accel%context%cl_context, accel%device%cl_device, &
      cl_queue_profiling_enable, cl_status)
    if (cl_status /= cl_success) call opencl_print_error(cl_status, "CreateCommandQueue")
 
    call clgetdeviceinfo(accel%device%cl_device, cl_device_type, device_type, cl_status)
 
    select case (device_type)
    case (cl_device_type_gpu)
      accel%shared_mem = .true.
    case (cl_device_type_cpu, cl_device_type_accelerator)
      accel%shared_mem = .false.
    case default
      accel%shared_mem = .false.
    end select
 
#ifdef HAVE_CLBLAS
    call clblassetup(cl_status)
    if (cl_status /= clblassuccess) call clblas_print_error(cl_status, 'clblasSetup')
#endif
 
#ifdef HAVE_CLFFT
    call clfftsetup(cl_status)
    if (cl_status /= clfft_success) call clfft_print_error(cl_status, 'clfftSetup')
#endif
 
    call profiling_out('CL_INIT')
#endif
 
    ! Get some device information that we will need later
 
    ! total memory
#ifdef HAVE_OPENCL
    call clgetdeviceinfo(accel%device%cl_device, cl_device_global_mem_size, accel%global_memory_size, cl_status)
    call clgetdeviceinfo(accel%device%cl_device, cl_device_local_mem_size, accel%local_memory_size, cl_status)
    call clgetdeviceinfo(accel%device%cl_device, cl_device_max_work_group_size, accel%max_workgroup_size, cl_status)
    accel%warp_size = 1
#endif
#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%local_memory_size)
    call cuda_device_max_threads_per_block(accel%device%cuda_device, accel%max_workgroup_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 (mpi_grp_is_root(base_grp)) 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"
#elif defined(HAVE_OPENCL)
    accel%debug_flag = "-g"
#endif
 
    call accel_kernel_start_call(set_zero, 'set_zero.cl', "set_zero")
    call accel_kernel_start_call(set_one, 'set_one.cl', "set_one")
    call accel_kernel_start_call(kernel_vpsi, 'vpsi.cl', "vpsi")
    call accel_kernel_start_call(kernel_vpsi_complex, 'vpsi.cl', "vpsi_complex")
    call accel_kernel_start_call(kernel_vpsi_spinors, 'vpsi.cl', "vpsi_spinors")
    call accel_kernel_start_call(kernel_vpsi_spinors_complex, 'vpsi.cl', "vpsi_spinors_complex")
    call accel_kernel_start_call(kernel_daxpy, 'axpy.cl', "daxpy", flags = '-DRTYPE_DOUBLE')
    call accel_kernel_start_call(kernel_zaxpy, 'axpy.cl', "zaxpy", flags = '-DRTYPE_COMPLEX')
    call accel_kernel_start_call(dkernel_batch_axpy, 'axpy.cl', "dbatch_axpy_function", &
      flags = ' -DRTYPE_DOUBLE')
    call accel_kernel_start_call(zkernel_batch_axpy, 'axpy.cl', "zbatch_axpy_function", &
      flags = '-DRTYPE_COMPLEX')
    call accel_kernel_start_call(dkernel_ax_function_py, 'axpy.cl', "dbatch_ax_function_py", &
      flags = '-DRTYPE_DOUBLE')
    call accel_kernel_start_call(zkernel_ax_function_py, 'axpy.cl', "zbatch_ax_function_py", &
      flags = '-DRTYPE_COMPLEX')
    call accel_kernel_start_call(dkernel_batch_dotp, 'mesh_batch_single.cl', "dbatch_mf_dotp")
    call accel_kernel_start_call(zkernel_batch_dotp, 'mesh_batch_single.cl', "zbatch_mf_dotp")
    call accel_kernel_start_call(dpack, 'pack.cl', "dpack")
    call accel_kernel_start_call(zpack, 'pack.cl', "zpack")
    call accel_kernel_start_call(dunpack, 'pack.cl', "dunpack")
    call accel_kernel_start_call(zunpack, 'pack.cl', "zunpack")
    call accel_kernel_start_call(kernel_copy, 'copy.cl', "copy")
    call accel_kernel_start_call(kernel_copy_complex_to_real, 'copy.cl', "copy_complex_to_real")
    call accel_kernel_start_call(kernel_copy_real_to_complex, 'copy.cl', "copy_real_to_complex")
    call accel_kernel_start_call(kernel_ghost_reorder, 'ghost.cl', "ghost_reorder")
    call accel_kernel_start_call(kernel_density_real, 'density.cl', "density_real")
    call accel_kernel_start_call(kernel_density_complex, 'density.cl', "density_complex")
    call accel_kernel_start_call(kernel_density_spinors, 'density.cl', "density_spinors")
    call accel_kernel_start_call(kernel_phase, 'phase.cl', "phase")
    call accel_kernel_start_call(dkernel_dot_matrix, 'mesh_batch.cl', "ddot_matrix")
    call accel_kernel_start_call(zkernel_dot_matrix, 'mesh_batch.cl', "zdot_matrix")
    call accel_kernel_start_call(zkernel_dot_matrix_spinors, 'mesh_batch.cl', "zdot_matrix_spinors")
 
 
    call accel_kernel_start_call(dzmul, 'mul.cl', "dzmul", flags = '-DRTYPE_DOUBLE')
    call accel_kernel_start_call(zzmul, 'mul.cl', "zzmul", flags = '-DRTYPE_COMPLEX')
 
    !%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 opencl_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 logical
    !%Section Execution::Accel
    !%Description
    !% Initialize new GPU buffers to zero on creation (use only for debugging, as it has a performance impact!).
    !%End
    call parse_variable(namespace, 'InitializeGPUBuffers', .false., accel%initialize_buffers)
 
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(accel_init)
 
  contains
 
#if defined(HAVE_OPENCL)
    subroutine select_device(idevice)
      integer, intent(inout) :: idevice
      integer :: irank
      character(len=256) :: device_name
 
      push_sub(accel_init.select_device)
 
      idevice = mod(base_grp%rank, ndevices)
 
      call base_grp%barrier()
      call messages_write('Info: CL device distribution:')
      call messages_info()
      do irank = 0, base_grp%size - 1
        if (irank == base_grp%rank) then
          call clgetdeviceinfo(alldevices(idevice + 1), cl_device_name, device_name, cl_status)
          call messages_write('      MPI node ')
          call messages_write(base_grp%rank)
          call messages_write(' -> CL device ')
          call messages_write(idevice)
          call messages_write(' : '//device_name)
          call messages_info(all_nodes = .true.)
        end if
        call base_grp%barrier()
      end do
 
      pop_sub(accel_init.select_device)
    end subroutine select_device
#endif
 
    subroutine device_info()
#ifdef HAVE_OPENCL
      integer(int64) :: val
#endif
#ifdef HAVE_CUDA
      integer :: version
#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_OPENCL
      call messages_write('      Framework              : OpenCL')
#endif
#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_OPENCL
      call clgetdeviceinfo(accel%device%cl_device, cl_device_type, val, cl_status)
      call messages_write('      Device type            :')
      select case (int(val, int32))
      case (cl_device_type_gpu)
        call messages_write(' GPU')
      case (cl_device_type_cpu)
        call messages_write(' CPU')
      case (cl_device_type_accelerator)
        call messages_write(' accelerator')
      end select
      call messages_new_line()
 
      call clgetdeviceinfo(accel%device%cl_device, cl_device_vendor, val_str, cl_status)
      call messages_write('      Device vendor          : '//trim(val_str))
      call messages_new_line()
#endif
 
#ifdef HAVE_OPENCL
      call clgetdeviceinfo(accel%device%cl_device, cl_device_name, val_str, cl_status)
#endif
#ifdef HAVE_CUDA
      call cuda_device_name(accel%device%cuda_device, 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_OPENCL
      call clgetdeviceinfo(accel%device%cl_device, cl_driver_version, val_str, cl_status)
      call messages_write('      Driver version         : '//trim(val_str))
#endif
#ifdef HAVE_CUDA
      call cuda_driver_version(version)
      call messages_write('      Driver version         : ')
      call messages_write(version)
#endif
      call messages_new_line()
 
 
#ifdef HAVE_OPENCL
      call clgetdeviceinfo(accel%device%cl_device, cl_device_max_compute_units, val, cl_status)
      call messages_write('      Compute units          :')
      call messages_write(val)
      call messages_new_line()
 
      call clgetdeviceinfo(accel%device%cl_device, cl_device_max_clock_frequency, val, cl_status)
      call messages_write('      Clock frequency        :')
      call messages_write(val)
      call messages_write(' GHz')
      call messages_new_line()
#endif
 
      call messages_write('      Device memory          :')
      call messages_write(accel%global_memory_size, units=unit_megabytes)
      call messages_new_line()
 
      call messages_write('      Local/shared memory    :')
      call messages_write(accel%local_memory_size, units=unit_kilobytes)
      call messages_new_line()
 
 
#ifdef HAVE_OPENCL
      call clgetdeviceinfo(accel%device%cl_device, cl_device_max_mem_alloc_size, val, cl_status)
      call messages_write('      Max alloc size         :')
      call messages_write(val, units = unit_megabytes)
      call messages_new_line()
 
      call clgetdeviceinfo(accel%device%cl_device, cl_device_global_mem_cache_size, val, cl_status)
      call messages_write('      Device cache           :')
      call messages_write(val, units = unit_kilobytes)
      call messages_new_line()
 
      call clgetdeviceinfo(accel%device%cl_device, cl_device_max_constant_buffer_size, val, cl_status)
      call messages_write('      Constant memory        :')
      call messages_write(val, units = unit_kilobytes)
      call messages_new_line()
#endif
 
      call messages_write('      Max. group/block size  :')
      call messages_write(accel%max_workgroup_size)
      call messages_new_line()
 
 
#ifdef HAVE_OPENCL
      call messages_write('      Extension cl_khr_fp64  :')
      call messages_write(f90_cl_device_has_extension(accel%device%cl_device, "cl_khr_fp64"))
      call messages_new_line()
 
      call messages_write('      Extension cl_amd_fp64  :')
      call messages_write(f90_cl_device_has_extension(accel%device%cl_device, "cl_amd_fp64"))
      call messages_new_line()
 
      call messages_write('      Extension cl_khr_int64_base_atomics  :')
      call messages_write(f90_cl_device_has_extension(accel%device%cl_device, "cl_khr_int64_base_atomics"))
      call messages_new_line()
 
#endif
 
      call messages_info()
 
 
      pop_sub(accel_init.device_info)
    end subroutine device_info
 
  end subroutine accel_init
 
  ! ------------------------------------------
#ifdef HAVE_OPENCL
  integer function get_platform_id(platform_name) result(platform_id)
    character(len=*), intent(in) :: platform_name
 
    platform_id = cl_plat_invalid
    if (index(platform_name, 'AMD') > 0)    platform_id = cl_plat_amd
    if (index(platform_name, 'ATI') > 0)    platform_id = cl_plat_ati
    if (index(platform_name, 'NVIDIA') > 0) platform_id = cl_plat_nvidia
    if (index(platform_name, 'Intel') > 0)  platform_id = cl_plat_intel
  end function get_platform_id
#endif
  ! ------------------------------------------
 
  subroutine accel_end(namespace)
    type(namespace_t), intent(in) :: namespace
 
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
    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
 
      do
        call alloc_cache_get(memcache, alloc_cache_any_size, found, tmp%mem)
        if (.not. found) exit
 
#ifdef HAVE_OPENCL
        call clreleasememobject(tmp%mem, ierr)
        if (ierr /= cl_success) call opencl_print_error(ierr, "clReleaseMemObject")
#endif
#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()
 
#ifdef HAVE_CLBLAS
    call clblasteardown()
#endif
 
#ifdef HAVE_CLFFT
    call clfftteardown()
#endif
 
    if (accel_is_enabled()) then
#ifdef HAVE_CUDA
      call cublas_end(accel%cublas_handle)
      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
 
#ifdef HAVE_OPENCL
      call clreleasecommandqueue(accel%command_queue, ierr)
 
      if (ierr /= cl_success) call opencl_print_error(ierr, "ReleaseCommandQueue")
      call clreleasecontext(accel%context%cl_context, cl_status)
#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
 
  ! ------------------------------------------
 
  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_workgroup_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)
    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
 
    call accel_create_buffer_8(this, flags, type, int(size, int64), set_zero)
  end subroutine accel_create_buffer_4
 
  ! ------------------------------------------
 
  subroutine accel_create_buffer_8(this, flags, type, size, set_zero)
    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
 
    integer(int64) :: fsize
    logical    :: found
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
 
    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_OPENCL
        this%mem = clcreatebuffer(accel%context%cl_context, flags, fsize, ierr)
        if (ierr /= cl_success) call opencl_print_error(ierr, "clCreateBuffer")
#endif
#ifdef HAVE_CUDA
        call cuda_mem_alloc(this%mem, fsize)
#endif
      end if
 
      buffer_alloc_count = buffer_alloc_count + 1
      allocated_mem = allocated_mem + fsize
 
    end if
 
    if(optional_default(set_zero,  accel%initialize_buffers)) then
      call accel_set_buffer_to_zero_i8(this, type, size)
    endif
 
    pop_sub(accel_create_buffer_8)
  end subroutine accel_create_buffer_8
 
  ! ------------------------------------------
 
  subroutine accel_release_buffer(this)
    type(accel_mem_t), intent(inout) :: this
 
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
    logical :: put
    integer(int64) :: fsize
 
    push_sub(accel_release_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_OPENCL
        call clreleasememobject(this%mem, ierr)
        if (ierr /= cl_success) call opencl_print_error(ierr, "clReleaseMemObject")
#endif
#ifdef HAVE_CUDA
        call cuda_mem_free(this%mem)
#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_release_buffer)
  end subroutine accel_release_buffer
 
  ! ------------------------------------------
 
  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()
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
 
    ! no push_sub, called too frequently
 
    if (accel_is_enabled()) then
#ifdef HAVE_OPENCL
      call clfinish(accel%command_queue, ierr)
      if (ierr /= cl_success) call opencl_print_error(ierr, 'clFinish')
#endif
#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
 
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
 
    assert(accel_buffer_is_allocated(buffer))
 
    ! no push_sub, called too frequently
#ifdef HAVE_OPENCL
    call clsetkernelarg(kernel%kernel, narg, buffer%mem, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "clSetKernelArg_buf")
#endif
 
#ifdef HAVE_CUDA
    call cuda_kernel_set_arg_buffer(kernel%arguments, buffer%mem, narg)
#endif
 
  end subroutine accel_set_kernel_arg_buffer
 
  ! ------------------------------------------
 
  subroutine accel_set_kernel_arg_local(kernel, narg, type, size)
    type(accel_kernel_t), intent(inout) :: kernel
    integer,              intent(in)    :: narg
    type(type_t),         intent(in)    :: type
    integer,              intent(in)    :: size
 
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
    integer(int64) :: size_in_bytes
 
    push_sub(accel_set_kernel_arg_local)
 
 
    size_in_bytes = int(size, int64)*types_get_size(type)
 
    if (size_in_bytes > accel%local_memory_size) then
      write(message(1), '(a,f12.6,a)') "CL Error: requested local memory: ", real(size_in_bytes, real64) /1024.0, " Kb"
      write(message(2), '(a,f12.6,a)') "          available local memory: ", real(accel%local_memory_size, real64) /1024.0, " Kb"
      call messages_fatal(2)
    else if (size_in_bytes <= 0) then
      write(message(1), '(a,i10)') "CL Error: invalid local memory size: ", size_in_bytes
      call messages_fatal(1)
    end if
 
#ifdef HAVE_CUDA
    kernel%cuda_shared_mem = size_in_bytes
#endif
 
#ifdef HAVE_OPENCL
    call clsetkernelarglocal(kernel%kernel, narg, size_in_bytes, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "set_kernel_arg_local")
#endif
 
    pop_sub(accel_set_kernel_arg_local)
  end subroutine accel_set_kernel_arg_local
 
  ! ------------------------------------------
 
  subroutine accel_kernel_run_8(kernel, globalsizes, localsizes)
    type(accel_kernel_t), intent(inout) :: kernel
    integer(int64),          intent(in)    :: globalsizes(:)
    integer(int64),          intent(in)    :: localsizes(:)
 
    integer :: dim
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
    integer(int64) :: gsizes(1:3)
    integer(int64) :: lsizes(1:3)
 
    ! no push_sub, called too frequently
 
    ! cuda needs all dimensions
    gsizes = 1
    lsizes = 1
 
    dim = ubound(globalsizes, dim=1)
 
    assert(dim == ubound(localsizes, dim=1))
 
    ! if one size is zero, there is nothing to do
    if (any(globalsizes == 0)) return
 
    assert(all(localsizes > 0))
    assert(all(localsizes <= accel_max_workgroup_size()))
    assert(all(mod(globalsizes, localsizes) == 0))
 
    gsizes(1:dim) = globalsizes(1:dim)
    lsizes(1:dim) = localsizes(1:dim)
 
#ifdef HAVE_OPENCL
    call clenqueuendrangekernel(accel%command_queue, kernel%kernel, gsizes(1:dim), lsizes(1:dim), ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "EnqueueNDRangeKernel")
#endif
 
#ifdef HAVE_CUDA
    ! Maximum dimension of a block
    if (any(lsizes(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: ", lsizes(1), " <= ", accel%max_block_dim(1)
      write(message(4), "(A, I8, A, I8)") "Dim 2: ", lsizes(2), " <= ", accel%max_block_dim(2)
      write(message(5), "(A, I8, A, I8)") "Dim 3: ", lsizes(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(lsizes) > accel_max_workgroup_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(lsizes), " <= ", accel_max_workgroup_size()
      message(4) = "This is an internal error, please contact the developers."
      call messages_fatal(4)
    end if
 
    gsizes(1:3) = gsizes(1:3)/lsizes(1:3)
 
    ! 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
 
    call cuda_launch_kernel(kernel%cuda_kernel, gsizes(1), lsizes(1), kernel%cuda_shared_mem, kernel%arguments)
 
    kernel%cuda_shared_mem = 0
#endif
 
  end subroutine accel_kernel_run_8
 
  ! -----------------------------------------------
 
  subroutine accel_kernel_run_4(kernel, globalsizes, localsizes)
    type(accel_kernel_t), intent(inout) :: kernel
    integer,              intent(in)    :: globalsizes(:)
    integer,              intent(in)    :: localsizes(:)
 
    call accel_kernel_run_8(kernel, int(globalsizes, int64), int(localsizes, int64))
 
  end subroutine accel_kernel_run_4
 
  ! -----------------------------------------------
 
  integer pure function accel_max_workgroup_size() result(max_workgroup_size)
    max_workgroup_size = accel%max_workgroup_size
  end function accel_max_workgroup_size
 
  ! -----------------------------------------------
 
  integer function accel_kernel_workgroup_size(kernel) result(workgroup_size)
    type(accel_kernel_t), intent(inout) :: kernel
 
#ifdef HAVE_OPENCL
    integer(int64) :: workgroup_size8
    integer :: ierr
#endif
#ifdef HAVE_CUDA
    integer :: max_workgroup_size
#endif
 
    workgroup_size = 0
 
#ifdef HAVE_OPENCL
    call clgetkernelworkgroupinfo(kernel%kernel, accel%device%cl_device, cl_kernel_work_group_size, workgroup_size8, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "EnqueueNDRangeKernel")
    workgroup_size = workgroup_size8
#endif
 
#ifdef HAVE_CUDA
    call cuda_kernel_max_threads_per_block(kernel%cuda_kernel, max_workgroup_size)
    if (debug%info .and. max_workgroup_size /= accel%max_workgroup_size) then
      write(message(1), "(A, I5, A)") "A kernel can use only less threads per block (", workgroup_size, ")", &
        "than available on the device (", accel%max_workgroup_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
    workgroup_size = 256
    ! make sure we do not use more threads per block than available for this kernel
    workgroup_size = min(workgroup_size, max_workgroup_size)
#endif
 
  end function accel_kernel_workgroup_size
 
  ! -----------------------------------------------
 
#ifdef HAVE_OPENCL
  subroutine opencl_build_program(prog, filename, flags)
    type(cl_program),           intent(inout) :: prog
    character(len=*),           intent(in)    :: filename
    character(len=*), optional, intent(in)    :: flags
 
    character(len = 1000) :: string
    character(len = 256) :: share_string
    integer :: ierr, ierrlog, iunit, irec, newlen
 
    push_sub(opencl_build_program)
 
    string = '#include "'//trim(filename)//'"'
 
    call messages_write("Building CL program '"//trim(filename)//"'.")
    call messages_info(debug_only=.true.)
 
    prog = clcreateprogramwithsource(accel%context%cl_context, trim(string), ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "clCreateProgramWithSource")
 
    ! build the compilation flags
    string='-w'
    ! full optimization
    string=trim(string)//' -cl-denorms-are-zero'
    ! The following flag gives an error with the Xeon Phi
    !    string=trim(string)//' -cl-strict-aliasing'
    string=trim(string)//' -cl-mad-enable'
    string=trim(string)//' -cl-unsafe-math-optimizations'
    string=trim(string)//' -cl-finite-math-only'
    string=trim(string)//' -cl-fast-relaxed-math'
 
    share_string='-I'//trim(conf%share)//'/opencl/'
 
    if (f90_cl_device_has_extension(accel%device%cl_device, "cl_khr_fp64")) then
      string = trim(string)//' -DEXT_KHR_FP64'
    else if (f90_cl_device_has_extension(accel%device%cl_device, "cl_amd_fp64")) then
      string = trim(string)//' -DEXT_AMD_FP64'
    else
      call messages_write('Octopus requires an OpenCL device with double-precision support.')
      call messages_fatal()
    end if
 
    if (accel_use_shared_mem()) then
      string = trim(string)//' -DSHARED_MEM'
    end if
 
    if (present(flags)) then
      string = trim(string)//' '//trim(flags)
    end if
 
    call messages_write("Debug info: compilation flags '"//trim(string), new_line = .true.)
    call messages_write('  '//trim(share_string)//"'.")
    call messages_info(debug_only=.true.)
 
    string = trim(string)//' '//trim(share_string)
 
    call clbuildprogram(prog, trim(string), ierr)
 
    if(ierr /= cl_success) then
      call clgetprogrambuildinfo(prog, accel%device%cl_device, cl_program_build_log, string, ierrlog)
      if (ierrlog /= cl_success) call opencl_print_error(ierrlog, "clGetProgramBuildInfo")
 
      ! CL_PROGRAM_BUILD_LOG seems to have a useless '\n' in it
      newlen = scan(string, achar(010), back = .true.) - 1
      if (newlen >= 0) string = string(1:newlen)
 
      if (len(trim(string)) > 0) write(stderr, '(a)') trim(string)
 
      call opencl_print_error(ierr, "clBuildProgram")
    end if
 
    pop_sub(opencl_build_program)
  end subroutine opencl_build_program
#endif
 
  ! -----------------------------------------------
#ifdef HAVE_OPENCL
  subroutine opencl_release_program(prog)
    type(cl_program),    intent(inout) :: prog
 
    integer :: ierr
 
    push_sub(opencl_release_program)
 
    call clreleaseprogram(prog, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "clReleaseProgram")
 
    pop_sub(opencl_release_program)
  end subroutine opencl_release_program
#endif
 
  ! -----------------------------------------------
 
#ifdef HAVE_OPENCL
  subroutine opencl_release_kernel(prog)
    type(cl_kernel),      intent(inout) :: prog
 
    integer :: ierr
 
    push_sub(opencl_release_kernel)
 
#ifdef HAVE_OPENCL
    call clreleasekernel(prog, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "clReleaseKernel")
#endif
 
    pop_sub(opencl_release_kernel)
  end subroutine opencl_release_kernel
#endif
 
#ifdef HAVE_OPENCL
  ! -----------------------------------------------
  subroutine opencl_create_kernel(kernel, prog, name)
    type(cl_kernel),  intent(inout) :: kernel
    type(cl_program), intent(inout) :: prog
    character(len=*), intent(in)    :: name
 
    integer :: ierr
 
    push_sub(opencl_create_kernel)
    call profiling_in("CL_BUILD_KERNEL", exclude = .true.)
 
#ifdef HAVE_OPENCL
    kernel = clcreatekernel(prog, name, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "clCreateKernel")
#endif
 
    call profiling_out("CL_BUILD_KERNEL")
    pop_sub(opencl_create_kernel)
  end subroutine opencl_create_kernel
#endif
 
  ! ------------------------------------------------
#ifdef HAVE_OPENCL
  subroutine opencl_print_error(ierr, name)
    integer,          intent(in) :: ierr
    character(len=*), intent(in) :: name
 
    character(len=40) :: errcode
 
    push_sub(opencl_print_error)
 
    select case (ierr)
    case (cl_success); errcode = 'CL_SUCCESS '
    case (cl_device_not_found); errcode = 'CL_DEVICE_NOT_FOUND '
    case (cl_device_not_available); errcode = 'CL_DEVICE_NOT_AVAILABLE '
    case (cl_compiler_not_available); errcode = 'CL_COMPILER_NOT_AVAILABLE '
    case (cl_mem_object_allocation_failure); errcode = 'CL_MEM_OBJECT_ALLOCATION_FAILURE '
    case (cl_out_of_resources); errcode = 'CL_OUT_OF_RESOURCES '
    case (cl_out_of_host_memory); errcode = 'CL_OUT_OF_HOST_MEMORY '
    case (cl_profiling_info_not_available); errcode = 'CL_PROFILING_INFO_NOT_AVAILABLE '
    case (cl_mem_copy_overlap); errcode = 'CL_MEM_COPY_OVERLAP '
    case (cl_image_format_mismatch); errcode = 'CL_IMAGE_FORMAT_MISMATCH '
    case (cl_image_format_not_supported); errcode = 'CL_IMAGE_FORMAT_NOT_SUPPORTED '
    case (cl_build_program_failure); errcode = 'CL_BUILD_PROGRAM_FAILURE '
    case (cl_map_failure); errcode = 'CL_MAP_FAILURE '
    case (cl_invalid_value); errcode = 'CL_INVALID_VALUE '
    case (cl_invalid_device_type); errcode = 'CL_INVALID_DEVICE_TYPE '
    case (cl_invalid_platform); errcode = 'CL_INVALID_PLATFORM '
    case (cl_invalid_device); errcode = 'CL_INVALID_DEVICE '
    case (cl_invalid_context); errcode = 'CL_INVALID_CONTEXT '
    case (cl_invalid_queue_properties); errcode = 'CL_INVALID_QUEUE_PROPERTIES '
    case (cl_invalid_command_queue); errcode = 'CL_INVALID_COMMAND_QUEUE '
    case (cl_invalid_host_ptr); errcode = 'CL_INVALID_HOST_PTR '
    case (cl_invalid_mem_object); errcode = 'CL_INVALID_MEM_OBJECT '
    case (cl_invalid_image_format_descriptor); errcode = 'CL_INVALID_IMAGE_FORMAT_DESCRIPTOR '
    case (cl_invalid_image_size); errcode = 'CL_INVALID_IMAGE_SIZE '
    case (cl_invalid_sampler); errcode = 'CL_INVALID_SAMPLER '
    case (cl_invalid_binary); errcode = 'CL_INVALID_BINARY '
    case (cl_invalid_build_options); errcode = 'CL_INVALID_BUILD_OPTIONS '
    case (cl_invalid_program); errcode = 'CL_INVALID_PROGRAM '
    case (cl_invalid_program_executable); errcode = 'CL_INVALID_PROGRAM_EXECUTABLE '
    case (cl_invalid_kernel_name); errcode = 'CL_INVALID_KERNEL_NAME '
    case (cl_invalid_kernel_definition); errcode = 'CL_INVALID_KERNEL_DEFINITION '
    case (cl_invalid_kernel); errcode = 'CL_INVALID_KERNEL '
    case (cl_invalid_arg_index); errcode = 'CL_INVALID_ARG_INDEX '
    case (cl_invalid_arg_value); errcode = 'CL_INVALID_ARG_VALUE '
    case (cl_invalid_arg_size); errcode = 'CL_INVALID_ARG_SIZE '
    case (cl_invalid_kernel_args); errcode = 'CL_INVALID_KERNEL_ARGS '
    case (cl_invalid_work_dimension); errcode = 'CL_INVALID_WORK_DIMENSION '
    case (cl_invalid_work_group_size); errcode = 'CL_INVALID_WORK_GROUP_SIZE '
    case (cl_invalid_work_item_size); errcode = 'CL_INVALID_WORK_ITEM_SIZE '
    case (cl_invalid_global_offset); errcode = 'CL_INVALID_GLOBAL_OFFSET '
    case (cl_invalid_event_wait_list); errcode = 'CL_INVALID_EVENT_WAIT_LIST '
    case (cl_invalid_event); errcode = 'CL_INVALID_EVENT '
    case (cl_invalid_operation); errcode = 'CL_INVALID_OPERATION '
    case (cl_invalid_gl_object); errcode = 'CL_INVALID_GL_OBJECT '
    case (cl_invalid_buffer_size); errcode = 'CL_INVALID_BUFFER_SIZE '
    case (cl_invalid_mip_level); errcode = 'CL_INVALID_MIP_LEVEL '
    case (cl_invalid_global_work_size); errcode = 'CL_INVALID_GLOBAL_WORK_SIZE '
    case (cl_platform_not_found_khr); errcode = 'CL_PLATFORM_NOT_FOUND_KHR'
    case default
      write(errcode, '(i10)') ierr
      errcode = 'UNKNOWN ERROR CODE ('//trim(adjustl(errcode))//')'
    end select
 
    message(1) = 'OpenCL '//trim(name)//' '//trim(errcode)
    call messages_fatal(1)
 
    pop_sub(opencl_print_error)
  end subroutine opencl_print_error
#endif
 
  ! ----------------------------------------------------
 
  subroutine clblas_print_error(ierr, name)
    integer,          intent(in) :: ierr
    character(len=*), intent(in) :: name
 
    character(len=40) :: errcode
 
    push_sub(clblas_print_error)
#if defined(HAVE_CLBLAS) || defined(HAVE_CLBLAST)
    select case (ierr)
    case (clblassuccess);                    errcode = 'clblasSuccess'
    case (clblasinvalidvalue);               errcode = 'clblasInvalidValue'
    case (clblasinvalidcommandqueue);        errcode = 'clblasInvalidCommandQueue'
    case (clblasinvalidcontext);             errcode = 'clblasInvalidContext'
    case (clblasinvalidmemobject);           errcode = 'clblasInvalidMemObject'
    case (clblasinvaliddevice);              errcode = 'clblasInvalidDevice'
    case (clblasinvalideventwaitlist);       errcode = 'clblasInvalidEventWaitList'
    case (clblasoutofresources);             errcode = 'clblasOutOfResources'
    case (clblasoutofhostmemory);            errcode = 'clblasOutOfHostMemory'
    case (clblasinvalidoperation);           errcode = 'clblasInvalidOperation'
    case (clblascompilernotavailable);       errcode = 'clblasCompilerNotAvailable'
    case (clblasbuildprogramfailure);        errcode = 'clblasBuildProgramFailure'
    case (clblasnotimplemented);             errcode = 'clblasNotImplemented'
    case (clblasnotinitialized);             errcode = 'clblasNotInitialized'
    case (clblasinvalidmata);                errcode = 'clblasInvalidMatA'
    case (clblasinvalidmatb);                errcode = 'clblasInvalidMatB'
    case (clblasinvalidmatc);                errcode = 'clblasInvalidMatC'
    case (clblasinvalidvecx);                errcode = 'clblasInvalidVecX'
    case (clblasinvalidvecy);                errcode = 'clblasInvalidVecY'
    case (clblasinvaliddim);                 errcode = 'clblasInvalidDim'
    case (clblasinvalidleaddima);            errcode = 'clblasInvalidLeadDimA'
    case (clblasinvalidleaddimb);            errcode = 'clblasInvalidLeadDimB'
    case (clblasinvalidleaddimc);            errcode = 'clblasInvalidLeadDimC'
    case (clblasinvalidincx);                errcode = 'clblasInvalidIncX'
    case (clblasinvalidincy);                errcode = 'clblasInvalidIncY'
    case (clblasinsufficientmemmata);        errcode = 'clblasInsufficientMemMatA'
    case (clblasinsufficientmemmatb);        errcode = 'clblasInsufficientMemMatB'
    case (clblasinsufficientmemmatc);        errcode = 'clblasInsufficientMemMatC'
    case (clblasinsufficientmemvecx);        errcode = 'clblasInsufficientMemVecX'
    case (clblasinsufficientmemvecy);        errcode = 'clblasInsufficientMemVecY'
#ifdef HAVE_CLBLAST
    case (clblastinsufficientmemorytemp);    errcode = 'clblastInsufficientMemoryTemp'
    case (clblastinvalidbatchcount);         errcode = 'clblastInvalidBatchCount'
    case (clblastinvalidoverridekernel);     errcode = 'clblastInvalidOverrideKernel'
    case (clblastmissingoverrideparameter);  errcode = 'clblastMissingOverrideParameter'
    case (clblastinvalidlocalmemusage);      errcode = 'clblastInvalidLocalMemUsage'
    case (clblastnohalfprecision);           errcode = 'clblastNoHalfPrecision'
    case (clblastnodoubleprecision);         errcode = 'clblastNoDoublePrecision'
    case (clblastinvalidvectorscalar);       errcode = 'clblastInvalidVectorScalar'
    case (clblastinsufficientmemoryscalar);  errcode = 'clblastInsufficientMemoryScalar'
    case (clblastdatabaseerror);             errcode = 'clblastDatabaseError'
    case (clblastunknownerror);              errcode = 'clblastUnknownError'
    case (clblastunexpectederror);           errcode = 'clblastUnexpectedError'
#endif
 
    case default
      write(errcode, '(i10)') ierr
      errcode = 'UNKNOWN ERROR CODE ('//trim(adjustl(errcode))//')'
    end select
#endif
 
    message(1) = 'Error in calling clblas routine '//trim(name)//' : '//trim(errcode)
    call messages_fatal(1)
 
    pop_sub(clblas_print_error)
  end subroutine clblas_print_error
 
  ! ----------------------------------------------------
  subroutine clfft_print_error(ierr, name)
    integer,          intent(in) :: ierr
    character(len=*), intent(in) :: name
 
    character(len=40) :: errcode
 
    push_sub(clfft_print_error)
#ifdef HAVE_CLFFT
    select case (ierr)
    case (clfft_invalid_global_work_size);          errcode = 'CLFFT_INVALID_GLOBAL_WORK_SIZE'
    case (clfft_invalid_mip_level);                 errcode = 'CLFFT_INVALID_MIP_LEVEL'
    case (clfft_invalid_buffer_size);               errcode = 'CLFFT_INVALID_BUFFER_SIZE'
    case (clfft_invalid_gl_object);                 errcode = 'CLFFT_INVALID_GL_OBJECT'
    case (clfft_invalid_operation);                 errcode = 'CLFFT_INVALID_OPERATION'
    case (clfft_invalid_event);                     errcode = 'CLFFT_INVALID_EVENT'
    case (clfft_invalid_event_wait_list);           errcode = 'CLFFT_INVALID_EVENT_WAIT_LIST'
    case (clfft_invalid_global_offset);             errcode = 'CLFFT_INVALID_GLOBAL_OFFSET'
    case (clfft_invalid_work_item_size);            errcode = 'CLFFT_INVALID_WORK_ITEM_SIZE'
    case (clfft_invalid_work_group_size);           errcode = 'CLFFT_INVALID_WORK_GROUP_SIZE'
    case (clfft_invalid_work_dimension);            errcode = 'CLFFT_INVALID_WORK_DIMENSION'
    case (clfft_invalid_kernel_args);               errcode = 'CLFFT_INVALID_KERNEL_ARGS'
    case (clfft_invalid_arg_size);                  errcode = 'CLFFT_INVALID_ARG_SIZE'
    case (clfft_invalid_arg_value);                 errcode = 'CLFFT_INVALID_ARG_VALUE'
    case (clfft_invalid_arg_index);                 errcode = 'CLFFT_INVALID_ARG_INDEX'
    case (clfft_invalid_kernel);                    errcode = 'CLFFT_INVALID_KERNEL'
    case (clfft_invalid_kernel_definition);         errcode = 'CLFFT_INVALID_KERNEL_DEFINITION'
    case (clfft_invalid_kernel_name);               errcode = 'CLFFT_INVALID_KERNEL_NAME'
    case (clfft_invalid_program_executable);        errcode = 'CLFFT_INVALID_PROGRAM_EXECUTABLE'
    case (clfft_invalid_program);                   errcode = 'CLFFT_INVALID_PROGRAM'
    case (clfft_invalid_build_options);             errcode = 'CLFFT_INVALID_BUILD_OPTIONS'
    case (clfft_invalid_binary);                    errcode = 'CLFFT_INVALID_BINARY'
    case (clfft_invalid_sampler);                   errcode = 'CLFFT_INVALID_SAMPLER'
    case (clfft_invalid_image_size);                errcode = 'CLFFT_INVALID_IMAGE_SIZE'
    case (clfft_invalid_image_format_descriptor);   errcode = 'CLFFT_INVALID_IMAGE_FORMAT_DESCRIPTOR'
    case (clfft_invalid_mem_object);                errcode = 'CLFFT_INVALID_MEM_OBJECT'
    case (clfft_invalid_host_ptr);                  errcode = 'CLFFT_INVALID_HOST_PTR'
    case (clfft_invalid_command_queue);             errcode = 'CLFFT_INVALID_COMMAND_QUEUE'
    case (clfft_invalid_queue_properties);          errcode = 'CLFFT_INVALID_QUEUE_PROPERTIES'
    case (clfft_invalid_context);                   errcode = 'CLFFT_INVALID_CONTEXT'
    case (clfft_invalid_device);                    errcode = 'CLFFT_INVALID_DEVICE'
    case (clfft_invalid_platform);                  errcode = 'CLFFT_INVALID_PLATFORM'
    case (clfft_invalid_device_type);               errcode = 'CLFFT_INVALID_DEVICE_TYPE'
    case (clfft_invalid_value);                     errcode = 'CLFFT_INVALID_VALUE'
    case (clfft_map_failure);                       errcode = 'CLFFT_MAP_FAILURE'
    case (clfft_build_program_failure);             errcode = 'CLFFT_BUILD_PROGRAM_FAILURE'
    case (clfft_image_format_not_supported);        errcode = 'CLFFT_IMAGE_FORMAT_NOT_SUPPORTED'
    case (clfft_image_format_mismatch);             errcode = 'CLFFT_IMAGE_FORMAT_MISMATCH'
    case (clfft_mem_copy_overlap);                  errcode = 'CLFFT_MEM_COPY_OVERLAP'
    case (clfft_profiling_info_not_available);      errcode = 'CLFFT_PROFILING_INFO_NOT_AVAILABLE'
    case (clfft_out_of_host_memory);                errcode = 'CLFFT_OUT_OF_HOST_MEMORY'
    case (clfft_out_of_resources);                  errcode = 'CLFFT_OUT_OF_RESOURCES'
    case (clfft_mem_object_allocation_failure);     errcode = 'CLFFT_MEM_OBJECT_ALLOCATION_FAILURE'
    case (clfft_compiler_not_available);            errcode = 'CLFFT_COMPILER_NOT_AVAILABLE'
    case (clfft_device_not_available);              errcode = 'CLFFT_DEVICE_NOT_AVAILABLE'
    case (clfft_device_not_found);                  errcode = 'CLFFT_DEVICE_NOT_FOUND'
    case (clfft_success);                           errcode = 'CLFFT_SUCCESS'
    case (clfft_bugcheck);                          errcode = 'CLFFT_BUGCHECK'
    case (clfft_notimplemented);                    errcode = 'CLFFT_NOTIMPLEMENTED'
    case (clfft_file_not_found);                    errcode = 'CLFFT_FILE_NOT_FOUND'
    case (clfft_file_create_failure);               errcode = 'CLFFT_FILE_CREATE_FAILURE'
    case (clfft_version_mismatch);                  errcode = 'CLFFT_VERSION_MISMATCH'
    case (clfft_invalid_plan);                      errcode = 'CLFFT_INVALID_PLAN'
    case (clfft_device_no_double);                  errcode = 'CLFFT_DEVICE_NO_DOUBLE'
    case (clfft_endstatus);                         errcode = 'CLFFT_ENDSTATUS'
    case default
      write(errcode, '(i10)') ierr
      errcode = 'UNKNOWN ERROR CODE ('//trim(adjustl(errcode))//')'
    end select
#endif
 
    message(1) = 'clfft '//trim(name)//' '//trim(errcode)
    call messages_fatal(1)
 
    pop_sub(clfft_print_error)
  end subroutine clfft_print_error
 
  ! ----------------------------------------------------
 
#ifdef HAVE_OPENCL
  logical function f90_cl_device_has_extension(device, extension) result(has)
    type(cl_device_id), intent(inout) :: device
    character(len=*),   intent(in)    :: extension
 
    integer :: cl_status
    character(len=2048) :: all_extensions
 
#ifdef HAVE_OPENCL
    call clgetdeviceinfo(device, cl_device_extensions, all_extensions, cl_status)
#endif
 
    has = index(all_extensions, extension) /= 0
 
  end function f90_cl_device_has_extension
#endif
 
  ! ---------------------------------------------------------
 
  integer(int64) pure function opencl_pad(size, blk) result(pad)
    integer(int64), intent(in) :: size
    integer, intent(in) :: blk
 
    integer(int64) :: mm
 
    mm = mod(size, int(blk, int64))
    if (mm == 0) then
      pad = size
    else
      pad = size + blk - mm
    end if
  end function opencl_pad
 
  ! ----------------------------------------------------
 
  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
 
    integer :: bsize
    integer(int64) :: nval_real, offset_real
 
    push_sub(accel_set_buffer_to_zero_i8)
 
    if (nval > 0) then
 
      nval_real = nval
      if (type == type_cmplx) nval_real = nval_real * 2
      if (present(offset)) then
        offset_real = offset
        if (type == type_cmplx) offset_real = offset_real * 2
      else
        offset_real = 0_int64
      end if
 
      assert(nval_real > 0)
 
      call accel_set_kernel_arg(set_zero, 0, nval_real)
      call accel_set_kernel_arg(set_zero, 1, offset_real)
      call accel_set_kernel_arg(set_zero, 2, buffer)
 
      bsize = accel_kernel_workgroup_size(set_zero)
 
      call accel_kernel_run(set_zero, (/ opencl_pad(nval_real, bsize) /), (/ int(bsize, int64) /))
 
      if(.not. optional_default(async, .false.)) call accel_finish()
 
    end if
 
    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 opencl_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_release_buffer(buff)
 
      safe_deallocate_a(data)
 
      size = int(size*2.0)
 
      if (size > 50000000) exit
    end do
  end subroutine opencl_check_bandwidth
 
  ! ----------------------------------------------------
 
  logical pure function accel_use_shared_mem() result(use_shared_mem)
 
    use_shared_mem = accel%shared_mem
 
  end function accel_use_shared_mem
 
  !------------------------------------------------------------
 
  subroutine accel_kernel_global_init()
 
    push_sub(accel_kernel_global_init)
 
    nullify(head)
 
    call cuda_module_map_init(accel%module_map)
 
    pop_sub(accel_kernel_global_init)
  end subroutine accel_kernel_global_init
 
  !------------------------------------------------------------
 
  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_module_map_end(accel%module_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_OPENCL
    type(cl_program) :: prog
#endif
#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)//'/opencl/'//" "//trim(accel%debug_flag)
 
    if (accel_use_shared_mem()) then
      all_flags = trim(all_flags)//' -DSHARED_MEM'
    end if
 
    if (present(flags)) then
      all_flags = trim(all_flags)//' '//trim(flags)
    end if
 
    call cuda_build_program(accel%module_map, this%cuda_module, accel%device%cuda_device, trim(file_name), trim(all_flags))
 
    call cuda_create_kernel(this%cuda_kernel, this%cuda_module, trim(kernel_name))
    call cuda_alloc_arg_array(this%arguments)
 
    this%cuda_shared_mem = 0
#endif
 
#ifdef HAVE_OPENCL
    call opencl_build_program(prog, trim(conf%share)//'/opencl/'//trim(file_name), flags = flags)
    call opencl_create_kernel(this%kernel, prog, trim(kernel_name))
    call opencl_release_program(prog)
#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
#ifdef HAVE_OPENCL
    integer :: ierr
#endif
 
    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
 
#ifdef HAVE_OPENCL
    call clreleasekernel(this%kernel, ierr)
    if (ierr /= cl_success) call opencl_print_error(ierr, "release_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_local_memory_size() result(size)
 
    size = accel%local_memory_size
 
  end function accel_local_memory_size
 
  !--------------------------------------------------------------
 
  integer pure function accel_max_size_per_dim(dim) result(size)
    integer, intent(in) :: dim
 
    size = 0
#ifdef HAVE_OPENCL
    size = 32768 ! Setting here arbitrarily higher dimensions to 32768, as 2**30 leads to a
    ! value of zero when multiplied by 2048 and converted to integer 4.
    if (dim == 1) size = 2**30
#endif
#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)
      call cublas_set_stream(accel%cublas_handle, accel%cuda_stream)
#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 = size
    thread_block_size = size
#else
    grid_size = size * accel%warp_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 "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: