Fix compiler warnings in OCL
[alexxy/gromacs.git] / src / gromacs / mdlib / nbnxn_ocl / nbnxn_ocl_data_mgmt.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013,2014,2015,2016,2017,2018, 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 /*! \internal \file
36  *  \brief Define OpenCL implementation of nbnxn_gpu_data_mgmt.h
37  *
38  *  \author Anca Hamuraru <anca@streamcomputing.eu>
39  *  \author Dimitrios Karkoulis <dimitris.karkoulis@gmail.com>
40  *  \author Teemu Virolainen <teemu@streamcomputing.eu>
41  *  \author Szilárd Páll <pall.szilard@gmail.com>
42  */
43 #include "gmxpre.h"
44
45 #include <assert.h>
46 #include <stdarg.h>
47 #include <stdio.h>
48 #include <stdlib.h>
49 #include <string.h>
50
51 #include <cmath>
52
53 #include "gromacs/gpu_utils/gpu_utils.h"
54 #include "gromacs/gpu_utils/oclutils.h"
55 #include "gromacs/hardware/gpu_hw_info.h"
56 #include "gromacs/math/vectypes.h"
57 #include "gromacs/mdlib/force_flags.h"
58 #include "gromacs/mdlib/nb_verlet.h"
59 #include "gromacs/mdlib/nbnxn_consts.h"
60 #include "gromacs/mdlib/nbnxn_gpu.h"
61 #include "gromacs/mdlib/nbnxn_gpu_common_utils.h"
62 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
63 #include "gromacs/mdlib/nbnxn_gpu_jit_support.h"
64 #include "gromacs/mdtypes/interaction_const.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/timing/gpu_timing.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/gmxassert.h"
71 #include "gromacs/utility/real.h"
72 #include "gromacs/utility/smalloc.h"
73
74 #include "nbnxn_ocl_internal.h"
75 #include "nbnxn_ocl_types.h"
76
77 /*! \brief This parameter should be determined heuristically from the
78  * kernel execution times
79  *
80  * This value is best for small systems on a single AMD Radeon R9 290X
81  * (and about 5% faster than 40, which is the default for CUDA
82  * devices). Larger simulation systems were quite insensitive to the
83  * value of this parameter.
84  */
85 static unsigned int gpu_min_ci_balanced_factor = 50;
86
87
88 /*! \brief Returns true if LJ combination rules are used in the non-bonded kernels.
89  *
90  * Full doc in nbnxn_ocl_internal.h */
91 bool useLjCombRule(int vdwType)
92 {
93     return (vdwType == evdwOclCUTCOMBGEOM ||
94             vdwType == evdwOclCUTCOMBLB);
95 }
96
97 /*! \brief Tabulates the Ewald Coulomb force and initializes the size/scale
98  * and the table GPU array.
99  *
100  * If called with an already allocated table, it just re-uploads the
101  * table.
102  */
103 static void init_ewald_coulomb_force_table(const interaction_const_t       *ic,
104                                            cl_nbparam_t                    *nbp,
105                                            const gmx_device_runtime_data_t *runData)
106 {
107     cl_mem       coul_tab;
108
109     cl_int       cl_error;
110
111     if (nbp->coulomb_tab_climg2d != nullptr)
112     {
113         freeDeviceBuffer(&(nbp->coulomb_tab_climg2d));
114     }
115
116     /* Switched from using textures to using buffers */
117     // TODO: decide which alternative is most efficient - textures or buffers.
118     /*
119        cl_image_format array_format;
120
121        array_format.image_channel_data_type = CL_FLOAT;
122        array_format.image_channel_order     = CL_R;
123
124        coul_tab = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
125        &array_format, tabsize, 1, 0, ftmp, &cl_error);
126      */
127
128     coul_tab = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, ic->tabq_size*sizeof(cl_float), ic->tabq_coul_F, &cl_error);
129     assert(cl_error == CL_SUCCESS);
130     // TODO: handle errors, check clCreateBuffer flags
131
132     nbp->coulomb_tab_climg2d  = coul_tab;
133     nbp->coulomb_tab_scale    = ic->tabq_scale;
134 }
135
136
137 /*! \brief Initializes the atomdata structure first time, it only gets filled at
138     pair-search.
139  */
140 static void init_atomdata_first(cl_atomdata_t *ad, int ntypes, gmx_device_runtime_data_t *runData)
141 {
142     cl_int cl_error;
143
144     ad->ntypes  = ntypes;
145
146     /* An element of the shift_vec device buffer has the same size as one element
147        of the host side shift_vec buffer. */
148     ad->shift_vec_elem_size = sizeof(*nbnxn_atomdata_t::shift_vec);
149
150     // TODO: handle errors, check clCreateBuffer flags
151     ad->shift_vec = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, SHIFTS * ad->shift_vec_elem_size, nullptr, &cl_error);
152     assert(cl_error == CL_SUCCESS);
153     ad->bShiftVecUploaded = false;
154
155     /* An element of the fshift device buffer has the same size as one element
156        of the host side fshift buffer. */
157     ad->fshift_elem_size = sizeof(*cl_nb_staging_t::fshift);
158
159     ad->fshift = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, SHIFTS * ad->fshift_elem_size, nullptr, &cl_error);
160     assert(cl_error == CL_SUCCESS);
161     // TODO: handle errors, check clCreateBuffer flags
162
163     ad->e_lj = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, sizeof(float), nullptr, &cl_error);
164     assert(cl_error == CL_SUCCESS);
165     // TODO: handle errors, check clCreateBuffer flags
166
167     ad->e_el = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, sizeof(float), nullptr, &cl_error);
168     assert(cl_error == CL_SUCCESS);
169     // TODO: handle errors, check clCreateBuffer flags
170
171     /* initialize to nullptr pointers to data that is not allocated here and will
172        need reallocation in nbnxn_gpu_init_atomdata */
173     ad->xq = nullptr;
174     ad->f  = nullptr;
175
176     /* size -1 indicates that the respective array hasn't been initialized yet */
177     ad->natoms = -1;
178     ad->nalloc = -1;
179 }
180
181 /*! \brief Copies all parameters related to the cut-off from ic to nbp
182  */
183 static void set_cutoff_parameters(cl_nbparam_t              *nbp,
184                                   const interaction_const_t *ic,
185                                   const NbnxnListParameters *listParams)
186 {
187     nbp->ewald_beta        = ic->ewaldcoeff_q;
188     nbp->sh_ewald          = ic->sh_ewald;
189     nbp->epsfac            = ic->epsfac;
190     nbp->two_k_rf          = 2.0 * ic->k_rf;
191     nbp->c_rf              = ic->c_rf;
192     nbp->rvdw_sq           = ic->rvdw * ic->rvdw;
193     nbp->rcoulomb_sq       = ic->rcoulomb * ic->rcoulomb;
194     nbp->rlistOuter_sq     = listParams->rlistOuter * listParams->rlistOuter;
195     nbp->rlistInner_sq     = listParams->rlistInner * listParams->rlistInner;
196     nbp->useDynamicPruning = listParams->useDynamicPruning;
197
198     nbp->sh_lj_ewald       = ic->sh_lj_ewald;
199     nbp->ewaldcoeff_lj     = ic->ewaldcoeff_lj;
200
201     nbp->rvdw_switch       = ic->rvdw_switch;
202     nbp->dispersion_shift  = ic->dispersion_shift;
203     nbp->repulsion_shift   = ic->repulsion_shift;
204     nbp->vdw_switch        = ic->vdw_switch;
205 }
206
207 /*! \brief Returns the kinds of electrostatics and Vdw OpenCL
208  *  kernels that will be used.
209  *
210  * Respectively, these values are from enum eelOcl and enum
211  * evdwOcl. */
212 static void
213 map_interaction_types_to_gpu_kernel_flavors(const interaction_const_t *ic,
214                                             int                        combRule,
215                                             int                       *gpu_eeltype,
216                                             int                       *gpu_vdwtype)
217 {
218     if (ic->vdwtype == evdwCUT)
219     {
220         switch (ic->vdw_modifier)
221         {
222             case eintmodNONE:
223             case eintmodPOTSHIFT:
224                 switch (combRule)
225                 {
226                     case ljcrNONE:
227                         *gpu_vdwtype = evdwOclCUT;
228                         break;
229                     case ljcrGEOM:
230                         *gpu_vdwtype = evdwOclCUTCOMBGEOM;
231                         break;
232                     case ljcrLB:
233                         *gpu_vdwtype = evdwOclCUTCOMBLB;
234                         break;
235                     default:
236                         gmx_incons("The requested LJ combination rule is not implemented in the OpenCL GPU accelerated kernels!");
237                         break;
238                 }
239                 break;
240             case eintmodFORCESWITCH:
241                 *gpu_vdwtype = evdwOclFSWITCH;
242                 break;
243             case eintmodPOTSWITCH:
244                 *gpu_vdwtype = evdwOclPSWITCH;
245                 break;
246             default:
247                 gmx_incons("The requested VdW interaction modifier is not implemented in the GPU accelerated kernels!");
248                 break;
249         }
250     }
251     else if (ic->vdwtype == evdwPME)
252     {
253         if (ic->ljpme_comb_rule == ljcrGEOM)
254         {
255             *gpu_vdwtype = evdwOclEWALDGEOM;
256         }
257         else
258         {
259             *gpu_vdwtype = evdwOclEWALDLB;
260         }
261     }
262     else
263     {
264         gmx_incons("The requested VdW type is not implemented in the GPU accelerated kernels!");
265     }
266
267     if (ic->eeltype == eelCUT)
268     {
269         *gpu_eeltype = eelOclCUT;
270     }
271     else if (EEL_RF(ic->eeltype))
272     {
273         *gpu_eeltype = eelOclRF;
274     }
275     else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
276     {
277         /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
278         *gpu_eeltype = nbnxn_gpu_pick_ewald_kernel_type(false);
279     }
280     else
281     {
282         /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
283         gmx_incons("The requested electrostatics type is not implemented in the GPU accelerated kernels!");
284     }
285 }
286
287 /*! \brief Initializes the nonbonded parameter data structure.
288  */
289 static void init_nbparam(cl_nbparam_t                    *nbp,
290                          const interaction_const_t       *ic,
291                          const NbnxnListParameters       *listParams,
292                          const nbnxn_atomdata_t          *nbat,
293                          const gmx_device_runtime_data_t *runData)
294 {
295     int         ntypes, nnbfp, nnbfp_comb;
296     cl_int      cl_error;
297
298
299     ntypes = nbat->ntype;
300
301     set_cutoff_parameters(nbp, ic, listParams);
302
303     map_interaction_types_to_gpu_kernel_flavors(ic,
304                                                 nbat->comb_rule,
305                                                 &(nbp->eeltype),
306                                                 &(nbp->vdwtype));
307
308     if (ic->vdwtype == evdwPME)
309     {
310         if (ic->ljpme_comb_rule == ljcrGEOM)
311         {
312             assert(nbat->comb_rule == ljcrGEOM);
313         }
314         else
315         {
316             assert(nbat->comb_rule == ljcrLB);
317         }
318     }
319     /* generate table for PME */
320     nbp->coulomb_tab_climg2d = nullptr;
321     if (nbp->eeltype == eelOclEWALD_TAB || nbp->eeltype == eelOclEWALD_TAB_TWIN)
322     {
323         init_ewald_coulomb_force_table(ic, nbp, runData);
324     }
325     else
326     // TODO: improvement needed.
327     // The image2d is created here even if eeltype is not eelCuEWALD_TAB or eelCuEWALD_TAB_TWIN because the OpenCL kernels
328     // don't accept nullptr values for image2D parameters.
329     {
330         /* Switched from using textures to using buffers */
331         // TODO: decide which alternative is most efficient - textures or buffers.
332         /*
333            cl_image_format array_format;
334
335            array_format.image_channel_data_type = CL_FLOAT;
336            array_format.image_channel_order     = CL_R;
337
338            nbp->coulomb_tab_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
339             &array_format, 1, 1, 0, nullptr, &cl_error);
340          */
341
342         nbp->coulomb_tab_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY, sizeof(cl_float), nullptr, &cl_error);
343         // TODO: handle errors
344     }
345
346     nnbfp      = 2*ntypes*ntypes;
347     nnbfp_comb = 2*ntypes;
348
349     {
350         /* Switched from using textures to using buffers */
351         // TODO: decide which alternative is most efficient - textures or buffers.
352         /*
353            cl_image_format array_format;
354
355            array_format.image_channel_data_type = CL_FLOAT;
356            array_format.image_channel_order     = CL_R;
357
358            nbp->nbfp_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
359             &array_format, nnbfp, 1, 0, nbat->nbfp, &cl_error);
360          */
361
362         nbp->nbfp_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, nnbfp*sizeof(cl_float), nbat->nbfp, &cl_error);
363         assert(cl_error == CL_SUCCESS);
364         // TODO: handle errors
365
366         if (ic->vdwtype == evdwPME)
367         {
368             /* Switched from using textures to using buffers */
369             // TODO: decide which alternative is most efficient - textures or buffers.
370             /*  nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
371                 &array_format, nnbfp_comb, 1, 0, nbat->nbfp_comb, &cl_error);*/
372             nbp->nbfp_comb_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, nnbfp_comb*sizeof(cl_float), nbat->nbfp_comb, &cl_error);
373
374
375             assert(cl_error == CL_SUCCESS);
376             // TODO: handle errors
377         }
378         else
379         {
380             // TODO: improvement needed.
381             // The image2d is created here even if vdwtype is not evdwPME because the OpenCL kernels
382             // don't accept nullptr values for image2D parameters.
383             /* Switched from using textures to using buffers */
384             // TODO: decide which alternative is most efficient - textures or buffers.
385             /* nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
386                 &array_format, 1, 1, 0, nullptr, &cl_error);*/
387             nbp->nbfp_comb_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY, sizeof(cl_float), nullptr, &cl_error);
388
389
390             assert(cl_error == CL_SUCCESS);
391             // TODO: handle errors
392         }
393     }
394 }
395
396 //! This function is documented in the header file
397 void nbnxn_gpu_pme_loadbal_update_param(const nonbonded_verlet_t    *nbv,
398                                         const interaction_const_t   *ic,
399                                         const NbnxnListParameters   *listParams)
400 {
401     if (!nbv || nbv->grp[0].kernel_type != nbnxnk8x8x8_GPU)
402     {
403         return;
404     }
405     gmx_nbnxn_ocl_t    *nb  = nbv->gpu_nbv;
406     cl_nbparam_t       *nbp = nb->nbparam;
407
408     set_cutoff_parameters(nbp, ic, listParams);
409
410     nbp->eeltype = nbnxn_gpu_pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw);
411
412     init_ewald_coulomb_force_table(ic, nb->nbparam, nb->dev_rundata);
413 }
414
415 /*! \brief Initializes the pair list data structure.
416  */
417 static void init_plist(cl_plist_t *pl)
418 {
419     /* initialize to nullptr pointers to data that is not allocated here and will
420        need reallocation in nbnxn_gpu_init_pairlist */
421     pl->sci     = nullptr;
422     pl->cj4     = nullptr;
423     pl->imask   = nullptr;
424     pl->excl    = nullptr;
425
426     /* size -1 indicates that the respective array hasn't been initialized yet */
427     pl->na_c           = -1;
428     pl->nsci           = -1;
429     pl->sci_nalloc     = -1;
430     pl->ncj4           = -1;
431     pl->cj4_nalloc     = -1;
432     pl->nimask         = -1;
433     pl->imask_nalloc   = -1;
434     pl->nexcl          = -1;
435     pl->excl_nalloc    = -1;
436     pl->haveFreshList  = false;
437 }
438
439 /*! \brief Initializes the timer data structure.
440  */
441 static void init_timers(cl_timers_t *t,
442                         bool         bUseTwoStreams)
443 {
444     for (int i = 0; i <= (bUseTwoStreams ? 1 : 0); i++)
445     {
446         t->didPairlistH2D[i]  = false;
447         t->didPrune[i]        = false;
448         t->didRollingPrune[i] = false;
449     }
450 }
451
452 /*! \brief Initializes the timings data structure.
453  */
454 static void init_timings(gmx_wallclock_gpu_nbnxn_t *t)
455 {
456     int i, j;
457
458     t->nb_h2d_t = 0.0;
459     t->nb_d2h_t = 0.0;
460     t->nb_c     = 0;
461     t->pl_h2d_t = 0.0;
462     t->pl_h2d_c = 0;
463     for (i = 0; i < 2; i++)
464     {
465         for (j = 0; j < 2; j++)
466         {
467             t->ktime[i][j].t = 0.0;
468             t->ktime[i][j].c = 0;
469         }
470     }
471
472     t->pruneTime.c        = 0;
473     t->pruneTime.t        = 0.0;
474     t->dynamicPruneTime.c = 0;
475     t->dynamicPruneTime.t = 0.0;
476 }
477
478 /*! \brief Creates context for OpenCL GPU given by \p mygpu
479  *
480  * A fatal error results if creation fails.
481  *
482  * \param[inout] runtimeData runtime data including program and context
483  * \param[in]    devInfo     device info struct
484  * \param[in]    rank        MPI rank (for error reporting)
485  */
486 static void
487 nbnxn_gpu_create_context(gmx_device_runtime_data_t *runtimeData,
488                          const gmx_device_info_t   *devInfo,
489                          int                        rank)
490 {
491     cl_context_properties     context_properties[3];
492     cl_platform_id            platform_id;
493     cl_device_id              device_id;
494     cl_context                context;
495     cl_int                    cl_error;
496
497     assert(runtimeData != nullptr);
498     assert(devInfo != nullptr);
499
500     platform_id      = devInfo->ocl_gpu_id.ocl_platform_id;
501     device_id        = devInfo->ocl_gpu_id.ocl_device_id;
502
503     context_properties[0] = CL_CONTEXT_PLATFORM;
504     context_properties[1] = (cl_context_properties) platform_id;
505     context_properties[2] = 0; /* Terminates the list of properties */
506
507     context = clCreateContext(context_properties, 1, &device_id, nullptr, nullptr, &cl_error);
508     if (CL_SUCCESS != cl_error)
509     {
510         gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s:\n OpenCL error %d: %s",
511                   rank,
512                   devInfo->device_name,
513                   cl_error, ocl_get_error_string(cl_error).c_str());
514         return;
515     }
516
517     runtimeData->context = context;
518 }
519
520 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
521 static cl_kernel nbnxn_gpu_create_kernel(gmx_nbnxn_ocl_t *nb,
522                                          const char      *kernel_name)
523 {
524     cl_kernel kernel;
525     cl_int    cl_error;
526
527     kernel = clCreateKernel(nb->dev_rundata->program, kernel_name, &cl_error);
528     if (CL_SUCCESS != cl_error)
529     {
530         gmx_fatal(FARGS, "Failed to create kernel '%s' for GPU #%s: OpenCL error %d",
531                   kernel_name,
532                   nb->dev_info->device_name,
533                   cl_error);
534     }
535
536     return kernel;
537 }
538
539 /*! \brief Clears nonbonded shift force output array and energy outputs on the GPU.
540  */
541 static void
542 nbnxn_ocl_clear_e_fshift(gmx_nbnxn_ocl_t *nb)
543 {
544
545     cl_int               cl_error;
546     cl_atomdata_t *      adat     = nb->atdat;
547     cl_command_queue     ls       = nb->stream[eintLocal];
548
549     size_t               local_work_size[3]   = {1, 1, 1};
550     size_t               global_work_size[3]  = {1, 1, 1};
551
552     cl_int               shifts   = SHIFTS*3;
553
554     cl_int               arg_no;
555
556     cl_kernel            zero_e_fshift = nb->kernel_zero_e_fshift;
557
558     local_work_size[0]   = 64;
559     // Round the total number of threads up from the array size
560     global_work_size[0]  = ((shifts + local_work_size[0] - 1)/local_work_size[0])*local_work_size[0];
561
562     arg_no    = 0;
563     cl_error  = clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->fshift));
564     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->e_lj));
565     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->e_el));
566     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_uint), &shifts);
567     GMX_ASSERT(cl_error == CL_SUCCESS, ocl_get_error_string(cl_error).c_str());
568
569     cl_error = clEnqueueNDRangeKernel(ls, zero_e_fshift, 3, nullptr, global_work_size, local_work_size, 0, nullptr, nullptr);
570     GMX_ASSERT(cl_error == CL_SUCCESS, ocl_get_error_string(cl_error).c_str());
571 }
572
573 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
574 static void nbnxn_gpu_init_kernels(gmx_nbnxn_ocl_t *nb)
575 {
576     /* Init to 0 main kernel arrays */
577     /* They will be later on initialized in select_nbnxn_kernel */
578     // TODO: consider always creating all variants of the kernels here so that there is no
579     // need for late call to clCreateKernel -- if that gives any advantage?
580     memset(nb->kernel_ener_noprune_ptr, 0, sizeof(nb->kernel_ener_noprune_ptr));
581     memset(nb->kernel_ener_prune_ptr, 0, sizeof(nb->kernel_ener_prune_ptr));
582     memset(nb->kernel_noener_noprune_ptr, 0, sizeof(nb->kernel_noener_noprune_ptr));
583     memset(nb->kernel_noener_prune_ptr, 0, sizeof(nb->kernel_noener_prune_ptr));
584
585     /* Init pruning kernels
586      *
587      * TODO: we could avoid creating kernels if dynamic pruning is turned off,
588      * but ATM that depends on force flags not passed into the initialization.
589      */
590     nb->kernel_pruneonly[epruneFirst]   = nbnxn_gpu_create_kernel(nb, "nbnxn_kernel_prune_opencl");
591     nb->kernel_pruneonly[epruneRolling] = nbnxn_gpu_create_kernel(nb, "nbnxn_kernel_prune_rolling_opencl");
592
593     /* Init auxiliary kernels */
594     nb->kernel_memset_f      = nbnxn_gpu_create_kernel(nb, "memset_f");
595     nb->kernel_memset_f2     = nbnxn_gpu_create_kernel(nb, "memset_f2");
596     nb->kernel_memset_f3     = nbnxn_gpu_create_kernel(nb, "memset_f3");
597     nb->kernel_zero_e_fshift = nbnxn_gpu_create_kernel(nb, "zero_e_fshift");
598 }
599
600 /*! \brief Initializes simulation constant data.
601  *
602  *  Initializes members of the atomdata and nbparam structs and
603  *  clears e/fshift output buffers.
604  */
605 static void nbnxn_ocl_init_const(gmx_nbnxn_ocl_t                *nb,
606                                  const interaction_const_t      *ic,
607                                  const NbnxnListParameters      *listParams,
608                                  const nbnxn_atomdata_t         *nbat)
609 {
610     init_atomdata_first(nb->atdat, nbat->ntype, nb->dev_rundata);
611     init_nbparam(nb->nbparam, ic, listParams, nbat, nb->dev_rundata);
612 }
613
614
615 //! This function is documented in the header file
616 void nbnxn_gpu_init(gmx_nbnxn_ocl_t          **p_nb,
617                     const gmx_device_info_t   *deviceInfo,
618                     const interaction_const_t *ic,
619                     const NbnxnListParameters *listParams,
620                     const nbnxn_atomdata_t    *nbat,
621                     int                        rank,
622                     gmx_bool                   bLocalAndNonlocal)
623 {
624     gmx_nbnxn_ocl_t            *nb;
625     cl_int                      cl_error;
626     cl_command_queue_properties queue_properties;
627
628     assert(ic);
629
630     if (p_nb == nullptr)
631     {
632         return;
633     }
634
635     snew(nb, 1);
636     snew(nb->atdat, 1);
637     snew(nb->nbparam, 1);
638     snew(nb->plist[eintLocal], 1);
639     if (bLocalAndNonlocal)
640     {
641         snew(nb->plist[eintNonlocal], 1);
642     }
643
644     nb->bUseTwoStreams = bLocalAndNonlocal;
645
646     nb->timers = new cl_timers_t();
647     snew(nb->timings, 1);
648
649     /* set device info, just point it to the right GPU among the detected ones */
650     nb->dev_info = deviceInfo;
651     snew(nb->dev_rundata, 1);
652
653     /* init nbst */
654     pmalloc((void**)&nb->nbst.e_lj, sizeof(*nb->nbst.e_lj));
655     pmalloc((void**)&nb->nbst.e_el, sizeof(*nb->nbst.e_el));
656     pmalloc((void**)&nb->nbst.fshift, SHIFTS * sizeof(*nb->nbst.fshift));
657
658     init_plist(nb->plist[eintLocal]);
659
660     /* OpenCL timing disabled if GMX_DISABLE_GPU_TIMING is defined. */
661     nb->bDoTime = (getenv("GMX_DISABLE_GPU_TIMING") == nullptr);
662
663     /* Create queues only after bDoTime has been initialized */
664     if (nb->bDoTime)
665     {
666         queue_properties = CL_QUEUE_PROFILING_ENABLE;
667     }
668     else
669     {
670         queue_properties = 0;
671     }
672
673     nbnxn_gpu_create_context(nb->dev_rundata, nb->dev_info, rank);
674
675     /* local/non-local GPU streams */
676     nb->stream[eintLocal] = clCreateCommandQueue(nb->dev_rundata->context, nb->dev_info->ocl_gpu_id.ocl_device_id, queue_properties, &cl_error);
677     if (CL_SUCCESS != cl_error)
678     {
679         gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
680                   rank,
681                   nb->dev_info->device_name,
682                   cl_error);
683         return;
684     }
685
686     if (nb->bUseTwoStreams)
687     {
688         init_plist(nb->plist[eintNonlocal]);
689
690         nb->stream[eintNonlocal] = clCreateCommandQueue(nb->dev_rundata->context, nb->dev_info->ocl_gpu_id.ocl_device_id, queue_properties, &cl_error);
691         if (CL_SUCCESS != cl_error)
692         {
693             gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
694                       rank,
695                       nb->dev_info->device_name,
696                       cl_error);
697             return;
698         }
699     }
700
701     if (nb->bDoTime)
702     {
703         init_timers(nb->timers, nb->bUseTwoStreams);
704         init_timings(nb->timings);
705     }
706
707     nbnxn_ocl_init_const(nb, ic, listParams, nbat);
708
709     /* Enable LJ param manual prefetch for AMD or if we request through env. var.
710      * TODO: decide about NVIDIA
711      */
712     nb->bPrefetchLjParam =
713         (getenv("GMX_OCL_DISABLE_I_PREFETCH") == nullptr) &&
714         ((nb->dev_info->vendor_e == OCL_VENDOR_AMD) || (getenv("GMX_OCL_ENABLE_I_PREFETCH") != nullptr));
715
716     /* NOTE: in CUDA we pick L1 cache configuration for the nbnxn kernels here,
717      * but sadly this is not supported in OpenCL (yet?). Consider adding it if
718      * it becomes supported.
719      */
720     nbnxn_gpu_compile_kernels(nb);
721     nbnxn_gpu_init_kernels(nb);
722
723     /* clear energy and shift force outputs */
724     nbnxn_ocl_clear_e_fshift(nb);
725
726     *p_nb = nb;
727
728     if (debug)
729     {
730         fprintf(debug, "Initialized OpenCL data structures.\n");
731     }
732 }
733
734 /*! \brief Clears the first natoms_clear elements of the GPU nonbonded force output array.
735  */
736 static void nbnxn_ocl_clear_f(gmx_nbnxn_ocl_t *nb, int natoms_clear)
737 {
738     if (natoms_clear == 0)
739     {
740         return;
741     }
742
743     cl_atomdata_t *      adat     = nb->atdat;
744     cl_command_queue     ls       = nb->stream[eintLocal];
745     cl_float             value    = 0.0f;
746
747     size_t               local_work_size[3]  = {1, 1, 1};
748     size_t               global_work_size[3] = {1, 1, 1};
749
750     cl_int               arg_no;
751
752     cl_kernel            memset_f = nb->kernel_memset_f;
753
754     cl_uint              natoms_flat = natoms_clear * (sizeof(rvec)/sizeof(real));
755
756     local_work_size[0]  = 64;
757     // Round the total number of threads up from the array size
758     global_work_size[0] = ((natoms_flat + local_work_size[0] - 1)/local_work_size[0])*local_work_size[0];
759
760
761     cl_int gmx_used_in_debug cl_error;
762     arg_no    = 0;
763     cl_error  = clSetKernelArg(memset_f, arg_no++, sizeof(cl_mem), &(adat->f));
764     cl_error |= clSetKernelArg(memset_f, arg_no++, sizeof(cl_float), &value);
765     cl_error |= clSetKernelArg(memset_f, arg_no++, sizeof(cl_uint), &natoms_flat);
766     assert(cl_error == CL_SUCCESS);
767
768     cl_error = clEnqueueNDRangeKernel(ls, memset_f, 3, nullptr, global_work_size, local_work_size, 0, nullptr, nullptr);
769     assert(cl_error == CL_SUCCESS);
770 }
771
772 //! This function is documented in the header file
773 void
774 nbnxn_gpu_clear_outputs(gmx_nbnxn_ocl_t   *nb,
775                         int                flags)
776 {
777     nbnxn_ocl_clear_f(nb, nb->atdat->natoms);
778     /* clear shift force array and energies if the outputs were
779        used in the current step */
780     if (flags & GMX_FORCE_VIRIAL)
781     {
782         nbnxn_ocl_clear_e_fshift(nb);
783     }
784
785     /* kick off buffer clearing kernel to ensure concurrency with constraints/update */
786     cl_int gmx_unused cl_error;
787     cl_error = clFlush(nb->stream[eintLocal]);
788     assert(CL_SUCCESS == cl_error);
789 }
790
791 //! This function is documented in the header file
792 void nbnxn_gpu_init_pairlist(gmx_nbnxn_ocl_t        *nb,
793                              const nbnxn_pairlist_t *h_plist,
794                              int                     iloc)
795 {
796     if (canSkipWork(nb, iloc))
797     {
798         return;
799     }
800
801     char             sbuf[STRLEN];
802     bool             bDoTime    = nb->bDoTime;
803     cl_command_queue stream     = nb->stream[iloc];
804     cl_plist_t      *d_plist    = nb->plist[iloc];
805
806     if (d_plist->na_c < 0)
807     {
808         d_plist->na_c = h_plist->na_ci;
809     }
810     else
811     {
812         if (d_plist->na_c != h_plist->na_ci)
813         {
814             sprintf(sbuf, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
815                     d_plist->na_c, h_plist->na_ci);
816             gmx_incons(sbuf);
817         }
818     }
819
820     if (bDoTime)
821     {
822         nb->timers->pl_h2d[iloc].openTimingRegion(stream);
823         nb->timers->didPairlistH2D[iloc] = true;
824     }
825
826     // TODO most of this function is same in CUDA and OpenCL, move into the header
827     Context context = nb->dev_rundata->context;
828
829     reallocateDeviceBuffer(&d_plist->sci, h_plist->nsci,
830                            &d_plist->nsci, &d_plist->sci_nalloc, context);
831     copyToDeviceBuffer(&d_plist->sci, h_plist->sci, 0, h_plist->nsci,
832                        stream, GpuApiCallBehavior::Async,
833                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
834
835     reallocateDeviceBuffer(&d_plist->cj4, h_plist->ncj4,
836                            &d_plist->ncj4, &d_plist->cj4_nalloc, context);
837     copyToDeviceBuffer(&d_plist->cj4, h_plist->cj4, 0, h_plist->ncj4,
838                        stream, GpuApiCallBehavior::Async,
839                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
840
841     reallocateDeviceBuffer(&d_plist->imask, h_plist->ncj4*c_nbnxnGpuClusterpairSplit,
842                            &d_plist->nimask, &d_plist->imask_nalloc, context);
843
844     reallocateDeviceBuffer(&d_plist->excl, h_plist->nexcl,
845                            &d_plist->nexcl, &d_plist->excl_nalloc, context);
846     copyToDeviceBuffer(&d_plist->excl, h_plist->excl, 0, h_plist->nexcl,
847                        stream, GpuApiCallBehavior::Async,
848                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
849
850     if (bDoTime)
851     {
852         nb->timers->pl_h2d[iloc].closeTimingRegion(stream);
853     }
854
855     /* need to prune the pair list during the next step */
856     d_plist->haveFreshList = true;
857 }
858
859 //! This function is documented in the header file
860 void nbnxn_gpu_upload_shiftvec(gmx_nbnxn_ocl_t        *nb,
861                                const nbnxn_atomdata_t *nbatom)
862 {
863     cl_atomdata_t   *adat  = nb->atdat;
864     cl_command_queue ls    = nb->stream[eintLocal];
865
866     /* only if we have a dynamic box */
867     if (nbatom->bDynamicBox || !adat->bShiftVecUploaded)
868     {
869         ocl_copy_H2D_async(adat->shift_vec, nbatom->shift_vec, 0,
870                            SHIFTS * adat->shift_vec_elem_size, ls, nullptr);
871         adat->bShiftVecUploaded = true;
872     }
873 }
874
875 //! This function is documented in the header file
876 void nbnxn_gpu_init_atomdata(gmx_nbnxn_ocl_t               *nb,
877                              const nbnxn_atomdata_t        *nbat)
878 {
879     cl_int           cl_error;
880     int              nalloc, natoms;
881     bool             realloced;
882     bool             bDoTime = nb->bDoTime;
883     cl_timers_t     *timers  = nb->timers;
884     cl_atomdata_t   *d_atdat = nb->atdat;
885     cl_command_queue ls      = nb->stream[eintLocal];
886
887     natoms    = nbat->natoms;
888     realloced = false;
889
890     if (bDoTime)
891     {
892         /* time async copy */
893         timers->atdat.openTimingRegion(ls);
894     }
895
896     /* need to reallocate if we have to copy more atoms than the amount of space
897        available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
898     if (natoms > d_atdat->nalloc)
899     {
900         nalloc = over_alloc_small(natoms);
901
902         /* free up first if the arrays have already been initialized */
903         if (d_atdat->nalloc != -1)
904         {
905             freeDeviceBuffer(&d_atdat->f);
906             freeDeviceBuffer(&d_atdat->xq);
907             freeDeviceBuffer(&d_atdat->lj_comb);
908             freeDeviceBuffer(&d_atdat->atom_types);
909         }
910
911         d_atdat->f_elem_size = sizeof(rvec);
912
913         // TODO: handle errors, check clCreateBuffer flags
914         d_atdat->f = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * d_atdat->f_elem_size, nullptr, &cl_error);
915         assert(CL_SUCCESS == cl_error);
916
917         // TODO: change the flag to read-only
918         d_atdat->xq = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(cl_float4), nullptr, &cl_error);
919         assert(CL_SUCCESS == cl_error);
920         // TODO: handle errors, check clCreateBuffer flags
921
922         if (useLjCombRule(nb->nbparam->vdwtype))
923         {
924             // TODO: change the flag to read-only
925             d_atdat->lj_comb = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(cl_float2), nullptr, &cl_error);
926             assert(CL_SUCCESS == cl_error);
927             // TODO: handle errors, check clCreateBuffer flags
928         }
929         else
930         {
931             // TODO: change the flag to read-only
932             d_atdat->atom_types = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(int), nullptr, &cl_error);
933             assert(CL_SUCCESS == cl_error);
934             // TODO: handle errors, check clCreateBuffer flags
935         }
936
937         d_atdat->nalloc = nalloc;
938         realloced       = true;
939     }
940
941     d_atdat->natoms       = natoms;
942     d_atdat->natoms_local = nbat->natoms_local;
943
944     /* need to clear GPU f output if realloc happened */
945     if (realloced)
946     {
947         nbnxn_ocl_clear_f(nb, nalloc);
948     }
949
950     if (useLjCombRule(nb->nbparam->vdwtype))
951     {
952         ocl_copy_H2D_async(d_atdat->lj_comb, nbat->lj_comb, 0,
953                            natoms*sizeof(cl_float2), ls, bDoTime ? timers->atdat.fetchNextEvent() : nullptr);
954     }
955     else
956     {
957         ocl_copy_H2D_async(d_atdat->atom_types, nbat->type, 0,
958                            natoms*sizeof(int), ls, bDoTime ? timers->atdat.fetchNextEvent() : nullptr);
959
960     }
961
962     if (bDoTime)
963     {
964         timers->atdat.closeTimingRegion(ls);
965     }
966
967     /* kick off the tasks enqueued above to ensure concurrency with the search */
968     cl_error = clFlush(ls);
969     assert(CL_SUCCESS == cl_error);
970 }
971
972 /*! \brief Releases an OpenCL kernel pointer */
973 static void free_kernel(cl_kernel *kernel_ptr)
974 {
975     cl_int gmx_unused cl_error;
976
977     assert(nullptr != kernel_ptr);
978
979     if (*kernel_ptr)
980     {
981         cl_error = clReleaseKernel(*kernel_ptr);
982         assert(cl_error == CL_SUCCESS);
983
984         *kernel_ptr = nullptr;
985     }
986 }
987
988 /*! \brief Releases a list of OpenCL kernel pointers */
989 static void free_kernels(cl_kernel *kernels, int count)
990 {
991     int i;
992
993     for (i = 0; i < count; i++)
994     {
995         free_kernel(kernels + i);
996     }
997 }
998
999 /*! \brief Free the OpenCL runtime data (context and program).
1000  *
1001  *  The function releases the OpenCL context and program assuciated with the
1002  *  device that the calling PP rank is running on.
1003  *
1004  *  \param runData [in]  porinter to the structure with runtime data.
1005  */
1006 static void free_gpu_device_runtime_data(gmx_device_runtime_data_t *runData)
1007 {
1008     if (runData == nullptr)
1009     {
1010         return;
1011     }
1012
1013     cl_int gmx_unused cl_error;
1014
1015     if (runData->context)
1016     {
1017         cl_error         = clReleaseContext(runData->context);
1018         runData->context = nullptr;
1019         assert(CL_SUCCESS == cl_error);
1020     }
1021
1022     if (runData->program)
1023     {
1024         cl_error         = clReleaseProgram(runData->program);
1025         runData->program = nullptr;
1026         assert(CL_SUCCESS == cl_error);
1027     }
1028
1029 }
1030
1031 //! This function is documented in the header file
1032 void nbnxn_gpu_free(gmx_nbnxn_ocl_t *nb)
1033 {
1034     if (nb == nullptr)
1035     {
1036         return;
1037     }
1038
1039     /* Free kernels */
1040     int kernel_count = sizeof(nb->kernel_ener_noprune_ptr) / sizeof(nb->kernel_ener_noprune_ptr[0][0]);
1041     free_kernels((cl_kernel*)nb->kernel_ener_noprune_ptr, kernel_count);
1042
1043     kernel_count = sizeof(nb->kernel_ener_prune_ptr) / sizeof(nb->kernel_ener_prune_ptr[0][0]);
1044     free_kernels((cl_kernel*)nb->kernel_ener_prune_ptr, kernel_count);
1045
1046     kernel_count = sizeof(nb->kernel_noener_noprune_ptr) / sizeof(nb->kernel_noener_noprune_ptr[0][0]);
1047     free_kernels((cl_kernel*)nb->kernel_noener_noprune_ptr, kernel_count);
1048
1049     kernel_count = sizeof(nb->kernel_noener_prune_ptr) / sizeof(nb->kernel_noener_prune_ptr[0][0]);
1050     free_kernels((cl_kernel*)nb->kernel_noener_prune_ptr, kernel_count);
1051
1052     free_kernel(&(nb->kernel_memset_f));
1053     free_kernel(&(nb->kernel_memset_f2));
1054     free_kernel(&(nb->kernel_memset_f3));
1055     free_kernel(&(nb->kernel_zero_e_fshift));
1056
1057     /* Free atdat */
1058     freeDeviceBuffer(&(nb->atdat->xq));
1059     freeDeviceBuffer(&(nb->atdat->f));
1060     freeDeviceBuffer(&(nb->atdat->e_lj));
1061     freeDeviceBuffer(&(nb->atdat->e_el));
1062     freeDeviceBuffer(&(nb->atdat->fshift));
1063     freeDeviceBuffer(&(nb->atdat->lj_comb));
1064     freeDeviceBuffer(&(nb->atdat->atom_types));
1065     freeDeviceBuffer(&(nb->atdat->shift_vec));
1066     sfree(nb->atdat);
1067
1068     /* Free nbparam */
1069     freeDeviceBuffer(&(nb->nbparam->nbfp_climg2d));
1070     freeDeviceBuffer(&(nb->nbparam->nbfp_comb_climg2d));
1071     freeDeviceBuffer(&(nb->nbparam->coulomb_tab_climg2d));
1072     sfree(nb->nbparam);
1073
1074     /* Free plist */
1075     auto *plist = nb->plist[eintLocal];
1076     freeDeviceBuffer(&plist->sci);
1077     freeDeviceBuffer(&plist->cj4);
1078     freeDeviceBuffer(&plist->imask);
1079     freeDeviceBuffer(&plist->excl);
1080     sfree(plist);
1081     if (nb->bUseTwoStreams)
1082     {
1083         auto *plist_nl = nb->plist[eintNonlocal];
1084         freeDeviceBuffer(&plist_nl->sci);
1085         freeDeviceBuffer(&plist_nl->cj4);
1086         freeDeviceBuffer(&plist_nl->imask);
1087         freeDeviceBuffer(&plist_nl->excl);
1088         sfree(plist_nl);
1089     }
1090
1091     /* Free nbst */
1092     pfree(nb->nbst.e_lj);
1093     nb->nbst.e_lj = nullptr;
1094
1095     pfree(nb->nbst.e_el);
1096     nb->nbst.e_el = nullptr;
1097
1098     pfree(nb->nbst.fshift);
1099     nb->nbst.fshift = nullptr;
1100
1101     /* Free command queues */
1102     clReleaseCommandQueue(nb->stream[eintLocal]);
1103     nb->stream[eintLocal] = nullptr;
1104     if (nb->bUseTwoStreams)
1105     {
1106         clReleaseCommandQueue(nb->stream[eintNonlocal]);
1107         nb->stream[eintNonlocal] = nullptr;
1108     }
1109     /* Free other events */
1110     if (nb->nonlocal_done)
1111     {
1112         clReleaseEvent(nb->nonlocal_done);
1113         nb->nonlocal_done = nullptr;
1114     }
1115     if (nb->misc_ops_and_local_H2D_done)
1116     {
1117         clReleaseEvent(nb->misc_ops_and_local_H2D_done);
1118         nb->misc_ops_and_local_H2D_done = nullptr;
1119     }
1120
1121     free_gpu_device_runtime_data(nb->dev_rundata);
1122     sfree(nb->dev_rundata);
1123
1124     /* Free timers and timings */
1125     delete nb->timers;
1126     sfree(nb->timings);
1127     sfree(nb);
1128
1129     if (debug)
1130     {
1131         fprintf(debug, "Cleaned up OpenCL data structures.\n");
1132     }
1133 }
1134
1135 //! This function is documented in the header file
1136 gmx_wallclock_gpu_nbnxn_t *nbnxn_gpu_get_timings(gmx_nbnxn_ocl_t *nb)
1137 {
1138     return (nb != nullptr && nb->bDoTime) ? nb->timings : nullptr;
1139 }
1140
1141 //! This function is documented in the header file
1142 void nbnxn_gpu_reset_timings(nonbonded_verlet_t* nbv)
1143 {
1144     if (nbv->gpu_nbv && nbv->gpu_nbv->bDoTime)
1145     {
1146         init_timings(nbv->gpu_nbv->timings);
1147     }
1148 }
1149
1150 //! This function is documented in the header file
1151 int nbnxn_gpu_min_ci_balanced(gmx_nbnxn_ocl_t *nb)
1152 {
1153     return nb != nullptr ?
1154            gpu_min_ci_balanced_factor * nb->dev_info->compute_units : 0;
1155 }
1156
1157 //! This function is documented in the header file
1158 gmx_bool nbnxn_gpu_is_kernel_ewald_analytical(const gmx_nbnxn_ocl_t *nb)
1159 {
1160     return ((nb->nbparam->eeltype == eelOclEWALD_ANA) ||
1161             (nb->nbparam->eeltype == eelOclEWALD_ANA_TWIN));
1162 }