Add hip support for fft on gpu (generic interface)

INSTALL.md

@@ -478,7 +478,11 @@ partially depending on installed libraries (see 2.)
 - `-D__LIBXC` use LIBXC
 - `-D__ELPA` use ELPA in place of SYEVD  to solve the eigenvalue problem
 - `-D__FFTW3` FFTW version 3 is recommended
-- `-D__PW_CUDA` CUDA FFT and associated gather/scatter on the GPU
+- `-D__PW_GPU` CUDA or hip FFT and associated gather/scatter on the GPU
+- `-D__PW_CUDA` CUDA FFT and associated gather/scatter on the GPU. `__PW_GPU`
+  needs to be set
+- `-D__PW_HIP` HIP FFT and associated gather/scatter on the GPU. `__PW_GPU`
+  needs to be set
 - `-D__MKL` link the MKL library for linear algebra and/or FFT
 - `-D__GRID_CORE=X` (with X=1..6) specific optimized core routines can be
   selected.  Reasonable defaults are [provided](./src/grid/collocate_fast.f90)

Makefile

@@ -61,39 +61,42 @@ ifeq ($(OFFLOAD_CC),)
 EXE_NAMES := $(basename $(notdir $(filter-out %.cu, $(ALL_EXE_FILES))))
 endif
 ifneq ($(LD_SHARED),)
- ARCHIVE_EXT := .so
+	ARCHIVE_EXT := .so
 else
- ARCHIVE_EXT := .a
+	ARCHIVE_EXT := .a
 endif
 include $(EXTSHOME)/Makefile.inc
 endif
 
 # Declare PHONY targets =====================================================
 .PHONY : $(VERSION) $(EXE_NAMES) \
-				 dirs makedep default_target all \
-				 toolversions exts extversions extclean \
-				 libcp2k cp2k_shell pkgconfig python-bindings \
-				 pre-commit pre-commit-clean \
-				 pretty precommit precommitclean doxygenclean doxygen \
-				 fpretty fprettyclean \
-				 doxify doxifyclean \
-				 install clean realclean distclean mrproper help \
-				 test testbg testclean testrealclean \
-				 data \
-	 $(EXTSPACKAGES)
+					dirs makedep default_target all \
+					toolversions exts extversions extclean \
+					libcp2k cp2k_shell pkgconfig python-bindings \
+					pre-commit pre-commit-clean \
+					pretty precommit precommitclean doxygenclean doxygen \
+					fpretty fprettyclean \
+					doxify doxifyclean \
+					install clean realclean distclean mrproper help \
+					test testbg testclean testrealclean \
+					data \
+		$(EXTSPACKAGES)
 
 # Discover files and directories ============================================
 ALL_SRC_DIRS := $(shell find $(SRCDIR) -type d ! -name preprettify | awk '{printf("%s:",$$1)}')
 ALL_PREPRETTY_DIRS = $(shell find $(SRCDIR) -type d -name preprettify)
 
 ALL_PKG_FILES  = $(shell find $(SRCDIR) ! -path "*/preprettify/*" -name "PACKAGE")
-OBJ_SRC_FILES  = $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -name "*.F")
+OBJ_SRC_FILES  = $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -name "*.F" ! -name "pw_gpu.F")
 OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" ! -path "*/python*" -name "*.c")
 OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -name "*.cpp")
 OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -name "*.cxx")
-OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -path "./grid/hip" -prune -type f -name "*.cc")
+OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" ! -path "./grid/hip"  ! -path "./pw/gpu" -prune -type f -name "pw_gpu.F" -prune -type f -name "*.cc")
 
 ifneq (,$(findstring hipcc,$(OFFLOAD_CC)))
+OBJ_SRC_FILES += ./pw/gpu/hip/pw_hip_z.cc
+OBJ_SRC_FILES += ./pw/gpu/pw_gpu_internal.cc
+OBJ_SRC_FILES += ./pw/pw_gpu.F
 OBJ_SRC_FILES += ./grid/hip/grid_hip_integrate.cc
 OBJ_SRC_FILES += ./grid/hip/grid_hip_collocate.cc
 OBJ_SRC_FILES += ./grid/hip/grid_hip_context.cc
@@ -104,6 +107,8 @@ endif
 # if we compile for cuda by directly calling nvcc then include all cuda files.
 ifneq (,$(findstring nvcc,$(OFFLOAD_CC)))
 OBJ_SRC_FILES += $(shell cd $(SRCDIR); find . ! -path "*/preprettify/*" -name "*.cu")
+OBJ_SRC_FILES += ./pw/gpu/pw_gpu_internal.cc
+OBJ_SRC_FILES += ./pw/pw_gpu.F
 endif
 
 
@@ -463,9 +468,9 @@ MAKEDEPMODE = "normal"
 ifeq ($(HACKDEP),yes)
 MAKEDEPMODE = "hackdep"
 else
- ifneq ($(MC),)
- MAKEDEPMODE = "mod_compiler"
- endif
+	ifneq ($(MC),)
+	MAKEDEPMODE = "mod_compiler"
+	endif
 endif
 
 # this happens on stage 3
@@ -516,10 +521,10 @@ cp2k_info.o: $(GIT_REF)
 
 # Add some practical metadata about the build.
 FCFLAGS += -D__COMPILE_ARCH="\"$(ARCH)\""\
-					 -D__COMPILE_DATE="\"$(shell date)\""\
-					 -D__COMPILE_HOST="\"$(shell hostname 2>/dev/null || hostnamectl --transient)\""\
-					 -D__COMPILE_REVISION="\"$(strip $(REVISION))\""\
-					 -D__DATA_DIR="\"$(DATA_DIR)\""
+						-D__COMPILE_DATE="\"$(shell date)\""\
+						-D__COMPILE_HOST="\"$(shell hostname 2>/dev/null || hostnamectl --transient)\""\
+						-D__COMPILE_REVISION="\"$(strip $(REVISION))\""\
+						-D__DATA_DIR="\"$(DATA_DIR)\""
 
 # $(FCLOGPIPE) can be used to store compiler output, e.g. warnings, for each F-file separately.
 # This is used e.g. by the convention checker.
@@ -558,6 +563,8 @@ ifneq (,$(findstring hipcc,$(OFFLOAD_CC)))
 
 #define specific rules for the hip backend since they need to be propcessed by hipcc
 
+pw_hip_z.o: pw_hip_z.cc
+	$(OFFLOAD_CC) -c $(OFFLOAD_FLAGS) $<
 grid_hip_collocate.o: grid_hip_collocate.cc
 	$(OFFLOAD_CC) -c $(OFFLOAD_FLAGS) $<
 grid_hip_integrate.o: grid_hip_integrate.cc

src/cp2k_info.F

@@ -143,18 +143,26 @@ CONTAINS
       flags = TRIM(flags)//" quip"
 #endif
 
+#if defined(__PW_GPU)
+      flags = TRIM(flags)//" pw_gpu"
+#if defined(__PW_HIP)
+      flags = TRIM(flags)//" pw_hip"
+#endif
 #if defined(__PW_CUDA)
       flags = TRIM(flags)//" pw_cuda"
 #endif
+#endif
+
+#if defined(__HAS_PATCHED_CUFFT_70)
+      flags = TRIM(flags)//" patched_cufft_70"
+#endif
+
 #if defined(__PW_FPGA)
       flags = TRIM(flags)//" pw_fpga"
 #endif
 #if defined(__PW_FPGA_SP)
       flags = TRIM(flags)//" pw_fpga_sp"
 #endif
-#if defined(__HAS_PATCHED_CUFFT_70)
-      flags = TRIM(flags)//" patched_cufft_70"
-#endif
 
 #if defined(__HAS_smm_vec)
       flags = TRIM(flags)//" smm_vec"

src/f77_interface.F

@@ -18,132 +18,134 @@
 !> \author fawzi & Johanna
 ! **************************************************************************************************
 MODULE f77_interface
-   USE base_hooks,                      ONLY: cp_abort_hook,&
-                                              cp_warn_hook,&
-                                              timeset_hook,&
-                                              timestop_hook
-   USE bibliography,                    ONLY: add_all_references
-   USE cell_types,                      ONLY: cell_type,&
-                                              init_cell
-   USE cp2k_info,                       ONLY: get_runtime_info
-   USE cp_error_handling,               ONLY: cp_error_handling_setup
-   USE cp_files,                        ONLY: init_preconnection_list,&
-                                              open_file
-   USE cp_log_handling,                 ONLY: &
-        cp_add_default_logger, cp_default_logger_stack_size, cp_failure_level, &
-        cp_get_default_logger, cp_logger_create, cp_logger_get_default_unit_nr, cp_logger_release, &
-        cp_logger_retain, cp_logger_type, cp_rm_default_logger, cp_to_string
-   USE cp_output_handling,              ONLY: cp_iterate
-   USE cp_para_env,                     ONLY: cp_para_env_create,&
-                                              cp_para_env_release,&
-                                              cp_para_env_retain
-   USE cp_para_types,                   ONLY: cp_para_env_type
-   USE cp_result_methods,               ONLY: get_results,&
-                                              test_for_result
-   USE cp_result_types,                 ONLY: cp_result_type
-   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
-                                              cp_subsys_set,&
-                                              cp_subsys_type,&
-                                              unpack_subsys_particles
-   USE dbcsr_api,                       ONLY: dbcsr_finalize_lib,&
-                                              dbcsr_init_lib
-   USE dbm_api,                         ONLY: dbm_library_finalize,&
-                                              dbm_library_init
-   USE eip_environment,                 ONLY: eip_init
-   USE eip_environment_types,           ONLY: eip_env_create,&
-                                              eip_env_release,&
-                                              eip_environment_type
-   USE embed_main,                      ONLY: embed_create_force_env
-   USE embed_types,                     ONLY: embed_env_release,&
-                                              embed_env_type
-   USE environment,                     ONLY: cp2k_finalize,&
-                                              cp2k_init,&
-                                              cp2k_read,&
-                                              cp2k_setup
-   USE fist_main,                       ONLY: fist_create_force_env
-   USE force_env_methods,               ONLY: force_env_calc_energy_force,&
-                                              force_env_create
-   USE force_env_types,                 ONLY: &
-        force_env_get, force_env_get_frc, force_env_get_natom, force_env_get_nparticle, &
-        force_env_get_pos, force_env_get_vel, force_env_release, force_env_retain, force_env_set, &
-        force_env_type, multiple_fe_list
-   USE fp_types,                        ONLY: fp_env_create,&
-                                              fp_env_read,&
-                                              fp_env_release,&
-                                              fp_env_write,&
-                                              fp_type
-   USE global_types,                    ONLY: global_environment_type,&
-                                              globenv_create,&
-                                              globenv_release,&
-                                              globenv_retain
-   USE grid_api,                        ONLY: grid_library_finalize,&
-                                              grid_library_init
-   USE input_constants,                 ONLY: &
-        do_eip, do_embed, do_fist, do_mixed, do_nnp, do_qmmm, do_qmmmx, do_qs, do_sirius
-   USE input_cp2k_check,                ONLY: check_cp2k_input
-   USE input_cp2k_force_eval,           ONLY: create_force_eval_section
-   USE input_cp2k_read,                 ONLY: empty_initial_variables,&
-                                              read_input
-   USE input_enumeration_types,         ONLY: enum_i2c,&
-                                              enumeration_type
-   USE input_keyword_types,             ONLY: keyword_get,&
-                                              keyword_type
-   USE input_section_types,             ONLY: &
-        section_get_keyword, section_release, section_type, section_vals_duplicate, &
-        section_vals_get, section_vals_get_subs_vals, section_vals_release, &
-        section_vals_remove_values, section_vals_retain, section_vals_type, section_vals_val_get, &
-        section_vals_write
-   USE kinds,                           ONLY: default_path_length,&
-                                              default_string_length,&
-                                              dp
-   USE machine,                         ONLY: default_output_unit,&
-                                              m_chdir,&
-                                              m_getcwd,&
-                                              m_memory
-   USE message_passing,                 ONLY: &
-        add_mp_perf_env, get_mp_perf_env, mp_comm_world, mp_environ, mp_max, mp_perf_env_release, &
-        mp_perf_env_retain, mp_perf_env_type, mp_world_finalize, mp_world_init, rm_mp_perf_env
-   USE metadynamics_types,              ONLY: meta_env_release,&
-                                              meta_env_type
-   USE metadynamics_utils,              ONLY: metadyn_read
-   USE mixed_environment_types,         ONLY: mixed_env_release,&
-                                              mixed_environment_type
-   USE mixed_main,                      ONLY: mixed_create_force_env
-   USE nnp_environment,                 ONLY: nnp_init
-   USE nnp_environment_types,           ONLY: nnp_env_create,&
-                                              nnp_env_release,&
-                                              nnp_type
-   USE offload_api,                     ONLY: offload_get_device_count,&
-                                              offload_get_device_id,&
-                                              offload_set_device_id
-   USE periodic_table,                  ONLY: init_periodic_table
-   USE pw_cuda,                         ONLY: pw_cuda_finalize,&
-                                              pw_cuda_init
-   USE pwdft_environment,               ONLY: pwdft_init
-   USE pwdft_environment_types,         ONLY: pwdft_env_create,&
-                                              pwdft_env_release,&
-                                              pwdft_environment_type
-   USE qmmm_create,                     ONLY: qmmm_env_create
-   USE qmmm_types,                      ONLY: qmmm_env_release,&
-                                              qmmm_env_type
-   USE qmmmx_create,                    ONLY: qmmmx_env_create
-   USE qmmmx_types,                     ONLY: qmmmx_env_release,&
-                                              qmmmx_env_type
-   USE qs_environment,                  ONLY: qs_init
-   USE qs_environment_types,            ONLY: qs_env_create,&
-                                              qs_env_release,&
-                                              qs_environment_type
-   USE reference_manager,               ONLY: remove_all_references
-   USE string_table,                    ONLY: string_table_allocate,&
-                                              string_table_deallocate
-   USE timings,                         ONLY: add_timer_env,&
-                                              get_timer_env,&
-                                              rm_timer_env,&
-                                              timer_env_release,&
-                                              timer_env_retain,&
-                                              timings_register_hooks
-   USE timings_types,                   ONLY: timer_env_type
-   USE virial_types,                    ONLY: virial_type
+   USE base_hooks, ONLY: cp_abort_hook, &
+                         cp_warn_hook, &
+                         timeset_hook, &
+                         timestop_hook
+   USE bibliography, ONLY: add_all_references
+   USE cell_types, ONLY: cell_type, &
+                         init_cell
+   USE cp2k_info, ONLY: get_runtime_info
+   USE cp_error_handling, ONLY: cp_error_handling_setup
+   USE cp_files, ONLY: init_preconnection_list, &
+                       open_file
+   USE cp_log_handling, ONLY: &
+      cp_add_default_logger, cp_default_logger_stack_size, cp_failure_level, &
+      cp_get_default_logger, cp_logger_create, cp_logger_get_default_unit_nr, cp_logger_release, &
+      cp_logger_retain, cp_logger_type, cp_rm_default_logger, cp_to_string
+   USE cp_output_handling, ONLY: cp_iterate
+   USE cp_para_env, ONLY: cp_para_env_create, &
+                          cp_para_env_release, &
+                          cp_para_env_retain
+   USE cp_para_types, ONLY: cp_para_env_type
+   USE cp_result_methods, ONLY: get_results, &
+                                test_for_result
+   USE cp_result_types, ONLY: cp_result_type
+   USE cp_subsys_types, ONLY: cp_subsys_get, &
+                              cp_subsys_set, &
+                              cp_subsys_type, &
+                              unpack_subsys_particles
+   USE dbcsr_api, ONLY: dbcsr_finalize_lib, &
+                        dbcsr_init_lib
+   USE dbm_api, ONLY: dbm_library_finalize, &
+                      dbm_library_init
+   USE eip_environment, ONLY: eip_init
+   USE eip_environment_types, ONLY: eip_env_create, &
+                                    eip_env_release, &
+                                    eip_environment_type
+   USE embed_main, ONLY: embed_create_force_env
+   USE embed_types, ONLY: embed_env_release, &
+                          embed_env_type
+   USE environment, ONLY: cp2k_finalize, &
+                          cp2k_init, &
+                          cp2k_read, &
+                          cp2k_setup
+   USE fist_main, ONLY: fist_create_force_env
+   USE force_env_methods, ONLY: force_env_calc_energy_force, &
+                                force_env_create
+   USE force_env_types, ONLY: &
+      force_env_get, force_env_get_frc, force_env_get_natom, force_env_get_nparticle, &
+      force_env_get_pos, force_env_get_vel, force_env_release, force_env_retain, force_env_set, &
+      force_env_type, multiple_fe_list
+   USE fp_types, ONLY: fp_env_create, &
+                       fp_env_read, &
+                       fp_env_release, &
+                       fp_env_write, &
+                       fp_type
+   USE global_types, ONLY: global_environment_type, &
+                           globenv_create, &
+                           globenv_release, &
+                           globenv_retain
+   USE grid_api, ONLY: grid_library_finalize, &
+                       grid_library_init
+   USE input_constants, ONLY: &
+      do_eip, do_embed, do_fist, do_mixed, do_nnp, do_qmmm, do_qmmmx, do_qs, do_sirius
+   USE input_cp2k_check, ONLY: check_cp2k_input
+   USE input_cp2k_force_eval, ONLY: create_force_eval_section
+   USE input_cp2k_read, ONLY: empty_initial_variables, &
+                              read_input
+   USE input_enumeration_types, ONLY: enum_i2c, &
+                                      enumeration_type
+   USE input_keyword_types, ONLY: keyword_get, &
+                                  keyword_type
+   USE input_section_types, ONLY: &
+      section_get_keyword, section_release, section_type, section_vals_duplicate, &
+      section_vals_get, section_vals_get_subs_vals, section_vals_release, &
+      section_vals_remove_values, section_vals_retain, section_vals_type, section_vals_val_get, &
+      section_vals_write
+   USE kinds, ONLY: default_path_length, &
+                    default_string_length, &
+                    dp
+   USE machine, ONLY: default_output_unit, &
+                      m_chdir, &
+                      m_getcwd, &
+                      m_memory
+   USE message_passing, ONLY: &
+      add_mp_perf_env, get_mp_perf_env, mp_comm_world, mp_environ, mp_max, mp_perf_env_release, &
+      mp_perf_env_retain, mp_perf_env_type, mp_world_finalize, mp_world_init, rm_mp_perf_env
+   USE metadynamics_types, ONLY: meta_env_release, &
+                                 meta_env_type
+   USE metadynamics_utils, ONLY: metadyn_read
+   USE mixed_environment_types, ONLY: mixed_env_release, &
+                                      mixed_environment_type
+   USE mixed_main, ONLY: mixed_create_force_env
+   USE nnp_environment, ONLY: nnp_init
+   USE nnp_environment_types, ONLY: nnp_env_create, &
+                                    nnp_env_release, &
+                                    nnp_type
+   USE offload_api, ONLY: offload_get_device_count, &
+                          offload_get_device_id, &
+                          offload_set_device_id
+   USE periodic_table, ONLY: init_periodic_table
+#if defined(__PW_GPU)
+   USE pw_gpu, ONLY: pw_gpu_finalize, &
+                     pw_gpu_init
+#endif
+   USE pwdft_environment, ONLY: pwdft_init
+   USE pwdft_environment_types, ONLY: pwdft_env_create, &
+                                      pwdft_env_release, &
+                                      pwdft_environment_type
+   USE qmmm_create, ONLY: qmmm_env_create
+   USE qmmm_types, ONLY: qmmm_env_release, &
+                         qmmm_env_type
+   USE qmmmx_create, ONLY: qmmmx_env_create
+   USE qmmmx_types, ONLY: qmmmx_env_release, &
+                          qmmmx_env_type
+   USE qs_environment, ONLY: qs_init
+   USE qs_environment_types, ONLY: qs_env_create, &
+                                   qs_env_release, &
+                                   qs_environment_type
+   USE reference_manager, ONLY: remove_all_references
+   USE string_table, ONLY: string_table_allocate, &
+                           string_table_deallocate
+   USE timings, ONLY: add_timer_env, &
+                      get_timer_env, &
+                      rm_timer_env, &
+                      timer_env_release, &
+                      timer_env_retain, &
+                      timings_register_hooks
+   USE timings_types, ONLY: timer_env_type
+   USE virial_types, ONLY: virial_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -292,8 +294,9 @@ CONTAINS
          CALL dbcsr_init_lib(default_para_env%group, timeset_hook, timestop_hook, &
                              cp_abort_hook, cp_warn_hook, io_unit=unit_nr)
 #endif
-         CALL pw_cuda_init()
-
+#if defined(__PW_CUDA) || defined(__PW_HIP)
+         CALL pw_gpu_init()
+#endif
          CALL grid_library_init()
 
          CALL dbm_library_init()
@@ -337,8 +340,9 @@ CONTAINS
 
          ! Finalize the DBCSR configuration
          CALL dbcsr_finalize_lib()
-         CALL pw_cuda_finalize()
-
+#if defined(__PW_CUDA) || defined(__PW_HIP)
+         CALL pw_gpu_finalize()
+#endif
          CALL cp_para_env_release(default_para_env)
          ierr = 0
          CALL cp_rm_default_logger()

src/grid/hip/grid_hip_context.hpp

@@ -23,7 +23,7 @@ extern "C" {
 }
 
 #include "../../offload/offload_library.h"
