DBM: introduced OpenCL support (#3321)

INSTALL.md

@@ -400,37 +400,38 @@ serving DBM/DBT and other libraries depends on DBCSR's OpenCL backend.
     necessary if CP2K/DBCSR is compiled from source. Please note, some implementations ship with
     outdated OpenCL headers which can prevent using latest features (if an application discovers
     such features only at compile-time). When building from source, for instance `libOpenCL.so` is
-    sufficient (ICD loader) at link-time. However, an Installable Client Driver (ICD) is finally
+    sufficient at link-time (ICD loader). However, an Installable Client Driver (ICD) is finally
     necessary at runtime.
   - Nvidia CUDA, AMD HIP, and Intel OneAPI are fully equipped with an OpenCL runtime (if
-    `opencl-headers` package is not installed, CPATH can be needed to point into such an
+    `opencl-headers` package is not installed, CPATH can be needed to point into the former
     installation, similarly `LIBRARY_PATH` for finding `libOpenCL.so` at link-time). Installing a
     minimal or stand-alone OpenCL is also possible, e.g., following the instructions for Debian (or
     Ubuntu) as given for every [release](https://github.com/intel/compute-runtime/releases) of the
     [Intel Compute Runtime](https://github.com/intel/compute-runtime).
-  - The environment variable `ACC_OPENCL_VERBOSE=2` prints information at runtime of CP2K about
-    kernels generated and used (and thereby checks the installation).
+  - The environment variable `ACC_OPENCL_VERBOSE` prints information at runtime of CP2K about
+    kernels generated (`ACC_OPENCL_VERBOSE=2`) or executed (`ACC_OPENCL_VERBOSE=3`) which can be
+    used to check an installation.
 - Building CP2K with OpenCL-based DBCSR
   - CP2K's toolchain supports `--enable-opencl` to select DBCSR's OpenCL backend. This can be
-    combined with `--enable-cuda` (`--gpu-ver` is then imposed) to use a GPU for CP2K's grid and
-    DBM/DBT components (no OpenCL support yet).
+    combined with `--enable-cuda` (`--gpu-ver` is then imposed) to use a GPU for CP2K's GRID and PW
+    components (no OpenCL support yet) with DBM's CUDA implementation to be preferred.
   - For manually writing an ARCH-file, add `-D__OPENCL` and `-D__DBCSR_ACC` to `CFLAGS` and add
     `-lOpenCL` to the `LIBS` variable, i.e., `OFFLOAD_CC` and `OFFLOAD_FLAGS` can duplicate `CC` and
     `CFLAGS` (no special offload compiler needed). Please also set `OFFLOAD_TARGET = opencl` to
     enable the OpenCL backend in DBCSR. For OpenCL, it is not necessary to specify a GPU version
     (e.g., `GPUVER = V100` would map/limit to
     `exts/dbcsr/src/acc/opencl/smm/params/tune_multiply_V100.csv`). In fact, `GPUVER` limits tuned
-    parameters to the specified GPU, and by default all tuned parameters are embedded
+    parameters to the specified GPU, whereas by default all tuned parameters are embedded
     (`exts/dbcsr/src/acc/opencl/smm/params/*.csv`) and applied at runtime. If auto-tuned parameters
     are not available for DBCSR, well-chosen defaults will be used to populate kernels at runtime.
-  - Refer to the toolchain method (above) for an ARCH-file that blends, e.g., OpenCL and CUDA.
   - Auto-tuned parameters are embedded into the binary, i.e., CP2K does not rely on a hard-coded
     location. Setting `OPENCL_LIBSMM_SMM_PARAMS=/path/to/csv-file` environment variable can supply
     parameters for an already built application, or `OPENCL_LIBSMM_SMM_PARAMS=0` can disable using
     tuned parameters. Refer to <https://cp2k.github.io/dbcsr/> on how to tune kernels (parameters).
-- Building CP2K with OpenCL-based DBM/DBT libraries
+- Building CP2K with OpenCL-based DBM library
   - For manually writing an ARCH-file, add `-D__OFFLOAD_OPENCL` to `CFLAGS` in addition to following
-    above instructions for "Building CP2K with OpenCL-based DBCSR".
+    above instructions for "Building CP2K with OpenCL-based DBCSR". An additional Makefile rule can
+    be necessary to transform OpenCL code into a ressource header file.
 
 ### 2y. matrix-matrix multiplication offloading on GPU using SPLA
 

src/CMakeLists.txt

@@ -1372,7 +1372,8 @@ set(CP2K_GRID_SRCS_HIP
     grid/hip/grid_hip_context.cu)
 
 set(CP2K_DBM_SRCS_GPU dbm/dbm_multiply_gpu_kernel.cu)
-set(CP2K_DBM_SRCS_GPU_C dbm/dbm_multiply_gpu.c)
+set(CP2K_DBM_SRCS_GPU_C dbm/dbm_multiply_opencl.c dbm/dbm_multiply_gpu.c)
+set(CP2K_DBM_SRCS_OPENCL dbm/dbm_multiply_opencl.cl)
 
 set(CP2K_PW_SRCS_GPU pw/gpu/pw_gpu_kernels.cu)
 

src/dbm/.gitignore

@@ -0,0 +1,2 @@
+dbm_multiply_opencl.cl.h
+dbm_miniapp.x

src/dbm/Makefile

@@ -1,44 +1,5 @@
-.PHONY : all clean
-
-all: dbm_miniapp.x
-
-clean:
-	rm -fv *.o */*.o ../offload/*.o
-
-realclean: clean
-	rm -fv *.x
-
-CFLAGS := -fopenmp -g -O3 -march=native -Wall -Wextra -Wno-vla-parameter
-
-NVCC := $(shell which nvcc)
-NVARCH := sm_70
-NVFLAGS := -g -O3 -lineinfo -arch $(NVARCH) -Wno-deprecated-gpu-targets -Xcompiler "$(CFLAGS)" -D__OFFLOAD_CUDA
-
-HIPCC := $(shell which hipcc)
-HIPARCH := gfx90a
-HIPFLAGS := -fPIE -g -O3 --offload-arch=$(HIPARCH) -Wall -Wextra -Werror -I${ROCM_PATH}/include -D__OFFLOAD_HIP -D__HIP_PLATFORM_AMD__ -D__HIP_PLATFORM_AMD__
-
-ifneq ($(LIBXSMMROOT),)
-LIBS += \
-        $(LIBXSMMROOT)/lib/libxsmm.a $(LIBXSMMROOT)/lib/libxsmmext.a \
-        -lpthread -lrt
-CFLAGS += -D__LIBXSMM -I$(LIBXSMMROOT)/include
-endif
-
-ifneq ($(MKLROOT),)
-LIBS += \
-        -L$(MKLROOT)/lib/intel64 \
-        -Wl,--start-group \
-        -lmkl_gf_lp64 \
-        -lmkl_core \
-        -lmkl_sequential \
-        -Wl,--end-group
-else
-LIBS += -lblas
-endif
-
-LIBS += -fopenmp -ldl -lstdc++ -lc -lm
-
+PROJHOME := $(abspath $(dir $(lastword $(MAKEFILE_LIST))))
+CP2KHOME := $(abspath $(PROJHOME)/../..)
 ALL_HEADERS := $(shell find . -name "*.h") $(shell find ../offload/ -name "*.h")
 ALL_OBJECTS := ../offload/offload_library.o \
         dbm_distribution.o \
@@ -51,12 +12,133 @@ ALL_OBJECTS := ../offload/offload_library.o \
         dbm_multiply_cpu.o \
         dbm_shard.o
 
+# Optimization level
+ifeq (,$(filter-out 0,$(DBG)))
+OPT ?= 3
+else
+OPT ?= 0
+endif
+
+# Optimization flag derived from OPT flag
+OPTFLAG ?= -O$(patsubst O%,%,$(OPT))
+
+CFLAGS := -fopenmp -g $(OPTFLAG) -march=native -Wall -Wextra -Wcast-qual
+
+# Some additional flags when relying on default-CC
+ifeq ($(CC),)
+CFLAGS += -Wno-vla-parameter
+endif
+
+NVCC := $(shell which nvcc 2>/dev/null)
+NVARCH := sm_70
+NVFLAGS := -g $(OPTFLAG) -lineinfo -arch $(NVARCH) -Wno-deprecated-gpu-targets -Xcompiler "$(CFLAGS)" -D__OFFLOAD_CUDA
+
+HIPCC := $(shell which hipcc 2>/dev/null)
+HIPARCH := gfx90a
+HIPFLAGS := -fPIE -g $(OPTFLAG) --offload-arch=$(HIPARCH) -Wall -Wextra -Werror -I${ROCM_PATH}/include -D__OFFLOAD_HIP -D__HIP_PLATFORM_AMD__ -D__HIP_PLATFORM_AMD__
+
+# prefer NVCC or HIPCC over OpenCL
+ifneq ($(NVCC)$(HIPCC),)
+NVCC_PATH := $(if $(NVCC),$(wildcard $(dir $(NVCC))/..))
+CUDA_FILE := $(wildcard $(NVCC_PATH)/../cuda/include/cuda.h)
+CUDA_PATH := $(if $(CUDA_FILE),$(NVCC_PATH)/../cuda,$(NVCC_PATH))
+CUDA_LIBS := $(if $(wildcard $(CUDA_PATH)/lib64),lib64,lib)
+OPENCL := 0
+endif
+ifneq ($(OPENCL),0)
+OPENCL_OFFLOAD := $(firstword $(wildcard $(CP2KHOME)/lib/*/*/exts/dbcsr/libdbcsr.a))
+OPENCL_GENKRNL := $(PROJHOME)/dbm_multiply_opencl.cl.h
+endif
+
+# Make foundational runtimes available
+LIBS += -fopenmp -ldl -lstdc++ -lc -lm
+
+
+# Make BLAS/LAPACK available
+ifneq ($(MKLROOT),)
+LIBS += \
+        -L$(MKLROOT)/lib/intel64 \
+        -Wl,--start-group \
+        -lmkl_gf_lp64 \
+        -lmkl_core \
+        -lmkl_sequential \
+        -Wl,--end-group
+else
+LIBS += -lblas
+endif
+
+.PHONY : all clean
+
+all: dbm_miniapp.x
+
+clean:
+	rm -fv *.o ../offload/*.o $(OPENCL_GENKRNL)
+
+realclean: clean
+	rm -fv *.x
+
+
+# Enable OpenCL when DBCSR library was prebuilt (assume __DBCSR_ACC).
+ifneq ($(OPENCL_OFFLOAD),)
+OPENCL_SRC := $(PROJHOME)/dbm_multiply_opencl.cl
+OPENCL_CMN := $(wildcard $(CP2KHOME)/exts/dbcsr/src/acc/opencl/common/*.h)
+OPENCL_GEN := $(CP2KHOME)/exts/dbcsr/src/acc/opencl/acc_opencl.sh
+ALL_HEADERS += $(OPENCL_GENKRNL)
+
+ALL_OBJECTS += dbm_multiply_gpu.o dbm_multiply_opencl.o
+CFLAGS += -I$(CP2KHOME)/exts/dbcsr -D__OFFLOAD_OPENCL -D__DBCSR_ACC
+LIBS += $(OPENCL_OFFLOAD) -lgfortran
+
+ifneq (Darwin,$(shell uname))
+OPENCL_LIB := $(shell ldconfig -p 2>/dev/null | grep -m1 OpenCL | rev | cut -d' ' -f1 | rev)
+ifeq (,$(OPENCL_LIB))
+OPENCL_LIB := $(wildcard /usr/lib/x86_64-linux-gnu/libOpenCL.so.1)
+endif
+ifneq (,$(CUDA_PATH))
+ifeq (,$(wildcard $(OPENCL_INC)))
+CFLAGS += -I$(CUDA_PATH)/include
+endif
+ifeq (,$(wildcard $(OPENCL_LIB)))
+LIBS += -L$(CUDA_PATH)/$(CUDA_LIBS)
+LIBS += -Wl,-rpath=$(CUDA_PATH)/$(CUDA_LIBS)
+endif
+else ifneq (,$(wildcard $(OPENCL_ROOT)/include/CL/cl.h))
+ifeq (,$(wildcard $(OPENCL_INC)))
+CFLAGS += -I$(OPENCL_ROOT)/include
+endif
+LIBS += -L$(OPENCL_ROOT)/lib64
+else ifneq (,$(wildcard $(OCL_ROOT)/include/CL/cl.h))
+ifeq (,$(wildcard $(OPENCL_INC)))
+CFLAGS += -I$(OCL_ROOT)/include
+endif
+LIBS += -L$(OCL_ROOT)/lib64
+endif
+# OPENCL_INC: directory containing CL/cl.h.
+ifneq (,$(wildcard $(OPENCL_INC)))
+CFLAGS += -I$(OPENCL_INC)
+endif
+# OPENCL_LIB: file/library to be linked
+ifneq (,$(wildcard $(OPENCL_LIB)))
+LIBS += $(OPENCL_LIB)
+else
+LIBS += -l:libOpenCL.so.1
+endif
+else # macOS
+LIBS += -framework OpenCL
+endif
+
+$(OPENCL_GENKRNL): $(OPENCL_GEN) $(OPENCL_SRC) $(OPENCL_CMN)
+	$(OPENCL_GEN) -b 6 -p "" $(OPENCL_SRC) $@
+
 
 # Enable Cuda when nvcc compiler is present.
-ifneq ($(NVCC),)
-LIBS += -lcudart -lcuda -lcublas -L${CUDA_PATH}/lib64
-CFLAGS += -I${CUDA_PATH}/include -D__OFFLOAD_CUDA
+else ifneq ($(NVCC),)
 ALL_OBJECTS += dbm_multiply_gpu.o dbm_multiply_gpu_kernel.o
+CFLAGS += -I${CUDA_PATH}/include -D__OFFLOAD_CUDA
+LIBS += -lcudart -lcuda -lcublas -L${CUDA_PATH}/$(CUDA_LIBS)
+ifneq ($(wildcard $(NVCC_PATH)/../math_libs/$(CUDA_LIBS)),)
+LIBS += -L$(NVCC_PATH)/../math_libs/$(CUDA_LIBS)
+endif
 
 %.o: %.cu $(ALL_HEADERS)
 	cd $(dir $<); $(NVCC) -c $(NVFLAGS) $(notdir $<)
@@ -64,15 +146,24 @@ ALL_OBJECTS += dbm_multiply_gpu.o dbm_multiply_gpu_kernel.o
 
 # Enable HIP/ROCm when hipcc compiler is present.
 else ifneq ($(HIPCC),)
-LIBS += -L${ROCM_PATH}/lib -lamdhip64 -lhipfft -lhipblas -lrocblas
-CFLAGS += -I${ROCM_PATH}/include -D__OFFLOAD_HIP -D__HIP_PLATFORM_AMD__
 ALL_OBJECTS += dbm_multiply_gpu.o
+CFLAGS += -I${ROCM_PATH}/include -D__OFFLOAD_HIP -D__HIP_PLATFORM_AMD__
+LIBS += -L${ROCM_PATH}/lib -lamdhip64 -lhipfft -lhipblas -lrocblas
 
 %.o: %.cu $(ALL_HEADERS)
 	cd $(dir $<); $(HIPCC) -c $(HIPFLAGS) $(notdir $<)
 endif
 
 
+# Make LIBXSMM available
+ifneq ($(LIBXSMMROOT),)
+CFLAGS += -D__LIBXSMM -I$(LIBXSMMROOT)/include
+LIBS += \
+        $(LIBXSMMROOT)/lib/libxsmm.a $(LIBXSMMROOT)/lib/libxsmmext.a \
+        -lpthread -lrt
+endif
+
+
 # Build with MPI when mpicc compiler is present.
 ifneq ($(MPICC),)
 %.o: %.c $(ALL_HEADERS)

src/dbm/dbm_multiply.c

@@ -6,6 +6,7 @@
 /*----------------------------------------------------------------------------*/
 
 #include <assert.h>
+#include <limits.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -19,11 +20,24 @@
 #include "dbm_multiply_internal.h"
 
 /*******************************************************************************
- * \brief Returns the larger of two given integer (missing from the C standard)
+ * \brief Returns the larger of two given integer (missing from the C standard).
  * \author Ole Schuett
  ******************************************************************************/
 static inline int imax(int x, int y) { return (x > y ? x : y); }
 
+/*******************************************************************************
+ * \brief Updates the min/max of a range of values (initially {INT_MAX, 0}).
+ * \author Hans Pabst
+ ******************************************************************************/
+static inline void min_max(int result[2], int value) {
+  if (value < result[0]) {
+    result[0] = value;
+  }
+  if (result[1] < value) {
+    result[1] = value;
+  }
+}
+
 /*******************************************************************************
  * \brief Private routine for computing the max filter threshold for each row.
  * \author Ole Schuett
@@ -110,7 +124,8 @@ static void backend_upload_packs(const dbm_pack_t *pack_a,
  * \author Ole Schuett
  ******************************************************************************/
 static void backend_process_batch(const int ntasks, dbm_task_t batch[ntasks],
-                                  const double alpha, const dbm_pack_t *pack_a,
+                                  const int mnk_range[3][2], const double alpha,
+                                  const dbm_pack_t *pack_a,
                                   const dbm_pack_t *pack_b, const int kshard,
                                   dbm_shard_t *shard_c,
                                   backend_context_t *ctx) {
@@ -118,9 +133,11 @@ static void backend_process_batch(const int ntasks, dbm_task_t batch[ntasks],
   (void)pack_a; // mark as used
   (void)pack_b;
   (void)shard_c;
-  dbm_multiply_gpu_process_batch(ntasks, batch, alpha, kshard, &ctx->gpu);
+  dbm_multiply_gpu_process_batch(ntasks, batch, mnk_range, alpha, kshard,
+                                 &ctx->gpu);
 #else
-  (void)kshard; // mark as used
+  (void)mnk_range; // mark as used
+  (void)kshard;
   (void)ctx;
   dbm_multiply_cpu_process_batch(ntasks, batch, alpha, pack_a, pack_b, shard_c);
 #endif
@@ -208,6 +225,7 @@ static void multiply_packs(const bool transa, const bool transb,
         const int ishard = shard_row * nshard_cols + shard_col;
         dbm_shard_t *shard_c = &matrix_c->shards[ishard];
         dbm_task_t batch[MAX_BATCH_SIZE];
+        int mnk_range[][2] = {{INT_MAX, 0}, {INT_MAX, 0}, {INT_MAX, 0}};
         int ntasks = 0;
 
         // Use a merge-join to find pairs of blocks with matching sum indices.
@@ -276,15 +294,22 @@ static void multiply_packs(const bool transa, const bool transb,
             batch[ntasks].offset_c = blk_c->offset;
             ntasks++;
 
+            // track MxN-shape covering an entire batch
+            min_max(mnk_range[0], m);
+            min_max(mnk_range[1], n);
+            min_max(mnk_range[2], k);
+
             if (ntasks == MAX_BATCH_SIZE) {
-              backend_process_batch(ntasks, batch, alpha, pack_a, pack_b,
-                                    ishard, shard_c, ctx);
+              backend_process_batch(ntasks, batch, mnk_range, alpha, pack_a,
+                                    pack_b, ishard, shard_c, ctx);
+              mnk_range[0][0] = mnk_range[1][0] = mnk_range[2][0] = INT_MAX;
+              mnk_range[0][1] = mnk_range[1][1] = mnk_range[2][1] = 0;
               ntasks = 0;
             }
           }
         }
-        backend_process_batch(ntasks, batch, alpha, pack_a, pack_b, ishard,
-                              shard_c, ctx);
+        backend_process_batch(ntasks, batch, mnk_range, alpha, pack_a, pack_b,
+                              ishard, shard_c, ctx);
       }
     }
   }

