Move smalloc.h to utility/
[alexxy/gromacs.git] / 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, 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 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <stdlib.h>
40 #include <stdio.h>
41 #include <assert.h>
42
43 #include <cuda.h>
44
45 #include "gmx_fatal.h"
46 #include "gromacs/utility/smalloc.h"
47 #include "tables.h"
48 #include "typedefs.h"
49 #include "types/enums.h"
50 #include "types/nb_verlet.h"
51 #include "types/interaction_const.h"
52 #include "types/force_flags.h"
53 #include "../nbnxn_consts.h"
54 #include "gmx_detect_hardware.h"
55
56 #include "nbnxn_cuda_types.h"
57 #include "../../gmxlib/cuda_tools/cudautils.cuh"
58 #include "nbnxn_cuda_data_mgmt.h"
59 #include "pmalloc_cuda.h"
60 #include "gpu_utils.h"
61
62 #include "gromacs/utility/common.h"
63
64 static bool bUseCudaEventBlockingSync = false; /* makes the CPU thread block */
65
66 /* This is a heuristically determined parameter for the Fermi architecture for
67  * the minimum size of ci lists by multiplying this constant with the # of
68  * multiprocessors on the current device.
69  */
70 static unsigned int gpu_min_ci_balanced_factor = 40;
71
72 /* Functions from nbnxn_cuda.cu */
73 extern void nbnxn_cuda_set_cacheconfig(cuda_dev_info_t *devinfo);
74 extern const struct texture<float, 1, cudaReadModeElementType> &nbnxn_cuda_get_nbfp_texref();
75 extern const struct texture<float, 1, cudaReadModeElementType> &nbnxn_cuda_get_nbfp_comb_texref();
76 extern const struct texture<float, 1, cudaReadModeElementType> &nbnxn_cuda_get_coulomb_tab_texref();
77
78 /* We should actually be using md_print_warn in md_logging.c,
79  * but we can't include mpi.h in CUDA code.
80  */
81 static void md_print_warn(FILE       *fplog,
82                           const char *fmt, ...)
83 {
84     va_list ap;
85
86     if (fplog != NULL)
87     {
88         /* We should only print to stderr on the master node,
89          * in most cases fplog is only set on the master node, so this works.
90          */
91         va_start(ap, fmt);
92         fprintf(stderr, "\n");
93         vfprintf(stderr, fmt, ap);
94         fprintf(stderr, "\n");
95         va_end(ap);
96
97         va_start(ap, fmt);
98         fprintf(fplog, "\n");
99         vfprintf(fplog, fmt, ap);
100         fprintf(fplog, "\n");
101         va_end(ap);
102     }
103 }
104
105
106 /* Fw. decl. */
107 static void nbnxn_cuda_clear_e_fshift(nbnxn_cuda_ptr_t cu_nb);
108
109
110 /*! Tabulates the Ewald Coulomb force and initializes the size/scale
111     and the table GPU array. If called with an already allocated table,
112     it just re-uploads the table.
113  */
114 static void init_ewald_coulomb_force_table(cu_nbparam_t          *nbp,
115                                            const cuda_dev_info_t *dev_info)
116 {
117     float       *ftmp, *coul_tab;
118     int          tabsize;
119     double       tabscale;
120     cudaError_t  stat;
121
122     tabsize     = GPU_EWALD_COULOMB_FORCE_TABLE_SIZE;
123     /* Subtract 2 iso 1 to avoid access out of range due to rounding */
124     tabscale    = (tabsize - 2) / sqrt(nbp->rcoulomb_sq);
125
126     pmalloc((void**)&ftmp, tabsize*sizeof(*ftmp));
127
128     table_spline3_fill_ewald_lr(ftmp, NULL, NULL, tabsize,
129                                 1/tabscale, nbp->ewald_beta, v_q_ewald_lr);
130
131     /* If the table pointer == NULL the table is generated the first time =>
132        the array pointer will be saved to nbparam and the texture is bound.
133      */
134     coul_tab = nbp->coulomb_tab;
135     if (coul_tab == NULL)
136     {
137         stat = cudaMalloc((void **)&coul_tab, tabsize*sizeof(*coul_tab));
138         CU_RET_ERR(stat, "cudaMalloc failed on coul_tab");
139
140         nbp->coulomb_tab = coul_tab;
141
142 #ifdef TEXOBJ_SUPPORTED
143         /* Only device CC >= 3.0 (Kepler and later) support texture objects */
144         if (dev_info->prop.major >= 3)
145         {
146             cudaResourceDesc rd;
147             memset(&rd, 0, sizeof(rd));
148             rd.resType                  = cudaResourceTypeLinear;
149             rd.res.linear.devPtr        = nbp->coulomb_tab;
150             rd.res.linear.desc.f        = cudaChannelFormatKindFloat;
151             rd.res.linear.desc.x        = 32;
152             rd.res.linear.sizeInBytes   = tabsize*sizeof(*coul_tab);
153
154             cudaTextureDesc td;
155             memset(&td, 0, sizeof(td));
156             td.readMode                 = cudaReadModeElementType;
157             stat = cudaCreateTextureObject(&nbp->coulomb_tab_texobj, &rd, &td, NULL);
158             CU_RET_ERR(stat, "cudaCreateTextureObject on coulomb_tab_texobj failed");
159         }
160         else
161 #endif
162         {
163             GMX_UNUSED_VALUE(dev_info);
164             cudaChannelFormatDesc cd   = cudaCreateChannelDesc<float>();
165             stat = cudaBindTexture(NULL, &nbnxn_cuda_get_coulomb_tab_texref(),
166                                    coul_tab, &cd, tabsize*sizeof(*coul_tab));
167             CU_RET_ERR(stat, "cudaBindTexture on coulomb_tab_texref failed");
168         }
169     }
170
171     cu_copy_H2D(coul_tab, ftmp, tabsize*sizeof(*coul_tab));
172
173     nbp->coulomb_tab_size     = tabsize;
174     nbp->coulomb_tab_scale    = tabscale;
175
176     pfree(ftmp);
177 }
178
179
180 /*! Initializes the atomdata structure first time, it only gets filled at
181     pair-search. */
182 static void init_atomdata_first(cu_atomdata_t *ad, int ntypes)
183 {
184     cudaError_t stat;
185
186     ad->ntypes  = ntypes;
187     stat        = cudaMalloc((void**)&ad->shift_vec, SHIFTS*sizeof(*ad->shift_vec));
188     CU_RET_ERR(stat, "cudaMalloc failed on ad->shift_vec");
189     ad->bShiftVecUploaded = false;
190
191     stat = cudaMalloc((void**)&ad->fshift, SHIFTS*sizeof(*ad->fshift));
192     CU_RET_ERR(stat, "cudaMalloc failed on ad->fshift");
193
194     stat = cudaMalloc((void**)&ad->e_lj, sizeof(*ad->e_lj));
195     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_lj");
196     stat = cudaMalloc((void**)&ad->e_el, sizeof(*ad->e_el));
197     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_el");
198
199     /* initialize to NULL poiters to data that is not allocated here and will
200        need reallocation in nbnxn_cuda_init_atomdata */
201     ad->xq = NULL;
202     ad->f  = NULL;
203
204     /* size -1 indicates that the respective array hasn't been initialized yet */
205     ad->natoms = -1;
206     ad->nalloc = -1;
207 }
208
209 /*! Selects the Ewald kernel type, analytical on SM 3.0 and later, tabulated on
210     earlier GPUs, single or twin cut-off. */
211 static int pick_ewald_kernel_type(bool                   bTwinCut,
212                                   const cuda_dev_info_t *dev_info)
213 {
214     bool bUseAnalyticalEwald, bForceAnalyticalEwald, bForceTabulatedEwald;
215     int  kernel_type;
216
217     /* Benchmarking/development environment variables to force the use of
218        analytical or tabulated Ewald kernel. */
219     bForceAnalyticalEwald = (getenv("GMX_CUDA_NB_ANA_EWALD") != NULL);
220     bForceTabulatedEwald  = (getenv("GMX_CUDA_NB_TAB_EWALD") != NULL);
221
222     if (bForceAnalyticalEwald && bForceTabulatedEwald)
223     {
224         gmx_incons("Both analytical and tabulated Ewald CUDA non-bonded kernels "
225                    "requested through environment variables.");
226     }
227
228     /* By default, on SM 3.0 and later use analytical Ewald, on earlier tabulated. */
229     if ((dev_info->prop.major >= 3 || bForceAnalyticalEwald) && !bForceTabulatedEwald)
230     {
231         bUseAnalyticalEwald = true;
232
233         if (debug)
234         {
235             fprintf(debug, "Using analytical Ewald CUDA kernels\n");
236         }
237     }
238     else
239     {
240         bUseAnalyticalEwald = false;
241
242         if (debug)
243         {
244             fprintf(debug, "Using tabulated Ewald CUDA kernels\n");
245         }
246     }
247
248     /* Use twin cut-off kernels if requested by bTwinCut or the env. var.
249        forces it (use it for debugging/benchmarking only). */
250     if (!bTwinCut && (getenv("GMX_CUDA_NB_EWALD_TWINCUT") == NULL))
251     {
252         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA : eelCuEWALD_TAB;
253     }
254     else
255     {
256         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA_TWIN : eelCuEWALD_TAB_TWIN;
257     }
258
259     return kernel_type;
260 }
261
262 /*! Copies all parameters related to the cut-off from ic to nbp */
263 static void set_cutoff_parameters(cu_nbparam_t              *nbp,
264                                   const interaction_const_t *ic)
265 {
266     nbp->ewald_beta       = ic->ewaldcoeff_q;
267     nbp->sh_ewald         = ic->sh_ewald;
268     nbp->epsfac           = ic->epsfac;
269     nbp->two_k_rf         = 2.0 * ic->k_rf;
270     nbp->c_rf             = ic->c_rf;
271     nbp->rvdw_sq          = ic->rvdw * ic->rvdw;
272     nbp->rcoulomb_sq      = ic->rcoulomb * ic->rcoulomb;
273     nbp->rlist_sq         = ic->rlist * ic->rlist;
274
275     nbp->sh_lj_ewald      = ic->sh_lj_ewald;
276     nbp->ewaldcoeff_lj    = ic->ewaldcoeff_lj;
277
278     nbp->rvdw_switch      = ic->rvdw_switch;
279     nbp->dispersion_shift = ic->dispersion_shift;
280     nbp->repulsion_shift  = ic->repulsion_shift;
281     nbp->vdw_switch       = ic->vdw_switch;
282 }
283
284 /*! Initializes the nonbonded parameter data structure. */
285 static void init_nbparam(cu_nbparam_t              *nbp,
286                          const interaction_const_t *ic,
287                          const nbnxn_atomdata_t    *nbat,
288                          const cuda_dev_info_t     *dev_info)
289 {
290     cudaError_t stat;
291     int         ntypes, nnbfp, nnbfp_comb;
292
293     ntypes  = nbat->ntype;
294
295     set_cutoff_parameters(nbp, ic);
296
297     if (ic->vdwtype == evdwCUT)
298     {
299         switch (ic->vdw_modifier)
300         {
301             case eintmodNONE:
302             case eintmodPOTSHIFT:
303                 nbp->vdwtype = evdwCuCUT;
304                 break;
305             case eintmodFORCESWITCH:
306                 nbp->vdwtype = evdwCuFSWITCH;
307                 break;
308             case eintmodPOTSWITCH:
309                 nbp->vdwtype = evdwCuPSWITCH;
310                 break;
311             default:
312                 gmx_incons("The requested VdW interaction modifier is not implemented in the CUDA GPU accelerated kernels!");
313                 break;
314         }
315     }
316     else if (ic->vdwtype == evdwPME)
317     {
318         if (ic->ljpme_comb_rule == ljcrGEOM)
319         {
320             assert(nbat->comb_rule == ljcrGEOM);
321             nbp->vdwtype = evdwCuEWALDGEOM;
322         }
323         else
324         {
325             assert(nbat->comb_rule == ljcrLB);
326             nbp->vdwtype = evdwCuEWALDLB;
327         }
328     }
329     else
330     {
331         gmx_incons("The requested VdW type is not implemented in the CUDA GPU accelerated kernels!");
332     }
333
334     if (ic->eeltype == eelCUT)
335     {
336         nbp->eeltype = eelCuCUT;
337     }
338     else if (EEL_RF(ic->eeltype))
339     {
340         nbp->eeltype = eelCuRF;
341     }
342     else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
343     {
344         /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
345         nbp->eeltype = pick_ewald_kernel_type(false, dev_info);
346     }
347     else
348     {
349         /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
350         gmx_incons("The requested electrostatics type is not implemented in the CUDA GPU accelerated kernels!");
351     }
352
353     /* generate table for PME */
354     nbp->coulomb_tab = NULL;
355     if (nbp->eeltype == eelCuEWALD_TAB || nbp->eeltype == eelCuEWALD_TAB_TWIN)
356     {
357         init_ewald_coulomb_force_table(nbp, dev_info);
358     }
359
360     nnbfp      = 2*ntypes*ntypes;
361     nnbfp_comb = 2*ntypes;
362
363     stat  = cudaMalloc((void **)&nbp->nbfp, nnbfp*sizeof(*nbp->nbfp));
364     CU_RET_ERR(stat, "cudaMalloc failed on nbp->nbfp");
365     cu_copy_H2D(nbp->nbfp, nbat->nbfp, nnbfp*sizeof(*nbp->nbfp));
366
367
368     if (ic->vdwtype == evdwPME)
369     {
370         stat  = cudaMalloc((void **)&nbp->nbfp_comb, nnbfp_comb*sizeof(*nbp->nbfp_comb));
371         CU_RET_ERR(stat, "cudaMalloc failed on nbp->nbfp_comb");
372         cu_copy_H2D(nbp->nbfp_comb, nbat->nbfp_comb, nnbfp_comb*sizeof(*nbp->nbfp_comb));
373     }
374
375 #ifdef TEXOBJ_SUPPORTED
376     /* Only device CC >= 3.0 (Kepler and later) support texture objects */
377     if (dev_info->prop.major >= 3)
378     {
379         cudaResourceDesc rd;
380         cudaTextureDesc  td;
381
382         memset(&rd, 0, sizeof(rd));
383         rd.resType                  = cudaResourceTypeLinear;
384         rd.res.linear.devPtr        = nbp->nbfp;
385         rd.res.linear.desc.f        = cudaChannelFormatKindFloat;
386         rd.res.linear.desc.x        = 32;
387         rd.res.linear.sizeInBytes   = nnbfp*sizeof(*nbp->nbfp);
388
389         memset(&td, 0, sizeof(td));
390         td.readMode                 = cudaReadModeElementType;
391         stat = cudaCreateTextureObject(&nbp->nbfp_texobj, &rd, &td, NULL);
392         CU_RET_ERR(stat, "cudaCreateTextureObject on nbfp_texobj failed");
393
394         if (ic->vdwtype == evdwPME)
395         {
396             memset(&rd, 0, sizeof(rd));
397             rd.resType                  = cudaResourceTypeLinear;
398             rd.res.linear.devPtr        = nbp->nbfp_comb;
399             rd.res.linear.desc.f        = cudaChannelFormatKindFloat;
400             rd.res.linear.desc.x        = 32;
401             rd.res.linear.sizeInBytes   = nnbfp_comb*sizeof(*nbp->nbfp_comb);
402
403             memset(&td, 0, sizeof(td));
404             td.readMode = cudaReadModeElementType;
405             stat        = cudaCreateTextureObject(&nbp->nbfp_comb_texobj, &rd, &td, NULL);
406             CU_RET_ERR(stat, "cudaCreateTextureObject on nbfp_comb_texobj failed");
407         }
408     }
409     else
410 #endif
411     {
412         cudaChannelFormatDesc cd = cudaCreateChannelDesc<float>();
413         stat = cudaBindTexture(NULL, &nbnxn_cuda_get_nbfp_texref(),
414                                nbp->nbfp, &cd, nnbfp*sizeof(*nbp->nbfp));
415         CU_RET_ERR(stat, "cudaBindTexture on nbfp_texref failed");
416
417         if (ic->vdwtype == evdwPME)
418         {
419             stat = cudaBindTexture(NULL, &nbnxn_cuda_get_nbfp_comb_texref(),
420                                    nbp->nbfp_comb, &cd, nnbfp_comb*sizeof(*nbp->nbfp_comb));
421             CU_RET_ERR(stat, "cudaBindTexture on nbfp_comb_texref failed");
422         }
423     }
424 }
425
426 /*! Re-generate the GPU Ewald force table, resets rlist, and update the
427  *  electrostatic type switching to twin cut-off (or back) if needed. */
428 void nbnxn_cuda_pme_loadbal_update_param(nbnxn_cuda_ptr_t           cu_nb,
429                                          const interaction_const_t *ic)
430 {
431     cu_nbparam_t *nbp = cu_nb->nbparam;
432
433     set_cutoff_parameters(nbp, ic);
434
435     nbp->eeltype        = pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw,
436                                                  cu_nb->dev_info);
437
438     init_ewald_coulomb_force_table(cu_nb->nbparam, cu_nb->dev_info);
439 }
440
441 /*! Initializes the pair list data structure. */
442 static void init_plist(cu_plist_t *pl)
443 {
444     /* initialize to NULL pointers to data that is not allocated here and will
445        need reallocation in nbnxn_cuda_init_pairlist */
446     pl->sci     = NULL;
447     pl->cj4     = NULL;
448     pl->excl    = NULL;
449
450     /* size -1 indicates that the respective array hasn't been initialized yet */
451     pl->na_c        = -1;
452     pl->nsci        = -1;
453     pl->sci_nalloc  = -1;
454     pl->ncj4        = -1;
455     pl->cj4_nalloc  = -1;
456     pl->nexcl       = -1;
457     pl->excl_nalloc = -1;
458     pl->bDoPrune    = false;
459 }
460
461 /*! Initializes the timer data structure. */
462 static void init_timers(cu_timers_t *t, bool bUseTwoStreams)
463 {
464     cudaError_t stat;
465     int         eventflags = ( bUseCudaEventBlockingSync ? cudaEventBlockingSync : cudaEventDefault );
466
467     stat = cudaEventCreateWithFlags(&(t->start_atdat), eventflags);
468     CU_RET_ERR(stat, "cudaEventCreate on start_atdat failed");
469     stat = cudaEventCreateWithFlags(&(t->stop_atdat), eventflags);
470     CU_RET_ERR(stat, "cudaEventCreate on stop_atdat failed");
471
472     /* The non-local counters/stream (second in the array) are needed only with DD. */
473     for (int i = 0; i <= (bUseTwoStreams ? 1 : 0); i++)
474     {
475         stat = cudaEventCreateWithFlags(&(t->start_nb_k[i]), eventflags);
476         CU_RET_ERR(stat, "cudaEventCreate on start_nb_k failed");
477         stat = cudaEventCreateWithFlags(&(t->stop_nb_k[i]), eventflags);
478         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_k failed");
479
480
481         stat = cudaEventCreateWithFlags(&(t->start_pl_h2d[i]), eventflags);
482         CU_RET_ERR(stat, "cudaEventCreate on start_pl_h2d failed");
483         stat = cudaEventCreateWithFlags(&(t->stop_pl_h2d[i]), eventflags);
484         CU_RET_ERR(stat, "cudaEventCreate on stop_pl_h2d failed");
485
486         stat = cudaEventCreateWithFlags(&(t->start_nb_h2d[i]), eventflags);
487         CU_RET_ERR(stat, "cudaEventCreate on start_nb_h2d failed");
488         stat = cudaEventCreateWithFlags(&(t->stop_nb_h2d[i]), eventflags);
489         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_h2d failed");
490
491         stat = cudaEventCreateWithFlags(&(t->start_nb_d2h[i]), eventflags);
492         CU_RET_ERR(stat, "cudaEventCreate on start_nb_d2h failed");
493         stat = cudaEventCreateWithFlags(&(t->stop_nb_d2h[i]), eventflags);
494         CU_RET_ERR(stat, "cudaEventCreate on stop_nb_d2h failed");
495     }
496 }
497
498 /*! Initializes the timings data structure. */
499 static void init_timings(wallclock_gpu_t *t)
500 {
501     int i, j;
502
503     t->nb_h2d_t = 0.0;
504     t->nb_d2h_t = 0.0;
505     t->nb_c     = 0;
506     t->pl_h2d_t = 0.0;
507     t->pl_h2d_c = 0;
508     for (i = 0; i < 2; i++)
509     {
510         for (j = 0; j < 2; j++)
511         {
512             t->ktime[i][j].t = 0.0;
513             t->ktime[i][j].c = 0;
514         }
515     }
516 }
517
518 void nbnxn_cuda_init(FILE                 *fplog,
519                      nbnxn_cuda_ptr_t     *p_cu_nb,
520                      const gmx_gpu_info_t *gpu_info,
521                      const gmx_gpu_opt_t  *gpu_opt,
522                      int                   my_gpu_index,
523                      gmx_bool              bLocalAndNonlocal)
524 {
525     cudaError_t       stat;
526     nbnxn_cuda_ptr_t  nb;
527     char              sbuf[STRLEN];
528     bool              bStreamSync, bNoStreamSync, bTMPIAtomics, bX86, bOldDriver;
529     int               cuda_drv_ver;
530
531     assert(gpu_info);
532
533     if (p_cu_nb == NULL)
534     {
535         return;
536     }
537
538     snew(nb, 1);
539     snew(nb->atdat, 1);
540     snew(nb->nbparam, 1);
541     snew(nb->plist[eintLocal], 1);
542     if (bLocalAndNonlocal)
543     {
544         snew(nb->plist[eintNonlocal], 1);
545     }
546
547     nb->bUseTwoStreams = bLocalAndNonlocal;
548
549     snew(nb->timers, 1);
550     snew(nb->timings, 1);
551
552     /* init nbst */
553     pmalloc((void**)&nb->nbst.e_lj, sizeof(*nb->nbst.e_lj));
554     pmalloc((void**)&nb->nbst.e_el, sizeof(*nb->nbst.e_el));
555     pmalloc((void**)&nb->nbst.fshift, SHIFTS * sizeof(*nb->nbst.fshift));
556
557     init_plist(nb->plist[eintLocal]);
558
559     /* set device info, just point it to the right GPU among the detected ones */
560     nb->dev_info = &gpu_info->cuda_dev[get_gpu_device_id(gpu_info, gpu_opt, my_gpu_index)];
561
562     /* local/non-local GPU streams */
563     stat = cudaStreamCreate(&nb->stream[eintLocal]);
564     CU_RET_ERR(stat, "cudaStreamCreate on stream[eintLocal] failed");
565     if (nb->bUseTwoStreams)
566     {
567         init_plist(nb->plist[eintNonlocal]);
568
569         /* CUDA stream priority available in the CUDA RT 5.5 API.
570          * Note that the device we're running on does not have to support
571          * priorities, because we are querying the priority range which in this
572          * case will be a single value.
573          */
574 #if CUDA_VERSION >= 5500
575         {
576             int highest_priority;
577             stat = cudaDeviceGetStreamPriorityRange(NULL, &highest_priority);
578             CU_RET_ERR(stat, "cudaDeviceGetStreamPriorityRange failed");
579
580             stat = cudaStreamCreateWithPriority(&nb->stream[eintNonlocal],
581                                                 cudaStreamDefault,
582                                                 highest_priority);
583             CU_RET_ERR(stat, "cudaStreamCreateWithPriority on stream[eintNonlocal] failed");
584         }
585 #else
586         stat = cudaStreamCreate(&nb->stream[eintNonlocal]);
587         CU_RET_ERR(stat, "cudaStreamCreate on stream[eintNonlocal] failed");
588 #endif
589     }
590
591     /* init events for sychronization (timing disabled for performance reasons!) */
592     stat = cudaEventCreateWithFlags(&nb->nonlocal_done, cudaEventDisableTiming);
593     CU_RET_ERR(stat, "cudaEventCreate on nonlocal_done failed");
594     stat = cudaEventCreateWithFlags(&nb->misc_ops_done, cudaEventDisableTiming);
595     CU_RET_ERR(stat, "cudaEventCreate on misc_ops_one failed");
596
597     /* On GPUs with ECC enabled, cudaStreamSynchronize shows a large overhead
598      * (which increases with shorter time/step) caused by a known CUDA driver bug.
599      * To work around the issue we'll use an (admittedly fragile) memory polling
600      * waiting to preserve performance. This requires support for atomic
601      * operations and only works on x86/x86_64.
602      * With polling wait event-timing also needs to be disabled.
603      *
604      * The overhead is greatly reduced in API v5.0 drivers and the improvement
605      * is independent of runtime version. Hence, with API v5.0 drivers and later
606      * we won't switch to polling.
607      *
608      * NOTE: Unfortunately, this is known to fail when GPUs are shared by (t)MPI,
609      * ranks so we will also disable it in that case.
610      */
611
612     bStreamSync    = getenv("GMX_CUDA_STREAMSYNC") != NULL;
613     bNoStreamSync  = getenv("GMX_NO_CUDA_STREAMSYNC") != NULL;
614
615 #ifdef TMPI_ATOMICS
616     bTMPIAtomics = true;
617 #else
618     bTMPIAtomics = false;
619 #endif
620
621 #ifdef GMX_TARGET_X86
622     bX86 = true;
623 #else
624     bX86 = false;
625 #endif
626
627     if (bStreamSync && bNoStreamSync)
628     {
629         gmx_fatal(FARGS, "Conflicting environment variables: both GMX_CUDA_STREAMSYNC and GMX_NO_CUDA_STREAMSYNC defined");
630     }
631
632     stat = cudaDriverGetVersion(&cuda_drv_ver);
633     CU_RET_ERR(stat, "cudaDriverGetVersion failed");
634
635     bOldDriver = (cuda_drv_ver < 5000);
636
637     if ((nb->dev_info->prop.ECCEnabled == 1) && bOldDriver)
638     {
639         /* Polling wait should be used instead of cudaStreamSynchronize only if:
640          *   - ECC is ON & driver is old (checked above),
641          *   - we're on x86/x86_64,
642          *   - atomics are available, and
643          *   - GPUs are not being shared.
644          */
645         bool bShouldUsePollSync = (bX86 && bTMPIAtomics &&
646                                    (gmx_count_gpu_dev_shared(gpu_opt) < 1));
647
648         if (bStreamSync)
649         {
650             nb->bUseStreamSync = true;
651
652             /* only warn if polling should be used */
653             if (bShouldUsePollSync)
654             {
655                 md_print_warn(fplog,
656                               "NOTE: Using a GPU with ECC enabled and CUDA driver API version <5.0, but\n"
657                               "      cudaStreamSynchronize waiting is forced by the GMX_CUDA_STREAMSYNC env. var.\n");
658             }
659         }
660         else
661         {
662             nb->bUseStreamSync = !bShouldUsePollSync;
663
664             if (bShouldUsePollSync)
665             {
666                 md_print_warn(fplog,
667                               "NOTE: Using a GPU with ECC enabled and CUDA driver API version <5.0, known to\n"
668                               "      cause performance loss. Switching to the alternative polling GPU wait.\n"
669                               "      If you encounter issues, switch back to standard GPU waiting by setting\n"
670                               "      the GMX_CUDA_STREAMSYNC environment variable.\n");
671             }
672             else
673             {
674                 /* Tell the user that the ECC+old driver combination can be bad */
675                 sprintf(sbuf,
676                         "NOTE: Using a GPU with ECC enabled and CUDA driver API version <5.0.\n"
677                         "      A known bug in this driver version can cause performance loss.\n"
678                         "      However, the polling wait workaround can not be used because\n%s\n"
679                         "      Consider updating the driver or turning ECC off.",
680                         (bX86 && bTMPIAtomics) ?
681                         "      GPU(s) are being oversubscribed." :
682                         "      atomic operations are not supported by the platform/CPU+compiler.");
683                 md_print_warn(fplog, sbuf);
684             }
685         }
686     }
687     else
688     {
689         if (bNoStreamSync)
690         {
691             nb->bUseStreamSync = false;
692
693             md_print_warn(fplog,
694                           "NOTE: Polling wait for GPU synchronization requested by GMX_NO_CUDA_STREAMSYNC\n");
695         }
696         else
697         {
698             /* no/off ECC, cudaStreamSynchronize not turned off by env. var. */
699             nb->bUseStreamSync = true;
700         }
701     }
702
703     /* CUDA timing disabled as event timers don't work:
704        - with multiple streams = domain-decomposition;
705        - with the polling waiting hack (without cudaStreamSynchronize);
706        - when turned off by GMX_DISABLE_CUDA_TIMING.
707      */
708     nb->bDoTime = (!nb->bUseTwoStreams && nb->bUseStreamSync &&
709                    (getenv("GMX_DISABLE_CUDA_TIMING") == NULL));
710
711     if (nb->bDoTime)
712     {
713         init_timers(nb->timers, nb->bUseTwoStreams);
714         init_timings(nb->timings);
715     }
716
717     /* set the kernel type for the current GPU */
718     /* pick L1 cache configuration */
719     nbnxn_cuda_set_cacheconfig(nb->dev_info);
720
721     *p_cu_nb = nb;
722
723     if (debug)
724     {
725         fprintf(debug, "Initialized CUDA data structures.\n");
726     }
727 }
728
729 void nbnxn_cuda_init_const(nbnxn_cuda_ptr_t                cu_nb,
730                            const interaction_const_t      *ic,
731                            const nonbonded_verlet_group_t *nbv_group)
732 {
733     init_atomdata_first(cu_nb->atdat, nbv_group[0].nbat->ntype);
734     init_nbparam(cu_nb->nbparam, ic, nbv_group[0].nbat, cu_nb->dev_info);
735
736     /* clear energy and shift force outputs */
737     nbnxn_cuda_clear_e_fshift(cu_nb);
738 }
739
740 void nbnxn_cuda_init_pairlist(nbnxn_cuda_ptr_t        cu_nb,
741                               const nbnxn_pairlist_t *h_plist,
742                               int                     iloc)
743 {
744     char          sbuf[STRLEN];
745     cudaError_t   stat;
746     bool          bDoTime    = cu_nb->bDoTime;
747     cudaStream_t  stream     = cu_nb->stream[iloc];
748     cu_plist_t   *d_plist    = cu_nb->plist[iloc];
749
750     if (d_plist->na_c < 0)
751     {
752         d_plist->na_c = h_plist->na_ci;
753     }
754     else
755     {
756         if (d_plist->na_c != h_plist->na_ci)
757         {
758             sprintf(sbuf, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
759                     d_plist->na_c, h_plist->na_ci);
760             gmx_incons(sbuf);
761         }
762     }
763
764     if (bDoTime)
765     {
766         stat = cudaEventRecord(cu_nb->timers->start_pl_h2d[iloc], stream);
767         CU_RET_ERR(stat, "cudaEventRecord failed");
768     }
769
770     cu_realloc_buffered((void **)&d_plist->sci, h_plist->sci, sizeof(*d_plist->sci),
771                         &d_plist->nsci, &d_plist->sci_nalloc,
772                         h_plist->nsci,
773                         stream, true);
774
775     cu_realloc_buffered((void **)&d_plist->cj4, h_plist->cj4, sizeof(*d_plist->cj4),
776                         &d_plist->ncj4, &d_plist->cj4_nalloc,
777                         h_plist->ncj4,
778                         stream, true);
779
780     cu_realloc_buffered((void **)&d_plist->excl, h_plist->excl, sizeof(*d_plist->excl),
781                         &d_plist->nexcl, &d_plist->excl_nalloc,
782                         h_plist->nexcl,
783                         stream, true);
784
785     if (bDoTime)
786     {
787         stat = cudaEventRecord(cu_nb->timers->stop_pl_h2d[iloc], stream);
788         CU_RET_ERR(stat, "cudaEventRecord failed");
789     }
790
791     /* need to prune the pair list during the next step */
792     d_plist->bDoPrune = true;
793 }
794
795 void nbnxn_cuda_upload_shiftvec(nbnxn_cuda_ptr_t        cu_nb,
796                                 const nbnxn_atomdata_t *nbatom)
797 {
798     cu_atomdata_t *adat  = cu_nb->atdat;
799     cudaStream_t   ls    = cu_nb->stream[eintLocal];
800
801     /* only if we have a dynamic box */
802     if (nbatom->bDynamicBox || !adat->bShiftVecUploaded)
803     {
804         cu_copy_H2D_async(adat->shift_vec, nbatom->shift_vec,
805                           SHIFTS * sizeof(*adat->shift_vec), ls);
806         adat->bShiftVecUploaded = true;
807     }
808 }
809
810 /*! Clears the first natoms_clear elements of the GPU nonbonded force output array. */
811 static void nbnxn_cuda_clear_f(nbnxn_cuda_ptr_t cu_nb, int natoms_clear)
812 {
813     cudaError_t    stat;
814     cu_atomdata_t *adat  = cu_nb->atdat;
815     cudaStream_t   ls    = cu_nb->stream[eintLocal];
816
817     stat = cudaMemsetAsync(adat->f, 0, natoms_clear * sizeof(*adat->f), ls);
818     CU_RET_ERR(stat, "cudaMemsetAsync on f falied");
819 }
820
821 /*! Clears nonbonded shift force output array and energy outputs on the GPU. */
822 static void nbnxn_cuda_clear_e_fshift(nbnxn_cuda_ptr_t cu_nb)
823 {
824     cudaError_t    stat;
825     cu_atomdata_t *adat  = cu_nb->atdat;
826     cudaStream_t   ls    = cu_nb->stream[eintLocal];
827
828     stat = cudaMemsetAsync(adat->fshift, 0, SHIFTS * sizeof(*adat->fshift), ls);
829     CU_RET_ERR(stat, "cudaMemsetAsync on fshift falied");
830     stat = cudaMemsetAsync(adat->e_lj, 0, sizeof(*adat->e_lj), ls);
831     CU_RET_ERR(stat, "cudaMemsetAsync on e_lj falied");
832     stat = cudaMemsetAsync(adat->e_el, 0, sizeof(*adat->e_el), ls);
833     CU_RET_ERR(stat, "cudaMemsetAsync on e_el falied");
834 }
835
836 void nbnxn_cuda_clear_outputs(nbnxn_cuda_ptr_t cu_nb, int flags)
837 {
838     nbnxn_cuda_clear_f(cu_nb, cu_nb->atdat->natoms);
839     /* clear shift force array and energies if the outputs were
840        used in the current step */
841     if (flags & GMX_FORCE_VIRIAL)
842     {
843         nbnxn_cuda_clear_e_fshift(cu_nb);
844     }
845 }
846
847 void nbnxn_cuda_init_atomdata(nbnxn_cuda_ptr_t        cu_nb,
848                               const nbnxn_atomdata_t *nbat)
849 {
850     cudaError_t    stat;
851     int            nalloc, natoms;
852     bool           realloced;
853     bool           bDoTime   = cu_nb->bDoTime;
854     cu_timers_t   *timers    = cu_nb->timers;
855     cu_atomdata_t *d_atdat   = cu_nb->atdat;
856     cudaStream_t   ls        = cu_nb->stream[eintLocal];
857
858     natoms    = nbat->natoms;
859     realloced = false;
860
861     if (bDoTime)
862     {
863         /* time async copy */
864         stat = cudaEventRecord(timers->start_atdat, ls);
865         CU_RET_ERR(stat, "cudaEventRecord failed");
866     }
867
868     /* need to reallocate if we have to copy more atoms than the amount of space
869        available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
870     if (natoms > d_atdat->nalloc)
871     {
872         nalloc = over_alloc_small(natoms);
873
874         /* free up first if the arrays have already been initialized */
875         if (d_atdat->nalloc != -1)
876         {
877             cu_free_buffered(d_atdat->f, &d_atdat->natoms, &d_atdat->nalloc);
878             cu_free_buffered(d_atdat->xq);
879             cu_free_buffered(d_atdat->atom_types);
880         }
881
882         stat = cudaMalloc((void **)&d_atdat->f, nalloc*sizeof(*d_atdat->f));
883         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->f");
884         stat = cudaMalloc((void **)&d_atdat->xq, nalloc*sizeof(*d_atdat->xq));
885         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->xq");
886
887         stat = cudaMalloc((void **)&d_atdat->atom_types, nalloc*sizeof(*d_atdat->atom_types));
888         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->atom_types");
889
890         d_atdat->nalloc = nalloc;
891         realloced       = true;
892     }
893
894     d_atdat->natoms       = natoms;
895     d_atdat->natoms_local = nbat->natoms_local;
896
897     /* need to clear GPU f output if realloc happened */
898     if (realloced)
899     {
900         nbnxn_cuda_clear_f(cu_nb, nalloc);
901     }
902
903     cu_copy_H2D_async(d_atdat->atom_types, nbat->type,
904                       natoms*sizeof(*d_atdat->atom_types), ls);
905
906     if (bDoTime)
907     {
908         stat = cudaEventRecord(timers->stop_atdat, ls);
909         CU_RET_ERR(stat, "cudaEventRecord failed");
910     }
911 }
912
913 void nbnxn_cuda_free(nbnxn_cuda_ptr_t cu_nb)
914 {
915     cudaError_t      stat;
916     cu_atomdata_t   *atdat;
917     cu_nbparam_t    *nbparam;
918     cu_plist_t      *plist, *plist_nl;
919     cu_timers_t     *timers;
920
921     if (cu_nb == NULL)
922     {
923         return;
924     }
925
926     atdat       = cu_nb->atdat;
927     nbparam     = cu_nb->nbparam;
928     plist       = cu_nb->plist[eintLocal];
929     plist_nl    = cu_nb->plist[eintNonlocal];
930     timers      = cu_nb->timers;
931
932     if (nbparam->eeltype == eelCuEWALD_TAB || nbparam->eeltype == eelCuEWALD_TAB_TWIN)
933     {
934
935 #ifdef TEXOBJ_SUPPORTED
936         /* Only device CC >= 3.0 (Kepler and later) support texture objects */
937         if (cu_nb->dev_info->prop.major >= 3)
938         {
939             stat = cudaDestroyTextureObject(nbparam->coulomb_tab_texobj);
940             CU_RET_ERR(stat, "cudaDestroyTextureObject on coulomb_tab_texobj failed");
941         }
942         else
943 #endif
944         {
945             stat = cudaUnbindTexture(nbnxn_cuda_get_coulomb_tab_texref());
946             CU_RET_ERR(stat, "cudaUnbindTexture on coulomb_tab_texref failed");
947         }
948         cu_free_buffered(nbparam->coulomb_tab, &nbparam->coulomb_tab_size);
949     }
950
951     stat = cudaEventDestroy(cu_nb->nonlocal_done);
952     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->nonlocal_done");
953     stat = cudaEventDestroy(cu_nb->misc_ops_done);
954     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->misc_ops_done");
955
956     if (cu_nb->bDoTime)
957     {
958         stat = cudaEventDestroy(timers->start_atdat);
959         CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_atdat");
960         stat = cudaEventDestroy(timers->stop_atdat);
961         CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_atdat");
962
963         /* The non-local counters/stream (second in the array) are needed only with DD. */
964         for (int i = 0; i <= (cu_nb->bUseTwoStreams ? 1 : 0); i++)
965         {
966             stat = cudaEventDestroy(timers->start_nb_k[i]);
967             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_k");
968             stat = cudaEventDestroy(timers->stop_nb_k[i]);
969             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_k");
970
971             stat = cudaEventDestroy(timers->start_pl_h2d[i]);
972             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_pl_h2d");
973             stat = cudaEventDestroy(timers->stop_pl_h2d[i]);
974             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_pl_h2d");
975
976             stat = cudaStreamDestroy(cu_nb->stream[i]);
977             CU_RET_ERR(stat, "cudaStreamDestroy failed on stream");
978
979             stat = cudaEventDestroy(timers->start_nb_h2d[i]);
980             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_h2d");
981             stat = cudaEventDestroy(timers->stop_nb_h2d[i]);
982             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_h2d");
983
984             stat = cudaEventDestroy(timers->start_nb_d2h[i]);
985             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->start_nb_d2h");
986             stat = cudaEventDestroy(timers->stop_nb_d2h[i]);
987             CU_RET_ERR(stat, "cudaEventDestroy failed on timers->stop_nb_d2h");
988         }
989     }
990
991 #ifdef TEXOBJ_SUPPORTED
992     /* Only device CC >= 3.0 (Kepler and later) support texture objects */
993     if (cu_nb->dev_info->prop.major >= 3)
994     {
995         stat = cudaDestroyTextureObject(nbparam->nbfp_texobj);
996         CU_RET_ERR(stat, "cudaDestroyTextureObject on nbfp_texobj failed");
997     }
998     else
999 #endif
1000     {
1001         stat = cudaUnbindTexture(nbnxn_cuda_get_nbfp_texref());
1002         CU_RET_ERR(stat, "cudaUnbindTexture on nbfp_texref failed");
1003     }
1004     cu_free_buffered(nbparam->nbfp);
1005
1006     if (nbparam->vdwtype == evdwCuEWALDGEOM || nbparam->vdwtype == evdwCuEWALDLB)
1007     {
1008 #ifdef TEXOBJ_SUPPORTED
1009         /* Only device CC >= 3.0 (Kepler and later) support texture objects */
1010         if (cu_nb->dev_info->prop.major >= 3)
1011         {
1012             stat = cudaDestroyTextureObject(nbparam->nbfp_comb_texobj);
1013             CU_RET_ERR(stat, "cudaDestroyTextureObject on nbfp_comb_texobj failed");
1014         }
1015         else
1016 #endif
1017         {
1018             stat = cudaUnbindTexture(nbnxn_cuda_get_nbfp_comb_texref());
1019             CU_RET_ERR(stat, "cudaUnbindTexture on nbfp_comb_texref failed");
1020         }
1021         cu_free_buffered(nbparam->nbfp_comb);
1022     }
1023
1024     stat = cudaFree(atdat->shift_vec);
1025     CU_RET_ERR(stat, "cudaFree failed on atdat->shift_vec");
1026     stat = cudaFree(atdat->fshift);
1027     CU_RET_ERR(stat, "cudaFree failed on atdat->fshift");
1028
1029     stat = cudaFree(atdat->e_lj);
1030     CU_RET_ERR(stat, "cudaFree failed on atdat->e_lj");
1031     stat = cudaFree(atdat->e_el);
1032     CU_RET_ERR(stat, "cudaFree failed on atdat->e_el");
1033
1034     cu_free_buffered(atdat->f, &atdat->natoms, &atdat->nalloc);
1035     cu_free_buffered(atdat->xq);
1036     cu_free_buffered(atdat->atom_types, &atdat->ntypes);
1037
1038     cu_free_buffered(plist->sci, &plist->nsci, &plist->sci_nalloc);
1039     cu_free_buffered(plist->cj4, &plist->ncj4, &plist->cj4_nalloc);
1040     cu_free_buffered(plist->excl, &plist->nexcl, &plist->excl_nalloc);
1041     if (cu_nb->bUseTwoStreams)
1042     {
1043         cu_free_buffered(plist_nl->sci, &plist_nl->nsci, &plist_nl->sci_nalloc);
1044         cu_free_buffered(plist_nl->cj4, &plist_nl->ncj4, &plist_nl->cj4_nalloc);
1045         cu_free_buffered(plist_nl->excl, &plist_nl->nexcl, &plist->excl_nalloc);
1046     }
1047
1048     sfree(atdat);
1049     sfree(nbparam);
1050     sfree(plist);
1051     if (cu_nb->bUseTwoStreams)
1052     {
1053         sfree(plist_nl);
1054     }
1055     sfree(timers);
1056     sfree(cu_nb->timings);
1057     sfree(cu_nb);
1058
1059     if (debug)
1060     {
1061         fprintf(debug, "Cleaned up CUDA data structures.\n");
1062     }
1063 }
1064
1065 void cu_synchstream_atdat(nbnxn_cuda_ptr_t cu_nb, int iloc)
1066 {
1067     cudaError_t  stat;
1068     cudaStream_t stream = cu_nb->stream[iloc];
1069
1070     stat = cudaStreamWaitEvent(stream, cu_nb->timers->stop_atdat, 0);
1071     CU_RET_ERR(stat, "cudaStreamWaitEvent failed");
1072 }
1073
1074 wallclock_gpu_t * nbnxn_cuda_get_timings(nbnxn_cuda_ptr_t cu_nb)
1075 {
1076     return (cu_nb != NULL && cu_nb->bDoTime) ? cu_nb->timings : NULL;
1077 }
1078
1079 void nbnxn_cuda_reset_timings(nbnxn_cuda_ptr_t cu_nb)
1080 {
1081     if (cu_nb->bDoTime)
1082     {
1083         init_timings(cu_nb->timings);
1084     }
1085 }
1086
1087 int nbnxn_cuda_min_ci_balanced(nbnxn_cuda_ptr_t cu_nb)
1088 {
1089     return cu_nb != NULL ?
1090            gpu_min_ci_balanced_factor*cu_nb->dev_info->prop.multiProcessorCount : 0;
1091
1092 }
1093
1094 gmx_bool nbnxn_cuda_is_kernel_ewald_analytical(const nbnxn_cuda_ptr_t cu_nb)
1095 {
1096     return ((cu_nb->nbparam->eeltype == eelCuEWALD_ANA) ||
1097             (cu_nb->nbparam->eeltype == eelCuEWALD_ANA_TWIN));
1098 }