-
+#include "../../offload/offload_operations.h"
 namespace rocm_backend {
 // a little helper class in the same spirit than std::vector. it must exist
 // somewhere. Maybe possible to get the same thing with std::vector and
@@ -369,7 +369,7 @@ public:
   ~context_info() { clear(); }
 
   void clear() {
-    offloadSetDevice(device_id_);
+    hipSetDevice(device_id_);
     tasks_dev.reset();
     block_offsets_dev.reset();
     coef_dev_.reset();
@@ -459,7 +459,7 @@ public:
     offloadDeviceSynchronize();
   }
 
-  void set_device() { offloadSetDevice(device_id_); }
+  void set_device() { hipSetDevice(device_id_); }
 
   void collocate_one_grid_level(const int level, const enum grid_func func,
                                 int *lp_diff);

src/offload/offload_api.F

@@ -32,6 +32,7 @@ MODULE offload_api
    PUBLIC :: offload_set_device_id, offload_get_device_id, offload_set_device
    PUBLIC :: offload_timeset, offload_timestop, offload_mem_info
    PUBLIC :: offload_buffer_type, offload_create_buffer, offload_free_buffer
+   PUBLIC :: offload_malloc_pinned_mem, offload_free_pinned_mem
 
    TYPE offload_buffer_type
       REAL(KIND=dp), DIMENSION(:), POINTER :: host_buffer => Null()
@@ -40,6 +41,51 @@ MODULE offload_api
 
 CONTAINS
 
+! **************************************************************************************************
+!> \brief allocate pinned memory.
+!> \param buffer address of the buffer
+!> \param length length of the buffer
+!> \return 0
+! **************************************************************************************************
+   FUNCTION offload_malloc_pinned_mem(buffer, length) RESULT(res)
+      TYPE(C_PTR)                                        :: buffer
+      INTEGER(C_SIZE_T), VALUE                           :: length
+      INTEGER                                            :: res
+
+      INTERFACE
+         FUNCTION offload_malloc_pinned_mem_c(buffer, length) &
+            BIND(C, name="offload_host_malloc")
+            IMPORT C_SIZE_T, C_PTR, C_INT
+            TYPE(C_PTR)              :: buffer
+            INTEGER(C_SIZE_T), VALUE :: length
+            INTEGER(KIND=C_INT)                :: offload_malloc_pinned_mem_c
+         END FUNCTION offload_malloc_pinned_mem_c
+      END INTERFACE
+
+      res = offload_malloc_pinned_mem_c(buffer, length)
+   END FUNCTION offload_malloc_pinned_mem
+
+! **************************************************************************************************
+!> \brief free pinned memory
+!> \param buffer address of the buffer
+!> \return 0
+! **************************************************************************************************
+   FUNCTION offload_free_pinned_mem(buffer) RESULT(res)
+      TYPE(C_PTR), VALUE                                 :: buffer
+      INTEGER                                            :: res
+
+      INTERFACE
+         FUNCTION offload_free_pinned_mem_c(buffer) &
+            BIND(C, name="offload_host_free")
+            IMPORT C_PTR, C_INT
+            INTEGER(KIND=C_INT)                :: offload_free_pinned_mem_c
+            TYPE(C_PTR), VALUE            :: buffer
+         END FUNCTION offload_free_pinned_mem_c
+      END INTERFACE
+
+      res = offload_free_pinned_mem_c(buffer)
+   END FUNCTION offload_free_pinned_mem
+
 ! **************************************************************************************************
 !> \brief Returns the number of available devices.
 !> \return ...
@@ -266,5 +312,4 @@ CONTAINS
 
       CALL timestop(handle)
    END SUBROUTINE offload_free_buffer
-
 END MODULE offload_api

src/offload/offload_buffer.c

@@ -11,7 +11,7 @@
 
 #include "offload_buffer.h"
 #include "offload_library.h"
-
+#include "offload_operations.h"
 /*******************************************************************************
  * \brief Allocates a buffer of given length, ie. number of elements.
  * \author Ole Schuett

src/offload/offload_cuda_internal.h

@@ -8,6 +8,8 @@
 #ifndef OFFLOAD_CUDA_INTERNAL_H
 #define OFFLOAD_CUDA_INTERNAL_H
 #include <cuda_runtime.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 #ifdef __cplusplus
 extern "C" {
@@ -26,21 +28,24 @@ typedef cudaEvent_t offloadEvent_t;
     abort();                                                                   \
   }
 
-static inline void offloadMemsetAsync(void *ptr__, int val__, size_t size__,
+static inline void offloadMemsetAsync(void *const ptr__, const int val__,
+                                      const size_t size__,
                                       offloadStream_t stream__) {
   OFFLOAD_CHECK(cudaMemsetAsync(ptr__, val__, size__, stream__));
 }
 
-static inline void offloadMemcpyAsyncHtoD(void *ptr1__, void *ptr2__,
-                                          size_t size__,
+static inline void offloadMemcpyAsyncHtoD(void *const ptr1__,
+                                          const void *ptr2__,
+                                          const size_t size__,
                                           offloadStream_t stream__) {
   OFFLOAD_CHECK(cudaMemcpyAsync(ptr1__, ptr2__, size__, cudaMemcpyHostToDevice,
                                 stream__));
 }
 
-static inline void offloadMemcpyAsyncDtoH(void *ptr1__, void *ptr2__,
-                                          size_t size__,
-                                          offloadStream_t stream__) {
+static inline void offloadMemcpyAsyncDtoH(void *const ptr1__,
+                                          const void *ptr2__,
+                                          const size_t size__,
+                                          const offloadStream_t stream__) {
   OFFLOAD_CHECK(cudaMemcpyAsync(ptr1__, ptr2__, size__, cudaMemcpyDeviceToHost,
                                 stream__));
 }
@@ -74,12 +79,20 @@ static inline void offloadEventRecord(offloadEvent_t event__,
   OFFLOAD_CHECK(cudaEventRecord(event__, stream__));
 }
 
-static inline void offloadMalloc(void *ptr__, size_t size__) {
-  OFFLOAD_CHECK(cudaMalloc((void **)ptr__, size__));
+static inline void offloadMallocHost(void **ptr__, size_t size__) {
+  OFFLOAD_CHECK(cudaMallocHost(ptr__, size__));
+}
+
+static inline void offloadMalloc(void **ptr__, size_t size__) {
+  OFFLOAD_CHECK(cudaMalloc(ptr__, size__));
 }
 
 static inline void offloadFree(void *ptr__) { OFFLOAD_CHECK(cudaFree(ptr__)); }
 
+static inline void offloadFreeHost(void *ptr__) {
+  OFFLOAD_CHECK(cudaFreeHost(ptr__));
+}
+
 static inline void offloadStreamWaitEvent(offloadStream_t stream__,
                                           offloadEvent_t event__,
                                           const int val__) {

src/offload/offload_hip_internal.h

@@ -8,6 +8,8 @@
 #ifndef OFFLOAD_HIP_INTERNAL_H
 #define OFFLOAD_HIP_INTERNAL_H
 #include <hip/hip_runtime_api.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 #ifdef __cplusplus
 extern "C" {
@@ -27,33 +29,37 @@ typedef hipEvent_t offloadEvent_t;
     abort();                                                                   \
   }
 
-static inline void offloadMemsetAsync(void *ptr__, int val__, size_t size__,
-                                      offloadStream_t stream__) {
-  OFFLOAD_CHECK(hipMemsetAsync(ptr__, val__, size__, stream__));
+static inline void offloadMemsetAsync(void *ptr_device__, int val__,
+                                      size_t size__, offloadStream_t stream__) {
+  OFFLOAD_CHECK(hipMemsetAsync(ptr_device__, val__, size__, stream__));
 }
 
-static inline void offloadMemcpyAsyncHtoD(void *ptr1__, const void *ptr2__,
+static inline void offloadMemcpyAsyncHtoD(void *ptr_device__,
+                                          const void *ptr_host__,
                                           const size_t size__,
                                           offloadStream_t stream__) {
-  OFFLOAD_CHECK(
-      hipMemcpyAsync(ptr1__, ptr2__, size__, hipMemcpyHostToDevice, stream__));
+  OFFLOAD_CHECK(hipMemcpyAsync(ptr_device__, ptr_host__, size__,
+                               hipMemcpyHostToDevice, stream__));
 }
 
-static inline void offloadMemcpyAsyncDtoH(void *ptr1__, const void *ptr2__,
+static inline void offloadMemcpyAsyncDtoH(void *ptr_host__,
+                                          const void *ptr_device__,
                                           const size_t size__,
                                           offloadStream_t stream__) {
+  OFFLOAD_CHECK(hipMemcpyAsync(ptr_host__, ptr_device__, size__,
+                               hipMemcpyDeviceToHost, stream__));
+}
+
+static inline void offloadMemcpyHtoD(void *ptr_device__, const void *ptr_host__,
+                                     const size_t size__) {
   OFFLOAD_CHECK(
-      hipMemcpyAsync(ptr1__, ptr2__, size__, hipMemcpyDeviceToHost, stream__));
+      hipMemcpy(ptr_device__, ptr_host__, size__, hipMemcpyHostToDevice));
 }
 
-static inline void offloadMemcpyHtoD(void *ptr1__, const void *ptr2__,
+static inline void offloadMemcpyDtoH(void *ptr_host__, const void *ptr_device__,
                                      const size_t size__) {
-  OFFLOAD_CHECK(hipMemcpy(ptr1__, ptr2__, size__, hipMemcpyHostToDevice));
-}
-
-static inline void offloadMemcpyDtoH(void *ptr1__, const void *ptr2__,
-                                     const size_t size__) {
-  OFFLOAD_CHECK(hipMemcpy(ptr1__, ptr2__, size__, hipMemcpyDeviceToHost));
+  OFFLOAD_CHECK(
+      hipMemcpy(ptr_host__, ptr_device__, size__, hipMemcpyDeviceToHost));
 }
 
 static inline void offloadEventCreate(offloadEvent_t *event__) {
@@ -85,8 +91,8 @@ static inline void offloadEventRecord(offloadEvent_t event__,
   OFFLOAD_CHECK(hipEventRecord(event__, stream__));
 }
 
-static inline void offloadMalloc(void *ptr__, size_t size__) {
-  OFFLOAD_CHECK(hipMalloc((void **)ptr__, size__));
+static inline void offloadMalloc(void **ptr__, size_t size__) {
+  OFFLOAD_CHECK(hipMalloc(ptr__, size__));
 }
 
 static inline void offloadFree(void *ptr__) { OFFLOAD_CHECK(hipFree(ptr__)); }
@@ -97,10 +103,6 @@ static inline void offloadStreamWaitEvent(offloadStream_t stream__,
   OFFLOAD_CHECK(hipStreamWaitEvent(stream__, event__, val__));
 }
 
-static inline void offloadSetDevice(const int dev_id__) {
-  OFFLOAD_CHECK(hipSetDevice(dev_id__));
-}
-
 static inline void offloadDeviceSynchronize() {
   OFFLOAD_CHECK(hipDeviceSynchronize());
 }
@@ -108,6 +110,14 @@ static inline void offloadDeviceSynchronize() {
 static inline void offloadMemset(void *ptr__, const int val__, size_t size__) {
   OFFLOAD_CHECK(hipMemset(ptr__, val__, size__));
 }
+
+static inline void offloadMallocHost(void **ptr__, size_t size__) {
+  OFFLOAD_CHECK(hipMallocHost(ptr__, size__));
+}
+static inline void offloadFreeHost(void *ptr__) {
+  OFFLOAD_CHECK(hipFreeHost(ptr__));
+}
+
 #ifdef __cplusplus
 }
 #endif

src/offload/offload_library.c

@@ -12,11 +12,15 @@
 #include <string.h>
 
 #include "offload_library.h"
+#include "offload_operations.h"
 
 #if defined(__OFFLOAD_PROFILING)
 #if defined(__OFFLOAD_CUDA)
+#include <cuda.h>
+#include <cuda_runtime.h>
 #include <nvToolsExt.h>
 #elif defined(__OFFLOAD_HIP) && defined(__HIP_PLATFORM_AMD__)
+#include <hip/hip_runtime_api.h>
 #include <roctracer/roctx.h>
 #endif
 #endif
@@ -135,4 +139,24 @@ void offload_mem_info(size_t *free, size_t *total) {
 #endif
 }
 
+int offload_host_malloc(void **ptr__, const size_t size__) {
+#if defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)
+  /* API checks are included in the overloading of the function */
+  offloadMallocHost(ptr__, size__);
+#else
+  *ptr__ = malloc(size__);
+  return 0;
+#endif
+  return 0;
+}
+
+int offload_host_free(void *ptr__) {
+#if defined(__OFFLOAD_CUDA) || defined(__OFFLOAD_HIP)
+  offloadFreeHost(ptr__);
+#else
+  free(ptr__);
+#endif
+  return 0;
+}
+
 // EOF

src/offload/offload_library.h

@@ -7,24 +7,15 @@
 #ifndef OFFLOAD_LIBRARY_H
 #define OFFLOAD_LIBRARY_H
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
 #if defined(__GRID_CUDA) || defined(__DBM_CUDA) || defined(__PW_CUDA)
 #define __OFFLOAD_CUDA
-#elif defined(__GRID_HIP)
+#elif defined(__GRID_HIP) || defined(__PW_HIP)
 #define __OFFLOAD_HIP
 #endif
 
-#if defined(__OFFLOAD_CUDA)
-#include <cuda.h>
-#include <cuda_runtime.h>
-#elif defined(__OFFLOAD_HIP)
-#include <hip/hip_runtime_api.h>
+#ifdef __cplusplus
+extern "C" {
 #endif
-
-#include "offload_operations.h"
 /*******************************************************************************
  * \brief Returns the number of available devices.
  * \author Ole Schuett
@@ -67,6 +58,16 @@ void offload_timestop(void);
  ******************************************************************************/
 void offload_mem_info(size_t *free, size_t *total);
 
+/*******************************************************************************
+ * \brief Allocate pinned memory (or simple malloc when there is no gpu)
+ ******************************************************************************/
+int offload_host_malloc(void **ptr__, const size_t size__);
+
+/*******************************************************************************
+ * \brief free pinned memory (or simple free when there is no gpu)
+ ******************************************************************************/
+int offload_host_free(void *ptr__);
+
 #ifdef __cplusplus
 }
 #endif

src/pw/PACKAGE

@@ -4,7 +4,7 @@
         "fpga",
         "fft",
         "../common",
-        "cuda",
+        "gpu",
         "../base",
         "../mpiwrap",
         "../offload",

src/pw/cuda/fft_cuda.h

@@ -1,40 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#ifndef FFT_CUDA_H
-#define FFT_CUDA_H
-
-/*******************************************************************************
- * \author Benjamin G Levine, Andreas Gloess
- * \par History
- *  2012/05/18                 Refacturing - original files:
- *                              - cuda_tools/cufft.h
- *                              - cuda_tools/cufft.cu
- ******************************************************************************/
-#if defined(__PW_CUDA)
-#include <cufft.h>
-
-// Double precision complex procedures
-extern "C" void fftcu_run_3d_z_(const int fsign, const int *n,
-                                const double scale, cufftDoubleComplex *data,
-                                const cudaStream_t cuda_stream);
-
-extern "C" void fftcu_run_2dm_z_(const int fsign, const int *n,
-                                 const double scale,
-                                 cufftDoubleComplex *data_in,
-                                 cufftDoubleComplex *data_out,
-                                 const cudaStream_t cuda_stream);
-
-extern "C" void fftcu_run_1dm_z_(const int fsign, const int n, const int m,
-                                 const double scale,
-                                 cufftDoubleComplex *data_in,
-                                 cufftDoubleComplex *data_out,
-                                 const cudaStream_t cuda_stream);
-
-extern "C" void fftcu_release_();
-#endif
-#endif

src/pw/cuda/fft_cuda_internal.h

@@ -1,39 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#ifndef FFT_CUDA_INTERNAL_H
-#define FFT_CUDA_INTERNAL_H
-/*******************************************************************************
- * \author Benjamin G Levine, Andreas Gloess
- * \par History
- *  2012/05/18                 Refacturing - original files:
- *                              - cuda_tools/cufft.h
- *                              - cuda_tools/cufft.cu
- ******************************************************************************/
-#if defined(__PW_CUDA)
-
-/*******************************************************************************
- * \brief   Sets up static data for FFT plan storage and reuse.
- * \author  Andreas Gloess
- * \date    2012-05-18
- * \version 0.01
- ******************************************************************************/
-static const int max_3d_plans = 30;
-static const int max_2d_plans = 0;
-static const int max_1d_plans = 30;
-
-static const int sum_plans = max_3d_plans + max_2d_plans + max_1d_plans;
-static const int max_plans = sum_plans > 1 ? sum_plans : 1;
-
-#if (__CUDACC_VER_MAJOR__ < 8)
-// configuration(s)
-#define FFT_ALIGNMENT CUFFT_COMPATIBILITY_NATIVE // potentially faster
-//#define FFT_ALIGNMENT CUFFT_COMPATIBILITY_FFTW_PADDING // the default
-#endif
-
-#endif
-#endif

src/pw/cuda/fft_cuda_utils.cu

@@ -1,67 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-/*******************************************************************************
- * \author Andreas Gloess
- ******************************************************************************/
-#if defined(__PW_CUDA)
-
-// global dependencies
-#include <cuda_runtime.h>
-#include <cufft.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-// debug flag
-#define CHECK 1
-#define VERBOSE 0
-
-// --- CODE -------------------------------------------------------------------
-
-extern void cufft_error_check(cufftResult_t cufftError, int line) {
-  int pid;
-  size_t free, total;
-  cudaError_t cErr2;
-
-  cErr2 = cudaGetLastError();
-  if (cufftError != CUFFT_SUCCESS || cErr2 != cudaSuccess) {
-    pid = getpid();
-    printf("%d CUDA FFT Error line: %d \n", pid, line);
-    switch (cufftError) {
-    case CUFFT_INVALID_PLAN:
-      printf("%d CUDA FFT1 Error (CUFFT_INVALID_PLAN)\n", pid);
-      break;
-    case CUFFT_ALLOC_FAILED:
-      printf("%d CUDA FFT1 Error (CUFFT_ALLOC_FAILED)\n", pid);
-      break;
-    case CUFFT_INVALID_VALUE:
-      printf("%d CUDA FFT1 Error (CUFFT_INVALID_VALUE)\n", pid);
-      break;
-    case CUFFT_INTERNAL_ERROR:
-      printf("%d CUDA FFT1 Error (CUFFT_INTERNAL_ERROR)\n", pid);
-      break;
-    case CUFFT_EXEC_FAILED:
-      printf("%d CUDA FFT1 Error (CUFFT_EXEC_FAILED)\n", pid);
-      break;
-    case CUFFT_INVALID_SIZE:
-      printf("%d CUDA FFT1 Error (CUFFT_INVALID_SIZE)\n", pid);
-      break;
-    default:
-      printf("%d CUDA FFT1 Error (--unimplemented--) %d %d\n", pid, cufftError,
-             cErr2);
-      break;
-    }
-    printf("%d CUDA FFT2 Error %s \n", pid, cudaGetErrorString(cErr2));
-    cudaMemGetInfo(&free, &total);
-    printf("%d Free: %zu , Total: %zu\n", pid, free, total);
-    fflush(stdout);
-    exit(-1);
-  }
-}
-
-#endif

src/pw/cuda/fft_cuda_utils.h

@@ -1,14 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#ifndef FFT_CUDA_UTILS_H
-#define FFT_CUDA_UTILS_H
-#include <cufft.h>
-
-extern void cufft_error_check(cufftResult_t cufftError, int line);
-
-#endif

src/pw/cuda/fft_cuda_z.cu

@@ -1,405 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-/*******************************************************************************
- * \author Benjamin G Levine, Andreas Gloess
- * \par History
- *  2012/05/18                 Refacturing - original files:
- *                              - cuda_tools/cufft.h
- *                              - cuda_tools/cufft.cu
- ******************************************************************************/
-#if defined(__PW_CUDA)
-
-// global dependencies
-#include <cublas.h>
-#include <cuda_runtime.h>
-#include <cufft.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-// local dependencies
-#include "fft_cuda_internal.h"
-#include "fft_cuda_utils.h"
-#include "pw_cuda_utils.h"
-
-// debug flag
-#define CHECK 1
-#define VERBOSE 0
-
-// dimensions
-static int n_plans = 0;
-static cufftHandle saved_plans[max_plans];
-static int iplandims[max_plans][5];
-
-// --- CODE -------------------------------------------------------------------
-
-/*******************************************************************************
- * \brief   Sets up and save a double precision complex 3D-FFT plan on the GPU.
- *          Saved plans are reused if they fit the requirements.
- * \author  Andreas Gloess
- * \date    2012-05-18
- * \version 0.01
- ******************************************************************************/
-void fftcu_plan3d_z(cufftHandle &plan, int &ioverflow, const int *n,
-                    const cudaStream_t cuda_stream) {
-
-  int i;
-  cufftResult_t cErr;
-
-  ioverflow = 0;
-  for (i = 0; i < n_plans; i++) {
-    if (iplandims[i][0] == 3 && iplandims[i][1] == n[0] &&
-        iplandims[i][2] == n[1] && iplandims[i][3] == n[2] &&
-        iplandims[i][4] == 0) {
-      plan = saved_plans[i];
-      return;
-    }
-  }
-
-  if (VERBOSE)
-    printf("FFT 3D (%d-%d-%d)\n", n[0], n[1], n[2]);
-  cErr = cufftPlan3d(&plan, n[2], n[1], n[0], CUFFT_Z2Z);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-  cErr = cufftSetStream(plan, cuda_stream);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#if (__CUDACC_VER_MAJOR__ < 8)
-  cErr = cufftSetCompatibilityMode(plan, FFT_ALIGNMENT);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#endif
-
-  if (n_plans < max_3d_plans) {
-    saved_plans[n_plans] = plan;
-    iplandims[n_plans][0] = 3;
-    iplandims[n_plans][1] = n[0];
-    iplandims[n_plans][2] = n[1];
-    iplandims[n_plans][3] = n[2];
-    iplandims[n_plans][4] = 0;
-    n_plans++;
-    return;
-  }
-  ioverflow = 1;
-}
-
-/*******************************************************************************
- * \brief   Sets up and save a double precision complex 2D-FFT plan on the GPU.
- *          Saved plans are reused if they fit the requirements.
- * \author  Andreas Gloess
- * \date    2012-07-16
- * \version 0.01
- ******************************************************************************/
-void fftcu_plan2dm_z(cufftHandle &plan, int &ioverflow, const int *n,
-                     const int fsign, const cudaStream_t cuda_stream) {
-
-  int i, istride, idist, ostride, odist, batch;
-  int nsize[2], inembed[2], onembed[2];
-  cufftResult_t cErr;
-
-  ioverflow = 0;
-  for (i = 0; i < n_plans; i++) {
-    if (iplandims[i][0] == 2 && iplandims[i][1] == n[0] &&
-        iplandims[i][2] == n[1] && iplandims[i][3] == n[2] &&
-        iplandims[i][4] == fsign) {
-      plan = saved_plans[i];
-      return;
-    }
-  }
-
-  nsize[0] = n[2];
-  nsize[1] = n[1];
-  inembed[0] = n[2];
-  inembed[1] = n[1];
-  onembed[0] = n[2];
-  onembed[1] = n[1];
-  batch = n[0];
-  if (fsign == +1) {
-    istride = n[0];
-    idist = 1;
-    ostride = 1;
-    odist = n[1] * n[2];
-  } else {
-    istride = 1;
-    idist = n[1] * n[2];
-    ostride = n[0];
-    odist = 1;
-  }
-
-  if (VERBOSE)
-    printf("FFT 2D (%d) (%d-%d-%d) %d %d %d %d\n", fsign, n[0], n[1], n[2],
-           istride, idist, ostride, odist);
-  cErr = cufftPlanMany(&plan, 2, nsize, inembed, istride, idist, onembed,
-                       ostride, odist, CUFFT_Z2Z, batch);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-  cErr = cufftSetStream(plan, cuda_stream);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#if (__CUDACC_VER_MAJOR__ < 8)
-  cErr = cufftSetCompatibilityMode(plan, FFT_ALIGNMENT);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#endif
-
-  if (n_plans < max_2d_plans) {
-    saved_plans[n_plans] = plan;
-    iplandims[n_plans][0] = 2;
-    iplandims[n_plans][1] = n[0];
-    iplandims[n_plans][2] = n[1];
-    iplandims[n_plans][3] = n[2];
-    iplandims[n_plans][4] = fsign;
-    n_plans++;
-    return;
-  }
-
-  ioverflow = 1;
-}
-
-/*******************************************************************************
- * \brief   Sets up and save a double precision complex 1D-FFT plan on the GPU.
- *          Saved plans are reused if they fit the requirements.
- * \author  Andreas Gloess
- * \date    2012-07-04
- * \version 0.01
- ******************************************************************************/
-void fftcu_plan1dm_z(cufftHandle &plan, int &ioverflow, const int n,
-                     const int m, const int fsign,
-                     const cudaStream_t cuda_stream) {
-
-  int i, istride, idist, ostride, odist, batch;
-  int nsize[1], inembed[1], onembed[1];
-  cufftResult_t cErr;
-
-  ioverflow = 0;
-  for (i = 0; i < n_plans; i++) {
-    if (iplandims[i][0] == 1 && iplandims[i][1] == n && iplandims[i][2] == m &&
-        iplandims[i][3] == 0 && iplandims[i][4] == fsign) {
-      plan = saved_plans[i];
-      return;
-    }
-  }
-
-  nsize[0] = n;
-  inembed[0] = 0; // is ignored, but is not allowed to be NULL pointer (for adv.
-                  // strided I/O)
-  onembed[0] = 0; // is ignored, but is not allowed to be NULL pointer (for adv.
-                  // strided I/O)
-  batch = m;
-  if (fsign == +1) {
-    istride = m;
-    idist = 1;
-    ostride = 1;
-    odist = n;
-  } else {
-    istride = 1;
-    idist = n;
-    ostride = m;
-    odist = 1;
-  }
-
-  if (VERBOSE)
-    printf("FFT 1D (%d) (%d-%d) %d %d %d %d\n", fsign, n, m, istride, idist,
-           ostride, odist);
-  cErr = cufftPlanMany(&plan, 1, nsize, inembed, istride, idist, onembed,
-                       ostride, odist, CUFFT_Z2Z, batch);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-  cErr = cufftSetStream(plan, cuda_stream);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#if (__CUDACC_VER_MAJOR__ < 8)
-  cErr = cufftSetCompatibilityMode(plan, FFT_ALIGNMENT);
-  if (CHECK)
-    cufft_error_check(cErr, __LINE__);
-#endif
-
-  if (n_plans < max_1d_plans) {
-    saved_plans[n_plans] = plan;
-    iplandims[n_plans][0] = 1;
-    iplandims[n_plans][1] = n;
-    iplandims[n_plans][2] = m;
-    iplandims[n_plans][3] = 0;
-    iplandims[n_plans][4] = fsign;
-    n_plans++;
-    return;
-  }
-
-  ioverflow = 1;
-}
-
-/*******************************************************************************
- * \brief   Performs a scaled double precision complex 3D-FFT on the GPU.
- *          Input/output is a DEVICE pointer (data).
- * \author  Andreas Gloess
- * \date    2012-05-18
- * \version 0.01
- ******************************************************************************/
-extern "C" void fftcu_run_3d_z_(const int fsign, const int *n,
-                                const double scale, cufftDoubleComplex *data,
-                                const cudaStream_t cuda_stream) {
-
-  int ioverflow, lmem;
-  cufftHandle plan;
-  cufftResult_t cErr;
-  cudaError_t cuErr;
-
-  lmem = n[0] * n[1] * n[2];
-
-  fftcu_plan3d_z(plan, ioverflow, n, cuda_stream);
-  if (fsign < 0) {
-    cErr = cufftExecZ2Z(plan, data, data, CUFFT_INVERSE);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  } else {
-    cErr = cufftExecZ2Z(plan, data, data, CUFFT_FORWARD);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-
-  if (scale != 1.0e0) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cublasDscal(2 * lmem, scale, (double *)data, 1);
-  }
-
-  if (ioverflow) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cErr = cufftDestroy(plan);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a scaled double precision complex 2D-FFT many times on
- *          the GPU.
- *          Input/output are DEVICE pointers (data_in, date_out).
- * \author  Andreas Gloess
- * \date    2012-07-16
- * \version 0.01
- ******************************************************************************/
-extern "C" void fftcu_run_2dm_z_(const int fsign, const int *n,
-                                 const double scale,
-                                 cufftDoubleComplex *data_in,
-                                 cufftDoubleComplex *data_out,
-                                 const cudaStream_t cuda_stream) {
-
-  int ioverflow, lmem;
-  cufftHandle plan;
-  cufftResult_t cErr;
-  cudaError_t cuErr;
-
-  lmem = n[0] * n[1] * n[2];
-
-  fftcu_plan2dm_z(plan, ioverflow, n, fsign, cuda_stream);
-  if (fsign < 0) {
-    cErr = cufftExecZ2Z(plan, data_in, data_out, CUFFT_INVERSE);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  } else {
-    cErr = cufftExecZ2Z(plan, data_in, data_out, CUFFT_FORWARD);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-
-  if (scale != 1.0e0) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cublasDscal(2 * lmem, scale, (double *)data_out, 1);
-  }
-
-  if (ioverflow) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cErr = cufftDestroy(plan);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a scaled double precision complex 1D-FFT many times on
- *          the GPU.
- *          Input/output are DEVICE pointers (data_in, date_out).
- * \author  Andreas Gloess
- * \date    2012-05-18
- * \version 0.01
- ******************************************************************************/
-extern "C" void fftcu_run_1dm_z_(const int fsign, const int n, const int m,
-                                 const double scale,
-                                 cufftDoubleComplex *data_in,
-                                 cufftDoubleComplex *data_out,
-                                 const cudaStream_t cuda_stream) {
-
-  int ioverflow, lmem;
-  cufftHandle plan;
-  cufftResult_t cErr;
-  cudaError_t cuErr;
-
-  lmem = n * m;
-
-  fftcu_plan1dm_z(plan, ioverflow, n, m, fsign, cuda_stream);
-  if (fsign < 0) {
-    cErr = cufftExecZ2Z(plan, data_in, data_out, CUFFT_INVERSE);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  } else {
-    cErr = cufftExecZ2Z(plan, data_in, data_out, CUFFT_FORWARD);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-
-  if (scale != 1.0e0) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cublasDscal(2 * lmem, scale, (double *)data_out, 1);
-  }
-
-  if (ioverflow) {
-    cuErr = cudaStreamSynchronize(cuda_stream);
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cErr = cufftDestroy(plan);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-}
-
-/*******************************************************************************
- * \brief   Release all stored plans.
- *
- * \author  Andreas Gloess
- * \date    2013-04-23
- * \version 0.01
- ******************************************************************************/
-extern "C" void fftcu_release_() {
-  int i;
-  cufftHandle plan;
-  cufftResult_t cErr;
-  cudaError_t cuErr;
-
-  for (i = 0; i < n_plans; i++) {
-    plan = saved_plans[i];
-    cuErr = cudaDeviceSynchronize();
-    if (CHECK)
-      pw_cuda_error_check(cuErr, __LINE__);
-    cErr = cufftDestroy(plan);
-    if (CHECK)
-      cufft_error_check(cErr, __LINE__);
-  }
-  n_plans = 0;
-}
-
-#endif

src/pw/cuda/pw_cuda.h

@@ -1,54 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#ifndef PW_CUDA_H
-#define PW_CUDA_H
-
-/*******************************************************************************
- * \author Benjamin G Levine, Andreas Gloess
- * \par History
- *  2012/05/18                 Refacturing - original files:
- *                              - cuda_tools/cuda_pw_cu.cu
- ******************************************************************************/
-#if defined(__PW_CUDA)
-
-/* Double precision complex procedures */
-extern "C" void pw_cuda_cfffg_z_(const double *din, cuDoubleComplex *zout,
-                                 const int *ghatmap, const int *npts,
-                                 const int ngpts, const double scale);
-
-extern "C" void pw_cuda_sfffc_z_(const cuDoubleComplex *zin, double *dout,
-                                 const int *ghatmap, const int *npts,
-                                 const int ngpts, const int nmaps,
-                                 const double scale);
-
-extern "C" void pw_cuda_cff_z_(const double *din, cuDoubleComplex *zout,
-                               const int *npts);
-
-extern "C" void pw_cuda_ffc_z_(const cuDoubleComplex *zin, double *dout,
-                               const int *npts);
-
-extern "C" void pw_cuda_cf_z_(const double *din, cuDoubleComplex *zout,
-                              const int *npts);
-
-extern "C" void pw_cuda_fc_z_(const cuDoubleComplex *zin, double *dout,
-                              const int *npts);
-
-extern "C" void pw_cuda_f_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                             const int dir, const int n, const int m);
-
-extern "C" void pw_cuda_fg_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                              const int *ghatmap, const int *npts,
-                              const int mmax, const int ngpts,
-                              const double scale);
-
-extern "C" void pw_cuda_sf_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                              const int *ghatmap, const int *npts,
-                              const int mmax, const int ngpts, const int nmaps,
-                              const double scale);
-#endif
-#endif