src/dbm/dbm_multiply_gpu.c

@@ -102,6 +102,7 @@ void dbm_multiply_gpu_upload_packs(const dbm_pack_t *pack_a,
  * \author Ole Schuett
  ******************************************************************************/
 void dbm_multiply_gpu_process_batch(const int ntasks, const dbm_task_t *batch,
+                                    const int mnk_range[3][2],
                                     const double alpha, const int kshard,
                                     dbm_multiply_gpu_context_t *ctx) {
 
@@ -123,7 +124,7 @@ void dbm_multiply_gpu_process_batch(const int ntasks, const dbm_task_t *batch,
   offloadEventCreate(&batch_uploaded);
   offloadEventRecord(batch_uploaded, shard_c_dev->stream);
 
-  // Reallocate shard_c_dev->data if nessecary.
+  // Reallocate shard_c_dev->data if necessary.
   if (shard_c_host->data_promised > shard_c_dev->data_allocated) {
     double *old_data_dev = shard_c_dev->data;
     shard_c_dev->data_allocated =
@@ -138,7 +139,7 @@ void dbm_multiply_gpu_process_batch(const int ntasks, const dbm_task_t *batch,
     dbm_mempool_free(old_data_dev);
   }
 
-  // Zero new blocks if nessecary.
+  // Zero new blocks if necessary.
   if (shard_c_host->data_promised > shard_c_dev->data_size) {
     const int tail = shard_c_host->data_promised - shard_c_dev->data_size;
     offloadMemsetAsync(&shard_c_dev->data[shard_c_dev->data_size], 0,
@@ -147,9 +148,9 @@ void dbm_multiply_gpu_process_batch(const int ntasks, const dbm_task_t *batch,
   }
 
   // Launch kernel.
-  dbm_multiply_gpu_launch_kernel(shard_c_dev->stream, alpha, ntasks, batch_dev,
-                                 ctx->pack_a_dev.data, ctx->pack_b_dev.data,
-                                 shard_c_dev->data);
+  dbm_multiply_gpu_launch_kernel(shard_c_dev->stream, mnk_range, alpha, ntasks,
+                                 batch_dev, ctx->pack_a_dev.data,
+                                 ctx->pack_b_dev.data, shard_c_dev->data);
   OFFLOAD_CHECK(offloadGetLastError());
 
   // Wait for batch to be uploaded before refilling it.
@@ -167,7 +168,7 @@ void dbm_multiply_gpu_download_results(dbm_multiply_gpu_context_t *ctx) {
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < ctx->nshards; i++) {
-    // Grow host buffer if nessecary.
+    // Grow host buffer if necessary.
     dbm_shard_t *shard_c_host = &ctx->shards_c_host[i];
     dbm_shard_allocate_promised_blocks(shard_c_host);
 

src/dbm/dbm_multiply_gpu.h

@@ -64,6 +64,7 @@ void dbm_multiply_gpu_upload_packs(const dbm_pack_t *pack_a,
  * \author Ole Schuett
  ******************************************************************************/
 void dbm_multiply_gpu_process_batch(const int ntasks, const dbm_task_t *batch,
+                                    const int mnk_range[3][2],
                                     const double alpha, const int kshard,
                                     dbm_multiply_gpu_context_t *ctx);
 

src/dbm/dbm_multiply_gpu_kernel.cu

@@ -199,15 +199,14 @@ __global__ static void process_batch_kernel(const double alpha,
  *        All arguments are assumed to be device pointers.
  * \author Ole Schuett
  ******************************************************************************/
-void dbm_multiply_gpu_launch_kernel(const offloadStream_t stream,
-                                    const double alpha, const int ntasks,
-                                    const dbm_task_t *batch,
-                                    const double *pack_a_data,
-                                    const double *pack_b_data,
-                                    double *shard_c_data) {
+void dbm_multiply_gpu_launch_kernel(
+    const offloadStream_t stream, const int mnk_range[3][2], const double alpha,
+    const int ntasks, const dbm_task_t *batch, const double *pack_a_data,
+    const double *pack_b_data, double *shard_c_data) {
   const int nblocks = ntasks; // TODO tune launch parameters.
   const int threads_per_block = NUM_THREADS;
   const size_t smem_per_block = 0;
+  (void)mnk_range; // mark used
   process_batch_kernel<<<nblocks, threads_per_block, smem_per_block, stream>>>(
       alpha, batch, pack_a_data, pack_b_data, shard_c_data);
 }

src/dbm/dbm_multiply_gpu_kernel.h

@@ -22,12 +22,10 @@ extern "C" {
  *        All arguments are assumed to be device pointers.
  * \author Ole Schuett
  ******************************************************************************/
-void dbm_multiply_gpu_launch_kernel(const offloadStream_t stream,
-                                    const double alpha, const int ntasks,
-                                    const dbm_task_t *batch,
-                                    const double *pack_a_data,
-                                    const double *pack_b_data,
-                                    double *shard_c_data);
+void dbm_multiply_gpu_launch_kernel(
+    const offloadStream_t stream, const int mnk_range[3][2], const double alpha,
+    const int ntasks, const dbm_task_t *batch, const double *pack_a_data,
+    const double *pack_b_data, double *shard_c_data);
 
 #ifdef __cplusplus
 }

src/dbm/dbm_multiply_opencl.c

@@ -0,0 +1,119 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2024 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: BSD-3-Clause                                     */
+/*----------------------------------------------------------------------------*/
+
+#include "../offload/offload_runtime.h"
+#if defined(__OFFLOAD_OPENCL) && !defined(__NO_OFFLOAD_DBM)
+
+#include "dbm_multiply_gpu_kernel.h"
+#include "dbm_multiply_opencl.cl.h"
+
+void dbm_multiply_gpu_launch_kernel(const offloadStream_t stream,
+                                    const int mnk_range[3][2], double alpha,
+                                    int ntasks, const dbm_task_t *tasks,
+                                    const double *pack_a_data,
+                                    const double *pack_b_data,
+                                    double *shard_c_data) {
+  static cl_kernel kernel = NULL;
+  static int task_split = 0;
+  int result = EXIT_SUCCESS, verbosity = c_dbcsr_acc_opencl_config.verbosity;
+  cl_event event, *const perf_event =
+                      ((0 <= verbosity && 2 >= verbosity) ? NULL : &event);
+  const c_dbcsr_acc_opencl_stream_t *const str = ACC_OPENCL_STREAM(stream);
+  const int max_m = mnk_range[0][1], max_n = mnk_range[1][1];
+  const size_t work_tasks = ntasks, wgsize = 0;
+  size_t ibatch = 0, iadata = 0, ibdata = 0, icdata = 0, work_size = 0;
+  c_dbcsr_acc_opencl_info_memptr_t adata, bdata, cdata, batch;
+  assert(NULL != pack_a_data && NULL != pack_b_data && NULL != shard_c_data);
+  assert(0 < mnk_range[0][0] && 0 < max_m && mnk_range[0][0] <= max_m);
+  assert(0 < mnk_range[1][0] && 0 < max_n && mnk_range[1][0] <= max_n);
+  assert(0 < mnk_range[2][0] && 0 < mnk_range[2][1] &&
+         mnk_range[2][0] <= mnk_range[2][1]);
+  assert(NULL != str && NULL != str->queue);
+  assert(0 < ntasks && NULL != tasks);
+  /* creating/calling kernel must be consistent across threads */
+  ACC_OPENCL_ACQUIRE(c_dbcsr_acc_opencl_config.lock_main);
+#if defined(OPENCL_DBM_SOURCE_MULTIPLY_OPENCL)
+  if (NULL == kernel) { /* first-time check if kernel is present */
+    char params[ACC_OPENCL_BUFFERSIZE] = "-DLU=0 -DBN=4";
+    const char *const flags = "-cl-fast-relaxed-math -cl-denorms-are-zero";
+    const char *extensions[] = {NULL, NULL};
+    const size_t nextensions = sizeof(extensions) / sizeof(*extensions);
+    const size_t offset = strlen(params);
+    const libxsmm_timer_tickint start = libxsmm_timer_tick();
+    const int nchar = c_dbcsr_acc_opencl_flags_atomics(
+        &c_dbcsr_acc_opencl_config.device, c_dbcsr_acc_opencl_atomic_fp_64,
+        extensions, nextensions, params + offset, sizeof(params) - offset);
+    if (0 < nchar && sizeof(params) > (offset + nchar)) {
+      const char *const task_split_env = getenv("LIBDBM_TASK_SPLIT");
+      task_split = (NULL == task_split_env ? 1 : atoi(task_split_env));
+      result |= c_dbcsr_acc_opencl_kernel(
+          0 /*source_is_file*/, OPENCL_DBM_SOURCE_MULTIPLY_OPENCL,
+          "dbm_multiply", params,
+          0 == c_dbcsr_acc_opencl_config.debug ? flags : NULL, NULL /*try*/,
+          NULL /*try_ok*/, extensions, nextensions, &kernel);
+      if (2 <= verbosity || 0 > verbosity) {
+        if (EXIT_SUCCESS == result) {
+          fprintf(stderr, "INFO ACC/LIBDBM: DBM-kernel ms=%.1f\n",
+                  1E3 * libxsmm_timer_duration(start, libxsmm_timer_tick()));
+        } else {
+          fprintf(stderr, "INFO ACC/LIBDBM: DBM-kernel failed to generate\n");
+        }
+      }
+    }
+  }
+#else
+#error "OpenCL kernel code not found!"
+#endif
+  result |= c_dbcsr_acc_opencl_info_devptr_lock(&adata, NULL /*lock*/,
+                                                pack_a_data, 1 /*esize*/,
+                                                NULL /*amount*/, &iadata);
+  result |= c_dbcsr_acc_opencl_info_devptr_lock(&bdata, NULL /*lock*/,
+                                                pack_b_data, 1 /*esize*/,
+                                                NULL /*amount*/, &ibdata);
+  result |= c_dbcsr_acc_opencl_info_devptr_lock(&cdata, NULL /*lock*/,
+                                                shard_c_data, 1 /*esize*/,
+                                                NULL /*amount*/, &icdata);
+  result |= c_dbcsr_acc_opencl_info_devptr_lock(
+      &batch, NULL /*lock*/, tasks /*batch*/, sizeof(dbm_task_t), &work_tasks,
+      &ibatch);
+  assert(0 == iadata && 0 == ibdata && 0 == icdata);
+  work_size = (0 != task_split ? (work_tasks * max_m) : work_tasks);
+  result |= clSetKernelArg(kernel, 0, sizeof(cl_double), &alpha);
+  result |= clSetKernelArg(kernel, 1, sizeof(cl_int), &ibatch);
+  result |= clSetKernelArg(kernel, 2, sizeof(cl_int), &ntasks);
+  result |= c_dbcsr_acc_opencl_set_kernel_ptr(kernel, 3, batch.memory);
+  result |= c_dbcsr_acc_opencl_set_kernel_ptr(kernel, 4, adata.memory);
+  result |= c_dbcsr_acc_opencl_set_kernel_ptr(kernel, 5, bdata.memory);
+  result |= c_dbcsr_acc_opencl_set_kernel_ptr(kernel, 6, cdata.memory);
+  result |= clEnqueueNDRangeKernel(str->queue, kernel, 1 /*work_dim*/,
+                                   NULL /*offset*/, &work_size,
+                                   0 != wgsize ? &wgsize : NULL, 0 /*num_wait*/,
+                                   NULL /*wait_list*/, perf_event);
+  if (NULL != perf_event && EXIT_SUCCESS == result) {
+    cl_ulong begin = 0, end = 0;
+    clWaitForEvents(1, perf_event);
+    result |= clGetEventProfilingInfo(*perf_event, CL_PROFILING_COMMAND_START,
+                                      sizeof(cl_ulong), &begin, NULL);
+    result |= clGetEventProfilingInfo(*perf_event, CL_PROFILING_COMMAND_END,
+                                      sizeof(cl_ulong), &end, NULL);
+    if (EXIT_SUCCESS == result) {
+      const double duration_ms = 1E-6 * LIBXSMM_DELTA(begin, end);
+      const double gflops =
+          (1E-6 * max_m * max_n * mnk_range[2][1] * ntasks) / duration_ms;
+      fprintf(stderr,
+              "INFO ACC/LIBDBM: DBM-kernel mnk=%ix%ix%i "
+              "ntasks=%i gflops=%.1f ms=%.2g\n",
+              max_m, max_n, mnk_range[2][1], ntasks, gflops, duration_ms);
+    }
+  }
+  ACC_OPENCL_RELEASE(c_dbcsr_acc_opencl_config.lock_main);
+  OFFLOAD_CHECK(result);
+}
+
+#endif // defined(__OFFLOAD_OPENCL) && !defined(__NO_OFFLOAD_DBM)
+
+// EOF

src/dbm/dbm_multiply_opencl.cl

@@ -0,0 +1,77 @@
+/*----------------------------------------------------------------------------*/
+/*  CP2K: A general program to perform molecular dynamics simulations         */
+/*  Copyright 2000-2024 CP2K developers group <https://cp2k.org>              */
+/*                                                                            */
+/*  SPDX-License-Identifier: BSD-3-Clause                                     */
+/*----------------------------------------------------------------------------*/
+#include "../../exts/dbcsr/src/acc/opencl/common/opencl_atomics.h"
+#include "dbm_multiply_internal.h"
+
+#define IDX(I, J, K, M, N) ((I) * (N) + (J) + (K))
+#define IDT(I, J, K, M, N) IDX(J, I, K, N, M)
+
+#define DBM_MULTIPLY_KERNEL(TASK, A, B, CV, M, N0, N1, UNROLL_N, UNROLL_K)     \
+  do {                                                                         \
+    UNROLL_K                                                                   \
+    for (int k = 0; k < (TASK).k; ++k) {                                       \
+      const int ia = IDT(M, k, (TASK).offset_a, (TASK).m, (TASK).k);           \
+      const double a = (A)[ia];                                                \
+      UNROLL_N                                                                 \
+      for (int n = 0; n < (N1); ++n) {                                         \
+        const int ib = IDX(k, n + (N0), (TASK).offset_b, (TASK).k, (TASK).n);  \
+        (CV)[n] = MAD(a, (B)[ib], (CV)[n]);                                    \
+      }                                                                        \
+    }                                                                          \
+  } while (0)
+
+#define DBM_MULTIPLY_ACCUMULATE(ALPHA, TASK, CV, C, M, BN, N0)                 \
+  do { /* flush private accumulator to global memory using atomics */          \
+    UNROLL_FORCE(BN)                                                           \
+    for (int n = 0; n < (BN); ++n) {                                           \
+      const int ic = IDT(M, n + (N0), (TASK).offset_c, (TASK).m, (TASK).n);    \
+      ACCUMULATE((C) + ic, (ALPHA) * (CV)[n]);                                 \
+      cv[n] = ZERO; /* reset */                                                \
+    }                                                                          \
+  } while (0)
+
+kernel void dbm_multiply(double alpha, int itask, int ntasks,
+                         global const dbm_task_t *tasks,
+                         global const double *restrict am,
+                         global const double *restrict bm,
+                         global double *restrict cm) {
+  double cv[BN] = {0}; /* private accumulator */
+  const int size = (int)get_global_size(0), i = (int)get_global_id(0);
+
+  if (size != ntasks) { /* LIBDBM_TASK_SPLIT */
+    const int max_m = size / ntasks, tid = i / max_m;
+    const dbm_task_t task = tasks[itask + min(tid, ntasks - 1)]; /* copy */
+    const int m = i - tid * max_m;
+    if (m < task.m) {
+      if ((BN) < task.n) {
+        UNROLL_AUTO
+        for (int n0 = 0; n0 < task.n; n0 += (BN)) {
+          const int n1 = min(BN, task.n - n0);
+          DBM_MULTIPLY_KERNEL(task, am, bm, cv, m, n0, n1, UNROLL_FORCE(BN),
+                              UNROLL_AUTO);
+          DBM_MULTIPLY_ACCUMULATE(alpha, task, cv, cm, m, BN, n0);
+        }
+      } else { /* small */
+        DBM_MULTIPLY_KERNEL(task, am, bm, cv, m, 0, task.n, UNROLL_FORCE(BN),
+                            UNROLL_FORCE(BN));
+        DBM_MULTIPLY_ACCUMULATE(alpha, task, cv, cm, m, BN, 0);
+      }
+    }
+  } else { /* full matrix multiplication */
+    const dbm_task_t task = tasks[itask + i];
+    UNROLL_OUTER(1)
+    for (int m = 0; m < task.m; ++m) {
+      UNROLL(1)
+      for (int n0 = 0; n0 < task.n; n0 += (BN)) {
+        const int n1 = min(BN, task.n - n0);
+        DBM_MULTIPLY_KERNEL(task, am, bm, cv, m, n0, n1, UNROLL_AUTO,
+                            UNROLL_AUTO);
+        DBM_MULTIPLY_ACCUMULATE(alpha, task, cv, cm, m, BN, n0);
+      }
+    }
+  }
+}

tools/precommit/README.md

@@ -6,7 +6,7 @@ The precommit system consists of the following tools that analyze and format the
 - [prettify](../prettify/) to format Fortran files.
 - [analyze_src](../conventions/analyze_src.py) to check copyright banners and a few other things.
 - [ast.parse](https://docs.python.org/3/library/ast.html) to check Python syntax.
-- [clang-format](https://clang.llvm.org/docs/ClangFormat.html) to format C and Cuda files.
+- [clang-format](https://clang.llvm.org/docs/ClangFormat.html) to format C, Cuda, and OpenCL files.
 - [black](https://github.com/psf/black) to format Python scripts.
 - [shfmt](https://github.com/mvdan/sh) to format Shell scripts.
 - [shellcheck](https://github.com/koalaman/shellcheck) to analyze Shell scripts.

tools/precommit/check_file_properties.py

@@ -21,7 +21,9 @@ FLAG_EXCEPTIONS = (
     r"\$\{.*\}\$",
     r"__.*__",
     r"CUDA_VERSION",
+    r"DBM_LIBXSMM_PREFETCH",
     r"DBM_VALIDATE_AGAINST_DBCSR",
+    r"OPENCL_DBM_SOURCE_MULTIPLY_OPENCL",
     r"FD_DEBUG",
     r"GRID_DO_COLLOCATE",
     r"INTEL_MKL_VERSION",
@@ -58,7 +60,6 @@ FLAG_EXCEPTIONS = (
     r"LIBXSMM_VERSION2",
     r"LIBXSMM_VERSION3",
     r"LIBXSMM_VERSION4",
-    r"DBM_LIBXSMM_PREFETCH",
 )
 
 FLAG_EXCEPTIONS_RE = re.compile(r"|".join(FLAG_EXCEPTIONS))

tools/precommit/precommit.py

@@ -194,7 +194,7 @@ def process_file(fn: str, allow_modifications: bool) -> None:
         run_local_tool("./tools/doxify/doxify.sh", fn)
         run_prettify(fn)
 
-    if re.match(r".*\.(c|cu|h)$", fn):
+    if re.match(r".*\.(c|cu|cl|h)$", fn):
         run_remote_tool("clangformat", fn)
 
     if re.match(r".*\.(cc|cpp|cxx|hcc|hpp|hxx)$", fn):