src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.cu

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016,2017, 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_data_mgmt.h
  37  *
  38  *  \author Szilard Pall <pall.szilard@gmail.com>
  39  */
  40 #include "gmxpre.h"
  41
  42 #include <assert.h>
  43 #include <stdarg.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46
  47 #include "gromacs/gpu_utils/cudautils.cuh"
  48 #include "gromacs/gpu_utils/gpu_utils.h"
  49 #include "gromacs/gpu_utils/pmalloc_cuda.h"
  50 #include "gromacs/hardware/gpu_hw_info.h"
  51 #include "gromacs/math/vectypes.h"
  52 #include "gromacs/mdlib/force_flags.h"
  53 #include "gromacs/mdlib/nb_verlet.h"
  54 #include "gromacs/mdlib/nbnxn_consts.h"
  55 #include "gromacs/mdlib/nbnxn_gpu_common.h"
  56 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  57 #include "gromacs/mdtypes/interaction_const.h"
  58 #include "gromacs/mdtypes/md_enums.h"
  59 #include "gromacs/pbcutil/ishift.h"
  60 #include "gromacs/timing/gpu_timing.h"
  61 #include "gromacs/utility/basedefinitions.h"
  62 #include "gromacs/utility/cstringutil.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/real.h"
  65 #include "gromacs/utility/smalloc.h"
  66
  67 #include "nbnxn_cuda.h"
  68 #include "nbnxn_cuda_types.h"
  69
  70 static bool bUseCudaEventBlockingSync = false; /* makes the CPU thread block */
  71
  72 /* This is a heuristically determined parameter for the Fermi, Kepler
  73  * and Maxwell architectures for the minimum size of ci lists by multiplying
  74  * this constant with the # of multiprocessors on the current device.
  75  * Since the maximum number of blocks per multiprocessor is 16, the ideal
  76  * count for small systems is 32 or 48 blocks per multiprocessor. Because
  77  * there is a bit of fluctuations in the generated block counts, we use
  78  * a target of 44 instead of the ideal value of 48.
  79  */
  80 static unsigned int gpu_min_ci_balanced_factor = 44;
  81
  82 /* Fw. decl. */
  83 static void nbnxn_cuda_clear_e_fshift(gmx_nbnxn_cuda_t *nb);
  84
  85 /* Fw. decl, */
  86 static void nbnxn_cuda_free_nbparam_table(cu_nbparam_t            *nbparam,
  87                                           const gmx_device_info_t *dev_info);
  88
  89 /*! \brief Return whether combination rules are used.
  90  *
  91  * \param[in]   pointer to nonbonded paramter struct
  92  * \return      true if combination rules are used in this run, false otherwise
  93  */
  94 static inline bool useLjCombRule(const cu_nbparam_t  *nbparam)
  95 {
  96     return (nbparam->vdwtype == evdwCuCUTCOMBGEOM ||
  97             nbparam->vdwtype == evdwCuCUTCOMBLB);
  98 }
  99
 100 /*! \brief Initialized the Ewald Coulomb correction GPU table.
 101
 102     Tabulates the Ewald Coulomb force and initializes the size/scale
 103     and the table GPU array. If called with an already allocated table,
 104     it just re-uploads the table.
 105  */
 106 static void init_ewald_coulomb_force_table(const interaction_const_t *ic,
 107                                            cu_nbparam_t              *nbp,
 108                                            const gmx_device_info_t   *dev_info)
 109 {
 110     if (nbp->coulomb_tab != NULL)
 111     {
 112         nbnxn_cuda_free_nbparam_table(nbp, dev_info);
 113     }
 114
 115     nbp->coulomb_tab_scale = ic->tabq_scale;
 116     initParamLookupTable(nbp->coulomb_tab, nbp->coulomb_tab_texobj,
 117                          &nbnxn_cuda_get_coulomb_tab_texref(),
 118                          ic->tabq_coul_F, ic->tabq_size, dev_info);
 119 }
 120
 121
 122 /*! Initializes the atomdata structure first time, it only gets filled at
 123     pair-search. */
 124 static void init_atomdata_first(cu_atomdata_t *ad, int ntypes)
 125 {
 126     cudaError_t stat;
 127
 128     ad->ntypes  = ntypes;
 129     stat        = cudaMalloc((void**)&ad->shift_vec, SHIFTS*sizeof(*ad->shift_vec));
 130     CU_RET_ERR(stat, "cudaMalloc failed on ad->shift_vec");
 131     ad->bShiftVecUploaded = false;
 132
 133     stat = cudaMalloc((void**)&ad->fshift, SHIFTS*sizeof(*ad->fshift));
 134     CU_RET_ERR(stat, "cudaMalloc failed on ad->fshift");
 135
 136     stat = cudaMalloc((void**)&ad->e_lj, sizeof(*ad->e_lj));
 137     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_lj");
 138     stat = cudaMalloc((void**)&ad->e_el, sizeof(*ad->e_el));
 139     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_el");
 140
 141     /* initialize to NULL poiters to data that is not allocated here and will
 142        need reallocation in nbnxn_cuda_init_atomdata */
 143     ad->xq = NULL;
 144     ad->f  = NULL;
 145
 146     /* size -1 indicates that the respective array hasn't been initialized yet */
 147     ad->natoms = -1;
 148     ad->nalloc = -1;
 149 }
 150
 151 /*! Selects the Ewald kernel type, analytical on SM 3.0 and later, tabulated on
 152     earlier GPUs, single or twin cut-off. */
 153 static int pick_ewald_kernel_type(bool                     bTwinCut,
 154                                   const gmx_device_info_t *dev_info)
 155 {
 156     bool bUseAnalyticalEwald, bForceAnalyticalEwald, bForceTabulatedEwald;
 157     int  kernel_type;
 158
 159     /* Benchmarking/development environment variables to force the use of
 160        analytical or tabulated Ewald kernel. */
 161     bForceAnalyticalEwald = (getenv("GMX_CUDA_NB_ANA_EWALD") != NULL);
 162     bForceTabulatedEwald  = (getenv("GMX_CUDA_NB_TAB_EWALD") != NULL);
 163
 164     if (bForceAnalyticalEwald && bForceTabulatedEwald)
 165     {
 166         gmx_incons("Both analytical and tabulated Ewald CUDA non-bonded kernels "
 167                    "requested through environment variables.");
 168     }
 169
 170     /* By default, on SM 3.0 and later use analytical Ewald, on earlier tabulated. */
 171     if ((dev_info->prop.major >= 3 || bForceAnalyticalEwald) && !bForceTabulatedEwald)
 172     {
 173         bUseAnalyticalEwald = true;
 174
 175         if (debug)
 176         {
 177             fprintf(debug, "Using analytical Ewald CUDA kernels\n");
 178         }
 179     }
 180     else
 181     {
 182         bUseAnalyticalEwald = false;
 183
 184         if (debug)
 185         {
 186             fprintf(debug, "Using tabulated Ewald CUDA kernels\n");
 187         }
 188     }
 189
 190     /* Use twin cut-off kernels if requested by bTwinCut or the env. var.
 191        forces it (use it for debugging/benchmarking only). */
 192     if (!bTwinCut && (getenv("GMX_CUDA_NB_EWALD_TWINCUT") == NULL))
 193     {
 194         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA : eelCuEWALD_TAB;
 195     }
 196     else
 197     {
 198         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA_TWIN : eelCuEWALD_TAB_TWIN;
 199     }
 200
 201     return kernel_type;
 202 }
 203
 204 /*! Copies all parameters related to the cut-off from ic to nbp */
 205 static void set_cutoff_parameters(cu_nbparam_t              *nbp,
 206                                   const interaction_const_t *ic,
 207                                   const NbnxnListParameters *listParams)
 208 {
 209     nbp->ewald_beta        = ic->ewaldcoeff_q;
 210     nbp->sh_ewald          = ic->sh_ewald;
 211     nbp->epsfac            = ic->epsfac;
 212     nbp->two_k_rf          = 2.0 * ic->k_rf;
 213     nbp->c_rf              = ic->c_rf;
 214     nbp->rvdw_sq           = ic->rvdw * ic->rvdw;
 215     nbp->rcoulomb_sq       = ic->rcoulomb * ic->rcoulomb;
 216     nbp->rlistOuter_sq     = listParams->rlistOuter * listParams->rlistOuter;
 217     nbp->rlistInner_sq     = listParams->rlistInner * listParams->rlistInner;
 218     nbp->useDynamicPruning = listParams->useDynamicPruning;
 219
 220     nbp->sh_lj_ewald       = ic->sh_lj_ewald;
 221     nbp->ewaldcoeff_lj     = ic->ewaldcoeff_lj;
 222
 223     nbp->rvdw_switch       = ic->rvdw_switch;
 224     nbp->dispersion_shift  = ic->dispersion_shift;
 225     nbp->repulsion_shift   = ic->repulsion_shift;
 226     nbp->vdw_switch        = ic->vdw_switch;
 227 }
 228
 229 /*! Initializes the nonbonded parameter data structure. */
 230 static void init_nbparam(cu_nbparam_t              *nbp,
 231                          const interaction_const_t *ic,
 232                          const NbnxnListParameters *listParams,
 233                          const nbnxn_atomdata_t    *nbat,
 234                          const gmx_device_info_t   *dev_info)
 235 {
 236     int         ntypes;
 237
 238     ntypes  = nbat->ntype;
 239
 240     set_cutoff_parameters(nbp, ic, listParams);
 241
 242     /* The kernel code supports LJ combination rules (geometric and LB) for
 243      * all kernel types, but we only generate useful combination rule kernels.
 244      * We currently only use LJ combination rule (geometric and LB) kernels
 245      * for plain cut-off LJ. On Maxwell the force only kernels speed up 15%
 246      * with PME and 20% with RF, the other kernels speed up about half as much.
 247      * For LJ force-switch the geometric rule would give 7% speed-up, but this
 248      * combination is rarely used. LJ force-switch with LB rule is more common,
 249      * but gives only 1% speed-up.
 250      */
 251     if (ic->vdwtype == evdwCUT)
 252     {
 253         switch (ic->vdw_modifier)
 254         {
 255             case eintmodNONE:
 256             case eintmodPOTSHIFT:
 257                 switch (nbat->comb_rule)
 258                 {
 259                     case ljcrNONE:
 260                         nbp->vdwtype = evdwCuCUT;
 261                         break;
 262                     case ljcrGEOM:
 263                         nbp->vdwtype = evdwCuCUTCOMBGEOM;
 264                         break;
 265                     case ljcrLB:
 266                         nbp->vdwtype = evdwCuCUTCOMBLB;
 267                         break;
 268                     default:
 269                         gmx_incons("The requested LJ combination rule is not implemented in the CUDA GPU accelerated kernels!");
 270                         break;
 271                 }
 272                 break;
 273             case eintmodFORCESWITCH:
 274                 nbp->vdwtype = evdwCuFSWITCH;
 275                 break;
 276             case eintmodPOTSWITCH:
 277                 nbp->vdwtype = evdwCuPSWITCH;
 278                 break;
 279             default:
 280                 gmx_incons("The requested VdW interaction modifier is not implemented in the CUDA GPU accelerated kernels!");
 281                 break;
 282         }
 283     }
 284     else if (ic->vdwtype == evdwPME)
 285     {
 286         if (ic->ljpme_comb_rule == ljcrGEOM)
 287         {
 288             assert(nbat->comb_rule == ljcrGEOM);
 289             nbp->vdwtype = evdwCuEWALDGEOM;
 290         }
 291         else
 292         {
 293             assert(nbat->comb_rule == ljcrLB);
 294             nbp->vdwtype = evdwCuEWALDLB;
 295         }
 296     }
 297     else
 298     {
 299         gmx_incons("The requested VdW type is not implemented in the CUDA GPU accelerated kernels!");
 300     }
 301
 302     if (ic->eeltype == eelCUT)
 303     {
 304         nbp->eeltype = eelCuCUT;
 305     }
 306     else if (EEL_RF(ic->eeltype))
 307     {
 308         nbp->eeltype = eelCuRF;
 309     }
 310     else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
 311     {
 312         /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
 313         nbp->eeltype = pick_ewald_kernel_type(false, dev_info);
 314     }
 315     else
 316     {
 317         /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
 318         gmx_incons("The requested electrostatics type is not implemented in the CUDA GPU accelerated kernels!");
 319     }
 320
 321     /* generate table for PME */
 322     nbp->coulomb_tab = NULL;
 323     if (nbp->eeltype == eelCuEWALD_TAB || nbp->eeltype == eelCuEWALD_TAB_TWIN)
 324     {
 325         init_ewald_coulomb_force_table(ic, nbp, dev_info);
 326     }
 327
 328     /* set up LJ parameter lookup table */
 329     if (!useLjCombRule(nbp))
 330     {
 331         initParamLookupTable(nbp->nbfp, nbp->nbfp_texobj,
 332                              &nbnxn_cuda_get_nbfp_texref(),
 333                              nbat->nbfp, 2*ntypes*ntypes, dev_info);
 334     }
 335
 336     /* set up LJ-PME parameter lookup table */
 337     if (ic->vdwtype == evdwPME)
 338     {
 339         initParamLookupTable(nbp->nbfp_comb, nbp->nbfp_comb_texobj,
 340                              &nbnxn_cuda_get_nbfp_comb_texref(),
 341                              nbat->nbfp_comb, 2*ntypes, dev_info);
 342     }
 343 }
 344
 345 /*! Re-generate the GPU Ewald force table, resets rlist, and update the
 346  *  electrostatic type switching to twin cut-off (or back) if needed. */
 347 void nbnxn_gpu_pme_loadbal_update_param(const nonbonded_verlet_t    *nbv,
 348                                         const interaction_const_t   *ic,
 349                                         const NbnxnListParameters   *listParams)
 350 {
 351     if (!nbv || nbv->grp[0].kernel_type != nbnxnk8x8x8_GPU)
 352     {
 353         return;
 354     }
 355     gmx_nbnxn_cuda_t *nb    = nbv->gpu_nbv;
 356     cu_nbparam_t     *nbp   = nb->nbparam;
 357
 358     set_cutoff_parameters(nbp, ic, listParams);
 359
 360     nbp->eeltype        = pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw,
 361                                                  nb->dev_info);
 362
 363     init_ewald_coulomb_force_table(ic, nb->nbparam, nb->dev_info);
 364 }
 365
 366 /*! Initializes the pair list data structure. */
 367 static void init_plist(cu_plist_t *pl)
 368 {
 369     /* initialize to NULL pointers to data that is not allocated here and will
 370        need reallocation in nbnxn_gpu_init_pairlist */
 371     pl->sci      = NULL;
 372     pl->cj4      = NULL;
 373     pl->imask    = NULL;
 374     pl->excl     = NULL;
 375
 376     /* size -1 indicates that the respective array hasn't been initialized yet */
 377     pl->na_c           = -1;
 378     pl->nsci           = -1;
 379     pl->sci_nalloc     = -1;
 380     pl->ncj4           = -1;
 381     pl->cj4_nalloc     = -1;
 382     pl->nimask         = -1;
 383     pl->imask_nalloc   = -1;
 384     pl->nexcl          = -1;
 385     pl->excl_nalloc    = -1;
 386     pl->haveFreshList  = false;
 387 }
 388
 389 /*! Initializes the timer data structure. */
 390 static void init_timers(cu_timers_t *t, bool bUseTwoStreams)
 391 {
 392     cudaError_t stat;
 393     int         eventflags = ( bUseCudaEventBlockingSync ? cudaEventBlockingSync : cudaEventDefault );
 394
 395     stat = cudaEventCreateWithFlags(&(t->start_atdat), eventflags);
 396     CU_RET_ERR(stat, "cudaEventCreate on start_atdat failed");
 397     stat = cudaEventCreateWithFlags(&(t->stop_atdat), eventflags);
 398     CU_RET_ERR(stat, "cudaEventCreate on stop_atdat failed");
 399
 400     /* The non-local counters/stream (second in the array) are needed only with DD. */
 401     for (int i = 0; i <= (bUseTwoStreams ? 1 : 0); i++)
 402     {
 403         stat = cudaEventCreateWithFlags(&(t->start_nb_k[i]), eventflags);
 404         CU_RET_ERR(stat, "cudaEventCreate on start_nb_k failed");
 405         stat = cudaEventCreateWithFlags(&(t->stop_nb_k[i]), eventflags);
 406         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_k failed");
 407
 408         stat = cudaEventCreateWithFlags(&(t->start_prune_k[i]), eventflags);
 409         CU_RET_ERR(stat, "cudaEventCreate on start_prune_k failed");
 410         stat = cudaEventCreateWithFlags(&(t->stop_prune_k[i]), eventflags);
 411         CU_RET_ERR(stat, "cudaEventCreate on stop_prune_k failed");
 412
 413         stat = cudaEventCreateWithFlags(&(t->start_rollingPrune_k[i]), eventflags);
 414         CU_RET_ERR(stat, "cudaEventCreate on start_rollingPrune_k failed");
 415         stat = cudaEventCreateWithFlags(&(t->stop_rollingPrune_k[i]), eventflags);
 416         CU_RET_ERR(stat, "cudaEventCreate on stop_rollingPrune_k failed");
 417
 418         stat = cudaEventCreateWithFlags(&(t->start_pl_h2d[i]), eventflags);
 419         CU_RET_ERR(stat, "cudaEventCreate on start_pl_h2d failed");
 420         stat = cudaEventCreateWithFlags(&(t->stop_pl_h2d[i]), eventflags);
 421         CU_RET_ERR(stat, "cudaEventCreate on stop_pl_h2d failed");
 422
 423         stat = cudaEventCreateWithFlags(&(t->start_nb_h2d[i]), eventflags);
 424         CU_RET_ERR(stat, "cudaEventCreate on start_nb_h2d failed");
 425         stat = cudaEventCreateWithFlags(&(t->stop_nb_h2d[i]), eventflags);
 426         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_h2d failed");
 427
 428         stat = cudaEventCreateWithFlags(&(t->start_nb_d2h[i]), eventflags);
 429         CU_RET_ERR(stat, "cudaEventCreate on start_nb_d2h failed");
 430         stat = cudaEventCreateWithFlags(&(t->stop_nb_d2h[i]), eventflags);
 431         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_d2h failed");
 432
 433         t->didPairlistH2D[i]  = false;
 434         t->didPrune[i]        = false;
 435         t->didRollingPrune[i] = false;
 436     }
 437 }
 438
 439 /*! Initializes the timings data structure. */
 440 static void init_timings(gmx_wallclock_gpu_t *t)
 441 {
 442     int i, j;
 443
 444     t->nb_h2d_t = 0.0;
 445     t->nb_d2h_t = 0.0;
 446     t->nb_c     = 0;
 447     t->pl_h2d_t = 0.0;
 448     t->pl_h2d_c = 0;
 449     for (i = 0; i < 2; i++)
 450     {
 451         for (j = 0; j < 2; j++)
 452         {
 453             t->ktime[i][j].t = 0.0;
 454             t->ktime[i][j].c = 0;
 455         }
 456     }
 457     t->pruneTime.c        = 0;
 458     t->pruneTime.t        = 0.0;
 459     t->dynamicPruneTime.c = 0;
 460     t->dynamicPruneTime.t = 0.0;
 461 }
 462
 463 /*! Initializes simulation constant data. */
 464 static void nbnxn_cuda_init_const(gmx_nbnxn_cuda_t               *nb,
 465                                   const interaction_const_t      *ic,
 466                                   const NbnxnListParameters      *listParams,
 467                                   const nonbonded_verlet_group_t *nbv_group)
 468 {
 469     init_atomdata_first(nb->atdat, nbv_group[0].nbat->ntype);
 470     init_nbparam(nb->nbparam, ic, listParams, nbv_group[0].nbat, nb->dev_info);
 471
 472     /* clear energy and shift force outputs */
 473     nbnxn_cuda_clear_e_fshift(nb);
 474 }
 475
 476 void nbnxn_gpu_init(gmx_nbnxn_cuda_t         **p_nb,
 477                     const gmx_device_info_t   *deviceInfo,
 478                     const interaction_const_t *ic,
 479                     const NbnxnListParameters *listParams,
 480                     nonbonded_verlet_group_t  *nbv_grp,
 481                     int                        /*rank*/,
 482                     gmx_bool                   bLocalAndNonlocal)
 483 {
 484     cudaError_t       stat;
 485     gmx_nbnxn_cuda_t *nb;
 486
 487     if (p_nb == NULL)
 488     {
 489         return;
 490     }
 491
 492     snew(nb, 1);
 493     snew(nb->atdat, 1);
 494     snew(nb->nbparam, 1);
 495     snew(nb->plist[eintLocal], 1);
 496     if (bLocalAndNonlocal)
 497     {
 498         snew(nb->plist[eintNonlocal], 1);
 499     }
 500
 501     nb->bUseTwoStreams = bLocalAndNonlocal;
 502
 503     snew(nb->timers, 1);
 504     snew(nb->timings, 1);
 505
 506     /* init nbst */
 507     pmalloc((void**)&nb->nbst.e_lj, sizeof(*nb->nbst.e_lj));
 508     pmalloc((void**)&nb->nbst.e_el, sizeof(*nb->nbst.e_el));
 509     pmalloc((void**)&nb->nbst.fshift, SHIFTS * sizeof(*nb->nbst.fshift));
 510
 511     init_plist(nb->plist[eintLocal]);
 512
 513     /* set device info, just point it to the right GPU among the detected ones */
 514     nb->dev_info = deviceInfo;
 515
 516     /* local/non-local GPU streams */
 517     stat = cudaStreamCreate(&nb->stream[eintLocal]);
 518     CU_RET_ERR(stat, "cudaStreamCreate on stream[eintLocal] failed");
 519     if (nb->bUseTwoStreams)
 520     {
 521         init_plist(nb->plist[eintNonlocal]);
 522
 523         /* Note that the device we're running on does not have to support
 524          * priorities, because we are querying the priority range which in this
 525          * case will be a single value.
 526          */
 527         int highest_priority;
 528         stat = cudaDeviceGetStreamPriorityRange(NULL, &highest_priority);
 529         CU_RET_ERR(stat, "cudaDeviceGetStreamPriorityRange failed");
 530
 531         stat = cudaStreamCreateWithPriority(&nb->stream[eintNonlocal],
 532                                             cudaStreamDefault,
 533                                             highest_priority);
 534         CU_RET_ERR(stat, "cudaStreamCreateWithPriority on stream[eintNonlocal] failed");
 535     }
 536
 537     /* init events for sychronization (timing disabled for performance reasons!) */
 538     stat = cudaEventCreateWithFlags(&nb->nonlocal_done, cudaEventDisableTiming);
 539     CU_RET_ERR(stat, "cudaEventCreate on nonlocal_done failed");
 540     stat = cudaEventCreateWithFlags(&nb->misc_ops_and_local_H2D_done, cudaEventDisableTiming);
 541     CU_RET_ERR(stat, "cudaEventCreate on misc_ops_and_local_H2D_done failed");
 542
 543     /* CUDA timing disabled as event timers don't work:
 544        - with multiple streams = domain-decomposition;
 545        - when turned off by GMX_DISABLE_CUDA_TIMING/GMX_DISABLE_GPU_TIMING.
 546      */
 547     nb->bDoTime = (!nb->bUseTwoStreams &&
 548                    (getenv("GMX_DISABLE_CUDA_TIMING") == NULL) &&
 549                    (getenv("GMX_DISABLE_GPU_TIMING") == NULL));
 550
 551     if (nb->bDoTime)
 552     {
 553         init_timers(nb->timers, nb->bUseTwoStreams);
 554         init_timings(nb->timings);
 555     }
 556
 557     /* set the kernel type for the current GPU */
 558     /* pick L1 cache configuration */
 559     nbnxn_cuda_set_cacheconfig(nb->dev_info);
 560
 561     nbnxn_cuda_init_const(nb, ic, listParams, nbv_grp);
 562
 563     *p_nb = nb;
 564
 565     if (debug)
 566     {
 567         fprintf(debug, "Initialized CUDA data structures.\n");
 568     }
 569 }
 570
 571 void nbnxn_gpu_init_pairlist(gmx_nbnxn_cuda_t       *nb,
 572                              const nbnxn_pairlist_t *h_plist,
 573                              int                     iloc)
 574 {
 575     if (canSkipWork(nb, iloc))
 576     {
 577         return;
 578     }
 579
 580     char          sbuf[STRLEN];
 581     cudaError_t   stat;
 582     bool          bDoTime    = nb->bDoTime;
 583     cudaStream_t  stream     = nb->stream[iloc];
 584     cu_plist_t   *d_plist    = nb->plist[iloc];
 585
 586     if (d_plist->na_c < 0)
 587     {
 588         d_plist->na_c = h_plist->na_ci;
 589     }
 590     else
 591     {
 592         if (d_plist->na_c != h_plist->na_ci)
 593         {
 594             sprintf(sbuf, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
 595                     d_plist->na_c, h_plist->na_ci);
 596             gmx_incons(sbuf);
 597         }
 598     }
 599
 600     if (bDoTime)
 601     {
 602         stat = cudaEventRecord(nb->timers->start_pl_h2d[iloc], stream);
 603         CU_RET_ERR(stat, "cudaEventRecord failed");
 604         nb->timers->didPairlistH2D[iloc] = true;
 605     }
 606
 607     cu_realloc_buffered((void **)&d_plist->sci, h_plist->sci, sizeof(*d_plist->sci),
 608                         &d_plist->nsci, &d_plist->sci_nalloc,
 609                         h_plist->nsci,
 610                         stream, true);
 611
 612     cu_realloc_buffered((void **)&d_plist->cj4, h_plist->cj4, sizeof(*d_plist->cj4),
 613                         &d_plist->ncj4, &d_plist->cj4_nalloc,
 614                         h_plist->ncj4,
 615                         stream, true);
 616
 617     /* this call only allocates space on the device (no data is transferred) */
 618     cu_realloc_buffered((void **)&d_plist->imask, NULL, sizeof(*d_plist->imask),
 619                         &d_plist->nimask, &d_plist->imask_nalloc,
 620                         h_plist->ncj4*c_nbnxnGpuClusterpairSplit,
 621                         stream, true);
 622
 623     cu_realloc_buffered((void **)&d_plist->excl, h_plist->excl, sizeof(*d_plist->excl),
 624                         &d_plist->nexcl, &d_plist->excl_nalloc,
 625                         h_plist->nexcl,
 626                         stream, true);
 627
 628     if (bDoTime)
 629     {
 630         stat = cudaEventRecord(nb->timers->stop_pl_h2d[iloc], stream);
 631         CU_RET_ERR(stat, "cudaEventRecord failed");
 632     }
 633
 634     /* the next use of thist list we be the first one, so we need to prune */
 635     d_plist->haveFreshList = true;
 636 }
 637
 638 void nbnxn_gpu_upload_shiftvec(gmx_nbnxn_cuda_t       *nb,
 639                                const nbnxn_atomdata_t *nbatom)
 640 {
 641     cu_atomdata_t *adat  = nb->atdat;
 642     cudaStream_t   ls    = nb->stream[eintLocal];
 643
 644     /* only if we have a dynamic box */
 645     if (nbatom->bDynamicBox || !adat->bShiftVecUploaded)
 646     {
 647         cu_copy_H2D_async(adat->shift_vec, nbatom->shift_vec,
 648                           SHIFTS * sizeof(*adat->shift_vec), ls);
 649         adat->bShiftVecUploaded = true;
 650     }
 651 }
 652
 653 /*! Clears the first natoms_clear elements of the GPU nonbonded force output array. */
 654 static void nbnxn_cuda_clear_f(gmx_nbnxn_cuda_t *nb, int natoms_clear)
 655 {
 656     cudaError_t    stat;
 657     cu_atomdata_t *adat  = nb->atdat;
 658     cudaStream_t   ls    = nb->stream[eintLocal];
 659
 660     stat = cudaMemsetAsync(adat->f, 0, natoms_clear * sizeof(*adat->f), ls);
 661     CU_RET_ERR(stat, "cudaMemsetAsync on f falied");
 662 }
 663
 664 /*! Clears nonbonded shift force output array and energy outputs on the GPU. */
 665 static void nbnxn_cuda_clear_e_fshift(gmx_nbnxn_cuda_t *nb)
 666 {
 667     cudaError_t    stat;
 668     cu_atomdata_t *adat  = nb->atdat;
 669     cudaStream_t   ls    = nb->stream[eintLocal];
 670
 671     stat = cudaMemsetAsync(adat->fshift, 0, SHIFTS * sizeof(*adat->fshift), ls);
 672     CU_RET_ERR(stat, "cudaMemsetAsync on fshift falied");
 673     stat = cudaMemsetAsync(adat->e_lj, 0, sizeof(*adat->e_lj), ls);
 674     CU_RET_ERR(stat, "cudaMemsetAsync on e_lj falied");
 675     stat = cudaMemsetAsync(adat->e_el, 0, sizeof(*adat->e_el), ls);
 676     CU_RET_ERR(stat, "cudaMemsetAsync on e_el falied");
 677 }
 678
 679 void nbnxn_gpu_clear_outputs(gmx_nbnxn_cuda_t *nb, int flags)
 680 {
 681     nbnxn_cuda_clear_f(nb, nb->atdat->natoms);
 682     /* clear shift force array and energies if the outputs were
 683        used in the current step */
 684     if (flags & GMX_FORCE_VIRIAL)
 685     {
 686         nbnxn_cuda_clear_e_fshift(nb);
 687     }
 688 }
 689
 690 void nbnxn_gpu_init_atomdata(gmx_nbnxn_cuda_t              *nb,
 691                              const struct nbnxn_atomdata_t *nbat)
 692 {
 693     cudaError_t    stat;
 694     int            nalloc, natoms;
 695     bool           realloced;
 696     bool           bDoTime   = nb->bDoTime;
 697     cu_timers_t   *timers    = nb->timers;
 698     cu_atomdata_t *d_atdat   = nb->atdat;
 699     cudaStream_t   ls        = nb->stream[eintLocal];
 700
 701     natoms    = nbat->natoms;
 702     realloced = false;
 703
 704     if (bDoTime)
 705     {
 706         /* time async copy */
 707         stat = cudaEventRecord(timers->start_atdat, ls);
 708         CU_RET_ERR(stat, "cudaEventRecord failed");
 709     }
 710
 711     /* need to reallocate if we have to copy more atoms than the amount of space
 712        available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
 713     if (natoms > d_atdat->nalloc)
 714     {
 715         nalloc = over_alloc_small(natoms);
 716
 717         /* free up first if the arrays have already been initialized */
 718         if (d_atdat->nalloc != -1)
 719         {
 720             cu_free_buffered(d_atdat->f, &d_atdat->natoms, &d_atdat->nalloc);
 721             cu_free_buffered(d_atdat->xq);
 722             cu_free_buffered(d_atdat->atom_types);
 723             cu_free_buffered(d_atdat->lj_comb);
 724         }
 725
 726         stat = cudaMalloc((void **)&d_atdat->f, nalloc*sizeof(*d_atdat->f));
 727         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->f");
 728         stat = cudaMalloc((void **)&d_atdat->xq, nalloc*sizeof(*d_atdat->xq));
 729         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->xq");
 730         if (useLjCombRule(nb->nbparam))
 731         {
 732             stat = cudaMalloc((void **)&d_atdat->lj_comb, nalloc*sizeof(*d_atdat->lj_comb));
 733             CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->lj_comb");
 734         }
 735         else
 736         {
 737             stat = cudaMalloc((void **)&d_atdat->atom_types, nalloc*sizeof(*d_atdat->atom_types));
 738             CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->atom_types");
 739         }
 740
 741         d_atdat->nalloc = nalloc;
 742         realloced       = true;
 743     }
 744
 745     d_atdat->natoms       = natoms;
 746     d_atdat->natoms_local = nbat->natoms_local;
 747
 748     /* need to clear GPU f output if realloc happened */
 749     if (realloced)
 750     {
 751         nbnxn_cuda_clear_f(nb, nalloc);
 752     }
 753
 754     if (useLjCombRule(nb->nbparam))
 755     {
 756         cu_copy_H2D_async(d_atdat->lj_comb, nbat->lj_comb,
 757                           natoms*sizeof(*d_atdat->lj_comb), ls);
 758     }
 759     else
 760     {
 761         cu_copy_H2D_async(d_atdat->atom_types, nbat->type,
 762                           natoms*sizeof(*d_atdat->atom_types), ls);
 763     }
 764
 765     if (bDoTime)
 766     {
 767         stat = cudaEventRecord(timers->stop_atdat, ls);
 768         CU_RET_ERR(stat, "cudaEventRecord failed");
 769     }
 770 }
 771
 772 static void nbnxn_cuda_free_nbparam_table(cu_nbparam_t            *nbparam,
 773                                           const gmx_device_info_t *dev_info)
 774 {
 775     cudaError_t stat;
 776
 777     if (nbparam->eeltype == eelCuEWALD_TAB || nbparam->eeltype == eelCuEWALD_TAB_TWIN)
 778     {
 779         if (!c_disableCudaTextures)
 780         {
 781             /* Only device CC >= 3.0 (Kepler and later) support texture objects */
 782             if (use_texobj(dev_info))
 783             {
 784                 stat = cudaDestroyTextureObject(nbparam->coulomb_tab_texobj);
 785                 CU_RET_ERR(stat, "cudaDestroyTextureObject on coulomb_tab_texobj failed");
 786             }
 787             else
 788             {
 789                 GMX_UNUSED_VALUE(dev_info);
 790                 stat = cudaUnbindTexture(nbnxn_cuda_get_coulomb_tab_texref());
 791                 CU_RET_ERR(stat, "cudaUnbindTexture on coulomb_tab_texref failed");
 792             }
 793         }
 794         cu_free_buffered(nbparam->coulomb_tab);
 795     }
 796 }
 797
 798 void nbnxn_gpu_free(gmx_nbnxn_cuda_t *nb)
 799 {
 800     cudaError_t      stat;
 801     cu_atomdata_t   *atdat;
 802     cu_nbparam_t    *nbparam;
 803     cu_plist_t      *plist, *plist_nl;
 804     cu_timers_t     *timers;
 805
 806     if (nb == NULL)
 807     {
 808         return;
 809     }
 810
 811     atdat       = nb->atdat;
 812     nbparam     = nb->nbparam;
 813     plist       = nb->plist[eintLocal];
 814     plist_nl    = nb->plist[eintNonlocal];
 815     timers      = nb->timers;
 816
 817     nbnxn_cuda_free_nbparam_table(nbparam, nb->dev_info);
 818
 819     stat = cudaEventDestroy(nb->nonlocal_done);
 820     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->nonlocal_done");
 821     stat = cudaEventDestroy(nb->misc_ops_and_local_H2D_done);
 822     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->misc_ops_and_local_H2D_done");
 823
 824     if (nb->bDoTime)
 825     {
 826         stat = cudaEventDestroy(timers->start_atdat);
 827         CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_atdat");
 828         stat = cudaEventDestroy(timers->stop_atdat);
 829         CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_atdat");
 830
 831         /* The non-local counters/stream (second in the array) are needed only with DD. */
 832         for (int i = 0; i <= (nb->bUseTwoStreams ? 1 : 0); i++)
 833         {
 834             stat = cudaEventDestroy(timers->start_nb_k[i]);
 835             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_k");
 836             stat = cudaEventDestroy(timers->stop_nb_k[i]);
 837             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_k");
 838
 839             stat = cudaEventDestroy(timers->start_prune_k[i]);
 840             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_prune_k");
 841             stat = cudaEventDestroy(timers->stop_prune_k[i]);
 842             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_prune_k");
 843
 844             stat = cudaEventDestroy(timers->start_rollingPrune_k[i]);
 845             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_rollingPrune_k");
 846             stat = cudaEventDestroy(timers->stop_rollingPrune_k[i]);
 847             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_rollingPrune_k");
 848
 849             stat = cudaEventDestroy(timers->start_pl_h2d[i]);
 850             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_pl_h2d");
 851             stat = cudaEventDestroy(timers->stop_pl_h2d[i]);
 852             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_pl_h2d");
 853
 854             stat = cudaStreamDestroy(nb->stream[i]);
 855             CU_RET_ERR(stat, "cudaStreamDestroy failed on stream");
 856
 857             stat = cudaEventDestroy(timers->start_nb_h2d[i]);
 858             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_h2d");
 859             stat = cudaEventDestroy(timers->stop_nb_h2d[i]);
 860             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_h2d");
 861
 862             stat = cudaEventDestroy(timers->start_nb_d2h[i]);
 863             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_d2h");
 864             stat = cudaEventDestroy(timers->stop_nb_d2h[i]);
 865             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_d2h");
 866         }
 867     }
 868
 869     if (!useLjCombRule(nb->nbparam))
 870     {
 871         if (!c_disableCudaTextures)
 872         {
 873             /* Only device CC >= 3.0 (Kepler and later) support texture objects */
 874             if (use_texobj(nb->dev_info))
 875             {
 876                 stat = cudaDestroyTextureObject(nbparam->nbfp_texobj);
 877                 CU_RET_ERR(stat, "cudaDestroyTextureObject on nbfp_texobj failed");
 878             }
 879             else
 880             {
 881                 stat = cudaUnbindTexture(nbnxn_cuda_get_nbfp_texref());
 882                 CU_RET_ERR(stat, "cudaUnbindTexture on nbfp_texref failed");
 883             }
 884         }
 885         cu_free_buffered(nbparam->nbfp);
 886     }
 887
 888     if (nbparam->vdwtype == evdwCuEWALDGEOM || nbparam->vdwtype == evdwCuEWALDLB)
 889     {
 890         if (!c_disableCudaTextures)
 891         {
 892             /* Only device CC >= 3.0 (Kepler and later) support texture objects */
 893             if (use_texobj(nb->dev_info))
 894             {
 895                 stat = cudaDestroyTextureObject(nbparam->nbfp_comb_texobj);
 896                 CU_RET_ERR(stat, "cudaDestroyTextureObject on nbfp_comb_texobj failed");
 897             }
 898             else
 899             {
 900                 stat = cudaUnbindTexture(nbnxn_cuda_get_nbfp_comb_texref());
 901                 CU_RET_ERR(stat, "cudaUnbindTexture on nbfp_comb_texref failed");
 902             }
 903         }
 904         cu_free_buffered(nbparam->nbfp_comb);
 905     }
 906
 907     stat = cudaFree(atdat->shift_vec);
 908     CU_RET_ERR(stat, "cudaFree failed on atdat->shift_vec");
 909     stat = cudaFree(atdat->fshift);
 910     CU_RET_ERR(stat, "cudaFree failed on atdat->fshift");
 911
 912     stat = cudaFree(atdat->e_lj);
 913     CU_RET_ERR(stat, "cudaFree failed on atdat->e_lj");
 914     stat = cudaFree(atdat->e_el);
 915     CU_RET_ERR(stat, "cudaFree failed on atdat->e_el");
 916
 917     cu_free_buffered(atdat->f, &atdat->natoms, &atdat->nalloc);
 918     cu_free_buffered(atdat->xq);
 919     cu_free_buffered(atdat->atom_types, &atdat->ntypes);
 920     cu_free_buffered(atdat->lj_comb);
 921
 922     cu_free_buffered(plist->sci, &plist->nsci, &plist->sci_nalloc);
 923     cu_free_buffered(plist->cj4, &plist->ncj4, &plist->cj4_nalloc);
 924     cu_free_buffered(plist->imask, &plist->nimask, &plist->imask_nalloc);
 925     cu_free_buffered(plist->excl, &plist->nexcl, &plist->excl_nalloc);
 926     if (nb->bUseTwoStreams)
 927     {
 928         cu_free_buffered(plist_nl->sci, &plist_nl->nsci, &plist_nl->sci_nalloc);
 929         cu_free_buffered(plist_nl->cj4, &plist_nl->ncj4, &plist_nl->cj4_nalloc);
 930         cu_free_buffered(plist_nl->imask, &plist_nl->nimask, &plist_nl->imask_nalloc);
 931         cu_free_buffered(plist_nl->excl, &plist_nl->nexcl, &plist->excl_nalloc);
 932     }
 933
 934     sfree(atdat);
 935     sfree(nbparam);
 936     sfree(plist);
 937     if (nb->bUseTwoStreams)
 938     {
 939         sfree(plist_nl);
 940     }
 941     sfree(timers);
 942     sfree(nb->timings);
 943     sfree(nb);
 944
 945     if (debug)
 946     {
 947         fprintf(debug, "Cleaned up CUDA data structures.\n");
 948     }
 949 }
 950
 951 gmx_wallclock_gpu_t * nbnxn_gpu_get_timings(gmx_nbnxn_cuda_t *nb)
 952 {
 953     return (nb != NULL && nb->bDoTime) ? nb->timings : NULL;
 954 }
 955
 956 void nbnxn_gpu_reset_timings(nonbonded_verlet_t* nbv)
 957 {
 958     if (nbv->gpu_nbv && nbv->gpu_nbv->bDoTime)
 959     {
 960         init_timings(nbv->gpu_nbv->timings);
 961     }
 962 }
 963
 964 int nbnxn_gpu_min_ci_balanced(gmx_nbnxn_cuda_t *nb)
 965 {
 966     return nb != NULL ?
 967            gpu_min_ci_balanced_factor*nb->dev_info->prop.multiProcessorCount : 0;
 968
 969 }
 970
 971 gmx_bool nbnxn_gpu_is_kernel_ewald_analytical(const gmx_nbnxn_cuda_t *nb)
 972 {
 973     return ((nb->nbparam->eeltype == eelCuEWALD_ANA) ||
 974             (nb->nbparam->eeltype == eelCuEWALD_ANA_TWIN));
 975 }