src/pw/cuda/pw_cuda_utils.cu

@@ -1,199 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#if defined(__PW_CUDA)
-
-/*******************************************************************************
- * \brief This file contains memory management routines for device memory.
- * \par History
- *      Revised: Sept. 2012, Andreas Gloess
- * \author Benjamin G Levine
- ******************************************************************************/
-
-// global dependencies
-#include <cuda_runtime.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include <unistd.h>
-
-#include "../../offload/offload_library.h"
-
-// local dependencies
-#include "fft_cuda.h"
-#include "fft_cuda_utils.h"
-
-// debug flag
-#define CHECK 1
-#define VERBOSE 0
-
-// --- CODE --------------------------------------------------------------------
-static const int nstreams = 3;
-static const int nevents = 2;
-
-cudaError_t cErr;
-static cudaStream_t *cuda_streams;
-static cudaEvent_t *cuda_events;
-static int is_configured = 0;
-
-extern void pw_cuda_error_check(cudaError_t cudaError, int line) {
-  int pid;
-  size_t free, total;
-  cudaError_t cErr2;
-
-  cErr2 = cudaGetLastError();
-  if (cudaError != cudaSuccess || cErr2 != cudaSuccess) {
-    pid = getpid();
-    printf("%d CUDA RT Error line %d\n", pid, line);
-    printf("%d CUDA RT1 Error: %s\n", pid, cudaGetErrorString(cudaError));
-    printf("%d CUDA RT2 Error: %s\n", pid, cudaGetErrorString(cErr2));
-    cudaMemGetInfo(&free, &total);
-    printf("%d Free: %zu , Total: %zu\n", pid, free, total);
-    fflush(stdout);
-    exit(-1);
-  }
-}
-
-// STREAMS INIT/GET/RELEASE
-void pw_cuda_device_streams_alloc(cudaStream_t **streams) {
-  cudaStream_t *cuda_streams_ptr;
-  cuda_streams_ptr = (cudaStream_t *)malloc(nstreams * sizeof(cudaStream_t));
-  for (int i = 0; i < nstreams; i++) {
-    cErr =
-        cudaStreamCreateWithFlags(&cuda_streams_ptr[i], cudaStreamNonBlocking);
-    if (CHECK)
-      pw_cuda_error_check(cErr, __LINE__);
-  }
-  *streams = cuda_streams_ptr;
-}
-
-extern void pw_cuda_get_streams(cudaStream_t **streams) {
-  *streams = cuda_streams;
-}
-
-void pw_cuda_device_streams_release(cudaStream_t **streams) {
-  cudaStream_t *cuda_streams_ptr;
-  cuda_streams_ptr = *streams;
-  for (int i = 0; i < nstreams; i++) {
-    cErr = cudaStreamDestroy(cuda_streams_ptr[i]);
-    if (CHECK)
-      pw_cuda_error_check(cErr, __LINE__);
-  }
-  free(cuda_streams_ptr);
-  cuda_streams_ptr = NULL;
-}
-
-// EVENTS INIT/GET/RELEASE
-void pw_cuda_device_events_alloc(cudaEvent_t **events) {
-  cudaEvent_t *cuda_events_ptr;
-  cuda_events_ptr = (cudaEvent_t *)malloc(nevents * sizeof(cudaEvent_t));
-  for (int i = 0; i < nevents; i++) {
-    cErr =
-        cudaEventCreateWithFlags(&cuda_events_ptr[i], cudaEventDisableTiming);
-    // cErr = cudaEventCreateWithFlags(&(cuda_events_ptr[i]), cudaEventDefault);
-    // cErr = cudaEventCreateWithFlags(&(cuda_events_ptr[i]),
-    // cudaEventBlockingSync);
-    if (CHECK)
-      pw_cuda_error_check(cErr, __LINE__);
-  }
-  *events = cuda_events_ptr;
-}
-
-extern void pw_cuda_get_events(cudaEvent_t **events) { *events = cuda_events; }
-
-void pw_cuda_device_events_release(cudaEvent_t **events) {
-  cudaEvent_t *cuda_events_ptr;
-  cuda_events_ptr = *events;
-  for (int i = 0; i < nevents; i++) {
-    cErr = cudaEventDestroy(cuda_events_ptr[i]);
-    if (CHECK)
-      pw_cuda_error_check(cErr, __LINE__);
-  }
-  free(cuda_events_ptr);
-  cuda_events_ptr = NULL;
-}
-
-// MEMORY ALLOC/RELEASE
-extern void pw_cuda_device_mem_alloc(int **ptr, int n) {
-  cErr = cudaMalloc((void **)ptr, (size_t)sizeof(int) * n);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-}
-
-extern void pw_cuda_device_mem_alloc(float **ptr, int n) {
-  cErr = cudaMalloc((void **)ptr, (size_t)sizeof(float) * n);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-}
-
-extern void pw_cuda_device_mem_alloc(double **ptr, int n) {
-  cErr = cudaMalloc((void **)ptr, (size_t)sizeof(double) * n);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-}
-
-extern void pw_cuda_device_mem_free(int **ptr) {
-  cErr = cudaFree((void *)*ptr);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  *ptr = NULL;
-}
-
-extern void pw_cuda_device_mem_free(float **ptr) {
-  cErr = cudaFree((void *)*ptr);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  *ptr = NULL;
-}
-
-extern void pw_cuda_device_mem_free(double **ptr) {
-  cErr = cudaFree((void *)*ptr);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  *ptr = NULL;
-}
-
-// INIT/RELEASE
-extern "C" int pw_cuda_init() {
-  if (is_configured == 0) {
-    int version;
-    cufftResult_t cufftErr;
-    offload_set_device();
-    pw_cuda_device_streams_alloc(&cuda_streams);
-    pw_cuda_device_events_alloc(&cuda_events);
-    is_configured = 1;
-    cufftErr = cufftGetVersion(&version);
-    if (CHECK)
-      cufft_error_check(cufftErr, __LINE__);
-    if (version == 7000) {
-#if defined(__HAS_PATCHED_CUFFT_70)
-//      message would be printed to all ranks, if the user has
-//      __HAS_PATCHED_CUFFT_70 she's supposed to know what she's doing
-//      printf("CUFFT 7.0 enabled on user request
-//      (-D__HAS_PATCHED_CUFFT_70).\n"); printf("Please ensure that CUFFT is
-//      patched (libcufft.so.x.y.z, libcufftw.so.x,y,z; x.y.z >= 7.0.35).\n");
-#else
-      printf("CUFFT 7.0 disabled due to an unresolved bug.\n");
-      printf("Please upgrade to CUDA 7.5 or later or apply CUFFT patch (see "
-             "-D__HAS_PATCHED_CUFFT_70).\n");
-      return -1;
-#endif
-    }
-  }
-  return 0;
-}
-
-extern "C" void pw_cuda_finalize() {
-  if (is_configured == 1) {
-    offload_set_device();
-    fftcu_release_();
-    pw_cuda_device_streams_release(&cuda_streams);
-    pw_cuda_device_events_release(&cuda_events);
-    is_configured = 0;
-  }
-}
-
-#endif

src/pw/cuda/pw_cuda_utils.h

@@ -1,35 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-#ifndef PW_CUDA_UTILS_H
-#define PW_CUDA_UTILS_H
-
-extern void pw_cuda_error_check(cudaError_t cudaError, int line);
-
-// STREAMS INIT/GET/RELEASE
-extern void pw_cuda_device_streams_alloc(cudaStream_t **streams);
-extern void pw_cuda_get_streams(cudaStream_t **streams);
-extern void pw_cuda_device_streams_release(cudaStream_t **streams);
-
-// EVENTS INIT/GET/RELEASE
-extern void pw_cuda_device_events_alloc(cudaEvent_t **events);
-extern void pw_cuda_get_events(cudaEvent_t **events);
-extern void pw_cuda_device_events_release(cudaEvent_t **events);
-
-// MEMORY ALLOC/RELEASE
-extern void pw_cuda_device_mem_alloc(int **ptr, int n);
-extern void pw_cuda_device_mem_alloc(float **ptr, int n);
-extern void pw_cuda_device_mem_alloc(double **ptr, int n);
-extern void pw_cuda_device_mem_free(int **ptr);
-extern void pw_cuda_device_mem_free(float **ptr);
-extern void pw_cuda_device_mem_free(double **ptr);
-
-// DEVICE INIT/RELEASE
-extern "C" int pw_cuda_init();
-extern "C" void pw_cuda_release();
-
-#endif

src/pw/cuda/pw_cuda_z.cu

