src/gromacs/nbnxm/cuda/nbnxm_cuda.cu

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \file
  36  *  \brief Define CUDA implementation of nbnxn_gpu.h
  37  *
  38  *  \author Szilard Pall <pall.szilard@gmail.com>
  39  */
  40 #include "gmxpre.h"
  41
  42 #include "config.h"
  43
  44 #include <assert.h>
  45 #include <stdlib.h>
  46
  47 #include "gromacs/nbnxm/nbnxm_gpu.h"
  48
  49 #if defined(_MSVC)
  50 #include <limits>
  51 #endif
  52
  53
  54 #include "nbnxm_cuda.h"
  55
  56 #include "gromacs/gpu_utils/cudautils.cuh"
  57 #include "gromacs/gpu_utils/gpueventsynchronizer.cuh"
  58 #include "gromacs/gpu_utils/vectype_ops.cuh"
  59 #include "gromacs/mdlib/ppforceworkload.h"
  60 #include "gromacs/nbnxm/atomdata.h"
  61 #include "gromacs/nbnxm/gpu_common.h"
  62 #include "gromacs/nbnxm/gpu_common_utils.h"
  63 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  64 #include "gromacs/nbnxm/grid.h"
  65 #include "gromacs/nbnxm/nbnxm.h"
  66 #include "gromacs/nbnxm/pairlist.h"
  67 #include "gromacs/nbnxm/cuda/nbnxm_buffer_ops_kernels.cuh"
  68 #include "gromacs/timing/gpu_timing.h"
  69 #include "gromacs/utility/cstringutil.h"
  70 #include "gromacs/utility/gmxassert.h"
  71
  72 #include "nbnxm_cuda_types.h"
  73
  74 /***** The kernel declarations/definitions come here *****/
  75
  76 /* Top-level kernel declaration generation: will generate through multiple
  77  * inclusion the following flavors for all kernel declarations:
  78  * - force-only output;
  79  * - force and energy output;
  80  * - force-only with pair list pruning;
  81  * - force and energy output with pair list pruning.
  82  */
  83 #define FUNCTION_DECLARATION_ONLY
  84 /** Force only **/
  85 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  86 /** Force & energy **/
  87 #define CALC_ENERGIES
  88 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  89 #undef CALC_ENERGIES
  90
  91 /*** Pair-list pruning kernels ***/
  92 /** Force only **/
  93 #define PRUNE_NBL
  94 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  95 /** Force & energy **/
  96 #define CALC_ENERGIES
  97 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  98 #undef CALC_ENERGIES
  99 #undef PRUNE_NBL
 100
 101 /* Prune-only kernels */
 102 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_pruneonly.cuh"
 103 #undef FUNCTION_DECLARATION_ONLY
 104
 105 /* Now generate the function definitions if we are using a single compilation unit. */
 106 #if GMX_CUDA_NB_SINGLE_COMPILATION_UNIT
 107 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_F_noprune.cu"
 108 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_F_prune.cu"
 109 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_VF_noprune.cu"
 110 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_VF_prune.cu"
 111 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_pruneonly.cu"
 112 #endif /* GMX_CUDA_NB_SINGLE_COMPILATION_UNIT */
 113
 114
 115 namespace Nbnxm
 116 {
 117
 118 //! Number of CUDA threads in a block
 119 //TODO Optimize this through experimentation
 120 constexpr static int c_bufOpsThreadsPerBlock = 128;
 121
 122 /*! Nonbonded kernel function pointer type */
 123 typedef void (*nbnxn_cu_kfunc_ptr_t)(const cu_atomdata_t,
 124                                      const cu_nbparam_t,
 125                                      const cu_plist_t,
 126                                      bool);
 127
 128 /*********************************/
 129
 130 /*! Returns the number of blocks to be used for the nonbonded GPU kernel. */
 131 static inline int calc_nb_kernel_nblock(int nwork_units, const gmx_device_info_t *dinfo)
 132 {
 133     int max_grid_x_size;
 134
 135     assert(dinfo);
 136     /* CUDA does not accept grid dimension of 0 (which can happen e.g. with an
 137        empty domain) and that case should be handled before this point. */
 138     assert(nwork_units > 0);
 139
 140     max_grid_x_size = dinfo->prop.maxGridSize[0];
 141
 142     /* do we exceed the grid x dimension limit? */
 143     if (nwork_units > max_grid_x_size)
 144     {
 145         gmx_fatal(FARGS, "Watch out, the input system is too large to simulate!\n"
 146                   "The number of nonbonded work units (=number of super-clusters) exceeds the"
 147                   "maximum grid size in x dimension (%d > %d)!", nwork_units, max_grid_x_size);
 148     }
 149
 150     return nwork_units;
 151 }
 152
 153
 154 /* Constant arrays listing all kernel function pointers and enabling selection
 155    of a kernel in an elegant manner. */
 156
 157 /*! Pointers to the non-bonded kernels organized in 2-dim arrays by:
 158  *  electrostatics and VDW type.
 159  *
 160  *  Note that the row- and column-order of function pointers has to match the
 161  *  order of corresponding enumerated electrostatics and vdw types, resp.,
 162  *  defined in nbnxn_cuda_types.h.
 163  */
 164
 165 /*! Force-only kernel function pointers. */
 166 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_noprune_ptr[eelCuNR][evdwCuNR] =
 167 {
 168     { nbnxn_kernel_ElecCut_VdwLJ_F_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_F_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_F_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_F_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_F_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_cuda            },
 169     { nbnxn_kernel_ElecRF_VdwLJ_F_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_F_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_F_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_F_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_F_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_cuda             },
 170     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_cuda        },
 171     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_cuda },
 172     { nbnxn_kernel_ElecEw_VdwLJ_F_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_F_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_F_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_F_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_F_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_cuda             },
 173     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_cuda      }
 174 };
 175
 176 /*! Force + energy kernel function pointers. */
 177 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_noprune_ptr[eelCuNR][evdwCuNR] =
 178 {
 179     { nbnxn_kernel_ElecCut_VdwLJ_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_cuda            },
 180     { nbnxn_kernel_ElecRF_VdwLJ_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_cuda             },
 181     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_cuda        },
 182     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_cuda },
 183     { nbnxn_kernel_ElecEw_VdwLJ_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_cuda             },
 184     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_cuda      }
 185 };
 186
 187 /*! Force + pruning kernel function pointers. */
 188 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_prune_ptr[eelCuNR][evdwCuNR] =
 189 {
 190     { nbnxn_kernel_ElecCut_VdwLJ_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_prune_cuda             },
 191     { nbnxn_kernel_ElecRF_VdwLJ_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_prune_cuda              },
 192     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_prune_cuda         },
 193     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_prune_cuda  },
 194     { nbnxn_kernel_ElecEw_VdwLJ_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_prune_cuda              },
 195     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_prune_cuda       }
 196 };
 197
 198 /*! Force + energy + pruning kernel function pointers. */
 199 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_prune_ptr[eelCuNR][evdwCuNR] =
 200 {
 201     { nbnxn_kernel_ElecCut_VdwLJ_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_prune_cuda            },
 202     { nbnxn_kernel_ElecRF_VdwLJ_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_prune_cuda             },
 203     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_prune_cuda        },
 204     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_prune_cuda },
 205     { nbnxn_kernel_ElecEw_VdwLJ_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_prune_cuda             },
 206     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_prune_cuda      }
 207 };
 208
 209 /*! Return a pointer to the kernel version to be executed at the current step. */
 210 static inline nbnxn_cu_kfunc_ptr_t select_nbnxn_kernel(int                                  eeltype,
 211                                                        int                                  evdwtype,
 212                                                        bool                                 bDoEne,
 213                                                        bool                                 bDoPrune,
 214                                                        const gmx_device_info_t gmx_unused  *devInfo)
 215 {
 216     nbnxn_cu_kfunc_ptr_t res;
 217
 218     GMX_ASSERT(eeltype < eelCuNR,
 219                "The electrostatics type requested is not implemented in the CUDA kernels.");
 220     GMX_ASSERT(evdwtype < evdwCuNR,
 221                "The VdW type requested is not implemented in the CUDA kernels.");
 222
 223     /* assert assumptions made by the kernels */
 224     GMX_ASSERT(c_nbnxnGpuClusterSize*c_nbnxnGpuClusterSize/c_nbnxnGpuClusterpairSplit == devInfo->prop.warpSize,
 225                "The CUDA kernels require the cluster_size_i*cluster_size_j/nbnxn_gpu_clusterpair_split to match the warp size of the architecture targeted.");
 226
 227     if (bDoEne)
 228     {
 229         if (bDoPrune)
 230         {
 231             res = nb_kfunc_ener_prune_ptr[eeltype][evdwtype];
 232         }
 233         else
 234         {
 235             res = nb_kfunc_ener_noprune_ptr[eeltype][evdwtype];
 236         }
 237     }
 238     else
 239     {
 240         if (bDoPrune)
 241         {
 242             res = nb_kfunc_noener_prune_ptr[eeltype][evdwtype];
 243         }
 244         else
 245         {
 246             res = nb_kfunc_noener_noprune_ptr[eeltype][evdwtype];
 247         }
 248     }
 249
 250     return res;
 251 }
 252
 253 /*! \brief Calculates the amount of shared memory required by the nonbonded kernel in use. */
 254 static inline int calc_shmem_required_nonbonded(const int num_threads_z, const gmx_device_info_t gmx_unused *dinfo, const cu_nbparam_t *nbp)
 255 {
 256     int shmem;
 257
 258     assert(dinfo);
 259
 260     /* size of shmem (force-buffers/xq/atom type preloading) */
 261     /* NOTE: with the default kernel on sm3.0 we need shmem only for pre-loading */
 262     /* i-atom x+q in shared memory */
 263     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 264     /* cj in shared memory, for each warp separately */
 265     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 266
 267     if (nbp->vdwtype == evdwCuCUTCOMBGEOM ||
 268         nbp->vdwtype == evdwCuCUTCOMBLB)
 269     {
 270         /* i-atom LJ combination parameters in shared memory */
 271         shmem += c_numClPerSupercl * c_clSize * sizeof(float2);
 272     }
 273     else
 274     {
 275         /* i-atom types in shared memory */
 276         shmem += c_numClPerSupercl * c_clSize * sizeof(int);
 277     }
 278
 279     return shmem;
 280 }
 281
 282 /*! \brief Sync the nonlocal stream with dependent tasks in the local queue.
 283  *
 284  *  As the point where the local stream tasks can be considered complete happens
 285  *  at the same call point where the nonlocal stream should be synced with the
 286  *  the local, this function records the event if called with the local stream as
 287  *  argument and inserts in the GPU stream a wait on the event on the nonlocal.
 288  */
 289 void nbnxnInsertNonlocalGpuDependency(const gmx_nbnxn_cuda_t   *nb,
 290                                       const InteractionLocality interactionLocality)
 291 {
 292     cudaStream_t stream  = nb->stream[interactionLocality];
 293
 294     /* When we get here all misc operations issued in the local stream as well as
 295        the local xq H2D are done,
 296        so we record that in the local stream and wait for it in the nonlocal one.
 297        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 298        compute on interactions between local and nonlocal atoms.
 299      */
 300     if (nb->bUseTwoStreams)
 301     {
 302         if (interactionLocality == InteractionLocality::Local)
 303         {
 304             cudaError_t stat = cudaEventRecord(nb->misc_ops_and_local_H2D_done, stream);
 305             CU_RET_ERR(stat, "cudaEventRecord on misc_ops_and_local_H2D_done failed");
 306         }
 307         else
 308         {
 309             cudaError_t stat = cudaStreamWaitEvent(stream, nb->misc_ops_and_local_H2D_done, 0);
 310             CU_RET_ERR(stat, "cudaStreamWaitEvent on misc_ops_and_local_H2D_done failed");
 311         }
 312     }
 313 }
 314
 315 /*! \brief Launch asynchronously the xq buffer host to device copy. */
 316 void gpu_copy_xq_to_gpu(gmx_nbnxn_cuda_t       *nb,
 317                         const nbnxn_atomdata_t *nbatom,
 318                         const AtomLocality      atomLocality)
 319 {
 320     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 321
 322     GMX_ASSERT(atomLocality == AtomLocality::Local || atomLocality == AtomLocality::NonLocal,
 323                "Only local and non-local xq transfers are supported");
 324
 325     const InteractionLocality iloc = gpuAtomToInteractionLocality(atomLocality);
 326
 327     int                       adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 328
 329     cu_atomdata_t            *adat    = nb->atdat;
 330     cu_plist_t               *plist   = nb->plist[iloc];
 331     cu_timers_t              *t       = nb->timers;
 332     cudaStream_t              stream  = nb->stream[iloc];
 333
 334     bool                      bDoTime     = nb->bDoTime;
 335
 336     /* Don't launch the non-local H2D copy if there is no dependent
 337        work to do: neither non-local nor other (e.g. bonded) work
 338        to do that has as input the nbnxn coordaintes.
 339        Doing the same for the local kernel is more complicated, since the
 340        local part of the force array also depends on the non-local kernel.
 341        So to avoid complicating the code and to reduce the risk of bugs,
 342        we always call the local local x+q copy (and the rest of the local
 343        work in nbnxn_gpu_launch_kernel().
 344      */
 345     if ((iloc == InteractionLocality::NonLocal) && !haveGpuShortRangeWork(*nb, iloc))
 346     {
 347         plist->haveFreshList = false;
 348
 349         return;
 350     }
 351
 352     /* calculate the atom data index range based on locality */
 353     if (atomLocality == AtomLocality::Local)
 354     {
 355         adat_begin  = 0;
 356         adat_len    = adat->natoms_local;
 357     }
 358     else
 359     {
 360         adat_begin  = adat->natoms_local;
 361         adat_len    = adat->natoms - adat->natoms_local;
 362     }
 363
 364     /* HtoD x, q */
 365     /* beginning of timed HtoD section */
 366     if (bDoTime)
 367     {
 368         t->xf[atomLocality].nb_h2d.openTimingRegion(stream);
 369     }
 370
 371     cu_copy_H2D_async(adat->xq + adat_begin, static_cast<const void *>(nbatom->x().data() + adat_begin * 4),
 372                       adat_len * sizeof(*adat->xq), stream);
 373
 374     if (bDoTime)
 375     {
 376         t->xf[atomLocality].nb_h2d.closeTimingRegion(stream);
 377     }
 378
 379     /* When we get here all misc operations issued in the local stream as well as
 380        the local xq H2D are done,
 381        so we record that in the local stream and wait for it in the nonlocal one.
 382        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 383        compute on interactions between local and nonlocal atoms.
 384      */
 385     nbnxnInsertNonlocalGpuDependency(nb, iloc);
 386 }
 387
 388 /*! As we execute nonbonded workload in separate streams, before launching
 389    the kernel we need to make sure that he following operations have completed:
 390    - atomdata allocation and related H2D transfers (every nstlist step);
 391    - pair list H2D transfer (every nstlist step);
 392    - shift vector H2D transfer (every nstlist step);
 393    - force (+shift force and energy) output clearing (every step).
 394
 395    These operations are issued in the local stream at the beginning of the step
 396    and therefore always complete before the local kernel launch. The non-local
 397    kernel is launched after the local on the same device/context hence it is
 398    inherently scheduled after the operations in the local stream (including the
 399    above "misc_ops") on pre-GK110 devices with single hardware queue, but on later
 400    devices with multiple hardware queues the dependency needs to be enforced.
 401    We use the misc_ops_and_local_H2D_done event to record the point where
 402    the local x+q H2D (and all preceding) tasks are complete and synchronize
 403    with this event in the non-local stream before launching the non-bonded kernel.
 404  */
 405 void gpu_launch_kernel(gmx_nbnxn_cuda_t          *nb,
 406                        const gmx::ForceFlags     &forceFlags,
 407                        const InteractionLocality  iloc)
 408 {
 409     cu_atomdata_t       *adat    = nb->atdat;
 410     cu_nbparam_t        *nbp     = nb->nbparam;
 411     cu_plist_t          *plist   = nb->plist[iloc];
 412     cu_timers_t         *t       = nb->timers;
 413     cudaStream_t         stream  = nb->stream[iloc];
 414
 415     bool                 bDoTime     = nb->bDoTime;
 416
 417     /* Don't launch the non-local kernel if there is no work to do.
 418        Doing the same for the local kernel is more complicated, since the
 419        local part of the force array also depends on the non-local kernel.
 420        So to avoid complicating the code and to reduce the risk of bugs,
 421        we always call the local kernel, and later (not in
 422        this function) the stream wait, local f copyback and the f buffer
 423        clearing. All these operations, except for the local interaction kernel,
 424        are needed for the non-local interactions. The skip of the local kernel
 425        call is taken care of later in this function. */
 426     if (canSkipNonbondedWork(*nb, iloc))
 427     {
 428         plist->haveFreshList = false;
 429
 430         return;
 431     }
 432
 433     if (nbp->useDynamicPruning && plist->haveFreshList)
 434     {
 435         /* Prunes for rlistOuter and rlistInner, sets plist->haveFreshList=false
 436            (TODO: ATM that's the way the timing accounting can distinguish between
 437            separate prune kernel and combined force+prune, maybe we need a better way?).
 438          */
 439         gpu_launch_kernel_pruneonly(nb, iloc, 1);
 440     }
 441
 442     if (plist->nsci == 0)
 443     {
 444         /* Don't launch an empty local kernel (not allowed with CUDA) */
 445         return;
 446     }
 447
 448     /* beginning of timed nonbonded calculation section */
 449     if (bDoTime)
 450     {
 451         t->interaction[iloc].nb_k.openTimingRegion(stream);
 452     }
 453
 454     /* Kernel launch config:
 455      * - The thread block dimensions match the size of i-clusters, j-clusters,
 456      *   and j-cluster concurrency, in x, y, and z, respectively.
 457      * - The 1D block-grid contains as many blocks as super-clusters.
 458      */
 459     int num_threads_z = 1;
 460     if (nb->dev_info->prop.major == 3 && nb->dev_info->prop.minor == 7)
 461     {
 462         num_threads_z = 2;
 463     }
 464     int nblock    = calc_nb_kernel_nblock(plist->nsci, nb->dev_info);
 465
 466
 467     KernelLaunchConfig config;
 468     config.blockSize[0]     = c_clSize;
 469     config.blockSize[1]     = c_clSize;
 470     config.blockSize[2]     = num_threads_z;
 471     config.gridSize[0]      = nblock;
 472     config.sharedMemorySize = calc_shmem_required_nonbonded(num_threads_z, nb->dev_info, nbp);
 473     config.stream           = stream;
 474
 475     if (debug)
 476     {
 477         fprintf(debug, "Non-bonded GPU launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 478                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 479                 "\tShMem: %zu\n",
 480                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 481                 config.gridSize[0], config.gridSize[1], plist->nsci*c_numClPerSupercl,
 482                 c_numClPerSupercl, plist->na_c,
 483                 config.sharedMemorySize);
 484     }
 485
 486     auto       *timingEvent = bDoTime ? t->interaction[iloc].nb_k.fetchNextEvent() : nullptr;
 487     const auto  kernel      = select_nbnxn_kernel(nbp->eeltype,
 488                                                   nbp->vdwtype,
 489                                                   forceFlags.computeEnergy,
 490                                                   (plist->haveFreshList && !nb->timers->interaction[iloc].didPrune),
 491                                                   nb->dev_info);
 492     const auto kernelArgs  = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &forceFlags.computeVirial);
 493     launchGpuKernel(kernel, config, timingEvent, "k_calc_nb", kernelArgs);
 494
 495     if (bDoTime)
 496     {
 497         t->interaction[iloc].nb_k.closeTimingRegion(stream);
 498     }
 499
 500     if (GMX_NATIVE_WINDOWS)
 501     {
 502         /* Windows: force flushing WDDM queue */
 503         cudaStreamQuery(stream);
 504     }
 505 }
 506
 507 /*! Calculates the amount of shared memory required by the CUDA kernel in use. */
 508 static inline int calc_shmem_required_prune(const int num_threads_z)
 509 {
 510     int shmem;
 511
 512     /* i-atom x in shared memory */
 513     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 514     /* cj in shared memory, for each warp separately */
 515     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 516
 517     return shmem;
 518 }
 519
 520 void gpu_launch_kernel_pruneonly(gmx_nbnxn_cuda_t          *nb,
 521                                  const InteractionLocality  iloc,
 522                                  const int                  numParts)
 523 {
 524     cu_atomdata_t       *adat    = nb->atdat;
 525     cu_nbparam_t        *nbp     = nb->nbparam;
 526     cu_plist_t          *plist   = nb->plist[iloc];
 527     cu_timers_t         *t       = nb->timers;
 528     cudaStream_t         stream  = nb->stream[iloc];
 529
 530     bool                 bDoTime = nb->bDoTime;
 531
 532     if (plist->haveFreshList)
 533     {
 534         GMX_ASSERT(numParts == 1, "With first pruning we expect 1 part");
 535
 536         /* Set rollingPruningNumParts to signal that it is not set */
 537         plist->rollingPruningNumParts = 0;
 538         plist->rollingPruningPart     = 0;
 539     }
 540     else
 541     {
 542         if (plist->rollingPruningNumParts == 0)
 543         {
 544             plist->rollingPruningNumParts = numParts;
 545         }
 546         else
 547         {
 548             GMX_ASSERT(numParts == plist->rollingPruningNumParts, "It is not allowed to change numParts in between list generation steps");
 549         }
 550     }
 551
 552     /* Use a local variable for part and update in plist, so we can return here
 553      * without duplicating the part increment code.
 554      */
 555     int part = plist->rollingPruningPart;
 556
 557     plist->rollingPruningPart++;
 558     if (plist->rollingPruningPart >= plist->rollingPruningNumParts)
 559     {
 560         plist->rollingPruningPart = 0;
 561     }
 562
 563     /* Compute the number of list entries to prune in this pass */
 564     int numSciInPart = (plist->nsci - part)/numParts;
 565
 566     /* Don't launch the kernel if there is no work to do (not allowed with CUDA) */
 567     if (numSciInPart <= 0)
 568     {
 569         plist->haveFreshList = false;
 570
 571         return;
 572     }
 573
 574     GpuRegionTimer *timer = nullptr;
 575     if (bDoTime)
 576     {
 577         timer = &(plist->haveFreshList ? t->interaction[iloc].prune_k : t->interaction[iloc].rollingPrune_k);
 578     }
 579
 580     /* beginning of timed prune calculation section */
 581     if (bDoTime)
 582     {
 583         timer->openTimingRegion(stream);
 584     }
 585
 586     /* Kernel launch config:
 587      * - The thread block dimensions match the size of i-clusters, j-clusters,
 588      *   and j-cluster concurrency, in x, y, and z, respectively.
 589      * - The 1D block-grid contains as many blocks as super-clusters.
 590      */
 591     int                num_threads_z  = c_cudaPruneKernelJ4Concurrency;
 592     int                nblock         = calc_nb_kernel_nblock(numSciInPart, nb->dev_info);
 593     KernelLaunchConfig config;
 594     config.blockSize[0]     = c_clSize;
 595     config.blockSize[1]     = c_clSize;
 596     config.blockSize[2]     = num_threads_z;
 597     config.gridSize[0]      = nblock;
 598     config.sharedMemorySize = calc_shmem_required_prune(num_threads_z);
 599     config.stream           = stream;
 600
 601     if (debug)
 602     {
 603         fprintf(debug, "Pruning GPU kernel launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 604                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 605                 "\tShMem: %zu\n",
 606                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 607                 config.gridSize[0], config.gridSize[1], numSciInPart*c_numClPerSupercl,
 608                 c_numClPerSupercl, plist->na_c,
 609                 config.sharedMemorySize);
 610     }
 611
 612     auto          *timingEvent  = bDoTime ? timer->fetchNextEvent() : nullptr;
 613     constexpr char kernelName[] = "k_pruneonly";
 614     const auto     kernel       = plist->haveFreshList ? nbnxn_kernel_prune_cuda<true> : nbnxn_kernel_prune_cuda<false>;
 615     const auto     kernelArgs   = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &numParts, &part);
 616     launchGpuKernel(kernel, config, timingEvent, kernelName, kernelArgs);
 617
 618     /* TODO: consider a more elegant way to track which kernel has been called
 619        (combined or separate 1st pass prune, rolling prune). */
 620     if (plist->haveFreshList)
 621     {
 622         plist->haveFreshList                   = false;
 623         /* Mark that pruning has been done */
 624         nb->timers->interaction[iloc].didPrune = true;
 625     }
 626     else
 627     {
 628         /* Mark that rolling pruning has been done */
 629         nb->timers->interaction[iloc].didRollingPrune = true;
 630     }
 631
 632     if (bDoTime)
 633     {
 634         timer->closeTimingRegion(stream);
 635     }
 636
 637     if (GMX_NATIVE_WINDOWS)
 638     {
 639         /* Windows: force flushing WDDM queue */
 640         cudaStreamQuery(stream);
 641     }
 642 }
 643
 644 void gpu_launch_cpyback(gmx_nbnxn_cuda_t       *nb,
 645                         nbnxn_atomdata_t       *nbatom,
 646                         const gmx::ForceFlags  &forceFlags,
 647                         const AtomLocality      atomLocality,
 648                         const bool              copyBackNbForce)
 649 {
 650     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 651
 652     cudaError_t stat;
 653     int         adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 654
 655     /* determine interaction locality from atom locality */
 656     const InteractionLocality iloc = gpuAtomToInteractionLocality(atomLocality);
 657
 658     /* extract the data */
 659     cu_atomdata_t   *adat    = nb->atdat;
 660     cu_timers_t     *t       = nb->timers;
 661     bool             bDoTime = nb->bDoTime;
 662     cudaStream_t     stream  = nb->stream[iloc];
 663
 664     /* don't launch non-local copy-back if there was no non-local work to do */
 665     if ((iloc == InteractionLocality::NonLocal) && !haveGpuShortRangeWork(*nb, iloc))
 666     {
 667         return;
 668     }
 669
 670     getGpuAtomRange(adat, atomLocality, &adat_begin, &adat_len);
 671
 672     /* beginning of timed D2H section */
 673     if (bDoTime)
 674     {
 675         t->xf[atomLocality].nb_d2h.openTimingRegion(stream);
 676     }
 677
 678     /* With DD the local D2H transfer can only start after the non-local
 679        kernel has finished. */
 680     if (iloc == InteractionLocality::Local && nb->bUseTwoStreams)
 681     {
 682         stat = cudaStreamWaitEvent(stream, nb->nonlocal_done, 0);
 683         CU_RET_ERR(stat, "cudaStreamWaitEvent on nonlocal_done failed");
 684     }
 685
 686     /* DtoH f */
 687     if (copyBackNbForce)
 688     {
 689         cu_copy_D2H_async(nbatom->out[0].f.data() + adat_begin * 3, adat->f + adat_begin,
 690                           (adat_len)*sizeof(*adat->f), stream);
 691     }
 692
 693     /* After the non-local D2H is launched the nonlocal_done event can be
 694        recorded which signals that the local D2H can proceed. This event is not
 695        placed after the non-local kernel because we want the non-local data
 696        back first. */
 697     if (iloc == InteractionLocality::NonLocal)
 698     {
 699         stat = cudaEventRecord(nb->nonlocal_done, stream);
 700         CU_RET_ERR(stat, "cudaEventRecord on nonlocal_done failed");
 701     }
 702
 703     /* only transfer energies in the local stream */
 704     if (iloc == InteractionLocality::Local)
 705     {
 706         /* DtoH fshift when virial is needed */
 707         if (forceFlags.computeVirial)
 708         {
 709             cu_copy_D2H_async(nb->nbst.fshift, adat->fshift,
 710                               SHIFTS * sizeof(*nb->nbst.fshift), stream);
 711         }
 712
 713         /* DtoH energies */
 714         if (forceFlags.computeEnergy)
 715         {
 716             cu_copy_D2H_async(nb->nbst.e_lj, adat->e_lj,
 717                               sizeof(*nb->nbst.e_lj), stream);
 718             cu_copy_D2H_async(nb->nbst.e_el, adat->e_el,
 719                               sizeof(*nb->nbst.e_el), stream);
 720         }
 721     }
 722
 723     if (bDoTime)
 724     {
 725         t->xf[atomLocality].nb_d2h.closeTimingRegion(stream);
 726     }
 727 }
 728
 729 void cuda_set_cacheconfig()
 730 {
 731     cudaError_t stat;
 732
 733     for (int i = 0; i < eelCuNR; i++)
 734     {
 735         for (int j = 0; j < evdwCuNR; j++)
 736         {
 737             /* Default kernel 32/32 kB Shared/L1 */
 738             cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 739             cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 740             cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 741             stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 742             CU_RET_ERR(stat, "cudaFuncSetCacheConfig failed");
 743         }
 744     }
 745 }
 746
 747 /* X buffer operations on GPU: copies coordinates to the GPU in rvec format. */
 748 void nbnxn_gpu_copy_x_to_gpu(const Nbnxm::Grid               &grid,
 749                              gmx_nbnxn_gpu_t                 *nb,
 750                              const Nbnxm::AtomLocality        locality,
 751                              const rvec                      *coordinatesHost)
 752 {
 753     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 754
 755     bool                       bDoTime = nb->bDoTime;
 756
 757     Nbnxm::InteractionLocality interactionLoc            = gpuAtomToInteractionLocality(locality);
 758     int                        numCopyAtoms              = grid.srcAtomEnd() - grid.srcAtomBegin();
 759     int                        copyAtomStart             = grid.srcAtomBegin();
 760
 761     cudaStream_t               stream  = nb->stream[interactionLoc];
 762
 763     // empty domain avoid launching zero-byte copy
 764     if (numCopyAtoms == 0)
 765     {
 766         return;
 767     }
 768     GMX_ASSERT(coordinatesHost,  "Need a valid host pointer");
 769
 770     if (bDoTime)
 771     {
 772         nb->timers->xf[locality].nb_h2d.openTimingRegion(stream);
 773     }
 774
 775     rvec       *devicePtrDest = reinterpret_cast<rvec *> (nb->xrvec[copyAtomStart]);
 776     const rvec *devicePtrSrc  = reinterpret_cast<const rvec *> (coordinatesHost[copyAtomStart]);
 777     copyToDeviceBuffer(&devicePtrDest, devicePtrSrc, 0, numCopyAtoms,
 778                        stream, GpuApiCallBehavior::Async, nullptr);
 779
 780     if (bDoTime)
 781     {
 782         nb->timers->xf[locality].nb_h2d.closeTimingRegion(stream);
 783     }
 784 }
 785
 786 DeviceBuffer<float> nbnxn_gpu_get_x_gpu(gmx_nbnxn_gpu_t *nb)
 787 {
 788     return reinterpret_cast< DeviceBuffer<float> >(nb->xrvec);
 789 }
 790
 791 /* X buffer operations on GPU: performs conversion from rvec to nb format. */
 792 void nbnxn_gpu_x_to_nbat_x(const Nbnxm::Grid               &grid,
 793                            bool                             setFillerCoords,
 794                            gmx_nbnxn_gpu_t                 *nb,
 795                            DeviceBuffer<float>              coordinatesDevice,
 796                            const Nbnxm::AtomLocality        locality,
 797                            int                              gridId,
 798                            int                              numColumnsMax)
 799 {
 800     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 801
 802     cu_atomdata_t             *adat    = nb->atdat;
 803
 804     const int                  numColumns                = grid.numColumns();
 805     const int                  cellOffset                = grid.cellOffset();
 806     const int                  numAtomsPerCell           = grid.numAtomsPerCell();
 807     Nbnxm::InteractionLocality interactionLoc            = gpuAtomToInteractionLocality(locality);
 808
 809     cudaStream_t               stream  = nb->stream[interactionLoc];
 810
 811     int numAtoms = grid.srcAtomEnd() - grid.srcAtomBegin();
 812     // avoid empty kernel launch, skip to inserting stream dependency
 813     if (numAtoms != 0)
 814     {
 815         // TODO: This will only work with CUDA
 816         GMX_ASSERT(coordinatesDevice, "Need a valid device pointer");
 817
 818         KernelLaunchConfig config;
 819         config.blockSize[0]     = c_bufOpsThreadsPerBlock;
 820         config.blockSize[1]     = 1;
 821         config.blockSize[2]     = 1;
 822         config.gridSize[0]      = (grid.numCellsColumnMax()*numAtomsPerCell + c_bufOpsThreadsPerBlock - 1)/c_bufOpsThreadsPerBlock;
 823         config.gridSize[1]      = numColumns;
 824         config.gridSize[2]      = 1;
 825         GMX_ASSERT(config.gridSize[0] > 0, "Can not have empty grid, early return above avoids this");
 826         config.sharedMemorySize = 0;
 827         config.stream           = stream;
 828
 829         auto       kernelFn       = nbnxn_gpu_x_to_nbat_x_kernel;
 830         float     *xqPtr          = &(adat->xq->x);
 831         const int *d_atomIndices  = nb->atomIndices;
 832         const int *d_cxy_na       = &nb->cxy_na[numColumnsMax*gridId];
 833         const int *d_cxy_ind      = &nb->cxy_ind[numColumnsMax*gridId];
 834         const auto kernelArgs     = prepareGpuKernelArguments(kernelFn, config,
 835                                                               &numColumns,
 836                                                               &xqPtr,
 837                                                               &setFillerCoords,
 838                                                               &coordinatesDevice,
 839                                                               &d_atomIndices,
 840                                                               &d_cxy_na,
 841                                                               &d_cxy_ind,
 842                                                               &cellOffset,
 843                                                               &numAtomsPerCell);
 844         launchGpuKernel(kernelFn, config, nullptr, "XbufferOps", kernelArgs);
 845     }
 846
 847     // TODO: note that this is not necessary when there are no local atoms, that is:
 848     // (numAtoms == 0 && interactionLoc == InteractionLocality::Local)
 849     // but for now we avoid that optimization
 850     nbnxnInsertNonlocalGpuDependency(nb, interactionLoc);
 851 }
 852
 853 /* F buffer operations on GPU: performs force summations and conversion from nb to rvec format. */
 854 void nbnxn_gpu_add_nbat_f_to_f(const AtomLocality               atomLocality,
 855                                DeviceBuffer<float>              totalForcesDevice,
 856                                gmx_nbnxn_gpu_t                 *nb,
 857                                void                            *pmeForcesDevice,
 858                                GpuEventSynchronizer            *pmeForcesReady,
 859                                int                              atomStart,
 860                                int                              numAtoms,
 861                                bool                             useGpuFPmeReduction,
 862                                bool                             accumulateForce)
 863 {
 864     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 865     GMX_ASSERT(numAtoms != 0, "Cannot call function with no atoms");
 866     GMX_ASSERT(totalForcesDevice, "Need a valid totalForcesDevice pointer");
 867
 868     const InteractionLocality iLocality     = gpuAtomToInteractionLocality(atomLocality);
 869     cudaStream_t              stream        = nb->stream[iLocality];
 870     cu_atomdata_t            *adat          = nb->atdat;
 871
 872     if (useGpuFPmeReduction)
 873     {
 874         //Stream must wait for PME force completion
 875         pmeForcesReady->enqueueWaitEvent(stream);
 876     }
 877
 878     /* launch kernel */
 879
 880     KernelLaunchConfig config;
 881     config.blockSize[0]     = c_bufOpsThreadsPerBlock;
 882     config.blockSize[1]     = 1;
 883     config.blockSize[2]     = 1;
 884     config.gridSize[0]      = ((numAtoms+1)+c_bufOpsThreadsPerBlock-1)/c_bufOpsThreadsPerBlock;
 885     config.gridSize[1]      = 1;
 886     config.gridSize[2]      = 1;
 887     config.sharedMemorySize = 0;
 888     config.stream           = stream;
 889
 890     auto  kernelFn = accumulateForce ?
 891         nbnxn_gpu_add_nbat_f_to_f_kernel<true, false> :
 892         nbnxn_gpu_add_nbat_f_to_f_kernel<false, false>;
 893
 894     if (useGpuFPmeReduction)
 895     {
 896         GMX_ASSERT(pmeForcesDevice, "Need a valid pmeForcesDevice pointer");
 897         kernelFn = accumulateForce ?
 898             nbnxn_gpu_add_nbat_f_to_f_kernel<true, true> :
 899             nbnxn_gpu_add_nbat_f_to_f_kernel<false, true>;
 900     }
 901
 902     const float3     *d_fNB    = adat->f;
 903     const float3     *d_fPme   = (float3*) pmeForcesDevice;
 904     float3           *d_fTotal = (float3*) totalForcesDevice;
 905     const int        *d_cell   = nb->cell;
 906
 907     const auto        kernelArgs   = prepareGpuKernelArguments(kernelFn, config,
 908                                                                &d_fNB,
 909                                                                &d_fPme,
 910                                                                &d_fTotal,
 911                                                                &d_cell,
 912                                                                &atomStart,
 913                                                                &numAtoms);
 914
 915     launchGpuKernel(kernelFn, config, nullptr, "FbufferOps", kernelArgs);
 916
 917 }
 918
 919 DeviceBuffer<float> nbnxn_gpu_get_f_gpu(gmx_nbnxn_gpu_t *nb)
 920 {
 921     return reinterpret_cast< DeviceBuffer<float> >(nb->frvec);
 922 }
 923
 924 void nbnxn_launch_copy_f_to_gpu(const AtomLocality               atomLocality,
 925                                 const Nbnxm::GridSet            &gridSet,
 926                                 gmx_nbnxn_gpu_t                 *nb,
 927                                 rvec                            *f)
 928 {
 929     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 930
 931     const InteractionLocality iLocality = gpuAtomToInteractionLocality(atomLocality);
 932     cudaStream_t              stream    = nb->stream[iLocality];
 933
 934     bool                      bDoTime = nb->bDoTime;
 935     cu_timers_t              *t       = nb->timers;
 936
 937     int                       atomStart = 0, numCopyAtoms = 0;
 938
 939     nbnxn_get_atom_range(atomLocality, gridSet, &atomStart, &numCopyAtoms);
 940
 941     // Avoiding launching copy with no work
 942     if (numCopyAtoms == 0)
 943     {
 944         return;
 945     }
 946     GMX_ASSERT(f, "Need a valid f pointer");
 947
 948     if (bDoTime)
 949     {
 950         t->xf[atomLocality].nb_h2d.openTimingRegion(stream);
 951     }
 952
 953     rvec       *ptrDest  = reinterpret_cast<rvec *> (nb->frvec[atomStart]);
 954     rvec       *ptrSrc   = reinterpret_cast<rvec *> (f[atomStart]);
 955     //copyToDeviceBuffer(&ptrDest, ptrSrc, 0, numCopyAtoms,
 956     //                   stream, GpuApiCallBehavior::Async, nullptr);
 957     //TODO use above API call rather than direct memcpy when force has been implemented in a hostvector
 958     cudaMemcpyAsync(ptrDest, ptrSrc, numCopyAtoms*sizeof(rvec), cudaMemcpyHostToDevice,
 959                     stream);
 960
 961     if (bDoTime)
 962     {
 963         t->xf[atomLocality].nb_h2d.closeTimingRegion(stream);
 964     }
 965
 966     return;
 967 }
 968
 969 void nbnxn_launch_copy_f_from_gpu(const AtomLocality               atomLocality,
 970                                   const Nbnxm::GridSet            &gridSet,
 971                                   gmx_nbnxn_gpu_t                 *nb,
 972                                   rvec                            *f)
 973 {
 974     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 975
 976     const InteractionLocality iLocality = gpuAtomToInteractionLocality(atomLocality);
 977     cudaStream_t              stream    = nb->stream[iLocality];
 978
 979     bool                      bDoTime = nb->bDoTime;
 980     cu_timers_t              *t       = nb->timers;
 981     int                       atomStart, numCopyAtoms;
 982
 983     nbnxn_get_atom_range(atomLocality, gridSet, &atomStart, &numCopyAtoms);
 984
 985     // Avoiding launching copy with no work
 986     if (numCopyAtoms == 0)
 987     {
 988         return;
 989     }
 990     GMX_ASSERT(f, "Need a valid f pointer");
 991
 992     if (bDoTime)
 993     {
 994         t->xf[atomLocality].nb_d2h.openTimingRegion(stream);
 995     }
 996
 997     GMX_ASSERT(nb->frvec,  "Need a valid nb->frvec pointer");
 998     rvec       *ptrDest = reinterpret_cast<rvec *> (f[atomStart]);
 999     rvec       *ptrSrc  = reinterpret_cast<rvec *> (nb->frvec[atomStart]);
1000     //copyFromDeviceBuffer(ptrDest, &ptrSrc, 0, numCopyAtoms,
1001     //                   stream, GpuApiCallBehavior::Async, nullptr);
1002     //TODO use above API call rather than direct memcpy when force has been implemented in a hostvector
1003     cudaMemcpyAsync(ptrDest, ptrSrc, numCopyAtoms*sizeof(rvec), cudaMemcpyDeviceToHost,
1004                     stream);
1005
1006     if (bDoTime)
1007     {
1008         t->xf[atomLocality].nb_d2h.closeTimingRegion(stream);
1009     }
1010
1011     return;
1012 }
1013
1014 void nbnxn_wait_for_gpu_force_reduction(const AtomLocality      gmx_unused atomLocality,
1015                                         gmx_nbnxn_gpu_t                   *nb)
1016 {
1017     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
1018
1019     const InteractionLocality iLocality = gpuAtomToInteractionLocality(atomLocality);
1020
1021     cudaStream_t              stream    = nb->stream[iLocality];
1022
1023     cudaStreamSynchronize(stream);
1024
1025 }
1026
1027 } // namespace Nbnxm