@@ -1,962 +0,0 @@
-/*----------------------------------------------------------------------------*/
-/*  CP2K: A general program to perform molecular dynamics simulations         */
-/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
-/*                                                                            */
-/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
-/*----------------------------------------------------------------------------*/
-
-/*******************************************************************************
- * \author Benjamin G Levine, Andreas Gloess
- * \par History
- *  2012/05/18                 Refacturing - original files:
- *                              - cuda_tools/cuda_pw_cu.cu
- ******************************************************************************/
-#if defined(__PW_CUDA)
-
-// global dependencies
-#include <cublas.h>
-#include <cuda_runtime.h>
-#include <cufft.h>
-#include <stdio.h>
-
-#include "../../offload/offload_library.h"
-
-// local dependencies
-#include "fft_cuda.h"
-#include "pw_cuda_utils.h"
-
-// debug flag
-#define CHECK 1
-
-// configuration(s)
-#define NTHREADS 32
-#define MAXTHREADS 1024
-#define MAXGRIDX 65535
-
-// helper routine(s)
-void get_grid_params(const int ngpts, const int blocksize,
-                     dim3 &threadsPerBlock, dim3 &blocksPerGrid) {
-  int blocks;
-  if (blocksize <= MAXTHREADS) {
-    threadsPerBlock.x = blocksize;
-  } else {
-    threadsPerBlock.x = MAXTHREADS;
-    printf("WARNING: Number of threads per block (x) is too large!\n");
-    printf("WARNING: Number of threads per block (x) is set to: %d.\n",
-           MAXTHREADS);
-  }
-  threadsPerBlock.y = 1;
-  threadsPerBlock.z = 1;
-  blocks = (ngpts + threadsPerBlock.x - 1) / threadsPerBlock.x;
-  blocksPerGrid.x = (int)ceil(sqrt((double)blocks));
-  if (blocksPerGrid.x > MAXGRIDX) {
-    printf("CUDA: Not allowed grid dimensions!\n");
-    exit(1);
-  }
-  blocksPerGrid.y = (int)rint(sqrt((double)blocks));
-  blocksPerGrid.z = 1;
-}
-
-void cudaStreamBarrier(cudaStream_t cuda_stream) {
-  cudaError_t cErr;
-
-  // stop on previous errors
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  // might result in endless loop
-  do {
-    cErr = cudaStreamQuery(cuda_stream);
-  } while (cErr != cudaSuccess);
-}
-
-// --- CODE -------------------------------------------------------------------
-
-/*******************************************************************************
- * \brief   Performs a out-of-place copy of a double precision vector (first
- *          half filled) into a double precision complex vector on the GPU.
- *          It requires a global double precision vector 'zout' of length '2n'.
- *          [memory (shared):  none Byte
- *           memory (private): 4 Byte
- *           memory (global):  16*n Byte]
- *          n - size of double precision input vector
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-__global__ void pw_copy_rc_cu_z(const double *din, double *zout, const int n) {
-  const int igpt =
-      (gridDim.x * blockIdx.y + blockIdx.x) * blockDim.x + threadIdx.x;
-
-  if (igpt < n) {
-    zout[2 * igpt] = din[igpt];
-    zout[2 * igpt + 1] = 0.0e0;
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a out-of-place copy of a double precision complex vector
- *          (real part) into a double precision vector on the GPU.
- *          It requires a global double precision vector 'dout' of length 'n'.
- *          [memory (shared):  none Byte
- *           memory (private): 4 Byte
- *           memory (global):  16*n Byte]
- *          n - size of double precision output vector
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-__global__ void pw_copy_cr_cu_z(const double *zin, double *dout, const int n) {
-  const int igpt =
-      (gridDim.x * blockIdx.y + blockIdx.x) * blockDim.x + threadIdx.x;
-
-  if (igpt < n) {
-    dout[igpt] = zin[2 * igpt];
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) gather and scale on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-__global__ void pw_gather_cu_z(double *pwcc, const double *c,
-                               const double scale, const int ngpts,
-                               const int *ghatmap) {
-
-  const int igpt =
-      (gridDim.x * blockIdx.y + blockIdx.x) * blockDim.x + threadIdx.x;
-
-  if (igpt < ngpts) {
-    pwcc[2 * igpt] = scale * c[2 * ghatmap[igpt]];
-    pwcc[2 * igpt + 1] = scale * c[2 * ghatmap[igpt] + 1];
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) scatter and scale on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-__global__ void pw_scatter_cu_z(double *c, const double *pwcc,
-                                const double scale, const int ngpts,
-                                const int nmaps, const int *ghatmap) {
-
-  const int igpt =
-      (gridDim.x * blockIdx.y + blockIdx.x) * blockDim.x + threadIdx.x;
-
-  if (igpt < ngpts) {
-    c[2 * ghatmap[igpt]] = scale * pwcc[2 * igpt];
-    c[2 * ghatmap[igpt] + 1] = scale * pwcc[2 * igpt + 1];
-    if (nmaps == 2) {
-      c[2 * ghatmap[igpt + ngpts]] = scale * pwcc[2 * igpt];
-      c[2 * ghatmap[igpt + ngpts] + 1] = -scale * pwcc[2 * igpt + 1];
-    }
-  }
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) FFT, followed by a (double
- *          precision complex) gather, on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_cfffg_z_(const double *din, cuDoubleComplex *zout,
-                                 const int *ghatmap, const int *npts,
-                                 const int ngpts, const double scale) {
-  double *ptr_1, *ptr_2;
-  int *ghatmap_dev;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0 || ngpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts); //? (ngpts)
-  pw_cuda_device_mem_alloc(&ghatmap_dev, ngpts);
-
-  // convert the real (host) pointer 'din' into a complex (device)
-  // pointer 'ptr_1'
-  // copy to device (NOTE: only first half of ptr_1 is written!)
-  cErr = cudaMemcpyAsync(ptr_1, din, sizeof(double) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // CUDA blocking for pw_copy_rc (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // real to complex blow-up
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  pw_copy_rc_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // copy gather map array from host to the device
-  cErr = cudaMemcpyAsync(ghatmap_dev, ghatmap, sizeof(int) * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  fftcu_run_3d_z_(+1, npts, 1.0e0, (cufftDoubleComplex *)ptr_2,
-                  cuda_streams[1]);
-
-  // CUDA blocking for gather (currently only 2-D grid)
-  get_grid_params(ngpts, NTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // gather on the GPU
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  pw_gather_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, scale, ngpts, ghatmap_dev);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(zout, ptr_1, sizeof(cuDoubleComplex) * ngpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-  pw_cuda_device_mem_free(&ghatmap_dev);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) scatter, followed by a
- *          (double precision complex) FFT, on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_sfffc_z_(const cuDoubleComplex *zin, double *dout,
-                                 const int *ghatmap, const int *npts,
-                                 const int ngpts, const int nmaps,
-                                 const double scale) {
-  double *ptr_1, *ptr_2;
-  int *ghatmap_dev;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0 || ngpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ghatmap_dev, nmaps * ngpts);
-
-  // copy all arrays from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(ghatmap_dev, ghatmap, sizeof(int) * nmaps * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // CUDA blocking for scatter (currently only 2-D grid)
-  get_grid_params(ngpts, NTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // scatter on the GPU
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemsetAsync(
-      ptr_2, 0, sizeof(double) * 2 * nrpts,
-      cuda_streams[1]); // we need to do this only if spherical cut-off is used!
-  if (CHECK)
-    pw_cuda_error_check(
-        cErr, __LINE__); // but it turns out to be performance irrelevant
-  pw_scatter_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_2, ptr_1, scale, ngpts, nmaps, ghatmap_dev);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  fftcu_run_3d_z_(-1, npts, 1.0e0, (cufftDoubleComplex *)ptr_2,
-                  cuda_streams[1]);
-
-  // CUDA blocking for pw_copy_cr (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // convert the complex (device) pointer 'ptr_2' into a real (host)
-  // pointer 'dout' (NOTE: Only first half of ptr_1 is written!)
-  pw_copy_cr_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_2, ptr_1, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(dout, ptr_1, sizeof(double) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-  pw_cuda_device_mem_free(&ghatmap_dev);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double to complex double) blow-up and a (double
- *          precision complex) 2D-FFT on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_cff_z_(const double *din, cuDoubleComplex *zout,
-                               const int *npts) {
-  double *ptr_1, *ptr_2;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers for:
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-
-  // convert the real (host) pointer 'din' into a complex (device)
-  // pointer 'ptr_in'
-  // copy all arrays from host to the device (NOTE: Only first half of ptr_1 is
-  // written!)
-  cErr = cudaMemcpyAsync(ptr_1, din, sizeof(double) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // CUDA blocking for pw_copy_rc (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // real to complex blow-up
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  pw_copy_rc_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  // NOTE: the following works, but CUDA does 2D-FFT in C-shaped (not optimal)
-  // order fftcu_run_2dm_z_(1, npts, 1.0e0, (cufftDoubleComplex *) ptr_2,
-  // (cufftDoubleComplex *) ptr_1, cuda_streams[1]);
-  fftcu_run_1dm_z_(1, npts[2], npts[0] * npts[1], 1.0e0,
-                   (cufftDoubleComplex *)ptr_2, (cufftDoubleComplex *)ptr_1,
-                   cuda_streams[1]);
-  fftcu_run_1dm_z_(1, npts[1], npts[0] * npts[2], 1.0e0,
-                   (cufftDoubleComplex *)ptr_1, (cufftDoubleComplex *)ptr_2,
-                   cuda_streams[1]);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  // cErr = cudaMemcpyAsync(zout, ptr_1, sizeof(cuDoubleComplex) * nrpts,
-  // cudaMemcpyDeviceToHost, cuda_streams[2]);
-  cErr = cudaMemcpyAsync(zout, ptr_2, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) 2D-FFT and a (double complex
- *          to double) shrink-down on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_ffc_z_(const cuDoubleComplex *zin, double *dout,
-                               const int *npts) {
-  double *ptr_1, *ptr_2;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers for:
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-
-  // copy input data from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_stream[1])
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  // fftcu_run_2dm_z_(-1, npts, 1.0e0, (cufftDoubleComplex *) ptr_1,
-  // (cufftDoubleComplex *) ptr_2, cuda_streams[1]);
-  fftcu_run_1dm_z_(-1, npts[1], npts[0] * npts[2], 1.0e0,
-                   (cufftDoubleComplex *)ptr_1, (cufftDoubleComplex *)ptr_2,
-                   cuda_streams[1]);
-  fftcu_run_1dm_z_(-1, npts[2], npts[0] * npts[1], 1.0e0,
-                   (cufftDoubleComplex *)ptr_2, (cufftDoubleComplex *)ptr_1,
-                   cuda_streams[1]);
-
-  // CUDA blocking for pw_copy_cr (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // convert the complex (device) pointer 'ptr_1' into a real (host)
-  // pointer 'dout' (NOTE: Only first half of ptr_2 is written!)
-  // pw_copy_cr_cu_z<<<blocksPerGrid, threadsPerBlock, 0,
-  // cuda_streams[1]>>>(ptr_2, ptr_1, nrpts);
-  pw_copy_cr_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  // cErr = cudaMemcpyAsync(dout, ptr_1, sizeof(double) * nrpts,
-  // cudaMemcpyDeviceToHost, cuda_streams[2]);
-  cErr = cudaMemcpyAsync(dout, ptr_2, sizeof(double) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double to complex double) blow-up and a (double
- *          precision complex) 1D-FFT on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_cf_z_(const double *din, cuDoubleComplex *zout,
-                              const int *npts) {
-  double *ptr_1, *ptr_2;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers for:
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-
-  // convert the real (host) pointer 'din' into a complex (device)
-  // pointer 'ptr_2' (NOTE: Only first half of ptr_1 is written!)
-  // copy all arrays from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, din, sizeof(double) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // CUDA blocking for pw_copy_rc (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // real to complex blow-up
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  pw_copy_rc_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  fftcu_run_1dm_z_(1, npts[2], npts[0] * npts[1], 1.0e0,
-                   (cufftDoubleComplex *)ptr_2, (cufftDoubleComplex *)ptr_1,
-                   cuda_streams[1]);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(zout, ptr_1, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) 1D-FFT and a (double complex
- *          to double) shrink-down on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_fc_z_(const cuDoubleComplex *zin, double *dout,
-                              const int *npts) {
-  double *ptr_1, *ptr_2;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * npts[1] * npts[2];
-  if (nrpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers for:
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-
-  // copy input data from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_stream[1])
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  fftcu_run_1dm_z_(-1, npts[2], npts[0] * npts[1], 1.0e0,
-                   (cufftDoubleComplex *)ptr_1, (cufftDoubleComplex *)ptr_2,
-                   cuda_streams[1]);
-
-  // CUDA blocking for pw_copy_cr (currently only 2-D grid)
-  get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // convert the complex (device) pointer 'ptr_2' into a real (host)
-  // pointer 'dout' (NOTE: Only first half of ptr_1 is written!)
-  pw_copy_cr_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_2, ptr_1, nrpts);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(dout, ptr_1, sizeof(double) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) 1D-FFT on the GPU.
- * \author  Andreas Gloess
- * \date    2013-05-01
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_f_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                             const int dir, const int n, const int m) {
-  double *ptr_1, *ptr_2;
-  int nrpts;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of complex arrays
-  nrpts = n * m;
-  if (nrpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers for:
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-
-  // copy input data from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_stream[1])
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  fftcu_run_1dm_z_(dir, n, m, 1.0e0, (cufftDoubleComplex *)ptr_1,
-                   (cufftDoubleComplex *)ptr_2, cuda_streams[1]);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(zout, ptr_2, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) 1D-FFT, followed by a (double
- *          precision complex) gather, on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_fg_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                              const int *ghatmap, const int *npts,
-                              const int mmax, const int ngpts,
-                              const double scale) {
-  double *ptr_1, *ptr_2;
-  int *ghatmap_dev;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * mmax;
-  if (nrpts == 0 || ngpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ghatmap_dev, ngpts);
-
-  // transfer gather data from host to device
-  cErr = cudaMemcpyAsync(ghatmap_dev, ghatmap, sizeof(int) * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // transfer input data from host to device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  fftcu_run_1dm_z_(1, npts[0], mmax, 1.0e0, (cufftDoubleComplex *)ptr_1,
-                   (cufftDoubleComplex *)ptr_2, cuda_streams[1]);
-
-  // CUDA blocking for gather (currently only 2-D grid)
-  get_grid_params(ngpts, NTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // gather on the GPU
-  pw_gather_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_1, ptr_2, scale, ngpts, ghatmap_dev);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // transfer results from device to host
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(zout, ptr_1, sizeof(cuDoubleComplex) * ngpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-  pw_cuda_device_mem_free(&ghatmap_dev);
-}
-
-/*******************************************************************************
- * \brief   Performs a (double precision complex) scatter, followed by a
- *          (double precision complex) 1D-FFT, on the GPU.
- * \author  Andreas Gloess
- * \date    2013-03-07
- * \version 0.01
- ******************************************************************************/
-extern "C" void pw_cuda_sf_z_(const cuDoubleComplex *zin, cuDoubleComplex *zout,
-                              const int *ghatmap, const int *npts,
-                              const int mmax, const int ngpts, const int nmaps,
-                              const double scale) {
-  double *ptr_1, *ptr_2;
-  int *ghatmap_dev;
-  int nrpts;
-  dim3 blocksPerGrid, threadsPerBlock;
-  cudaStream_t *cuda_streams;
-  cudaEvent_t *cuda_events;
-  cudaError_t cErr;
-
-  // dimensions of double and complex arrays
-  nrpts = npts[0] * mmax;
-  if (nrpts == 0 || ngpts == 0)
-    return;
-
-  // get streams and events
-  offload_set_device();
-  pw_cuda_get_streams(&cuda_streams);
-  pw_cuda_get_events(&cuda_events);
-
-  // get device memory pointers
-  pw_cuda_device_mem_alloc(&ptr_1, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ptr_2, 2 * nrpts);
-  pw_cuda_device_mem_alloc(&ghatmap_dev, nmaps * ngpts);
-
-  // transfer input data from host to the device
-  cErr = cudaMemcpyAsync(ptr_1, zin, sizeof(cuDoubleComplex) * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // transfer scatter data from host to device
-  cErr = cudaMemcpyAsync(ghatmap_dev, ghatmap, sizeof(int) * nmaps * ngpts,
-                         cudaMemcpyHostToDevice, cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[0], cuda_streams[0]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // CUDA blocking for scatter (currently only 2-D grid)
-  get_grid_params(ngpts, NTHREADS, threadsPerBlock, blocksPerGrid);
-
-  // scatter on the GPU
-  cErr = cudaStreamWaitEvent(cuda_streams[1], cuda_events[0], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemsetAsync(
-      ptr_2, 0, sizeof(double) * 2 * nrpts,
-      cuda_streams[1]); // we need to do this only if spherical cut-off is used!
-  if (CHECK)
-    pw_cuda_error_check(
-        cErr, __LINE__); // but it turns out to be performance irrelevant
-  pw_scatter_cu_z<<<blocksPerGrid, threadsPerBlock, 0, cuda_streams[1]>>>(
-      ptr_2, ptr_1, scale, ngpts, nmaps, ghatmap_dev);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // fft on the GPU (cuda_streams[1])
-  fftcu_run_1dm_z_(-1, npts[0], mmax, 1.0e0, (cufftDoubleComplex *)ptr_2,
-                   (cufftDoubleComplex *)ptr_1, cuda_streams[1]);
-  cErr = cudaGetLastError();
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaEventRecord(cuda_events[1], cuda_streams[1]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // get results from device memory
-  cErr = cudaStreamWaitEvent(cuda_streams[2], cuda_events[1], 0);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-  cErr = cudaMemcpyAsync(zout, ptr_1, sizeof(cuDoubleComplex) * nrpts,
-                         cudaMemcpyDeviceToHost, cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // synchronize with respect to host
-  cErr = cudaStreamSynchronize(cuda_streams[2]);
-  if (CHECK)
-    pw_cuda_error_check(cErr, __LINE__);
-
-  // release memory stack
-  pw_cuda_device_mem_free(&ptr_1);
-  pw_cuda_device_mem_free(&ptr_2);
-  pw_cuda_device_mem_free(&ghatmap_dev);
-}
-#endif

src/pw/fft_tools.F

@@ -44,6 +44,8 @@ MODULE fft_tools
    USE message_passing,                 ONLY: &
         mp_alltoall, mp_cart_coords, mp_cart_rank, mp_cart_sub, mp_comm_compare, mp_comm_free, &
         mp_comm_null, mp_environ, mp_irecv, mp_isend, mp_rank_compare, mp_sum, mp_sync, mp_waitall
+   USE offload_api,                     ONLY: offload_free_pinned_mem,&
+                                              offload_malloc_pinned_mem
 
 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
 
@@ -147,7 +149,7 @@ MODULE fft_tools
    LOGICAL, SAVE :: use_fftsg_sizes = .TRUE.
    INTEGER, SAVE :: fft_plan_style = 1
 
-   ! these are only needed for pw_methods_cuda (-D__PW_CUDA)
+   ! these are only needed for pw_methods_cuda (-D__PW_GPU)
    PUBLIC :: get_fft_scratch, release_fft_scratch
    PUBLIC :: cube_transpose_1, cube_transpose_2
    PUBLIC :: yz_to_x, x_to_yz, xz_to_yz, yz_to_xz
@@ -158,26 +160,6 @@ MODULE fft_tools
       MODULE PROCEDURE fft3d_s, fft3d_ps, fft3d_pb
    END INTERFACE
 
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
-   INTERFACE
-      INTEGER(KIND=C_INT) FUNCTION cudaFreeHost(buffer) BIND(C, name="cudaFreeHost")
-         IMPORT
-         IMPLICIT NONE
-         TYPE(C_PTR), VALUE :: buffer
-      END FUNCTION cudaFreeHost
-   END INTERFACE
-
-   INTERFACE
-      INTEGER(KIND=C_INT) FUNCTION cudaHostAlloc(buffer, length, flag) BIND(C, name="cudaHostAlloc")
-         IMPORT
-         IMPLICIT NONE
-         TYPE(C_PTR)                   :: buffer
-         INTEGER(KIND=C_SIZE_T), VALUE :: length
-         INTEGER(KIND=C_INT), VALUE    :: flag
-      END FUNCTION cudaHostAlloc
-   END INTERFACE
-#endif
-
 ! **************************************************************************************************
 
 CONTAINS
@@ -2609,7 +2591,7 @@ CONTAINS
    SUBROUTINE deallocate_fft_scratch_type(fft_scratch)
       TYPE(fft_scratch_type), INTENT(INOUT)    :: fft_scratch
 
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU )
       INTEGER                   :: ierr
       COMPLEX(KIND=dp), POINTER :: dummy_ptr_z
 #endif
@@ -2625,33 +2607,33 @@ CONTAINS
          DEALLOCATE (fft_scratch%p1buf)
       END IF
       IF (ASSOCIATED(fft_scratch%p2buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU )
          dummy_ptr_z => fft_scratch%p2buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%p2buf)
 #endif
       END IF
       IF (ASSOCIATED(fft_scratch%p3buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%p3buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%p3buf)
 #endif
       END IF
       IF (ASSOCIATED(fft_scratch%p4buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%p4buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%p4buf)
 #endif
       END IF
       IF (ASSOCIATED(fft_scratch%p5buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%p5buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%p5buf)
 #endif
@@ -2660,17 +2642,17 @@ CONTAINS
          DEALLOCATE (fft_scratch%p6buf)
       END IF
       IF (ASSOCIATED(fft_scratch%p7buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%p7buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%p7buf)
 #endif
       END IF
       IF (ASSOCIATED(fft_scratch%r1buf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%r1buf(1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%r1buf)
 #endif
@@ -2679,9 +2661,9 @@ CONTAINS
          DEALLOCATE (fft_scratch%r2buf)
       END IF
       IF (ASSOCIATED(fft_scratch%tbuf)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
          dummy_ptr_z => fft_scratch%tbuf(1, 1, 1)
-         ierr = cudaFreeHost(c_loc(dummy_ptr_z))
+         ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
 #else
          DEALLOCATE (fft_scratch%tbuf)
 #endif
@@ -2900,7 +2882,7 @@ CONTAINS
       TYPE(fft_scratch_pool_type), POINTER     :: fft_scratch_current, &
                                                   fft_scratch_last, &
                                                   fft_scratch_new
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
       INTEGER                   :: ierr
       INTEGER(KIND=C_SIZE_T)    :: length
       TYPE(C_PTR)               :: cptr_r1buf, cptr_tbuf, &
@@ -3113,13 +3095,13 @@ CONTAINS
                np = fft_sizes%numtask
                nmray = fft_sizes%nmray
                nyzray = fft_sizes%nyzray
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
                length = INT(2*dp_size*MAX(mmax, 1)*MAX(lmax, 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_r1buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_r1buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_r1buf, fft_scratch_new%fft_scratch%r1buf, (/MAX(mmax, 1), MAX(lmax, 1)/))
                length = INT(2*dp_size*MAX(ny, 1)*MAX(nz, 1)*MAX(nx, 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_tbuf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_tbuf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_tbuf, fft_scratch_new%fft_scratch%tbuf, (/MAX(ny, 1), MAX(nz, 1), MAX(nx, 1)/))
 #else
@@ -3167,25 +3149,25 @@ CONTAINS
                nbz = fft_sizes%nbz
                ALLOCATE (fft_scratch_new%fft_scratch%p1buf(mx1*my1, n(3)))
                ALLOCATE (fft_scratch_new%fft_scratch%p6buf(lg, mg))
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU )
                length = INT(2*dp_size*MAX(n(3), 1)*MAX(mx1*my1, 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_p2buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_p2buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_p2buf, fft_scratch_new%fft_scratch%p2buf, (/MAX(n(3), 1), MAX(mx1*my1, 1)/))
                length = INT(2*dp_size*MAX(mx2*mz2, 1)*MAX(n(2), 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_p3buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_p3buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_p3buf, fft_scratch_new%fft_scratch%p3buf, (/MAX(mx2*mz2, 1), MAX(n(2), 1)/))
                length = INT(2*dp_size*MAX(n(2), 1)*MAX(mx2*mz2, 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_p4buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_p4buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_p4buf, fft_scratch_new%fft_scratch%p4buf, (/MAX(n(2), 1), MAX(mx2*mz2, 1)/))
                length = INT(2*dp_size*MAX(nyzray, 1)*MAX(n(1), 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_p5buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_p5buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_p5buf, fft_scratch_new%fft_scratch%p5buf, (/MAX(nyzray, 1), MAX(n(1), 1)/))
                length = INT(2*dp_size*MAX(mg, 1)*MAX(lg, 1), KIND=C_SIZE_T)
-               ierr = cudaHostAlloc(cptr_p7buf, length, cudaHostAllocDefault)
+               ierr = offload_malloc_pinned_mem(cptr_p7buf, length)
                CPASSERT(ierr == 0)
                CALL c_f_pointer(cptr_p7buf, fft_scratch_new%fft_scratch%p7buf, (/MAX(mg, 1), MAX(lg, 1)/))
 #else
@@ -3438,4 +3420,3 @@ CONTAINS
    END SUBROUTINE is_equal
 
 END MODULE fft_tools
-

src/pw/gpu/PACKAGE

@@ -0,0 +1,5 @@
+{
+    "description": "fft on gpu",
+    "archive": "libcp2kpwgpu",
+    "requires": ["../../offload", "hip", "cuda",],
+}

src/pw/gpu/cuda/PACKAGE

@@ -1,5 +1,5 @@
 {
     "description": "CUDA FFT acceleration",
-    "requires": ["../../offload"],
+    "requires": ["../../../offload", "../kernels"],
     "archive": "libcp2kpwcuda",
 }

src/pw/gpu/cuda/cuda_fft_private_header.hpp

@@ -0,0 +1,284 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+#ifndef CUDA_FFT_PRIVATE_HEADER_HPP
+#define CUDA_FFT_PRIVATE_HEADER_HPP
+
+#include <array>
+#include <cstdio>
+#include <cublas_v2.h>
+#include <cuda.h>
+#include <cufft.h>
+#include <map>
+#include <string>
+#include <unistd.h>
+#include <vector>
+
+#include "../../../offload/offload_library.h"
+#include "../../../offload/offload_operations.h"
+
+typedef cufftDoubleComplex pw_complex_type;
+typedef cublasHandle_t blasHandle_t;
+
+enum fft_direction {
+  FFT_FORWARD = CUFFT_FORWARD,
+  FFT_BACKWARD = CUFFT_INVERSE,
+  FFT_UNKOWN
+};
+class fft_plan;
+
+static void fft_error_check(cufftResult cufftError, int line,
+                            const std::string filename) {
+  int pid;
+
+  std::map<cufftResult, std::string> error_messages_ = {
+      {CUFFT_SUCCESS, "CUFFT_SUCCESS"},
+      {CUFFT_INVALID_PLAN, "CUFFT_INVALID_PLAN"},
+      {CUFFT_ALLOC_FAILED, "CUFFT_ALLOC_FAILED"},
+      {CUFFT_INVALID_TYPE, "CUFFT_INVALID_TYPE"},
+      {CUFFT_INVALID_VALUE, "CUFFT_INVALID_VALUE"},
+      {CUFFT_INTERNAL_ERROR, "CUFFT_INTERNAL_ERROR"},
+      {CUFFT_EXEC_FAILED, "CUFFT_EXEC_FAILED"},
+      {CUFFT_SETUP_FAILED, "CUFFT_SETUP_FAILED"},
+      {CUFFT_INVALID_SIZE, "CUFFT_INVALID_SIZE"},
+      {CUFFT_INCOMPLETE_PARAMETER_LIST, "CUFFT_INCOMPLETE_PARAMETER_LIST"},
+      {CUFFT_INVALID_DEVICE, "CUFFT_INVALID_DEVICE"},
+      {CUFFT_PARSE_ERROR, "CUFFT_PARSE_ERROR"},
+      {CUFFT_NO_WORKSPACE, "CUFFT_NO_WORKSPACE"},
+      {CUFFT_NOT_IMPLEMENTED, "CUFFT_NOT_IMPLEMENTED"},
+      {CUFFT_NOT_SUPPORTED, "CUFFT_NOT_SUPPORTED"}};
+
+  if (cufftError != CUFFT_SUCCESS) {
+    pid = getpid();
+    printf("%d CUDA FFT Error line: %s %d \n", pid, filename.c_str(), line);
+    printf("%d CUDA FFT Error (%s)\n", pid,
+           error_messages_.at(cufftError).c_str());
+    fflush(stdout);
+    abort();
+  }
+}
+
+static void blas_error_check(cublasStatus_t cublasError, int line,
+                             const std::string filename) {
+  int pid;
+
+  std::map<cublasStatus_t, std::string> error_messages_ = {
+      {CUBLAS_STATUS_SUCCESS, "CUBLAS_STATUS_SUCCESS"},
+      {CUBLAS_STATUS_NOT_INITIALIZED, "CUBLAS_STATUS_NOT_INITIALIZED"},
+      {CUBLAS_STATUS_ALLOC_FAILED, "CUBLAS_STATUS_ALLOC_FAILED"},
+      {CUBLAS_STATUS_INVALID_VALUE, "CUBLAS_STATUS_INVALID_VALUE"},
+      {CUBLAS_STATUS_MAPPING_ERROR, "CUBLAS_STATUS_MAPPING_ERROR"},
+      {CUBLAS_STATUS_EXECUTION_FAILED, "CUBLAS_STATUS_EXECUTION_FAILED"},
+      {CUBLAS_STATUS_INTERNAL_ERROR, "CUBLAS_STATUS_INTERNAL_ERROR"},
+      {CUBLAS_STATUS_NOT_SUPPORTED, "CUBLAS_STATUS_NOT_SUPPORTED"},
+      {CUBLAS_STATUS_ARCH_MISMATCH, "CUBLAS_STATUS_ARCH_MISMATCH"}};
+
+  if (cublasError != CUBLAS_STATUS_SUCCESS) {
+    pid = getpid();
+    printf("%d CUDA BLAS Error file: %s line: %d \n", pid, filename.c_str(),
+           line);
+    printf("%d CUDA BLAS Error (%s)\n", pid,
+           error_messages_.at(cublasError).c_str());
+    fflush(stdout);
+    abort();
+  }
+}
+
+static void error_check(cudaError_t cudaError, int line,
+                        const std::string filename) {
+  int pid;
+  if (!(cudaError == cudaSuccess || cudaError == cudaErrorNotReady)) {
+    pid = getpid();
+    printf("%d CUDA RT Error line %s %d\n", pid, filename.c_str(), line);
+    printf("%d CUDA RT1 Error: %s\n", pid, cudaGetErrorString(cudaError));
+    fflush(stdout);
+    abort();
+  }
+}
+
+class fft_plan {
+private:
+  cufftHandle plan_;
+  int size_[3] = {0, 0, 0};
+  int dim_{3};
+  int batch_size_ = -1;
+  // forward or backward / reverse
+  enum fft_direction direction_ { FFT_UNKOWN };
+  cudaStream_t stream_;
+  bool is_initialized_{false};
+
+public:
+  fft_plan() {}
+
+  fft_plan(const std::vector<int> &fft_size__, const int dim__,
+           const int batch_size__, const enum fft_direction direction__)
+      : dim_(dim__), batch_size_(batch_size__) {
+    int inembed[2] = {0, 0};
+    int onembed[2] = {0, 0};
+    int istride, idist, ostride, odist;
+
+    direction_ = direction__;
+
+    switch (dim__) {
+    case 3: {
+      size_[0] = fft_size__[0];
+      size_[1] = fft_size__[1];
+      size_[2] = fft_size__[2];
+    } break;
+    case 2: {
+      size_[0] = fft_size__[1];
+      size_[1] = fft_size__[0];
+      size_[2] = batch_size__;
+      inembed[0] = fft_size__[1];
+      inembed[1] = fft_size__[0];
+      onembed[0] = fft_size__[1];
+      onembed[1] = fft_size__[0];
+      batch_size_ = batch_size__;
+
+      if (direction_ == CUFFT_FORWARD) {
+        istride = batch_size__;
+        idist = 1;
+        ostride = 1;
+        odist = size_[0] * size_[1];
+      } else {
+        istride = 1;
+        idist = size_[0] * size_[1];
+        ostride = batch_size__;
+        odist = 1;
+      }
+    } break;
+    case 1: {
+      size_[0] = fft_size__[0];
+      size_[1] = 1;
+      size_[2] = batch_size__;
+      batch_size_ = batch_size__;
+      if (direction_ == CUFFT_FORWARD) {
+        istride = batch_size__;
+        idist = 1;
+        ostride = 1;
+        odist = fft_size__[0];
+      } else {
+        istride = 1;
+        idist = fft_size__[0];
+        ostride = batch_size__;
+        odist = 1;
+      }
+      break;
+    }
+    default:
+      abort();
+      break;
+    }
+    if (dim_ == 3) {
+      fft_error_check(
+          cufftPlan3d(&plan_, size_[2], size_[1], size_[0], CUFFT_Z2Z),
+          __LINE__, __FILE__);
+    } else {
+      fft_error_check(cufftPlanMany(&plan_, dim_, &size_[0], inembed, istride,
+                                    idist, onembed, ostride, odist, CUFFT_Z2Z,
+                                    batch_size_),
+                      __LINE__, __FILE__);
+    }
+
+    is_initialized_ = true;
+  }
+
+  void set_stream(const cudaStream_t &cuda_stream) {
+    stream_ = cuda_stream;
+    fft_error_check(cufftSetStream(plan_, stream_), __LINE__, __FILE__);
+  }
+
+  ~fft_plan() {
+    if (is_initialized_)
+      destroy();
+  }
+
+  void destroy() {
+    error_check(cudaStreamSynchronize(stream_), __LINE__, __FILE__);
+    fft_error_check(cufftDestroy(plan_), __LINE__, __FILE__);
+    is_initialized_ = false;
+  }
+
+  /// run the fft on the data inplace
+  void execute_fft(const enum fft_direction direction__,
+                   cufftDoubleComplex *data__) {
+    fft_error_check(cufftExecZ2Z(plan_, data__, data__, direction__), __LINE__,
+                    __FILE__);
+  }
+
+  /// run the fft on the data out of place
+  void execute_fft(const enum fft_direction direction__,
+                   cufftDoubleComplex *dataIn__,
+                   cufftDoubleComplex *dataOut__) {
+    cufftResult_t cErr;
+    // set the stream
+
+    fft_error_check(cufftExecZ2Z(plan_, dataIn__, dataOut__, direction__),
+                    __LINE__, __FILE__);
+  }
+
+  /// check if this plane can be used to execute the fft
+  bool is_it_valid(std::vector<int> size__, const int dim__, const int batch__,
+                   const enum fft_direction direction__) const {
+    if (dim_ != dim__)
+      return false;
+    if (batch__ != batch_size_)
+      return false;
+    switch (dim__) {
+    case 3:
+      return ((size_[0] != size__[0]) || (size_[1] != size__[1]) ||
+              (size_[2] == size__[2]));
+      break;
+    case 2:
+      return ((size_[0] != size__[0]) || (size_[1] != size__[1]));
+      break;
+    case 1:
+      return (size_[0] != size__[0]);
+      break;
+    default:
+      return false;
+      break;
+    }
+
+    // check for the direction
+
+    if ((direction_ != direction__) && (dim_ != 3)) {
+      return false;
+    } else {
+      return true;
+    }
+  }
+
+  bool is_initialized() const { return is_initialized_; }
+};
+
+static void blasCreate(blasHandle_t *handle__) {
+  blas_error_check(cublasCreate(handle__), __LINE__, __FILE__);
+}
+
+static void blasDestroy(blasHandle_t &handle__) {
+  blas_error_check(cublasDestroy(handle__), __LINE__, __FILE__);
+}
+
+static void blasSetStream(blasHandle_t &handle__, offloadStream_t &stream__) {
+  blas_error_check(cublasSetStream(handle__, stream__), __LINE__, __FILE__);
+}
+
+inline void gpuDcopy(blasHandle_t &handle__, int n__, const double *x__,
+                     int incx__, double *y__, int incy__) {
+  blas_error_check(cublasDcopy(handle__, n__, x__, incx__, y__, incy__),
+                   __LINE__, __FILE__);
+}
+
+void gpu_gather(offloadStream_t &stream__, const double scale__,
+                int num_points__, const int *map_index__,
+                const double *dataIn__, double *dataOut__);
+
+void gpu_scatter(offloadStream_t &stream__, const double scale__,
+                 int num_points__, const int nmaps__, const int *map_index__,
+                 const double *dataIn__, double *dataOut__);
+#endif

src/pw/gpu/cuda/pw_cuda_z.cu

@@ -0,0 +1,40 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+#include "../../../offload/offload_library.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+#define __PW_CUDA_HIP_KERNELS
+
+#include "../kernels/pw_kernels.hpp"
+
+void gpu_scatter(cudaStream_t &stream__, const double scale__, int num_points__,
+                 const int nmaps__, const int *map_index__,
+                 const double *dataIn__, double *dataOut__) {
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  blocksPerGrid.x = num_points__ / 512 + ((num_points__ % 512) != 0);
+  threadsPerBlock.x = 512;
+
+  pw_scatter_z<double><<<blocksPerGrid, threadsPerBlock, 0, stream__>>>(
+      scale__, nmaps__, num_points__, map_index__, dataIn__, dataOut__);
+}
+
+void gpu_gather(cudaStream_t &stream__, const double scale__, int num_points__,
+                const int *map_index__, const double *dataIn__,
+                double *dataOut__) {
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  blocksPerGrid.x = num_points__ / 512 + ((num_points__ % 512) != 0);
+  threadsPerBlock.x = 512;
+
+  pw_gather_z<double><<<blocksPerGrid, threadsPerBlock, 0, stream__>>>(
+      scale__, num_points__, map_index__, dataIn__, dataOut__);
+}

src/pw/gpu/hip/PACKAGE

@@ -0,0 +1,5 @@
+{
+    "description": "hip FFT acceleration",
+    "requires": ["../../../offload", "../kernels"],
+    "archive": "libcp2kpwhip",
+}

src/pw/gpu/hip/hip_fft_private_header.hpp

@@ -0,0 +1,274 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+#ifndef HIP_FFT_PRIVATE_HEADER_HPP
+#define HIP_FFT_PRIVATE_HEADER_HPP
+#include <array>
+#include <cstdio>
+#include <hipblas.h>
+#include <hipfft.h>
+#include <map>
+#include <string>
+#include <unistd.h>
+#include <vector>
+
+#include "../../../offload/offload_library.h"
+#include "../../../offload/offload_operations.h"
+
+typedef hipfftDoubleComplex pw_complex_type;
+typedef hipblasHandle_t blasHandle_t;
+
+class fft_plan;
+
+enum fft_direction {
+  FFT_FORWARD = HIPFFT_FORWARD,
+  FFT_BACKWARD = HIPFFT_BACKWARD,
+  FFT_UNKOWN
+};
+
+static void fft_error_check(hipfftResult hipfftError, int line,
+                            const std::string filename) {
+  int pid;
+
+  std::map<hipfftResult, std::string> error_messages_ = {
+      {HIPFFT_SUCCESS, "HIPFFT_SUCCESS"},
+      {HIPFFT_INVALID_PLAN, "HIPFFT_INVALID_PLAN"},
+      {HIPFFT_ALLOC_FAILED, "HIPFFT_ALLOC_FAILED"},
+      {HIPFFT_INVALID_TYPE, "HIPFFT_INVALID_TYPE"},
+      {HIPFFT_INVALID_VALUE, "HIPFFT_INVALID_VALUE"},
+      {HIPFFT_INTERNAL_ERROR, "HIPFFT_INTERNAL_ERROR"},
+      {HIPFFT_EXEC_FAILED, "HIPFFT_EXEC_FAILED"},
+      {HIPFFT_SETUP_FAILED, "HIPFFT_SETUP_FAILED"},
+      {HIPFFT_INVALID_SIZE, "HIPFFT_INVALID_SIZE"},
+      {HIPFFT_INCOMPLETE_PARAMETER_LIST, "HIPFFT_INCOMPLETE_PARAMETER_LIST"},
+      {HIPFFT_INVALID_DEVICE, "HIPFFT_INVALID_DEVICE"},
+      {HIPFFT_PARSE_ERROR, "HIPFFT_PARSE_ERROR"},
+      {HIPFFT_NO_WORKSPACE, "HIPFFT_NO_WORKSPACE"},
+      {HIPFFT_NOT_IMPLEMENTED, "HIPFFT_NOT_IMPLEMENTED"},
+      {HIPFFT_NOT_SUPPORTED, "HIPFFT_NOT_SUPPORTED"}};
+
+  if (hipfftError != HIPFFT_SUCCESS) {
+    pid = getpid();
+    printf("%d HIP FFT Error line: %s %d \n", pid, filename.c_str(), line);
+    printf("%d HIP FFT Error (%s)\n", pid,
+           error_messages_.at(hipfftError).c_str());
+    fflush(stdout);
+    exit(-1);
+  }
+}
+
+static void blas_error_check(hipblasStatus_t hipblasError, int line,
+                             const std::string filename) {
+  int pid;
+
+  std::map<hipblasStatus_t, std::string> error_messages_ = {
+      {HIPBLAS_STATUS_SUCCESS, "HIPBLAS_STATUS_SUCCESS"},
+      {HIPBLAS_STATUS_NOT_INITIALIZED, "HIPBLAS_STATUS_NOT_INITIALIZED"},
+      {HIPBLAS_STATUS_ALLOC_FAILED, "HIPBLAS_STATUS_ALLOC_FAILED"},
+      {HIPBLAS_STATUS_INVALID_VALUE, "HIPBLAS_STATUS_INVALID_VALUE"},
+      {HIPBLAS_STATUS_MAPPING_ERROR, "HIPBLAS_STATUS_MAPPING_ERROR"},
+      {HIPBLAS_STATUS_EXECUTION_FAILED, "HIPBLAS_STATUS_EXECUTION_FAILED"},
+      {HIPBLAS_STATUS_INTERNAL_ERROR, "HIPBLAS_STATUS_INTERNAL_ERROR"},
+      {HIPBLAS_STATUS_NOT_SUPPORTED, "HIPBLAS_STATUS_NOT_SUPPORTED"},
+      {HIPBLAS_STATUS_ARCH_MISMATCH, "HIPBLAS_STATUS_ARCH_MISMATCH"},
+      {HIPBLAS_STATUS_HANDLE_IS_NULLPTR, "HIPBLAS_STATUS_HANDLE_IS_NULLPTR"},
+      {HIPBLAS_STATUS_INVALID_ENUM, "HIPBLAS_STATUS_INVALID_ENUM"},
+      {HIPBLAS_STATUS_UNKNOWN, "HIPBLAS_STATUS_UNKNOWN"}};
+
+  if (hipblasError != HIPBLAS_STATUS_SUCCESS) {
+    pid = getpid();
+    printf("%d HIP BLAS Error file: %s line: %d \n", pid, filename.c_str(),
+           line);
+    printf("%d HIP BLAS Error (%s)\n", pid,
+           error_messages_.at(hipblasError).c_str());
+    fflush(stdout);
+    exit(-1);
+  }
+}
+
+class fft_plan {
+private:
+  hipfftHandle plan_;
+  int size_[3] = {0, 0, 0};
+  int sign_{0};
+  int dim_{3};
+  enum fft_direction direction_ { FFT_UNKOWN };
+  int batch_size_ = -1;
+  offloadStream_t stream_;
+  bool is_initialized_{false};
+
+public:
+  fft_plan() {}
+
+  fft_plan(const std::vector<int> &fft_size__, const int dim__,
+           const int batch_size__, const enum fft_direction direction__)
+      : dim_(dim__), batch_size_(batch_size__) {
+    int inembed[2] = {-1, -1};
+    int onembed[2] = {-1, -1};
+    int istride, idist, ostride, odist;
+
+    direction_ = direction__;
+
+    switch (dim__) {
+    case 3: {
+      size_[0] = fft_size__[0];
+      size_[1] = fft_size__[1];
+      size_[2] = fft_size__[2];
+    } break;
+    case 2: {
+      size_[0] = fft_size__[1];
+      size_[1] = fft_size__[0];
+      size_[2] = batch_size__;
+      inembed[0] = fft_size__[1];
+      inembed[1] = fft_size__[0];
+      onembed[0] = fft_size__[1];
+      onembed[1] = fft_size__[0];
+      batch_size_ = batch_size__;
+
+      if (direction_ == FFT_FORWARD) {
+        istride = batch_size__;
+        idist = 1;
+        ostride = 1;
+        odist = size_[0] * size_[1];
+      } else {
+        istride = 1;
+        idist = size_[0] * size_[1];
+        ostride = batch_size__;
+        odist = 1;
+      }
+    } break;
+    case 1: {
+      size_[0] = fft_size__[0];
+      size_[1] = 1;
+      size_[2] = batch_size__;
+      batch_size_ = batch_size__;
+      if (direction_ == FFT_FORWARD) {
+        istride = batch_size__;
+        idist = 1;
+        ostride = 1;
+        odist = fft_size__[0];
+      } else {
+        istride = 1;
+        idist = fft_size__[0];
+        ostride = batch_size__;
+        odist = 1;
+      }
+      break;
+    }
+    default:
+      abort();
+      break;
+    }
+
+    if (dim_ == 3) {
+      fft_error_check(
+          hipfftPlan3d(&plan_, size_[2], size_[1], size_[0], HIPFFT_Z2Z),
+          __LINE__, __FILE__);
+    } else {
+      fft_error_check(hipfftPlanMany(&plan_, dim_, &size_[0], inembed, istride,
+                                     idist, onembed, ostride, odist, HIPFFT_Z2Z,
+                                     batch_size_),
+                      __LINE__, __FILE__);
+    }
+
+    fft_error_check(hipfftSetAutoAllocation(plan_, 1), __LINE__, __FILE__);
+    is_initialized_ = true;
+  }
+
+  void set_stream(const offloadStream_t &hip_stream) {
+    stream_ = hip_stream;
+    fft_error_check(hipfftSetStream(plan_, stream_), __LINE__, __FILE__);
+  }
+
+  ~fft_plan() {
+    if (is_initialized_)
+      destroy();
+  }
+
+  void destroy() {
+    offloadStreamSynchronize(stream_);
+    fft_error_check(hipfftDestroy(plan_), __LINE__, __FILE__);
+    is_initialized_ = false;
+  }
+
+  /// run the fft on the data inplace
+  void execute_fft(const enum fft_direction direction__,
+                   hipfftDoubleComplex *data__) {
+    fft_error_check(hipfftExecZ2Z(plan_, data__, data__, direction__), __LINE__,
+                    __FILE__);
+  }
+
+  /// run the fft on the data out of place
+  void execute_fft(const enum fft_direction direction__,
+                   hipfftDoubleComplex *dataIn__,
+                   hipfftDoubleComplex *dataOut__) {
+    fft_error_check(hipfftExecZ2Z(plan_, dataIn__, dataOut__, direction__),
+                    __LINE__, __FILE__);
+  }
+
+  /// check if this plane can be used to execute the fft
+  bool is_it_valid(std::vector<int> size__, const int dim__, const int batch__,
+                   const enum fft_direction direction__) const {
+    if (dim_ != dim__)
+      return false;
+    if (batch__ != batch_size_)
+      return false;
+    switch (dim__) {
+    case 3:
+      return ((size_[0] != size__[0]) || (size_[1] != size__[1]) ||
+              (size_[2] == size__[2]));
+      break;
+    case 2:
+      return ((size_[0] != size__[0]) || (size_[1] != size__[1]));
+      break;
+    case 1:
+      return (size_[0] != size__[0]);
+      break;
+    default:
+      return false;
+      break;
+    }
+
+    // check for the direction
+
+    if ((direction__ != direction_) && (dim_ != 3)) {
+      return false;
+    } else {
+      return true;
+    }
+  }
+
+  bool is_initialized() const { return is_initialized_; }
+};
+
+static void blasCreate(blasHandle_t *handle__) {
+  blas_error_check(hipblasCreate(handle__), __LINE__, __FILE__);
+}
+
+static void blasDestroy(blasHandle_t &handle__) {
+  blas_error_check(hipblasDestroy(handle__), __LINE__, __FILE__);
+}
+
+static void blasSetStream(blasHandle_t &handle__, offloadStream_t &stream__) {
+  blas_error_check(hipblasSetStream(handle__, stream__), __LINE__, __FILE__);
+}
+
+static void gpuDcopy(blasHandle_t &handle__, int n__, const double *x__,
+                     int incx__, double *y__, int incy__) {
+  blas_error_check(hipblasDcopy(handle__, n__, x__, incx__, y__, incy__),
+                   __LINE__, __FILE__);
+}
+
+void gpu_gather(hipStream_t &stream__, const double scale__, int num_points__,
+                const int *map_index__, const double *dataIn__,
+                double *dataOut__);
+
+void gpu_scatter(hipStream_t &stream__, const double scale__, int num_points__,
+                 const int nmaps__, const int *map_index__,
+                 const double *dataIn__, double *dataOut__);
+
+#endif

src/pw/gpu/hip/pw_hip_z.cc

@@ -0,0 +1,40 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+// global dependencies
+#include <cstdio>
+#include <hip/hip_runtime.h>
+#include <hip/hip_runtime_api.h>
+#include <hipblas.h>
+#include <hipfft.h>
+
+#define __PW_CUDA_HIP_KERNELS
+#include "../kernels/pw_kernels.hpp"
+
+void gpu_scatter(hipStream_t &stream__, const double scale__, int num_points__,
+                 const int nmaps__, const int *map_index__,
+                 const double *dataIn__, double *dataOut__) {
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  blocksPerGrid.x = num_points__ / 512 + ((num_points__ % 512) != 0);
+  threadsPerBlock.x = 512;
+
+  pw_scatter_z<double><<<blocksPerGrid, threadsPerBlock, 0, stream__>>>(
+      scale__, nmaps__, num_points__, map_index__, dataIn__, dataOut__);
+}
+
+void gpu_gather(hipStream_t &stream__, const double scale__, int num_points__,
+                const int *map_index__, const double *dataIn__,
+                double *dataOut__) {
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  blocksPerGrid.x = num_points__ / 512 + ((num_points__ % 512) != 0);
+  threadsPerBlock.x = 512;
+
+  pw_gather_z<double><<<blocksPerGrid, threadsPerBlock, 0, stream__>>>(
+      scale__, num_points__, map_index__, dataIn__, dataOut__);
+}

src/pw/gpu/kernels/PACKAGE

@@ -0,0 +1,4 @@
+{
+    "description": "kernels for fft on gpu",
+    "requires": [],
+}

src/pw/gpu/kernels/pw_kernels.hpp

@@ -0,0 +1,53 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+#ifndef PW_KERNELS_HPP
+#define PW_KERNELS_HPP
+
+#ifndef __PW_CUDA_HIP_KERNELS
+#error "this header file can only be included in pw_hip_z.cc or pw_cuda_z.cu"
+#endif
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) gather and scale on the GPU.
+ * \version 0.01
+ ******************************************************************************/
+template <typename T>
+__global__ void pw_gather_z(const T scale, const int ngpts,
+                            const int *__restrict__ const ghatmap,
+                            const T *__restrict__ const c,
+                            T *__restrict__ pwcc) {
+
+  const int igpt = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if (igpt >= ngpts)
+    return;
+
+  pwcc[2 * igpt] = scale * c[2 * ghatmap[igpt]];
+  pwcc[2 * igpt + 1] = scale * c[2 * ghatmap[igpt] + 1];
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) scatter and scale on the GPU.
+ * \version 0.01
+ ******************************************************************************/
+template <typename T>
+__global__ void pw_scatter_z(const T scale, const int nmaps, const int ngpts,
+                             const int *ghatmap, const T *pwcc, T *c) {
+  const int igpt = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if (igpt >= ngpts)
+    return;
+
+  c[2 * ghatmap[igpt]] = scale * pwcc[2 * igpt];
+  c[2 * ghatmap[igpt] + 1] = scale * pwcc[2 * igpt + 1];
+  if (nmaps == 2) {
+    c[2 * ghatmap[igpt + ngpts]] = scale * pwcc[2 * igpt];
+    c[2 * ghatmap[igpt + ngpts] + 1] = -scale * pwcc[2 * igpt + 1];
+  }
+}
+#endif

src/pw/gpu/pw_gpu_internal.cc

@@ -0,0 +1,582 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2022 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: GPL-2.0-or-later                                 */
+/*----------------------------------------------------------------------------*/
+
+#include <array>
+#include <cstdio>
+#include <map>
+#include <string>
+#include <unistd.h>
+#include <vector>
+
+#include "../../offload/offload_library.h"
+#include "../../offload/offload_operations.h"
+
+#ifdef __PW_CUDA
+#include "cuda/cuda_fft_private_header.hpp"
+#endif
+
+#ifdef __PW_HIP
+#include "hip/hip_fft_private_header.hpp"
+#endif
+
+static bool is_configured{false};
+
+static blasHandle_t handle_;
+static offloadStream_t streams_;
+
+// INIT/RELEASE
+extern "C" int pw_gpu_init() {
+  if (!is_configured) {
+    offload_set_device();
+    /* stream allocation  */
+    offloadStreamCreate(&streams_);
+    blasCreate(&handle_);
+    blasSetStream(handle_, streams_);
+    is_configured = true;
+  }
+  return 0;
+}
+
+extern "C" void pw_gpu_finalize() {
+  if (is_configured) {
+    offload_set_device();
+
+    offloadStreamDestroy(streams_);
+    blasDestroy(handle_);
+    is_configured = false;
+  }
+}
+
+template <typename T> void offloadMalloc(T **ptr__, size_t size__) {
+  offloadMalloc((void **)ptr__, size__ * sizeof(T));
+}
+
+// Double precision complex procedures
+inline void fft_gpu_run_3d_z_(const offloadStream_t &stream,
+                              const enum fft_direction fsign, const int *n,
+                              pw_complex_type *data_in__,
+                              pw_complex_type *data_out__) {
+  std::vector<int> size(3);
+  size[0] = n[0];
+  size[1] = n[1];
+  size[2] = n[2];
+  fft_plan plan = fft_plan(size, 3, 0, fsign);
+
+  plan.set_stream(stream);
+
+  plan.execute_fft(fsign, data_in__, data_out__);
+
+  offloadStreamSynchronize(stream);
+}
+
+inline void fft_gpu_run_2dm_z_(const offloadStream_t &stream,
+                               const enum fft_direction fsign, const int *n,
+                               pw_complex_type *data_in,
+                               pw_complex_type *data_out) {
+  std::vector<int> size(2);
+  size[0] = n[1];
+  size[1] = n[2];
+  const int batch = n[0];
+  fft_plan plan = fft_plan(size, 2, batch, fsign);
+
+  plan.set_stream(stream);
+
+  plan.execute_fft(fsign, data_in, data_out);
+
+  offloadStreamSynchronize(stream);
+}
+
+inline void fft_gpu_run_1dm_z_(const offloadStream_t &stream,
+                               const enum fft_direction fsign, const int n,
+                               const int m, pw_complex_type *data_in,
+                               pw_complex_type *data_out) {
+  fft_plan plan = fft_plan({n}, 1, m, fsign);
+  plan.set_stream(stream);
+
+  plan.execute_fft(fsign, data_in, data_out);
+
+  offloadStreamSynchronize(stream);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) FFT, followed by a (double
+ *          precision complex) gather, on the GPU.
+ ******************************************************************************/
+
+extern "C" void pw_gpu_cfffg_z_(const double *din, double *zout,
+                                const int *ghatmap, const int *npts,
+                                const int ngpts, const double scale) {
+
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  int *ghatmap_dev;
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0 || ngpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts); //? (ngpts)
+  offloadMemsetAsync(ptr_2, 0, nrpts * sizeof(pw_complex_type), streams_);
+
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  offloadMalloc<int>(&ghatmap_dev, ngpts);
+  // copy gather map array from host to the device
+  offloadMemcpyAsyncHtoD(ghatmap_dev, ghatmap, sizeof(int) * ngpts, streams_);
+
+  // convert the real (host) pointer 'din' into a complex (device)
+  // pointer 'ptr_1'
+  // copy to device (NOTE: only first half of ptr_1 is written!)
+  offloadMemcpyAsyncHtoD(ptr_1, din, sizeof(double) * nrpts, streams_);
+
+  // real to complex blow-up
+  gpuDcopy(handle_, nrpts, ptr_1, 1, ptr_2, 2);
+
+  // fft on the GPU
+  fft_gpu_run_3d_z_(streams_, fft_direction::FFT_FORWARD, npts,
+                    (pw_complex_type *)ptr_2, (pw_complex_type *)ptr_1);
+
+  // gather on the GPU
+  gpu_gather(streams_, scale, ngpts, ghatmap_dev, ptr_1, ptr_2);
+  offloadMemcpyAsyncDtoH(zout, ptr_2, sizeof(pw_complex_type) * ngpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+  offloadFree(ghatmap_dev);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) scatter, followed by a inverse
+ *          (double precision complex) FFT, on the GPU.
+ ******************************************************************************/
+
+extern "C" void pw_gpu_sfffc_z_(const double *zin, double *dout,
+                                const int *ghatmap, const int *npts,
+                                const int ngpts, const int nmaps,
+                                const double scale) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  int *ghatmap_dev;
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0 || ngpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+  offloadMemsetAsync(
+      ptr_2, 0, sizeof(pw_complex_type) * nrpts,
+      streams_); // we need to do this only if spherical cut-off is used!
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  offloadMalloc<int>(&ghatmap_dev, nmaps * ngpts);
+
+  // copy all arrays from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * ngpts, streams_);
+  offloadMemcpyAsyncHtoD(ghatmap_dev, ghatmap, sizeof(int) * nmaps * ngpts,
+                         streams_);
+  gpu_scatter(streams_, scale, ngpts, nmaps, ghatmap_dev, ptr_1, ptr_2);
+
+  // fft on the GPU (streams_[1])
+  fft_gpu_run_3d_z_(streams_, fft_direction::FFT_BACKWARD, npts,
+                    (pw_complex_type *)ptr_2, (pw_complex_type *)ptr_1);
+
+  // convert the complex (device) pointer 'ptr_2' into a real (host)
+  // pointer 'dout' (NOTE: Only first half of ptr_1 is written!)
+  // pw_copy_cr_cu_z<<<blocksPerGrid, threadsPerBlock, 0, streams_>>>(ptr_2,
+  // ptr_1, nrpts);
+  gpuDcopy(handle_, nrpts, ptr_1, 2, ptr_2, 1);
+
+  offloadMemcpyAsyncDtoH(dout, ptr_2, sizeof(double) * nrpts, streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+  offloadFree(ghatmap_dev);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double to complex double) blow-up and a (double
+ *          precision complex) 2D-FFT on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_cff_z_(const double *din, double *zout,
+                              const int *npts) {
+  double *ptr_1, *ptr_2;
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers for:
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+  offloadMemsetAsync(ptr_2, 0, sizeof(pw_complex_type) * nrpts, streams_);
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+
+  // convert the real (host) pointer 'din' into a complex (device)
+  // pointer 'ptr_in'
+  // copy all arrays from host to the device (NOTE: Only first half of ptr_1 is
+  // written!)
+  offloadMemcpyAsyncHtoD(ptr_1, din, sizeof(double) * nrpts, streams_);
+  gpuDcopy(handle_, nrpts, ptr_1, 1, ptr_2, 2);
+
+  // pw_copy_rc_cu_z<<<blocksPerGrid, threadsPerBlock, 0, streams_>>>(
+  //     ptr_1, ptr_2, nrpts);
+
+  // fft on the GPU (streams_[1])
+  // NOTE: the following works, but CUDA does 2D-FFT in C-shaped (not optimal)
+  // order fftcu_run_2dm_z_(1, npts, 1.0e0, (cufftDoubleComplex *) ptr_2,
+  // (cufftDoubleComplex *) ptr_1, streams_[1]);
+
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_FORWARD, npts[2],
+                     npts[0] * npts[1], (pw_complex_type *)ptr_2,
+                     (pw_complex_type *)ptr_1);
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_FORWARD, npts[1],
+                     npts[0] * npts[2], (pw_complex_type *)ptr_1,
+                     (pw_complex_type *)ptr_2);
+
+  offloadMemcpyAsyncDtoH(zout, ptr_2, sizeof(pw_complex_type) * nrpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) 2D-FFT and a (double complex
+ *          to double) shrink-down on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_ffc_z_(const double *zin, double *dout,
+                              const int *npts) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers for:
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // copy input data from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * nrpts, streams_);
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+
+  // fftcu_run_2dm_z_(-1, npts, 1.0e0, (cufftDoubleComplex *) ptr_1,
+  // (cufftDoubleComplex *) ptr_2, streams_[1]);
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_BACKWARD, npts[1],
+                     npts[0] * npts[2], (pw_complex_type *)ptr_1,
+                     (pw_complex_type *)ptr_2);
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_BACKWARD, npts[2],
+                     npts[0] * npts[1], (pw_complex_type *)ptr_2,
+                     (pw_complex_type *)ptr_1);
+
+  // CUDA blocking for pw_copy_cr (currently only 2-D grid)
+  // get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
+
+  // take the real part of the FFT and store the result in dout
+  gpuDcopy(handle_, nrpts, ptr_1, 2, ptr_2, 1);
+
+  offloadMemcpyAsyncDtoH(dout, ptr_2, sizeof(double) * nrpts, streams_);
+
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double to complex double) blow-up and a (double
+ *          precision complex) 1D-FFT on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_cf_z_(const double *din, double *zout, const int *npts) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // convert the real (host) pointer 'din' into a complex (device)
+  // pointer 'ptr_2' (NOTE: Only first half of ptr_1 is written!)
+  // copy all arrays from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, din, sizeof(double) * nrpts, streams_);
+
+  // get device memory pointers for:
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+  offloadMemsetAsync(ptr_2, 0, sizeof(pw_complex_type) * nrpts, streams_);
+
+  gpuDcopy(handle_, nrpts, ptr_1, 1, ptr_2, 2);
+
+  // fft on the GPU (streams_[1])
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_FORWARD, npts[2],
+                     npts[0] * npts[1], (pw_complex_type *)ptr_2,
+                     (pw_complex_type *)ptr_1);
+
+  offloadMemcpyAsyncDtoH(zout, ptr_1, sizeof(pw_complex_type) * nrpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) 1D-FFT and a (double complex
+ *          to double) shrink-down on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_fc_z_(const double *zin, double *dout, const int *npts) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * npts[1] * npts[2];
+  if (nrpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers for:
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // copy input data from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * nrpts, streams_);
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_BACKWARD, npts[2],
+                     npts[0] * npts[1], (pw_complex_type *)ptr_1,
+                     (pw_complex_type *)ptr_2);
+
+  // CUDA blocking for pw_copy_cr (currently only 2-D grid)
+  // get_grid_params(nrpts, MAXTHREADS, threadsPerBlock, blocksPerGrid);
+
+  // convert the complex (device) pointer 'ptr_2' into a real (host)
+  // pointer 'dout' (NOTE: Only first half of ptr_1 is written!)
+  gpuDcopy(handle_, nrpts, ptr_2, 2, ptr_1, 1);
+
+  offloadMemcpyAsyncDtoH(dout, ptr_1, sizeof(double) * nrpts, streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) 1D-FFT on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_f_z_(const double *zin, double *zout, const int dir,
+                            const int n, const int m) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  const enum fft_direction direction =
+      (dir > 0) ? (fft_direction::FFT_FORWARD) : (fft_direction::FFT_BACKWARD);
+  // dimensions of complex arrays
+  const int nrpts = n * m;
+  if (nrpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers for:
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // copy input data from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * nrpts, streams_);
+
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+
+  fft_gpu_run_1dm_z_(streams_, direction, n, m, (pw_complex_type *)ptr_1,
+                     (pw_complex_type *)ptr_2);
+
+  offloadMemcpyAsyncDtoH(zout, ptr_2, sizeof(pw_complex_type) * nrpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) 1D-FFT, followed by a (double
+ *          precision complex) gather, on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_fg_z_(const double *zin, double *zout,
+                             const int *ghatmap, const int *npts,
+                             const int mmax, const int ngpts,
+                             const double scale) {
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  int *ghatmap_dev;
+  //  int nrpts;
+  dim3 blocksPerGrid, threadsPerBlock;
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * mmax;
+  if (nrpts == 0 || ngpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // transfer input data from host to device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * nrpts, streams_);
+  offloadMalloc<int>(&ghatmap_dev, ngpts);
+  // transfer gather data from host to device
+  offloadMemcpyAsyncHtoD(ghatmap_dev, ghatmap, sizeof(int) * ngpts, streams_);
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_FORWARD, npts[0], mmax,
+                     (pw_complex_type *)ptr_1, (pw_complex_type *)ptr_2);
+
+  // gather on the GPU
+  gpu_gather(streams_, scale, ngpts, ghatmap_dev, ptr_2, ptr_1);
+  offloadMemcpyAsyncDtoH(zout, ptr_1, sizeof(pw_complex_type) * ngpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+  offloadFree(ghatmap_dev);
+}
+
+/*******************************************************************************
+ * \brief   Performs a (double precision complex) scatter, followed by a
+ *          (double precision complex) 1D-FFT, on the GPU.
+ ******************************************************************************/
+extern "C" void pw_gpu_sf_z_(const double *zin, double *zout,
+                             const int *ghatmap, const int *npts,
+                             const int mmax, const int ngpts, const int nmaps,
+                             const double scale) {
+
+  if (!is_configured) {
+    fprintf(stderr, "call to pw_gpu_init is missing\n");
+    abort();
+  }
+
+  double *ptr_1, *ptr_2;
+  int *ghatmap_dev;
+
+  // dimensions of double and complex arrays
+  const int nrpts = npts[0] * mmax;
+  if (nrpts == 0 || ngpts == 0)
+    return;
+
+  // get streams and events
+  offload_set_device();
+
+  // get device memory pointers
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_1, nrpts);
+  // transfer input data from host to the device
+  offloadMemcpyAsyncHtoD(ptr_1, zin, sizeof(pw_complex_type) * ngpts, streams_);
+
+  offloadMalloc<int>(&ghatmap_dev, nmaps * ngpts);
+  // transfer scatter data from host to device
+  offloadMemcpyAsyncHtoD(ghatmap_dev, ghatmap, sizeof(int) * nmaps * ngpts,
+                         streams_);
+
+  offloadMalloc<pw_complex_type>((pw_complex_type **)&ptr_2, nrpts);
+  offloadMemsetAsync(
+      ptr_2, 0, sizeof(pw_complex_type) * nrpts,
+      streams_); // we need to do this only if spherical cut-off is used!
+
+  // scatter on the GPU
+  gpu_scatter(streams_, scale, ngpts, nmaps, ghatmap_dev, ptr_1, ptr_2);
+
+  // fft on the GPU
+  fft_gpu_run_1dm_z_(streams_, fft_direction::FFT_BACKWARD, npts[0], mmax,
+                     (pw_complex_type *)ptr_2, (pw_complex_type *)ptr_1);
+
+  offloadMemcpyAsyncDtoH(zout, ptr_1, sizeof(pw_complex_type) * nrpts,
+                         streams_);
+
+  // synchronize with respect to host
+  offloadStreamSynchronize(streams_);
+
+  // release memory stack
+  offloadFree(ptr_1);
+  offloadFree(ptr_2);
+  offloadFree(ghatmap_dev);
+}

src/pw/pw_gpu.F

@@ -7,8 +7,8 @@
 
 ! **************************************************************************************************
 !> \note
-!> This module contains routines necessary to operate on plane waves on NVIDIA
-!> GPUs using CUDA. It depends at execution time on NVIDIA's CUFFT library.
+!> This module contains routines necessary to operate on plane waves on GPUs
+!     >  independently of the GPU platform.
 !> \par History
 !>      BGL (06-Mar-2008)  : Created
 !>      AG  (18-May-2012)  : Refacturing:
@@ -17,9 +17,14 @@
 !>      AG  (11-Sept-2012) : Modifications:
 !>                          - use pointers if precision mapping is not required
 !>                          - use OMP for mapping
+!>      MT  (Jan 2022)     : Modifications
+!>                          - use a generic interface for fft calls to GPUs
+!>                          - Support both Nvidia and AMD GPUs. Other GPUs manufacturers
+!>                            can be added easily.
 !> \author Benjamin G. Levine
 ! **************************************************************************************************
-MODULE pw_cuda
+#if defined (__PW_GPU)
+MODULE pw_gpu
    USE ISO_C_BINDING,                   ONLY: C_DOUBLE,&
                                               C_INT,&
                                               C_LOC,&
@@ -42,15 +47,15 @@ MODULE pw_cuda
 
    PRIVATE
 
-   PUBLIC :: pw_cuda_r3dc1d_3d
-   PUBLIC :: pw_cuda_c1dr3d_3d
-   PUBLIC :: pw_cuda_r3dc1d_3d_ps
-   PUBLIC :: pw_cuda_c1dr3d_3d_ps
-   PUBLIC :: pw_cuda_init, pw_cuda_finalize
+   PUBLIC :: pw_gpu_r3dc1d_3d
+   PUBLIC :: pw_gpu_c1dr3d_3d
+   PUBLIC :: pw_gpu_r3dc1d_3d_ps
+   PUBLIC :: pw_gpu_c1dr3d_3d_ps
+   PUBLIC :: pw_gpu_init, pw_gpu_finalize
 
    ! Explicit interfaces to double precision complex transformation
-   ! routines (C+CUDA). For more details, see: ../cuda_tools!
-   INTERFACE pw_cuda_cfffg_cu
+   ! routines (C+GPU). For more details, see: gpu
+   INTERFACE pw_gpu_cfffg_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param din ...
@@ -60,8 +65,8 @@ MODULE pw_cuda
 !> \param ngpts ...
 !> \param scale ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_cfffg_z(din, zout, ghatmap, npts, ngpts, scale) &
-         BIND(C, name="pw_cuda_cfffg_z_")
+      SUBROUTINE pw_gpu_cfffg_z(din, zout, ghatmap, npts, ngpts, scale) &
+         BIND(C, name="pw_gpu_cfffg_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: din
          TYPE(C_PTR), VALUE                       :: zout
@@ -71,10 +76,10 @@ MODULE pw_cuda
          INTEGER(KIND=C_INT), INTENT(IN), VALUE   :: ngpts
          REAL(KIND=C_DOUBLE), INTENT(IN), VALUE   :: scale
 
-      END SUBROUTINE pw_cuda_cfffg_z
+      END SUBROUTINE pw_gpu_cfffg_z
    END INTERFACE
 
-   INTERFACE pw_cuda_sfffc_cu
+   INTERFACE pw_gpu_sfffc_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
@@ -85,8 +90,8 @@ MODULE pw_cuda
 !> \param nmaps ...
 !> \param scale ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_sfffc_z(zin, dout, ghatmap, npts, ngpts, nmaps, scale) &
-         BIND(C, name="pw_cuda_sfffc_z_")
+      SUBROUTINE pw_gpu_sfffc_z(zin, dout, ghatmap, npts, ngpts, nmaps, scale) &
+         BIND(C, name="pw_gpu_sfffc_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: dout
@@ -96,82 +101,82 @@ MODULE pw_cuda
          INTEGER(KIND=C_INT), INTENT(IN), VALUE   :: ngpts, nmaps
          REAL(KIND=C_DOUBLE), INTENT(IN), VALUE   :: scale
 
-      END SUBROUTINE pw_cuda_sfffc_z
+      END SUBROUTINE pw_gpu_sfffc_z
    END INTERFACE
 
-   INTERFACE pw_cuda_cff_cu
+   INTERFACE pw_gpu_cff_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param din ...
 !> \param zout ...
 !> \param npts ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_cff_z(din, zout, npts) &
-         BIND(C, name="pw_cuda_cff_z_")
+      SUBROUTINE pw_gpu_cff_z(din, zout, npts) &
+         BIND(C, name="pw_gpu_cff_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: din
          TYPE(C_PTR), VALUE                       :: zout
          INTEGER(KIND=C_INT), DIMENSION(*), &
             INTENT(IN)                             :: npts
 
-      END SUBROUTINE pw_cuda_cff_z
+      END SUBROUTINE pw_gpu_cff_z
    END INTERFACE
 
-   INTERFACE pw_cuda_ffc_cu
+   INTERFACE pw_gpu_ffc_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
 !> \param dout ...
 !> \param npts ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_ffc_z(zin, dout, npts) &
-         BIND(C, name="pw_cuda_ffc_z_")
+      SUBROUTINE pw_gpu_ffc_z(zin, dout, npts) &
+         BIND(C, name="pw_gpu_ffc_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: dout
          INTEGER(KIND=C_INT), DIMENSION(*), &
             INTENT(IN)                             :: npts
 
-      END SUBROUTINE pw_cuda_ffc_z
+      END SUBROUTINE pw_gpu_ffc_z
    END INTERFACE
 
-   INTERFACE pw_cuda_cf_cu
+   INTERFACE pw_gpu_cf_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param din ...
 !> \param zout ...
 !> \param npts ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_cf_z(din, zout, npts) &
-         BIND(C, name="pw_cuda_cf_z_")
+      SUBROUTINE pw_gpu_cf_z(din, zout, npts) &
+         BIND(C, name="pw_gpu_cf_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: din
          TYPE(C_PTR), VALUE                       :: zout
          INTEGER(KIND=C_INT), DIMENSION(*), &
             INTENT(IN)                             :: npts
 
-      END SUBROUTINE pw_cuda_cf_z
+      END SUBROUTINE pw_gpu_cf_z
    END INTERFACE
 
-   INTERFACE pw_cuda_fc_cu
+   INTERFACE pw_gpu_fc_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
 !> \param dout ...
 !> \param npts ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_fc_z(zin, dout, npts) &
-         BIND(C, name="pw_cuda_fc_z_")
+      SUBROUTINE pw_gpu_fc_z(zin, dout, npts) &
+         BIND(C, name="pw_gpu_fc_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: dout
          INTEGER(KIND=C_INT), DIMENSION(*), &
             INTENT(IN)                             :: npts
 
-      END SUBROUTINE pw_cuda_fc_z
+      END SUBROUTINE pw_gpu_fc_z
    END INTERFACE
 
-   INTERFACE pw_cuda_f_cu
+   INTERFACE pw_gpu_f_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
@@ -180,17 +185,17 @@ MODULE pw_cuda
 !> \param n ...
 !> \param m ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_f_z(zin, zout, dir, n, m) &
-         BIND(C, name="pw_cuda_f_z_")
+      SUBROUTINE pw_gpu_f_z(zin, zout, dir, n, m) &
+         BIND(C, name="pw_gpu_f_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: zout
          INTEGER(KIND=C_INT), INTENT(IN), VALUE   :: dir, n, m
 
-      END SUBROUTINE pw_cuda_f_z
+      END SUBROUTINE pw_gpu_f_z
    END INTERFACE
 
-   INTERFACE pw_cuda_fg_cu
+   INTERFACE pw_gpu_fg_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
@@ -201,8 +206,8 @@ MODULE pw_cuda
 !> \param ngpts ...
 !> \param scale ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_fg_z(zin, zout, ghatmap, npts, mmax, ngpts, scale) &
-         BIND(C, name="pw_cuda_fg_z_")
+      SUBROUTINE pw_gpu_fg_z(zin, zout, ghatmap, npts, mmax, ngpts, scale) &
+         BIND(C, name="pw_gpu_fg_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: zout
@@ -212,10 +217,10 @@ MODULE pw_cuda
          INTEGER(KIND=C_INT), INTENT(IN), VALUE   :: mmax, ngpts
          REAL(KIND=C_DOUBLE), INTENT(IN), VALUE   :: scale
 
-      END SUBROUTINE pw_cuda_fg_z
+      END SUBROUTINE pw_gpu_fg_z
    END INTERFACE
 
-   INTERFACE pw_cuda_sf_cu
+   INTERFACE pw_gpu_sf_cu
 ! **************************************************************************************************
 !> \brief ...
 !> \param zin ...
@@ -227,8 +232,8 @@ MODULE pw_cuda
 !> \param nmaps ...
 !> \param scale ...
 ! **************************************************************************************************
-      SUBROUTINE pw_cuda_sf_z(zin, zout, ghatmap, npts, mmax, ngpts, nmaps, scale) &
-         BIND(C, name="pw_cuda_sf_z_")
+      SUBROUTINE pw_gpu_sf_z(zin, zout, ghatmap, npts, mmax, ngpts, nmaps, scale) &
+         BIND(C, name="pw_gpu_sf_z_")
          IMPORT
          TYPE(C_PTR), INTENT(IN), VALUE           :: zin
          TYPE(C_PTR), VALUE                       :: zout
@@ -238,45 +243,42 @@ MODULE pw_cuda
          INTEGER(KIND=C_INT), INTENT(IN), VALUE   :: mmax, ngpts, nmaps
          REAL(KIND=C_DOUBLE), INTENT(IN), VALUE   :: scale
 
-      END SUBROUTINE pw_cuda_sf_z
+      END SUBROUTINE pw_gpu_sf_z
    END INTERFACE
 
    INTERFACE
-      FUNCTION pw_cuda_init_cu() RESULT(istat) BIND(C, name="pw_cuda_init")
+      FUNCTION pw_gpu_init_cu() RESULT(istat) BIND(C, name="pw_gpu_init")
          IMPORT
          INTEGER(KIND=C_INT)                    :: istat
-      END FUNCTION pw_cuda_init_cu
-      SUBROUTINE pw_cuda_finalize_cu() BIND(C, name="pw_cuda_finalize")
-      END SUBROUTINE pw_cuda_finalize_cu
+      END FUNCTION pw_gpu_init_cu
+      SUBROUTINE pw_gpu_finalize_cu() BIND(C, name="pw_gpu_finalize")
+      END SUBROUTINE pw_gpu_finalize_cu
    END INTERFACE
 
-   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_methods_cuda'
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_methods_gpu'
    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
 
 CONTAINS
 
 ! **************************************************************************************************
-!> \brief Allocates resources on the cuda device for cuda fft acceleration
+!> \brief Allocates resources on the gpu device for gpu fft acceleration
 !> \author Ole Schuett
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_init()
-#if defined (__PW_CUDA)
-      INTEGER :: istat
-      istat = pw_cuda_init_cu()
+   SUBROUTINE pw_gpu_init()
+      INTEGER                                            :: istat
+
+      istat = pw_gpu_init_cu()
       IF (istat /= 0) &
-         CPABORT("pw_cuda_init: failed")
-#endif
-   END SUBROUTINE pw_cuda_init
+         CPABORT("pw_gpu_init: failed")
+   END SUBROUTINE pw_gpu_init
 
 ! **************************************************************************************************
-!> \brief Releases resources on the cuda device for cuda fft acceleration
+!> \brief Releases resources on the gpu device for gpu fft acceleration
 !> \author Ole Schuett
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_finalize()
-#if defined (__PW_CUDA)
-      CALL pw_cuda_finalize_cu()
-#endif
-   END SUBROUTINE pw_cuda_finalize
+   SUBROUTINE pw_gpu_finalize()
+      CALL pw_gpu_finalize_cu()
+   END SUBROUTINE pw_gpu_finalize
 
 ! **************************************************************************************************
 !> \brief perform an fft followed by a gather on the gpu
@@ -285,25 +287,18 @@ CONTAINS
 !> \param scale ...
 !> \author Benjamin G Levine
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_r3dc1d_3d(pw1, pw2, scale)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
-      REAL(KIND=dp)                            :: scale
+   SUBROUTINE pw_gpu_r3dc1d_3d(pw1, pw2, scale)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
+      REAL(KIND=dp)                                      :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pw2)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_r3dc1d_3d'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_r3dc1d_3d'
 
-      INTEGER                                  :: handle, ngpts
-      INTEGER                                  :: l1, l2, l3
-      INTEGER, DIMENSION(:), POINTER           :: npts
-
-      REAL(KIND=dp), POINTER                   :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
-      INTEGER, POINTER                         :: ptr_ghatmap
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwout
+      INTEGER                                            :: handle, l1, l2, l3, ngpts
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      INTEGER, POINTER                                   :: ptr_ghatmap
+      REAL(KIND=dp), POINTER                             :: ptr_pwin
 
       CALL timeset(routineN, handle)
 
@@ -321,13 +316,12 @@ CONTAINS
       ptr_ghatmap => pw2%pw_grid%g_hatmap(1, 1)
 
       ! invoke the combined transformation
-      CALL pw_cuda_cfffg_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, ngpts, scale)
+      CALL pw_gpu_cfffg_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, ngpts, scale)
 
       pw2%in_space = RECIPROCALSPACE
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_r3dc1d_3d
+   END SUBROUTINE pw_gpu_r3dc1d_3d
 
 ! **************************************************************************************************
 !> \brief perform an scatter followed by a fft on the gpu
@@ -336,25 +330,18 @@ CONTAINS
 !> \param scale ...
 !> \author Benjamin G Levine
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_c1dr3d_3d(pw1, pw2, scale)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
-      REAL(KIND=dp)                            :: scale
+   SUBROUTINE pw_gpu_c1dr3d_3d(pw1, pw2, scale)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
+      REAL(KIND=dp)                                      :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pw2)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_c1dr3d_3d'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_c1dr3d_3d'
 
-      INTEGER                                  :: handle, ngpts, nmaps
-      INTEGER                                  :: l1, l2, l3
-      INTEGER, DIMENSION(:), POINTER           :: npts
-
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      REAL(KIND=dp), POINTER                   :: ptr_pwout
-      INTEGER, POINTER                         :: ptr_ghatmap
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin
+      INTEGER                                            :: handle, l1, l2, l3, ngpts, nmaps
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      INTEGER, POINTER                                   :: ptr_ghatmap
+      REAL(KIND=dp), POINTER                             :: ptr_pwout
 
       CALL timeset(routineN, handle)
 
@@ -373,13 +360,12 @@ CONTAINS
       ptr_ghatmap => pw1%pw_grid%g_hatmap(1, 1)
 
       ! invoke the combined transformation
-      CALL pw_cuda_sfffc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, ngpts, nmaps, scale)
+      CALL pw_gpu_sfffc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, ngpts, nmaps, scale)
 
       pw2%in_space = REALSPACE
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_c1dr3d_3d
+   END SUBROUTINE pw_gpu_c1dr3d_3d
 
 ! **************************************************************************************************
 !> \brief perform an parallel fft followed by a gather on the gpu
@@ -388,37 +374,30 @@ CONTAINS
 !> \param scale ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_r3dc1d_3d_ps(pw1, pw2, scale)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
-      REAL(KIND=dp)                            :: scale
+   SUBROUTINE pw_gpu_r3dc1d_3d_ps(pw1, pw2, scale)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
+      REAL(KIND=dp)                                      :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pw2)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_r3dc1d_3d_ps'
+      CHARACTER(len=*), PARAMETER :: routineN = 'pw_gpu_r3dc1d_3d_ps'
 
-      INTEGER                                  :: handle
-      INTEGER                                  :: ngpts, iout
-      INTEGER                                  :: lg, mg, mmax, lmax, rp, mx2, mz2, n1, n2, nmax
-      INTEGER                                  :: g_pos, numtask, numtask_r, numtask_g
-      INTEGER, DIMENSION(2)                    :: r_pos, r_dim
-      INTEGER                                  :: gs_group, rs_group
-      INTEGER, DIMENSION(:), POINTER           :: n, nloc
-      INTEGER, DIMENSION(:), POINTER           :: nyzray !nyzray(0:)
-      INTEGER, DIMENSION(:, :, :), POINTER       :: yzp !yzp(:,:,0:)
-      INTEGER, DIMENSION(:, :, :, :), POINTER     :: bo !bo(:,:,0:,:)
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER                                :: grays
-      COMPLEX(KIND=dp), DIMENSION(:, :, :), &
-         POINTER                                :: tbuf
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER                                :: pbuf, rbuf, sbuf, qbuf
-      INTEGER, ALLOCATABLE, DIMENSION(:)       :: p2p
-      TYPE(fft_scratch_sizes)                  :: fft_scratch_size
-      TYPE(fft_scratch_type), POINTER          :: fft_scratch
+      COMPLEX(KIND=dp), DIMENSION(:, :), POINTER         :: grays, pbuf, qbuf, rbuf, sbuf
+      COMPLEX(KIND=dp), DIMENSION(:, :, :), POINTER      :: tbuf
+      INTEGER                                            :: g_pos, gs_group, handle, iout, lg, lmax, &
+                                                            mg, mmax, mx2, mz2, n1, n2, ngpts, &
+                                                            nmax, numtask, numtask_g, numtask_r, &
+                                                            rp, rs_group
+      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: p2p
+      INTEGER, DIMENSION(2)                              :: r_dim, r_pos
+      INTEGER, DIMENSION(:), POINTER                     :: n, nloc, nyzray
+      INTEGER, DIMENSION(:, :, :), POINTER               :: yzp
+      INTEGER, DIMENSION(:, :, :, :), POINTER            :: bo
+      TYPE(fft_scratch_sizes)                            :: fft_scratch_size
+      TYPE(fft_scratch_type), POINTER                    :: fft_scratch
+
+!nyzray(0:)
+!yzp(:,:,0:)
+!bo(:,:,0:,:)
 
       CALL timeset(routineN, handle)
 
@@ -507,7 +486,7 @@ CONTAINS
             sbuf => fft_scratch%p5buf
 
             ! FFT along z
-            CALL pw_cuda_cf(pw1, qbuf)
+            CALL pw_gpu_cf(pw1, qbuf)
 
             ! Exchange data ( transpose of matrix )
             CALL cube_transpose_2(qbuf, bo(:, :, :, 1), bo(:, :, :, 2), rbuf, fft_scratch)
@@ -516,14 +495,14 @@ CONTAINS
             ! use the inbuild fft-lib
             ! CALL fft_1dm(fft_scratch%fft_plan(2), rbuf, pbuf, 1.0_dp, stat)
             ! or cufft (works faster, but is only faster if plans are stored)
-            CALL pw_cuda_f(rbuf, pbuf, +1, n(2), mx2*mz2)
+            CALL pw_gpu_f(rbuf, pbuf, +1, n(2), mx2*mz2)
 
             ! Exchange data ( transpose of matrix ) and sort
             CALL xz_to_yz(pbuf, rs_group, r_dim, g_pos, p2p, yzp, nyzray, &
                           bo(:, :, :, 2), sbuf, fft_scratch)
 
             ! FFT along x
-            CALL pw_cuda_fg(sbuf, pw2, scale)
+            CALL pw_gpu_fg(sbuf, pw2, scale)
 
             CALL release_fft_scratch(fft_scratch)
 
@@ -541,14 +520,14 @@ CONTAINS
             sbuf => fft_scratch%r1buf
 
             ! FFT along y and z
-            CALL pw_cuda_cff(pw1, tbuf)
+            CALL pw_gpu_cff(pw1, tbuf)
 
             ! Exchange data ( transpose of matrix ) and sort
             CALL yz_to_x(tbuf, gs_group, g_pos, p2p, yzp, nyzray, &
                          bo(:, :, :, 2), sbuf, fft_scratch)
 
             ! FFT along x
-            CALL pw_cuda_fg(sbuf, pw2, scale)
+            CALL pw_gpu_fg(sbuf, pw2, scale)
 
             CALL release_fft_scratch(fft_scratch)
 
@@ -558,7 +537,7 @@ CONTAINS
 
 !--------------------------------------------------------------------------
       ELSE
-         CPABORT("Not implemented (no ray_distr.) in: pw_cuda_r3dc1d_3d_ps.")
+         CPABORT("Not implemented (no ray_distr.) in: pw_gpu_r3dc1d_3d_ps.")
          !CALL fft3d ( dir, n, pwin, grays, pw1%pw_grid%para%rs_group, &
          !     pw1%pw_grid%para%bo, scale = scale, debug=test )
       END IF
@@ -566,8 +545,7 @@ CONTAINS
       pw2%in_space = RECIPROCALSPACE
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_r3dc1d_3d_ps
+   END SUBROUTINE pw_gpu_r3dc1d_3d_ps
 
 ! **************************************************************************************************
 !> \brief perform an parallel scatter followed by a fft on the gpu
@@ -576,37 +554,30 @@ CONTAINS
 !> \param scale ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_c1dr3d_3d_ps(pw1, pw2, scale)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
-      REAL(KIND=dp)                            :: scale
+   SUBROUTINE pw_gpu_c1dr3d_3d_ps(pw1, pw2, scale)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
+      REAL(KIND=dp)                                      :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pw2)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_c1dr3d_3d_ps'
+      CHARACTER(len=*), PARAMETER :: routineN = 'pw_gpu_c1dr3d_3d_ps'
 
-      INTEGER                                  :: handle
-      INTEGER                                  :: ngpts, iout
-      INTEGER                                  :: lg, mg, mmax, lmax, rp, mx2, mz2, n1, n2, nmax
-      INTEGER                                  :: g_pos, numtask, numtask_r, numtask_g
-      INTEGER, DIMENSION(2)                    :: r_pos, r_dim
-      INTEGER                                  :: gs_group, rs_group
-      INTEGER, DIMENSION(:), POINTER           :: n, nloc
-      INTEGER, DIMENSION(:), POINTER           :: nyzray !nyzray(0:)
-      INTEGER, DIMENSION(:, :, :), POINTER       :: yzp !yzp(:,:,0:)
-      INTEGER, DIMENSION(:, :, :, :), POINTER     :: bo !bo(:,:,0:,:)
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER                                :: grays
-      COMPLEX(KIND=dp), DIMENSION(:, :, :), &
-         POINTER                                :: tbuf
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER                                :: pbuf, rbuf, sbuf, qbuf
-      INTEGER, ALLOCATABLE, DIMENSION(:)       :: p2p
-      TYPE(fft_scratch_sizes)                  :: fft_scratch_size
-      TYPE(fft_scratch_type), POINTER          :: fft_scratch
+      COMPLEX(KIND=dp), DIMENSION(:, :), POINTER         :: grays, pbuf, qbuf, rbuf, sbuf
+      COMPLEX(KIND=dp), DIMENSION(:, :, :), POINTER      :: tbuf
+      INTEGER                                            :: g_pos, gs_group, handle, iout, lg, lmax, &
+                                                            mg, mmax, mx2, mz2, n1, n2, ngpts, &
+                                                            nmax, numtask, numtask_g, numtask_r, &
+                                                            rp, rs_group
+      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: p2p
+      INTEGER, DIMENSION(2)                              :: r_dim, r_pos
+      INTEGER, DIMENSION(:), POINTER                     :: n, nloc, nyzray
+      INTEGER, DIMENSION(:, :, :), POINTER               :: yzp
+      INTEGER, DIMENSION(:, :, :, :), POINTER            :: bo
+      TYPE(fft_scratch_sizes)                            :: fft_scratch_size
+      TYPE(fft_scratch_type), POINTER                    :: fft_scratch
+
+!nyzray(0:)
+!yzp(:,:,0:)
+!bo(:,:,0:,:)
 
       CALL timeset(routineN, handle)
 
@@ -695,7 +666,7 @@ CONTAINS
             sbuf => fft_scratch%p2buf
 
             ! FFT along x
-            CALL pw_cuda_sf(pw1, pbuf, scale)
+            CALL pw_gpu_sf(pw1, pbuf, scale)
 
             ! Exchange data ( transpose of matrix ) and sort
             IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(qbuf)
@@ -706,7 +677,7 @@ CONTAINS
             ! use the inbuild fft-lib
             ! CALL fft_1dm(fft_scratch%fft_plan(5), qbuf, rbuf, 1.0_dp, stat)
             ! or cufft (works faster, but is only faster if plans are stored)
-            CALL pw_cuda_f(qbuf, rbuf, -1, n(2), mx2*mz2)
+            CALL pw_gpu_f(qbuf, rbuf, -1, n(2), mx2*mz2)
 
             ! Exchange data ( transpose of matrix )
             IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(sbuf)
@@ -714,7 +685,7 @@ CONTAINS
             CALL cube_transpose_1(rbuf, bo(:, :, :, 2), bo(:, :, :, 1), sbuf, fft_scratch)
 
             ! FFT along z
-            CALL pw_cuda_fc(sbuf, pw2)
+            CALL pw_gpu_fc(sbuf, pw2)
 
             CALL release_fft_scratch(fft_scratch)
 
@@ -732,7 +703,7 @@ CONTAINS
             tbuf => fft_scratch%tbuf
 
             ! FFT along x
-            CALL pw_cuda_sf(pw1, sbuf, scale)
+            CALL pw_gpu_sf(pw1, sbuf, scale)
 
             ! Exchange data ( transpose of matrix ) and sort
             IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(tbuf)
@@ -740,7 +711,7 @@ CONTAINS
                          bo(:, :, :, 2), tbuf, fft_scratch)
 
             ! FFT along y and z
-            CALL pw_cuda_ffc(tbuf, pw2)
+            CALL pw_gpu_ffc(tbuf, pw2)
 
             CALL release_fft_scratch(fft_scratch)
 
@@ -750,7 +721,7 @@ CONTAINS
 
 !--------------------------------------------------------------------------
       ELSE
-         CPABORT("Not implemented (no ray_distr.) in: pw_cuda_c1dr3d_3d_ps.")
+         CPABORT("Not implemented (no ray_distr.) in: pw_gpu_c1dr3d_3d_ps.")
          !CALL fft3d ( dir, n, pwin, grays, pw1%pw_grid%para%rs_group, &
          !     pw1%pw_grid%para%bo, scale = scale, debug=test )
       END IF
@@ -758,8 +729,7 @@ CONTAINS
       pw2%in_space = REALSPACE
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_c1dr3d_3d_ps
+   END SUBROUTINE pw_gpu_c1dr3d_3d_ps
 
 ! **************************************************************************************************
 !> \brief perform a parallel real_to_complex copy followed by a 2D-FFT on the gpu
@@ -767,21 +737,17 @@ CONTAINS
 !> \param pwbuf ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_cff(pw1, pwbuf)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
+   SUBROUTINE pw_gpu_cff(pw1, pwbuf)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
       COMPLEX(KIND=dp), DIMENSION(:, :, :), &
-         POINTER, INTENT(INOUT)                 :: pwbuf
+         INTENT(INOUT), POINTER                          :: pwbuf
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pwbuf)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_cff'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_cff'
 
-      INTEGER                                  :: l1, l2, l3, handle
-      INTEGER, DIMENSION(:), POINTER           :: npts
-      REAL(KIND=dp), POINTER                   :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwout
+      INTEGER                                            :: handle, l1, l2, l3
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      REAL(KIND=dp), POINTER                             :: ptr_pwin
 
       CALL timeset(routineN, handle)
 
@@ -796,11 +762,10 @@ CONTAINS
       ptr_pwout => pwbuf(1, 1, 1)
 
       ! invoke the combined transformation
-      CALL pw_cuda_cff_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
+      CALL pw_gpu_cff_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_cff
+   END SUBROUTINE pw_gpu_cff
 
 ! **************************************************************************************************
 !> \brief perform a parallel 2D-FFT followed by a complex_to_real copy on the gpu
@@ -808,21 +773,17 @@ CONTAINS
 !> \param pw2 ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_ffc(pwbuf, pw2)
-      COMPLEX(KIND=dp), DIMENSION(:, :, :), &
-         POINTER, INTENT(IN)                    :: pwbuf
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
+   SUBROUTINE pw_gpu_ffc(pwbuf, pw2)
+      COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
+         POINTER                                         :: pwbuf
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pwbuf)
-      MARK_USED(pw2)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_ffc'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_ffc'
 
-      INTEGER                                  :: l1, l2, l3, handle
-      INTEGER, DIMENSION(:), POINTER           :: npts
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      REAL(KIND=dp), POINTER                   :: ptr_pwout
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin
+      INTEGER                                            :: handle, l1, l2, l3
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      REAL(KIND=dp), POINTER                             :: ptr_pwout
 
       CALL timeset(routineN, handle)
 
@@ -837,11 +798,10 @@ CONTAINS
       ptr_pwout => pw2%cr3d(l1, l2, l3)
 
       ! invoke the combined transformation
-      CALL pw_cuda_ffc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
+      CALL pw_gpu_ffc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_ffc
+   END SUBROUTINE pw_gpu_ffc
 
 ! **************************************************************************************************
 !> \brief perform a parallel real_to_complex copy followed by a 1D-FFT on the gpu
@@ -849,21 +809,17 @@ CONTAINS
 !> \param pwbuf ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_cf(pw1, pwbuf)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(INOUT)                 :: pwbuf
+   SUBROUTINE pw_gpu_cf(pw1, pwbuf)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
+         POINTER                                         :: pwbuf
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pwbuf)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_cf'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_cf'
 
-      INTEGER                                  :: l1, l2, l3, handle
-      INTEGER, DIMENSION(:), POINTER           :: npts
-      REAL(KIND=dp), POINTER                   :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwout
+      INTEGER                                            :: handle, l1, l2, l3
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      REAL(KIND=dp), POINTER                             :: ptr_pwin
 
       CALL timeset(routineN, handle)
 
@@ -878,11 +834,10 @@ CONTAINS
       ptr_pwout => pwbuf(1, 1)
 
       ! invoke the combined transformation
-      CALL pw_cuda_cf_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
+      CALL pw_gpu_cf_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_cf
+   END SUBROUTINE pw_gpu_cf
 
 ! **************************************************************************************************
 !> \brief perform a parallel 1D-FFT followed by a complex_to_real copy on the gpu
@@ -890,21 +845,17 @@ CONTAINS
 !> \param pw2 ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_fc(pwbuf, pw2)
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(IN)                    :: pwbuf
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
+   SUBROUTINE pw_gpu_fc(pwbuf, pw2)
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN), &
+         POINTER                                         :: pwbuf
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pwbuf)
-      MARK_USED(pw2)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_fc'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_fc'
 
-      INTEGER                                  :: l1, l2, l3, handle
-      INTEGER, DIMENSION(:), POINTER           :: npts
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      REAL(KIND=dp), POINTER                   :: ptr_pwout
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin
+      INTEGER                                            :: handle, l1, l2, l3
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      REAL(KIND=dp), POINTER                             :: ptr_pwout
 
       CALL timeset(routineN, handle)
 
@@ -918,11 +869,10 @@ CONTAINS
       ptr_pwout => pw2%cr3d(l1, l2, l3)
 
       ! invoke the combined transformation
-      CALL pw_cuda_fc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
+      CALL pw_gpu_fc_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), npts)
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_fc
+   END SUBROUTINE pw_gpu_fc
 
 ! **************************************************************************************************
 !> \brief perform a parallel 1D-FFT on the gpu
@@ -933,27 +883,17 @@ CONTAINS
 !> \param m ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_f(pwbuf1, pwbuf2, dir, n, m)
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(IN)                    :: pwbuf1
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(INOUT)                 :: pwbuf2
-      INTEGER, INTENT(IN)                      :: dir
-      INTEGER, INTENT(IN)                      :: n
-      INTEGER, INTENT(IN)                      :: m
+   SUBROUTINE pw_gpu_f(pwbuf1, pwbuf2, dir, n, m)
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN), &
+         POINTER                                         :: pwbuf1
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
+         POINTER                                         :: pwbuf2
+      INTEGER, INTENT(IN)                                :: dir, n, m
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pwbuf1)
-      MARK_USED(pwbuf2)
-      MARK_USED(dir)
-      MARK_USED(n)
-      MARK_USED(m)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_f'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_f'
 
-      INTEGER                                  :: handle
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin, ptr_pwout
+      INTEGER                                            :: handle
 
       CALL timeset(routineN, handle)
 
@@ -963,12 +903,11 @@ CONTAINS
          ptr_pwout => pwbuf2(1, 1)
 
          ! invoke the combined transformation
-         CALL pw_cuda_f_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), dir, n, m)
+         CALL pw_gpu_f_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), dir, n, m)
       END IF
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_f
+   END SUBROUTINE pw_gpu_f
 ! **************************************************************************************************
 !> \brief perform a parallel 1D-FFT followed by a gather on the gpu
 !> \param pwbuf ...
@@ -976,24 +915,18 @@ CONTAINS
 !> \param scale ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_fg(pwbuf, pw2, scale)
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(IN)                    :: pwbuf
-      TYPE(pw_type), TARGET, INTENT(INOUT)     :: pw2
-      REAL(KIND=dp), INTENT(IN)                :: scale
+   SUBROUTINE pw_gpu_fg(pwbuf, pw2, scale)
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN), &
+         POINTER                                         :: pwbuf
+      TYPE(pw_type), INTENT(INOUT), TARGET               :: pw2
+      REAL(KIND=dp), INTENT(IN)                          :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pwbuf)
-      MARK_USED(pw2)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_fg'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_fg'
 
-      INTEGER                                  :: mmax, mg, handle, ngpts
-      INTEGER, DIMENSION(:), POINTER           :: npts
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
-      INTEGER, POINTER                         :: ptr_ghatmap
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin, ptr_pwout
+      INTEGER                                            :: handle, mg, mmax, ngpts
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      INTEGER, POINTER                                   :: ptr_ghatmap
 
       CALL timeset(routineN, handle)
 
@@ -1012,12 +945,11 @@ CONTAINS
          ptr_ghatmap => pw2%pw_grid%g_hatmap(1, 1)
 
          ! invoke the combined transformation
-         CALL pw_cuda_fg_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, mmax, ngpts, scale)
+         CALL pw_gpu_fg_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, mmax, ngpts, scale)
       END IF
 
       CALL timestop(handle)
-#endif
-   END SUBROUTINE pw_cuda_fg
+   END SUBROUTINE pw_gpu_fg
 
 ! **************************************************************************************************
 !> \brief perform a parallel scatter followed by a 1D-FFT on the gpu
@@ -1026,25 +958,18 @@ CONTAINS
 !> \param scale ...
 !> \author Andreas Gloess
 ! **************************************************************************************************
-   SUBROUTINE pw_cuda_sf(pw1, pwbuf, scale)
-      TYPE(pw_type), TARGET, INTENT(IN)        :: pw1
-      COMPLEX(KIND=dp), DIMENSION(:, :), &
-         POINTER, INTENT(INOUT)                 :: pwbuf
-      REAL(KIND=dp), INTENT(IN)                :: scale
+   SUBROUTINE pw_gpu_sf(pw1, pwbuf, scale)
+      TYPE(pw_type), INTENT(IN), TARGET                  :: pw1
+      COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
+         POINTER                                         :: pwbuf
+      REAL(KIND=dp), INTENT(IN)                          :: scale
 
-#if ! defined (__PW_CUDA)
-      MARK_USED(pw1)
-      MARK_USED(pwbuf)
-      MARK_USED(scale)
-#else
-      CHARACTER(len=*), PARAMETER              :: routineN = 'pw_cuda_sf'
+      CHARACTER(len=*), PARAMETER                        :: routineN = 'pw_gpu_sf'
 
-      INTEGER                                  :: mmax, mg, ngpts, handle, nmaps
-      INTEGER, DIMENSION(:), POINTER           :: npts
-
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwin
-      COMPLEX(KIND=dp), POINTER                :: ptr_pwout
-      INTEGER, POINTER                         :: ptr_ghatmap
+      COMPLEX(KIND=dp), POINTER                          :: ptr_pwin, ptr_pwout
+      INTEGER                                            :: handle, mg, mmax, ngpts, nmaps
+      INTEGER, DIMENSION(:), POINTER                     :: npts
+      INTEGER, POINTER                                   :: ptr_ghatmap
 
       CALL timeset(routineN, handle)
 
@@ -1064,10 +989,11 @@ CONTAINS
          ptr_ghatmap => pw1%pw_grid%g_hatmap(1, 1)
 
          ! invoke the combined transformation
-         CALL pw_cuda_sf_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, mmax, ngpts, nmaps, scale)
+         CALL pw_gpu_sf_cu(c_loc(ptr_pwin), c_loc(ptr_pwout), c_loc(ptr_ghatmap), npts, mmax, ngpts, nmaps, scale)
       END IF
 
       CALL timestop(handle)
+   END SUBROUTINE pw_gpu_sf
+END MODULE pw_gpu
+
 #endif
-   END SUBROUTINE pw_cuda_sf
-END MODULE pw_cuda

src/pw/pw_grids.F

@@ -35,7 +35,6 @@
 ! **************************************************************************************************
 MODULE pw_grids
    USE ISO_C_BINDING,                   ONLY: C_F_POINTER,&
-                                              C_INT,&
                                               C_LOC,&
                                               C_PTR,&
                                               C_SIZE_T
@@ -48,7 +47,9 @@ MODULE pw_grids
    USE message_passing,                 ONLY: &
         mp_allgather, mp_cart_coords, mp_cart_create, mp_cart_rank, mp_comm_compare, mp_comm_dup, &
         mp_comm_free, mp_comm_self, mp_dims_create, mp_environ, mp_max, mp_min, mp_sum
-   USE offload_api,                     ONLY: offload_set_device
+   USE offload_api,                     ONLY: offload_free_pinned_mem,&
+                                              offload_malloc_pinned_mem,&
+                                              offload_set_device
    USE pw_grid_info,                    ONLY: pw_find_cutoff,&
                                               pw_grid_bounds_from_n,&
                                               pw_grid_init_setup
@@ -73,27 +74,6 @@ MODULE pw_grids
    INTEGER :: grid_tag = 0
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_grids'
 
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
-   INTERFACE
-      INTEGER(C_INT) FUNCTION cudaHostAlloc(buffer, length, flag) &
-         BIND(C, name="cudaHostAlloc")
-         IMPORT
-         IMPLICIT NONE
-         TYPE(C_PTR)              :: buffer
-         INTEGER(C_SIZE_T), VALUE :: length
-         INTEGER(C_INT), VALUE    :: flag
-      END FUNCTION cudaHostAlloc
-   END INTERFACE
-
-   INTERFACE
-      INTEGER(C_INT) FUNCTION cudaFreeHost(buffer) BIND(C, name="cudaFreeHost")
-         IMPORT
-         IMPLICIT NONE
-         TYPE(C_PTR), VALUE :: buffer
-      END FUNCTION cudaFreeHost
-   END INTERFACE
-#endif
-
    ! Distribution in g-space can be
    INTEGER, PARAMETER, PUBLIC               :: do_pw_grid_blocked_false = 0, &
                                                do_pw_grid_blocked_true = 1, &
@@ -422,7 +402,7 @@ CONTAINS
       CALL pw_grid_setup_internal(cell_hmat, cell_h_inv, cell_deth, pw_grid, bounds_local=bounds_local, &
                                   blocked=blocked, ref_grid=ref_grid, rs_dims=rs_dims, iounit=iounit)
 
-#if defined ( __PW_CUDA )
+#if defined ( __PW_GPU)
       CALL pw_grid_create_ghatmap(pw_grid)
 #endif
 
@@ -430,8 +410,8 @@ CONTAINS
 
    END SUBROUTINE pw_grid_setup
 
-#if defined ( __PW_CUDA )
-! **************************************************************************************************
+#if defined ( __PW_GPU)
+   ! **************************************************************************************************
 !> \brief sets up a combined index for CUDA gather and scatter
 !> \param pw_grid ...
 !> \author Gloess Andreas (xx-Dec-2012)
@@ -1716,10 +1696,9 @@ CONTAINS
       CHARACTER(len=*), PARAMETER :: routineN = 'pw_grid_allocate'
 
       INTEGER                                    :: nmaps
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
       INTEGER(KIND=C_SIZE_T)                     :: length
       TYPE(C_PTR)                                :: cptr_g_hatmap
-      INTEGER(C_INT), PARAMETER                  :: cudaHostAllocDefault = 0
       INTEGER                                    :: stat
 #endif
 
@@ -1733,14 +1712,13 @@ CONTAINS
 
       nmaps = 1
       IF (pw_grid%grid_span == HALFSPACE) nmaps = 2
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU )
       CALL offload_set_device()
+
       length = INT(int_size*MAX(ng, 1)*MAX(nmaps, 1), KIND=C_SIZE_T)
-      stat = cudaHostAlloc(cptr_g_hatmap, length, cudaHostAllocDefault)
+      stat = offload_malloc_pinned_mem(cptr_g_hatmap, length)
       CPASSERT(stat == 0)
       CALL c_f_pointer(cptr_g_hatmap, pw_grid%g_hatmap, (/MAX(ng, 1), MAX(nmaps, 1)/))
-#elif defined ( __PW_CUDA ) && defined ( __PW_CUDA_NO_HOSTALLOC )
-      ALLOCATE (pw_grid%g_hatmap(ng, nmaps))
 #else
       ALLOCATE (pw_grid%g_hatmap(1, 1))
 #endif
@@ -2192,7 +2170,7 @@ CONTAINS
 
       TYPE(pw_grid_type), POINTER              :: pw_grid
 
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
       INTEGER, POINTER :: dummy_ptr
       INTEGER          :: stat
 #endif
@@ -2214,9 +2192,9 @@ CONTAINS
                DEALLOCATE (pw_grid%g_hat)
             END IF
             IF (ASSOCIATED(pw_grid%g_hatmap)) THEN
-#if defined ( __PW_CUDA ) && !defined ( __PW_CUDA_NO_HOSTALLOC )
+#if defined ( __PW_GPU)
                dummy_ptr => pw_grid%g_hatmap(1, 1)
-               stat = cudaFreeHost(c_loc(dummy_ptr))
+               stat = offload_free_pinned_mem(c_loc(dummy_ptr))
                CPASSERT(stat == 0)
 #else
                DEALLOCATE (pw_grid%g_hatmap)
@@ -2274,4 +2252,3 @@ CONTAINS
    END SUBROUTINE pw_grid_release
 
 END MODULE pw_grids
-

src/pw/pw_methods.F

@@ -23,40 +23,42 @@
 !> \author apsi
 ! **************************************************************************************************
 MODULE pw_methods
-   USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit,&
-                                              cp_to_string
-   USE fast,                            ONLY: copy_cr,&
-                                              copy_rc,&
-                                              zero_c
-   USE fft_tools,                       ONLY: BWFFT,&
-                                              FWFFT,&
-                                              fft3d
-   USE kahan_sum,                       ONLY: accurate_sum
-   USE kinds,                           ONLY: dp
-   USE machine,                         ONLY: m_memory
-   USE message_passing,                 ONLY: mp_sum
-   USE pw_copy_all,                     ONLY: pw_copy_match
-   USE pw_cuda,                         ONLY: pw_cuda_c1dr3d_3d,&
-                                              pw_cuda_c1dr3d_3d_ps,&
-                                              pw_cuda_r3dc1d_3d,&
-                                              pw_cuda_r3dc1d_3d_ps
-   USE pw_fpga,                         ONLY: pw_fpga_c1dr3d_3d_dp,&
-                                              pw_fpga_c1dr3d_3d_sp,&
-                                              pw_fpga_init_bitstream,&
-                                              pw_fpga_r3dc1d_3d_dp,&
-                                              pw_fpga_r3dc1d_3d_sp
-   USE pw_grid_types,                   ONLY: HALFSPACE,&
-                                              PW_MODE_DISTRIBUTED,&
-                                              PW_MODE_LOCAL,&
-                                              pw_grid_type
-   USE pw_types,                        ONLY: COMPLEXDATA1D,&
-                                              COMPLEXDATA3D,&
-                                              NOSPACE,&
-                                              REALDATA1D,&
-                                              REALDATA3D,&
-                                              REALSPACE,&
-                                              RECIPROCALSPACE,&
-                                              pw_type
+   USE cp_log_handling, ONLY: cp_logger_get_default_io_unit, &
+                              cp_to_string
+   USE fast, ONLY: copy_cr, &
+                   copy_rc, &
+                   zero_c
+   USE fft_tools, ONLY: BWFFT, &
+                        FWFFT, &
+                        fft3d
+   USE kahan_sum, ONLY: accurate_sum
+   USE kinds, ONLY: dp
+   USE machine, ONLY: m_memory
+   USE message_passing, ONLY: mp_sum
+   USE pw_copy_all, ONLY: pw_copy_match
+   USE pw_fpga, ONLY: pw_fpga_c1dr3d_3d_dp, &
+                      pw_fpga_c1dr3d_3d_sp, &
+                      pw_fpga_init_bitstream, &
+                      pw_fpga_r3dc1d_3d_dp, &
+                      pw_fpga_r3dc1d_3d_sp
+#ifdef __PW_GPU
+   USE pw_gpu, ONLY: pw_gpu_c1dr3d_3d, &
+                     pw_gpu_c1dr3d_3d_ps, &
+                     pw_gpu_r3dc1d_3d, &
+                     pw_gpu_r3dc1d_3d_ps
+#endif
+   USE pw_grid_types, ONLY: HALFSPACE, &
+                            PW_MODE_DISTRIBUTED, &
+                            PW_MODE_LOCAL, &
+                            pw_grid_type
+   USE pw_types, ONLY: COMPLEXDATA1D, &
+                       COMPLEXDATA3D, &
+                       NOSPACE, &
+                       REALDATA1D, &
+                       REALDATA3D, &
+                       REALSPACE, &
+                       RECIPROCALSPACE, &
+                       pw_type
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -1657,7 +1659,10 @@ CONTAINS
                                                             out_space, out_unit
       INTEGER, DIMENSION(3)                              :: nloc
       INTEGER, DIMENSION(:), POINTER                     :: n
-      LOGICAL                                            :: test, use_pw_cuda
+      LOGICAL                                            :: test
+#if defined(__PW_GPU)
+      LOGICAL                                            :: use_pw_gpu
+#endif
       REAL(KIND=dp)                                      :: norm
 
       CALL timeset(routineN, handle2)
@@ -1745,8 +1750,8 @@ CONTAINS
             CALL pw_gather(pw2, c_out)
             DEALLOCATE (c_out)
          CASE ("FW_R3DC1D")
-#if defined (__PW_CUDA)
-            CALL pw_cuda_r3dc1d_3d(pw1, pw2, scale=norm)
+#if defined (__PW_GPU)
+            CALL pw_gpu_r3dc1d_3d(pw1, pw2, scale=norm)
 #elif defined (__PW_FPGA)
             ALLOCATE (c_out(n(1), n(2), n(3)))
             ! check if bitstream for the fft size is present
@@ -1791,8 +1796,8 @@ CONTAINS
             CALL pw_scatter(pw1, c_out)
             CALL fft3d(dir, n, c_out, scale=norm, debug=test)
          CASE ("BW_C1DR3D")
-#if defined (__PW_CUDA)
-            CALL pw_cuda_c1dr3d_3d(pw1, pw2, scale=norm)
+#if defined (__PW_GPU)
+            CALL pw_gpu_c1dr3d_3d(pw1, pw2, scale=norm)
 #elif defined (__PW_FPGA)
             ALLOCATE (c_out(n(1), n(2), n(3)))
             ! check if bitstream for the fft size is present
@@ -1884,16 +1889,14 @@ CONTAINS
                WRITE (out_unit, '(A)') "  PW_GATHER : 2d -> 1d "
             CALL pw_gather(pw2, grays)
          CASE ("FW_R3DC1D")
-#if defined (__PW_CUDA)
+#if defined (__PW_GPU)
             ! (no ray dist. is not efficient in CUDA)
-            use_pw_cuda = pw1%pw_grid%para%ray_distribution
-#else
-            use_pw_cuda = .FALSE.
-#endif
-            IF (use_pw_cuda) THEN
-               CALL pw_cuda_r3dc1d_3d_ps(pw1, pw2, scale=norm)
+            use_pw_gpu = pw1%pw_grid%para%ray_distribution
+            IF (use_pw_gpu) THEN
+               CALL pw_gpu_r3dc1d_3d_ps(pw1, pw2, scale=norm)
             ELSE
-               !..prepare input
+#endif
+!..   prepare input
                nloc = pw1%pw_grid%npts_local
                ALLOCATE (c_in(nloc(1), nloc(2), nloc(3)))
                CALL copy_rc(pw1%cr3d, c_in)
@@ -1913,7 +1916,10 @@ CONTAINS
                   WRITE (out_unit, '(A)') "  PW_GATHER : 2d -> 1d "
                CALL pw_gather(pw2, grays)
                DEALLOCATE (c_in)
+
+#if defined (__PW_GPU)
             END IF
+#endif
          CASE ("BW_C1DC3D")
             !..prepare input
             IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
@@ -1933,16 +1939,14 @@ CONTAINS
             END IF
             !..prepare output (nothing to do)
          CASE ("BW_C1DR3D")
-#if defined (__PW_CUDA)
+#if defined (__PW_GPU)
             ! (no ray dist. is not efficient in CUDA)
-            use_pw_cuda = pw1%pw_grid%para%ray_distribution
-#else
-            use_pw_cuda = .FALSE.
-#endif
-            IF (use_pw_cuda) THEN
-               CALL pw_cuda_c1dr3d_3d_ps(pw1, pw2, scale=norm)
+            use_pw_gpu = pw1%pw_grid%para%ray_distribution
+            IF (use_pw_gpu) THEN
+               CALL pw_gpu_c1dr3d_3d_ps(pw1, pw2, scale=norm)
             ELSE
-               !..prepare input
+#endif
+!..   prepare input
                IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) &
                   WRITE (out_unit, '(A)') "  PW_SCATTER : 2d -> 1d "
                CALL zero_c(grays)
@@ -1964,7 +1968,9 @@ CONTAINS
                   WRITE (out_unit, '(A)') "  Real part "
                CALL copy_cr(c_in, pw2%cr3d)
                DEALLOCATE (c_in)
+#if defined (__PW_GPU)
             END IF
+#endif
          END SELECT
 
       END IF
@@ -2461,4 +2467,3 @@ CONTAINS
    END SUBROUTINE pw_set
 
 END MODULE pw_methods
-

src/start/cp2k_runs.F

@@ -7,138 +7,140 @@
 
 ! **************************************************************************************************
 MODULE cp2k_runs
-   USE atom,                            ONLY: atom_code
-   USE bibliography,                    ONLY: Hutter2014,&
-                                              cite_reference
-   USE bsse,                            ONLY: do_bsse_calculation
-   USE cell_opt,                        ONLY: cp_cell_opt
-   USE cp2k_debug,                      ONLY: cp2k_debug_energy_and_forces
-   USE cp2k_info,                       ONLY: compile_date,&
-                                              compile_revision,&
-                                              cp2k_version,&
-                                              cp2k_year
-   USE cp_control_types,                ONLY: dft_control_type
-   USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_config
-   USE cp_files,                        ONLY: close_file,&
-                                              open_file
-   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
-                                              cp_logger_get_default_io_unit,&
-                                              cp_logger_type,&
-                                              cp_logger_would_log,&
-                                              cp_note_level
-   USE cp_output_handling,              ONLY: cp_add_iter_level,&
-                                              cp_print_key_finished_output,&
-                                              cp_print_key_unit_nr,&
-                                              cp_rm_iter_level
-   USE cp_para_env,                     ONLY: cp_para_env_create,&
-                                              cp_para_env_release,&
-                                              cp_para_env_retain
-   USE cp_para_types,                   ONLY: cp_para_env_type
-   USE cp_parser_methods,               ONLY: parser_search_string
-   USE cp_parser_types,                 ONLY: cp_parser_type,&
-                                              parser_create,&
-                                              parser_release
-   USE cp_units,                        ONLY: cp_unit_set_create,&
-                                              cp_unit_set_release,&
-                                              cp_unit_set_type,&
-                                              print_all_units
-   USE dbcsr_api,                       ONLY: dbcsr_finalize_lib,&
-                                              dbcsr_init_lib,&
-                                              dbcsr_print_config,&
-                                              dbcsr_print_statistics
-   USE environment,                     ONLY: cp2k_finalize,&
-                                              cp2k_init,&
-                                              cp2k_read,&
-                                              cp2k_setup
-   USE f77_interface,                   ONLY: create_force_env,&
-                                              destroy_force_env,&
-                                              f77_default_para_env => default_para_env,&
-                                              f_env_add_defaults,&
-                                              f_env_rm_defaults,&
-                                              f_env_type
-   USE farming_methods,                 ONLY: do_deadlock,&
-                                              do_nothing,&
-                                              do_wait,&
-                                              farming_parse_input,&
-                                              get_next_job
-   USE farming_types,                   ONLY: deallocate_farming_env,&
-                                              farming_env_type,&
-                                              init_farming_env,&
-                                              job_finished,&
-                                              job_running
-   USE force_env_methods,               ONLY: force_env_calc_energy_force
-   USE force_env_types,                 ONLY: force_env_get,&
-                                              force_env_type
-   USE geo_opt,                         ONLY: cp_geo_opt
-   USE global_types,                    ONLY: global_environment_type,&
-                                              globenv_create,&
-                                              globenv_release,&
-                                              globenv_retain
-   USE grid_api,                        ONLY: grid_library_print_stats,&
-                                              grid_library_set_config
-   USE input_constants,                 ONLY: &
-        bsse_run, cell_opt_run, debug_run, do_atom, do_band, do_cp2k, do_embed, do_farming, &
-        do_fist, do_mixed, do_nnp, do_opt_basis, do_optimize_input, do_qmmm, do_qs, do_sirius, &
-        do_swarm, do_tamc, do_test, do_tree_mc, do_tree_mc_ana, driver_run, ehrenfest, &
-        electronic_spectra_run, energy_force_run, energy_run, geo_opt_run, linear_response_run, &
-        mol_dyn_run, mon_car_run, negf_run, none_run, pint_run, real_time_propagation, &
-        tree_mc_run, vib_anal
-   USE input_cp2k,                      ONLY: create_cp2k_root_section
-   USE input_cp2k_check,                ONLY: check_cp2k_input
-   USE input_cp2k_global,               ONLY: create_global_section
-   USE input_cp2k_read,                 ONLY: read_input
-   USE input_keyword_types,             ONLY: keyword_release
-   USE input_parsing,                   ONLY: section_vals_parse
-   USE input_section_types,             ONLY: &
-        section_release, section_type, section_vals_create, section_vals_get_subs_vals, &
-        section_vals_release, section_vals_retain, section_vals_type, section_vals_val_get, &
-        section_vals_write, write_section_xml
-   USE ipi_driver,                      ONLY: run_driver
-   USE kinds,                           ONLY: default_path_length,&
-                                              default_string_length,&
-                                              dp,&
-                                              int_8
-   USE library_tests,                   ONLY: lib_test
-   USE machine,                         ONLY: default_output_unit,&
-                                              m_chdir,&
-                                              m_flush,&
-                                              m_getcwd,&
-                                              m_memory,&
-                                              m_memory_max,&
-                                              m_walltime
-   USE mc_run,                          ONLY: do_mon_car
-   USE md_run,                          ONLY: qs_mol_dyn
-   USE message_passing,                 ONLY: &
-        mp_any_source, mp_bcast, mp_comm_dup, mp_comm_free, mp_comm_split, mp_environ, mp_max, &
-        mp_recv, mp_send, mp_sum, mp_sync
-   USE mscfg_methods,                   ONLY: do_mol_loop,&
-                                              loop_over_molecules
-   USE neb_methods,                     ONLY: neb
-   USE negf_methods,                    ONLY: do_negf
-   USE offload_api,                     ONLY: offload_get_device_count,&
-                                              offload_get_device_id
-   USE optimize_basis,                  ONLY: run_optimize_basis
-   USE optimize_input,                  ONLY: run_optimize_input
-   USE pint_methods,                    ONLY: do_pint_run
-   USE pw_cuda,                         ONLY: pw_cuda_finalize,&
-                                              pw_cuda_init
-   USE pw_fpga,                         ONLY: pw_fpga_finalize,&
-                                              pw_fpga_init
-   USE qs_environment_types,            ONLY: get_qs_env
-   USE qs_linres_module,                ONLY: linres_calculation
-   USE qs_tddfpt_module,                ONLY: tddfpt_calculation
-   USE reference_manager,               ONLY: print_all_references,&
-                                              print_format_html
-   USE rt_propagation,                  ONLY: rt_prop_setup
-   USE sirius_interface,                ONLY: cp_sirius_finalize,&
-                                              cp_sirius_init
-   USE string_utilities,                ONLY: html_entity_table,&
-                                              write_html_tables
-   USE swarm,                           ONLY: run_swarm
-   USE tamc_run,                        ONLY: qs_tamc
-   USE tmc_setup,                       ONLY: do_analyze_files,&
-                                              do_tmc
-   USE vibrational_analysis,            ONLY: vb_anal
+   USE atom, ONLY: atom_code
+   USE bibliography, ONLY: Hutter2014, &
+                           cite_reference
+   USE bsse, ONLY: do_bsse_calculation
+   USE cell_opt, ONLY: cp_cell_opt
+   USE cp2k_debug, ONLY: cp2k_debug_energy_and_forces
+   USE cp2k_info, ONLY: compile_date, &
+                        compile_revision, &
+                        cp2k_version, &
+                        cp2k_year
+   USE cp_control_types, ONLY: dft_control_type
+   USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_config
+   USE cp_files, ONLY: close_file, &
+                       open_file
+   USE cp_log_handling, ONLY: cp_get_default_logger, &
+                              cp_logger_get_default_io_unit, &
+                              cp_logger_type, &
+                              cp_logger_would_log, &
+                              cp_note_level
+   USE cp_output_handling, ONLY: cp_add_iter_level, &
+                                 cp_print_key_finished_output, &
+                                 cp_print_key_unit_nr, &
+                                 cp_rm_iter_level
+   USE cp_para_env, ONLY: cp_para_env_create, &
+                          cp_para_env_release, &
+                          cp_para_env_retain
+   USE cp_para_types, ONLY: cp_para_env_type
+   USE cp_parser_methods, ONLY: parser_search_string
+   USE cp_parser_types, ONLY: cp_parser_type, &
+                              parser_create, &
+                              parser_release
+   USE cp_units, ONLY: cp_unit_set_create, &
+                       cp_unit_set_release, &
+                       cp_unit_set_type, &
+                       print_all_units
+   USE dbcsr_api, ONLY: dbcsr_finalize_lib, &
+                        dbcsr_init_lib, &
+                        dbcsr_print_config, &
+                        dbcsr_print_statistics
+   USE environment, ONLY: cp2k_finalize, &
+                          cp2k_init, &
+                          cp2k_read, &
+                          cp2k_setup
+   USE f77_interface, ONLY: create_force_env, &
+                            destroy_force_env, &
+                            f77_default_para_env => default_para_env, &
+                            f_env_add_defaults, &
+                            f_env_rm_defaults, &
+                            f_env_type
+   USE farming_methods, ONLY: do_deadlock, &
+                              do_nothing, &
+                              do_wait, &
+                              farming_parse_input, &
+                              get_next_job
+   USE farming_types, ONLY: deallocate_farming_env, &
+                            farming_env_type, &
+                            init_farming_env, &
+                            job_finished, &
+                            job_running
+   USE force_env_methods, ONLY: force_env_calc_energy_force
+   USE force_env_types, ONLY: force_env_get, &
+                              force_env_type
+   USE geo_opt, ONLY: cp_geo_opt
+   USE global_types, ONLY: global_environment_type, &
+                           globenv_create, &
+                           globenv_release, &
+                           globenv_retain
+   USE grid_api, ONLY: grid_library_print_stats, &
+                       grid_library_set_config
+   USE input_constants, ONLY: &
+      bsse_run, cell_opt_run, debug_run, do_atom, do_band, do_cp2k, do_embed, do_farming, &
+      do_fist, do_mixed, do_nnp, do_opt_basis, do_optimize_input, do_qmmm, do_qs, do_sirius, &
+      do_swarm, do_tamc, do_test, do_tree_mc, do_tree_mc_ana, driver_run, ehrenfest, &
+      electronic_spectra_run, energy_force_run, energy_run, geo_opt_run, linear_response_run, &
+      mol_dyn_run, mon_car_run, negf_run, none_run, pint_run, real_time_propagation, &
+      tree_mc_run, vib_anal
+   USE input_cp2k, ONLY: create_cp2k_root_section
+   USE input_cp2k_check, ONLY: check_cp2k_input
+   USE input_cp2k_global, ONLY: create_global_section
+   USE input_cp2k_read, ONLY: read_input
+   USE input_keyword_types, ONLY: keyword_release
+   USE input_parsing, ONLY: section_vals_parse
+   USE input_section_types, ONLY: &
+      section_release, section_type, section_vals_create, section_vals_get_subs_vals, &
+      section_vals_release, section_vals_retain, section_vals_type, section_vals_val_get, &
+      section_vals_write, write_section_xml
+   USE ipi_driver, ONLY: run_driver
+   USE kinds, ONLY: default_path_length, &
+                    default_string_length, &
+                    dp, &
+                    int_8
+   USE library_tests, ONLY: lib_test
+   USE machine, ONLY: default_output_unit, &
+                      m_chdir, &
+                      m_flush, &
+                      m_getcwd, &
+                      m_memory, &
+                      m_memory_max, &
+                      m_walltime
+   USE mc_run, ONLY: do_mon_car
+   USE md_run, ONLY: qs_mol_dyn
+   USE message_passing, ONLY: &
+      mp_any_source, mp_bcast, mp_comm_dup, mp_comm_free, mp_comm_split, mp_environ, mp_max, &
+      mp_recv, mp_send, mp_sum, mp_sync
+   USE mscfg_methods, ONLY: do_mol_loop, &
+                            loop_over_molecules
+   USE neb_methods, ONLY: neb
+   USE negf_methods, ONLY: do_negf
+   USE offload_api, ONLY: offload_get_device_count, &
+                          offload_get_device_id
+   USE optimize_basis, ONLY: run_optimize_basis
+   USE optimize_input, ONLY: run_optimize_input
+   USE pint_methods, ONLY: do_pint_run
+   USE pw_fpga, ONLY: pw_fpga_finalize, &
+                      pw_fpga_init
+#ifdef __PW_GPU
+   USE pw_gpu, ONLY: pw_gpu_finalize, &
+                     pw_gpu_init
+#endif
+   USE qs_environment_types, ONLY: get_qs_env
+   USE qs_linres_module, ONLY: linres_calculation
+   USE qs_tddfpt_module, ONLY: tddfpt_calculation
+   USE reference_manager, ONLY: print_all_references, &
+                                print_format_html
+   USE rt_propagation, ONLY: rt_prop_setup
+   USE sirius_interface, ONLY: cp_sirius_finalize, &
+                               cp_sirius_init
+   USE string_utilities, ONLY: html_entity_table, &
+                               write_html_tables
+   USE swarm, ONLY: run_swarm
+   USE tamc_run, ONLY: qs_tamc
+   USE tmc_setup, ONLY: do_analyze_files, &
+                        do_tmc
+   USE vibrational_analysis, ONLY: vb_anal
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -209,7 +211,9 @@ CONTAINS
       CALL dbcsr_init_lib(mpi_comm, io_unit=output_unit)
 #endif
 
-      CALL pw_cuda_init()
+#ifdef __PW_GPU
+      CALL pw_gpu_init()
+#endif
 
       CALL pw_fpga_init()
 
@@ -279,7 +283,11 @@ CONTAINS
          ! TODO: refactor cp2k's startup code
          CALL dbcsr_finalize_lib()
          IF (method_name_id == do_sirius) CALL cp_sirius_finalize()
-         CALL pw_cuda_finalize()
+
+#if defined(__PW_GPU)
+         CALL pw_gpu_finalize()
+#endif
+
          CALL pw_fpga_finalize()
          CALL farming_run(input_declaration, root_section, para_env, initial_variables)
 #if defined(__DBCSR_ACC)
@@ -292,7 +300,10 @@ CONTAINS
 #else
          CALL dbcsr_init_lib(mpi_comm, io_unit=output_unit)
 #endif
-         CALL pw_cuda_init()
+
+#if defined(__PW_GPU)
+         CALL pw_gpu_init()
+#endif
          CALL pw_fpga_init()
          IF (method_name_id == do_sirius) CALL cp_sirius_init()
       CASE (do_opt_basis)
@@ -417,7 +428,9 @@ CONTAINS
 
       IF (method_name_id == do_sirius) CALL cp_sirius_finalize()
 
-      CALL pw_cuda_finalize()
+#if defined(__PW_GPU)
+      CALL pw_gpu_finalize()
+#endif
 
       CALL pw_fpga_finalize()
 

tools/precommit/check_file_properties.py

@@ -62,6 +62,7 @@ FLAG_EXCEPTIONS = (
     r"LIBXSMM_VERSION_MINOR",
     r"LIBXSMM_VERSION_PATCH",
     r"LIBXSMM_VERSION4",
+    r"__PW_CUDA_HIP_KERNELS",
 )
 
 FLAG_EXCEPTIONS_RE = re.compile(r"|".join(FLAG_EXCEPTIONS))

tools/toolchain/Dockerfile.hip_cuda

@@ -8,40 +8,90 @@ FROM nvidia/cuda:11.3.1-devel-ubuntu20.04
 # Setup CUDA environment.
 ENV CUDA_PATH /usr/local/cuda
 ENV LD_LIBRARY_PATH /usr/local/cuda/lib64
+ENV HIP_PLATFORM nvidia
+ENV ROCM_VER 4.5.2
+ENV HIP_DIR /opt/HIP-rocm-4.5.2
+ENV HIPAMD_DIR /opt/hipamd-rocm-4.5.2
 
 RUN export DEBIAN_FRONTEND=noninteractive DEBCONF_NONINTERACTIVE_SEEN=true \
     && apt-get update -qq && apt-get install -qq --no-install-recommends \
     ca-certificates \
     build-essential \
     cmake \
+    git \
     gfortran \
     mpich \
     libmpich-dev \
     wget \
+    libssl-dev \
     && rm -rf /var/lib/apt/lists/*
 
 # Install HIP from source because the hip-nvcc package drags in 10GB of unnecessary dependencies.
-WORKDIR /opt/
-RUN wget -q https://github.com/ROCm-Developer-Tools/HIP/archive/refs/tags/rocm-4.2.0.tar.gz \
-    && echo "ecb929e0fc2eaaf7bbd16a1446a876a15baf72419c723734f456ee62e70b4c24 rocm-4.2.0.tar.gz" |  sha256sum --check \
-    && tar -xzf rocm-*.tar.gz \
-    && cd HIP-rocm-* \
-    && mkdir build \
-    && cd build \
-    && cmake .. >/dev/null 2>&1 \
-    && make -j >/dev/null 2>&1 \
-    && make install >/dev/null 2>&1 \
-    && cd ../../ \
-    && rm -rf rocm-*.tar.gz HIP-rocm-*
+WORKDIR /opt
+
+RUN wget -q https://github.com/Kitware/CMake/releases/download/v3.20.6/cmake-3.20.6-Linux-x86_64.sh \
+    && echo "4772100c2578927eed5aa9e1a80694c0d64410448c0fda73d31b0eae18645784  cmake-3.20.6-Linux-x86_64.sh" | sha256sum --check \
+    && sh cmake-3.20.6-Linux-x86_64.sh --prefix=/usr/local --skip-license \
+    && rm -f cmake-3.20.6-Linux-x86_64.sh
+
+#RUN rm -Rf cmake*.gz \
+#    && wget -q https://github.com/Kitware/CMake/archive/refs/tags/v3.20.6.tar.gz -O cmake-3.20.6.tar.gz \
+#    &&  echo "997278dd849979e2391854a5e4c1e0be3a61961fe1af849dbc68f953579ffc5f  cmake-3.20.6.tar.gz" | sha256sum --check \
+#    && tar -xzf cmake-3.20.6.tar.gz \
+#    && cd CMake-3.20.6 \
+#    && ./bootstrap \
+#    && make -j \
+#    && make install \
+#    && cd .. \
+#    && rm -Rf cmake-*
+
+RUN cmake --version
+
+RUN rm -Rf *rocm*
+
+RUN wget -q https://github.com/ROCm-Developer-Tools/HIP/archive/refs/tags/rocm-${ROCM_VER}.tar.gz -O HIP-rocm-${ROCM_VER}.tar.gz\
+    && echo "c2113dc3c421b8084cd507d91b6fbc0170765a464b71fb0d96bb875df368f160  HIP-rocm-${ROCM_VER}.tar.gz" |  sha256sum --check \
+    && tar -xzf HIP-rocm-*.tar.gz \
+    && wget -q https://github.com/ROCm-Developer-Tools/hipamd/archive/refs/tags/rocm-${ROCM_VER}.tar.gz -O hipamd-rocm-${ROCM_VER}.tar.gz \
+    && echo "b6f35b1a1d0c466b5af28e26baf646ae63267eccc4852204db1e0c7222a39ce2  hipamd-rocm-${ROCM_VER}.tar.gz" | sha256sum --check \
+    && tar -xzf hipamd-rocm-*.tar.gz \
+    && wget -q https://github.com/ROCmSoftwarePlatform/hipBLAS/archive/refs/tags/rocm-${ROCM_VER}.tar.gz -O hipBLAS-rocm-${ROCM_VER}.tar.gz \
+    && echo "82dd82a41bbadbb2a91a2a44a5d8e0d2e4f36d3078286ed4db3549b1fb6d6978  hipBLAS-rocm-${ROCM_VER}.tar.gz" | sha256sum --check \
+    && tar -xzf hipBLAS-rocm-*.tar.gz \
+    && wget -q https://github.com/ROCmSoftwarePlatform/hipFFT/archive/refs/tags/rocm-${ROCM_VER}.tar.gz -O hipFFT-rocm-${ROCM_VER}.tar.gz \
+    && echo "32ba6a5f50cfede3777a43794371ffb1363302131d8a0382d96df90ed7bc911a  hipFFT-rocm-${ROCM_VER}.tar.gz" | sha256sum --check \
+    && tar -xzf hipFFT-rocm-*.tar.gz
+
+RUN cd ${HIPAMD_DIR} \
+    && mkdir -p build \
+    && cd  build \
+    && mkdir /opt/rocm-${ROCM_VER} \
+    && cmake -DHIP_COMMON_DIR=${HIP_DIR} -DHIP_PLATFORM=nvidia -DCMAKE_INSTALL_PREFIX=/opt/rocm-${ROCM_VER}/hip .. > /dev/null 2>&1 \
+    && make -j > /dev/null 2>&1 \
+    && make install > /dev/null 2>&1 \
+    && cd ../..
+
+RUN ls /opt/rocm-${ROCM_VER}
 
 # Install hipBLAS from source.
-RUN wget -q https://github.com/ROCmSoftwarePlatform/hipBLAS/archive/refs/tags/rocm-4.2.0.tar.gz \
-    && echo "c7ce7f69c7596b5a54e666fb1373ef41d1f896dd29260a691e2eadfa863e2b1a  rocm-4.2.0.tar.gz" | sha256sum --check \
-    && tar -xzf rocm-*.tar.gz \
-    && cd hipBLAS-rocm-* \
-    && ./install.sh --cuda --install >/dev/null 2>&1 \
+RUN cd hipBLAS-rocm-* \
+    && mkdir build \
+    && cd build \
+    && cmake -DCMAKE_INSTALL_PREFIX=/opt/rocm-${ROCM_VER} -DUSE_CUDA=YES -DCMAKE_MODULE_PATH=/opt/rocm-${ROCM_VER} -DCMAKE_MODULE_PATH=/opt/rocm-${ROCM_VER}/hip/cmake .. > /dev/null 2>&1 \
+    && make -j > /dev/null 2>&1 \
+    && make install > /dev/null 2>&1 \
     && cd .. \
-    && rm -rf rocm-*.tar.gz hipBLAS-rocm-*
+    && rm -rf hipBLAS-rocm-*
+
+ENV CPATH ${CPATH}:/opt/rocm-${ROCM_VER}/hip/include
+# Install hipFFT from source.
+RUN cd hipFFT-rocm-* \
+    && mkdir build \
+    && cd build \
+    && cmake -DCMAKE_INSTALL_PREFIX=/opt/rocm-${ROCM_VER} -DBUILD_WITH_LIB=CUDA .. > /dev/null 2>&1 \
+    && make -j > /dev/null 2>&1 \
+    && make install > /dev/null 2>&1 \
+    && rm -rf hipFFT*
 
 # This is the alternative installation path via Ubuntu packages.
 ## https://rocmdocs.amd.com/en/latest/Installation_Guide/Installation-Guide.html#ubuntu
@@ -54,9 +104,9 @@ RUN wget -q https://github.com/ROCmSoftwarePlatform/hipBLAS/archive/refs/tags/ro
 #    && rm -rf /var/lib/apt/lists/*
 
 # Setup HIP environment.
-ENV ROCM_PATH /opt/rocm
+ENV ROCM_PATH /opt/rocm-${ROCM_VER}
 ENV PATH ${PATH}:${ROCM_PATH}/hip/bin
-ENV LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${ROCM_PATH}/hip/lib:${ROCM_PATH}/hipblas/lib
+ENV LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${ROCM_PATH}/hip/lib:${ROCM_PATH}/hipblas/lib:${ROCM_PATH}/hipfft/lib:${ROCM_PATH}/hipfft/lib64:${ROCM_PATH}/hipblas/lib64
 ENV HIP_PLATFORM nvidia
 RUN hipconfig
 

tools/toolchain/scripts/generate_arch_files.sh

@@ -102,7 +102,7 @@ FCFLAGS="$G_CFLAGS \$(FCDEBFLAGS) \$(WFLAGS) \$(DFLAGS)"
 
 # TODO: Remove -Wno-vla-parameter after upgrade to gcc 11.3.
 # https://gcc.gnu.org/bugzilla//show_bug.cgi?id=101289
-CFLAGS="$G_CFLAGS -std=c11 -Wall -Wextra -Werror -Wno-vla-parameter \$(DFLAGS)"
+CFLAGS="$G_CFLAGS -std=c11 -Wall -Wextra -Werror -Wno-vla-parameter -Wno-deprecated-declarations \$(DFLAGS)"
 
 # Linker flags
 # About --whole-archive see: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=52590
@@ -114,14 +114,14 @@ LDFLAGS="IF_STATIC(${STATIC_FLAGS}|) \$(FCFLAGS) ${CP_LDFLAGS}"
 # add standard libs
 LIBS="${CP_LIBS} -lstdc++"
 
+CXXFLAGS+=" --std=c++11 \$(DFLAGS) -Wno-deprecated-declarations"
 # CUDA handling
 CUDA_LIBS="-lcudart -lnvrtc -lcuda -lcufft -lcublas -lrt IF_DEBUG(-lnvToolsExt|)"
-CUDA_DFLAGS="-D__GRID_CUDA -D__DBM_CUDA -D__DBCSR_ACC -D__PW_CUDA IF_DEBUG(-D__OFFLOAD_PROFILING|)"
+CUDA_DFLAGS="-D__GRID_CUDA -D__DBM_CUDA -D__PW_GPU -D__DBCSR_ACC -D__PW_CUDA IF_DEBUG(-D__OFFLOAD_PROFILING|)"
 if [ "${ENABLE_CUDA}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
   LIBS="${LIBS} IF_CUDA(${CUDA_LIBS}|)"
   DFLAGS="IF_CUDA(${CUDA_DFLAGS}|) ${DFLAGS}"
   NVFLAGS="-g -arch sm_${ARCH_NUM} -O3 -allow-unsupported-compiler -Xcompiler='-fopenmp' --std=c++11 \$(DFLAGS)"
-  CXXFLAGS+=" --std=c++11"
   check_command nvcc "cuda"
   check_lib -lcudart "cuda"
   check_lib -lnvrtc "cuda"
@@ -157,6 +157,9 @@ if [ "${ENABLE_HIP}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
   check_command hipcc "hip"
   check_lib -lhipblas "hip"
   add_lib_from_paths HIP_LDFLAGS "libhipblas.*" $LIB_PATHS
+  check_lib -lhipfft "hip"
+  add_lib_from_paths HIP_LDFLAGS "libhipfft.*" $LIB_PATHS
+
   PLATFORM_FLAGS=''
   HIP_INCLUDES="-I${ROCM_PATH}/hip/include -I${ROCM_PATH}/hipblas/include -I${ROCM_PATH}/include"
   case "${GPUVER}" in
@@ -170,8 +173,8 @@ if [ "${ENABLE_HIP}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
       HIP_FLAGS+="-fPIE -D__HIP_PLATFORM_AMD__ -g --offload-arch=gfx906 -O3 --std=c++11 \$(DFLAGS)"
       LIBS+=" IF_HIP(-lhipblas -lamdhip64 IF_DEBUG(-lroctx64 -lroctracer64|)|)"
       PLATFORM_FLAGS='-D__HIP_PLATFORM_AMD__'
-      DFLAGS+=' IF_HIP(-D__GRID_HIP -D__HIP_PLATFORM_AMD__ IF_DEBUG(-D__OFFLOAD_PROFILING|)|) -D__DBCSR_ACC'
-      CXXFLAGS+=" \$(DFLAGS) -fopenmp -std=c++11"
+      DFLAGS+=' IF_HIP(-D__GRID_HIP -D__HIP_PLATFORM_AMD__ -D__PW_HIP -D__PW_GPU IF_DEBUG(-D__OFFLOAD_PROFILING|)|) -D__DBCSR_ACC'
+      CXXFLAGS+=" -fopenmp -std=c++11"
       ;;
     Mi100)
       check_lib -lamdhip64 "hip"
@@ -181,10 +184,10 @@ if [ "${ENABLE_HIP}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
       check_lib -lroctracer64 "hip"
       add_lib_from_paths HIP_LDFLAGS "libroctracer64.*" $LIB_PATHS
       HIP_FLAGS+="-fPIE -D__HIP_PLATFORM_AMD__ -g --offload-arch=gfx908 -O3 --std=c++11 \$(DFLAGS)"
-      LIBS+=" IF_HIP(-lhipblas -lamdhip64 IF_DEBUG(-lroctx64 -lroctracer64|)|)"
+      LIBS+=" IF_HIP( -lhipblas -lamdhip64 IF_DEBUG(-lroctx64 -lroctracer64|)|)"
       PLATFORM_FLAGS='-D__HIP_PLATFORM_AMD__ '
-      DFLAGS+=' IF_HIP(-D__GRID_HIP -D__HIP_PLATFORM_AMD__ IF_DEBUG(-D__OFFLOAD_PROFILING|)|) -D__DBCSR_ACC'
-      CXXFLAGS+=" \$(DFLAGS) -fopenmp -std=c++11"
+      DFLAGS+=' IF_HIP(-D__GRID_HIP -D__HIP_PLATFORM_AMD__ -D__PW_HIP -D__PW_GPU IF_DEBUG(-D__OFFLOAD_PROFILING|)|) -D__DBCSR_ACC'
+      CXXFLAGS+=" -fopenmp -std=c++11"
       ;;
     *)
       check_command nvcc "cuda"
@@ -193,7 +196,7 @@ if [ "${ENABLE_HIP}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
       check_lib -lcuda "cuda"
       check_lib -lcufft "cuda"
       check_lib -lcublas "cuda"
-      DFLAGS+=" IF_HIP(-D__HIP_PLATFORM_NVIDIA__ -D__GRID_HIP|)" # -D__DBCSR_ACC
+      DFLAGS+=" IF_HIP(-D__HIP_PLATFORM_NVIDIA__ -D__GRID_HIP -D__PW_HIP -D__PW_GPU |)" # -D__DBCSR_ACC
       HIP_FLAGS=" -g -arch sm_${ARCH_NUM} -O3 -Xcompiler='-fopenmp' --std=c++11 \$(DFLAGS)"
       add_include_from_paths CUDA_CFLAGS "cuda.h" $INCLUDE_PATHS
       HIP_INCLUDES+=" -I${CUDA_PATH:-${CUDA_HOME:-/CUDA_HOME-notset}}/include"
@@ -201,7 +204,7 @@ if [ "${ENABLE_HIP}" = __TRUE__ ] && [ "${GPUVER}" != no ]; then
       CFLAGS+=" -Wno-error ${CUDA_CFLAGS}"
       CXXFLAGS+=" -Wno-error ${CUDA_CFLAGS}"
       # Set LD-flags
-      LIBS+=' -lnvrtc -lcudart -lcufft -lcublas -lcuda'
+      LIBS+=' -lhipfft -lhipblas -lhipfft -lnvrtc -lcudart -lcufft -lcublas -lcuda'
       add_lib_from_paths HIP_LDFLAGS "libcudart.*" $LIB_PATHS
       add_lib_from_paths HIP_LDFLAGS "libnvrtc.*" $LIB_PATHS
       add_lib_from_paths HIP_LDFLAGS "libcuda.*" $LIB_PATHS