2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
36 * \brief Define OpenCL implementation of nbnxn_gpu_data_mgmt.h
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>
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"
74 #include "nbnxn_ocl_internal.h"
75 #include "nbnxn_ocl_types.h"
77 /*! \brief This parameter should be determined heuristically from the
78 * kernel execution times
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.
85 static unsigned int gpu_min_ci_balanced_factor
= 50;
88 /*! \brief Returns true if LJ combination rules are used in the non-bonded kernels.
90 * Full doc in nbnxn_ocl_internal.h */
91 bool useLjCombRule(int vdwType
)
93 return (vdwType
== evdwOclCUTCOMBGEOM
||
94 vdwType
== evdwOclCUTCOMBLB
);
97 /*! \brief Free device buffers
99 * If the pointers to the size variables are NULL no resetting happens.
101 static void ocl_free_buffered(cl_mem d_ptr
, int *n
, int *nalloc
)
103 cl_int gmx_unused cl_error
;
107 cl_error
= clReleaseMemObject(d_ptr
);
108 assert(cl_error
== CL_SUCCESS
);
109 // TODO: handle errors
123 /*! \brief Reallocation device buffers
125 * Reallocation of the memory pointed by d_ptr and copying of the data from
126 * the location pointed by h_src host-side pointer is done. Allocation is
127 * buffered and therefore freeing is only needed if the previously allocated
128 * space is not enough.
129 * The H2D copy is launched in command queue s and can be done synchronously or
130 * asynchronously (the default is the latter).
131 * If copy_event is not NULL, on return it will contain an event object
132 * identifying the H2D copy. The event can further be used to queue a wait
133 * for this operation or to query profiling information.
134 * OpenCL equivalent of cu_realloc_buffered.
136 static void ocl_realloc_buffered(cl_mem
*d_dest
, void *h_src
,
138 int *curr_size
, int *curr_alloc_size
,
143 cl_event
*copy_event
= NULL
)
145 if (d_dest
== NULL
|| req_size
< 0)
150 /* reallocate only if the data does not fit = allocation size is smaller
151 than the current requested size */
152 if (req_size
> *curr_alloc_size
)
154 cl_int gmx_unused cl_error
;
156 /* only free if the array has already been initialized */
157 if (*curr_alloc_size
>= 0)
159 ocl_free_buffered(*d_dest
, curr_size
, curr_alloc_size
);
162 *curr_alloc_size
= over_alloc_large(req_size
);
164 *d_dest
= clCreateBuffer(context
, CL_MEM_READ_WRITE
, *curr_alloc_size
* type_size
, NULL
, &cl_error
);
165 assert(cl_error
== CL_SUCCESS
);
166 // TODO: handle errors, check clCreateBuffer flags
169 /* size could have changed without actual reallocation */
170 *curr_size
= req_size
;
172 /* upload to device */
177 ocl_copy_H2D_async(*d_dest
, h_src
, 0, *curr_size
* type_size
, s
, copy_event
);
181 ocl_copy_H2D_sync(*d_dest
, h_src
, 0, *curr_size
* type_size
, s
);
186 /*! \brief Releases the input OpenCL buffer */
187 static void free_ocl_buffer(cl_mem
*buffer
)
189 cl_int gmx_unused cl_error
;
191 assert(NULL
!= buffer
);
195 cl_error
= clReleaseMemObject(*buffer
);
196 assert(CL_SUCCESS
== cl_error
);
201 /*! \brief Tabulates the Ewald Coulomb force and initializes the size/scale
202 * and the table GPU array.
204 * If called with an already allocated table, it just re-uploads the
207 static void init_ewald_coulomb_force_table(const interaction_const_t
*ic
,
209 const gmx_device_runtime_data_t
*runData
)
215 if (nbp
->coulomb_tab_climg2d
!= NULL
)
217 free_ocl_buffer(&(nbp
->coulomb_tab_climg2d
));
220 /* Switched from using textures to using buffers */
221 // TODO: decide which alternative is most efficient - textures or buffers.
223 cl_image_format array_format;
225 array_format.image_channel_data_type = CL_FLOAT;
226 array_format.image_channel_order = CL_R;
228 coul_tab = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
229 &array_format, tabsize, 1, 0, ftmp, &cl_error);
232 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
);
233 assert(cl_error
== CL_SUCCESS
);
234 // TODO: handle errors, check clCreateBuffer flags
236 nbp
->coulomb_tab_climg2d
= coul_tab
;
237 nbp
->coulomb_tab_scale
= ic
->tabq_scale
;
241 /*! \brief Initializes the atomdata structure first time, it only gets filled at
244 static void init_atomdata_first(cl_atomdata_t
*ad
, int ntypes
, gmx_device_runtime_data_t
*runData
)
250 /* An element of the shift_vec device buffer has the same size as one element
251 of the host side shift_vec buffer. */
252 ad
->shift_vec_elem_size
= sizeof(*(((nbnxn_atomdata_t
*)0)->shift_vec
));
254 // TODO: handle errors, check clCreateBuffer flags
255 ad
->shift_vec
= clCreateBuffer(runData
->context
, CL_MEM_READ_WRITE
, SHIFTS
* ad
->shift_vec_elem_size
, NULL
, &cl_error
);
256 assert(cl_error
== CL_SUCCESS
);
257 ad
->bShiftVecUploaded
= false;
259 /* An element of the fshift device buffer has the same size as one element
260 of the host side fshift buffer. */
261 ad
->fshift_elem_size
= sizeof(*(((cl_nb_staging_t
*)0)->fshift
));
263 ad
->fshift
= clCreateBuffer(runData
->context
, CL_MEM_READ_WRITE
, SHIFTS
* ad
->fshift_elem_size
, NULL
, &cl_error
);
264 assert(cl_error
== CL_SUCCESS
);
265 // TODO: handle errors, check clCreateBuffer flags
267 ad
->e_lj
= clCreateBuffer(runData
->context
, CL_MEM_READ_WRITE
, sizeof(float), NULL
, &cl_error
);
268 assert(cl_error
== CL_SUCCESS
);
269 // TODO: handle errors, check clCreateBuffer flags
271 ad
->e_el
= clCreateBuffer(runData
->context
, CL_MEM_READ_WRITE
, sizeof(float), NULL
, &cl_error
);
272 assert(cl_error
== CL_SUCCESS
);
273 // TODO: handle errors, check clCreateBuffer flags
275 /* initialize to NULL pointers to data that is not allocated here and will
276 need reallocation in nbnxn_gpu_init_atomdata */
280 /* size -1 indicates that the respective array hasn't been initialized yet */
285 /*! \brief Copies all parameters related to the cut-off from ic to nbp
287 static void set_cutoff_parameters(cl_nbparam_t
*nbp
,
288 const interaction_const_t
*ic
,
289 const NbnxnListParameters
*listParams
)
291 nbp
->ewald_beta
= ic
->ewaldcoeff_q
;
292 nbp
->sh_ewald
= ic
->sh_ewald
;
293 nbp
->epsfac
= ic
->epsfac
;
294 nbp
->two_k_rf
= 2.0 * ic
->k_rf
;
295 nbp
->c_rf
= ic
->c_rf
;
296 nbp
->rvdw_sq
= ic
->rvdw
* ic
->rvdw
;
297 nbp
->rcoulomb_sq
= ic
->rcoulomb
* ic
->rcoulomb
;
298 nbp
->rlistOuter_sq
= listParams
->rlistOuter
* listParams
->rlistOuter
;
299 nbp
->rlistInner_sq
= listParams
->rlistInner
* listParams
->rlistInner
;
300 nbp
->useDynamicPruning
= listParams
->useDynamicPruning
;
302 nbp
->sh_lj_ewald
= ic
->sh_lj_ewald
;
303 nbp
->ewaldcoeff_lj
= ic
->ewaldcoeff_lj
;
305 nbp
->rvdw_switch
= ic
->rvdw_switch
;
306 nbp
->dispersion_shift
= ic
->dispersion_shift
;
307 nbp
->repulsion_shift
= ic
->repulsion_shift
;
308 nbp
->vdw_switch
= ic
->vdw_switch
;
311 /*! \brief Returns the kinds of electrostatics and Vdw OpenCL
312 * kernels that will be used.
314 * Respectively, these values are from enum eelOcl and enum
317 map_interaction_types_to_gpu_kernel_flavors(const interaction_const_t
*ic
,
322 if (ic
->vdwtype
== evdwCUT
)
324 switch (ic
->vdw_modifier
)
327 case eintmodPOTSHIFT
:
331 *gpu_vdwtype
= evdwOclCUT
;
334 *gpu_vdwtype
= evdwOclCUTCOMBGEOM
;
337 *gpu_vdwtype
= evdwOclCUTCOMBLB
;
340 gmx_incons("The requested LJ combination rule is not implemented in the OpenCL GPU accelerated kernels!");
344 case eintmodFORCESWITCH
:
345 *gpu_vdwtype
= evdwOclFSWITCH
;
347 case eintmodPOTSWITCH
:
348 *gpu_vdwtype
= evdwOclPSWITCH
;
351 gmx_incons("The requested VdW interaction modifier is not implemented in the GPU accelerated kernels!");
355 else if (ic
->vdwtype
== evdwPME
)
357 if (ic
->ljpme_comb_rule
== ljcrGEOM
)
359 *gpu_vdwtype
= evdwOclEWALDGEOM
;
363 *gpu_vdwtype
= evdwOclEWALDLB
;
368 gmx_incons("The requested VdW type is not implemented in the GPU accelerated kernels!");
371 if (ic
->eeltype
== eelCUT
)
373 *gpu_eeltype
= eelOclCUT
;
375 else if (EEL_RF(ic
->eeltype
))
377 *gpu_eeltype
= eelOclRF
;
379 else if ((EEL_PME(ic
->eeltype
) || ic
->eeltype
== eelEWALD
))
381 /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
382 *gpu_eeltype
= nbnxn_gpu_pick_ewald_kernel_type(false);
386 /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
387 gmx_incons("The requested electrostatics type is not implemented in the GPU accelerated kernels!");
391 /*! \brief Initializes the nonbonded parameter data structure.
393 static void init_nbparam(cl_nbparam_t
*nbp
,
394 const interaction_const_t
*ic
,
395 const NbnxnListParameters
*listParams
,
396 const nbnxn_atomdata_t
*nbat
,
397 const gmx_device_runtime_data_t
*runData
)
399 int ntypes
, nnbfp
, nnbfp_comb
;
403 ntypes
= nbat
->ntype
;
405 set_cutoff_parameters(nbp
, ic
, listParams
);
407 map_interaction_types_to_gpu_kernel_flavors(ic
,
412 if (ic
->vdwtype
== evdwPME
)
414 if (ic
->ljpme_comb_rule
== ljcrGEOM
)
416 assert(nbat
->comb_rule
== ljcrGEOM
);
420 assert(nbat
->comb_rule
== ljcrLB
);
423 /* generate table for PME */
424 nbp
->coulomb_tab_climg2d
= NULL
;
425 if (nbp
->eeltype
== eelOclEWALD_TAB
|| nbp
->eeltype
== eelOclEWALD_TAB_TWIN
)
427 init_ewald_coulomb_force_table(ic
, nbp
, runData
);
430 // TODO: improvement needed.
431 // The image2d is created here even if eeltype is not eelCuEWALD_TAB or eelCuEWALD_TAB_TWIN because the OpenCL kernels
432 // don't accept NULL values for image2D parameters.
434 /* Switched from using textures to using buffers */
435 // TODO: decide which alternative is most efficient - textures or buffers.
437 cl_image_format array_format;
439 array_format.image_channel_data_type = CL_FLOAT;
440 array_format.image_channel_order = CL_R;
442 nbp->coulomb_tab_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
443 &array_format, 1, 1, 0, NULL, &cl_error);
446 nbp
->coulomb_tab_climg2d
= clCreateBuffer(runData
->context
, CL_MEM_READ_ONLY
, sizeof(cl_float
), NULL
, &cl_error
);
447 // TODO: handle errors
450 nnbfp
= 2*ntypes
*ntypes
;
451 nnbfp_comb
= 2*ntypes
;
454 /* Switched from using textures to using buffers */
455 // TODO: decide which alternative is most efficient - textures or buffers.
457 cl_image_format array_format;
459 array_format.image_channel_data_type = CL_FLOAT;
460 array_format.image_channel_order = CL_R;
462 nbp->nbfp_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
463 &array_format, nnbfp, 1, 0, nbat->nbfp, &cl_error);
466 nbp
->nbfp_climg2d
= clCreateBuffer(runData
->context
, CL_MEM_READ_ONLY
| CL_MEM_COPY_HOST_PTR
, nnbfp
*sizeof(cl_float
), nbat
->nbfp
, &cl_error
);
467 assert(cl_error
== CL_SUCCESS
);
468 // TODO: handle errors
470 if (ic
->vdwtype
== evdwPME
)
472 /* Switched from using textures to using buffers */
473 // TODO: decide which alternative is most efficient - textures or buffers.
474 /* nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
475 &array_format, nnbfp_comb, 1, 0, nbat->nbfp_comb, &cl_error);*/
476 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
);
479 assert(cl_error
== CL_SUCCESS
);
480 // TODO: handle errors
484 // TODO: improvement needed.
485 // The image2d is created here even if vdwtype is not evdwPME because the OpenCL kernels
486 // don't accept NULL values for image2D parameters.
487 /* Switched from using textures to using buffers */
488 // TODO: decide which alternative is most efficient - textures or buffers.
489 /* nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
490 &array_format, 1, 1, 0, NULL, &cl_error);*/
491 nbp
->nbfp_comb_climg2d
= clCreateBuffer(runData
->context
, CL_MEM_READ_ONLY
, sizeof(cl_float
), NULL
, &cl_error
);
494 assert(cl_error
== CL_SUCCESS
);
495 // TODO: handle errors
500 //! This function is documented in the header file
501 void nbnxn_gpu_pme_loadbal_update_param(const nonbonded_verlet_t
*nbv
,
502 const interaction_const_t
*ic
,
503 const NbnxnListParameters
*listParams
)
505 if (!nbv
|| nbv
->grp
[0].kernel_type
!= nbnxnk8x8x8_GPU
)
509 gmx_nbnxn_ocl_t
*nb
= nbv
->gpu_nbv
;
510 cl_nbparam_t
*nbp
= nb
->nbparam
;
512 set_cutoff_parameters(nbp
, ic
, listParams
);
514 nbp
->eeltype
= nbnxn_gpu_pick_ewald_kernel_type(ic
->rcoulomb
!= ic
->rvdw
);
516 init_ewald_coulomb_force_table(ic
, nb
->nbparam
, nb
->dev_rundata
);
519 /*! \brief Initializes the pair list data structure.
521 static void init_plist(cl_plist_t
*pl
)
523 /* initialize to NULL pointers to data that is not allocated here and will
524 need reallocation in nbnxn_gpu_init_pairlist */
530 /* size -1 indicates that the respective array hasn't been initialized yet */
537 pl
->imask_nalloc
= -1;
539 pl
->excl_nalloc
= -1;
540 pl
->haveFreshList
= false;
543 /*! \brief Initializes the timer data structure.
545 static void init_timers(cl_timers_t
*t
,
548 for (int i
= 0; i
<= (bUseTwoStreams
? 1 : 0); i
++)
550 t
->didPairlistH2D
[i
] = false;
551 t
->didPrune
[i
] = false;
552 t
->didRollingPrune
[i
] = false;
556 /*! \brief Initializes the timings data structure.
558 static void init_timings(gmx_wallclock_gpu_nbnxn_t
*t
)
567 for (i
= 0; i
< 2; i
++)
569 for (j
= 0; j
< 2; j
++)
571 t
->ktime
[i
][j
].t
= 0.0;
572 t
->ktime
[i
][j
].c
= 0;
577 t
->pruneTime
.t
= 0.0;
578 t
->dynamicPruneTime
.c
= 0;
579 t
->dynamicPruneTime
.t
= 0.0;
582 /*! \brief Creates context for OpenCL GPU given by \p mygpu
584 * A fatal error results if creation fails.
586 * \param[inout] runtimeData runtime data including program and context
587 * \param[in] devInfo device info struct
588 * \param[in] rank MPI rank (for error reporting)
591 nbnxn_gpu_create_context(gmx_device_runtime_data_t
*runtimeData
,
592 const gmx_device_info_t
*devInfo
,
595 cl_context_properties context_properties
[3];
596 cl_platform_id platform_id
;
597 cl_device_id device_id
;
601 assert(runtimeData
!= NULL
);
602 assert(devInfo
!= NULL
);
604 platform_id
= devInfo
->ocl_gpu_id
.ocl_platform_id
;
605 device_id
= devInfo
->ocl_gpu_id
.ocl_device_id
;
607 context_properties
[0] = CL_CONTEXT_PLATFORM
;
608 context_properties
[1] = (cl_context_properties
) platform_id
;
609 context_properties
[2] = 0; /* Terminates the list of properties */
611 context
= clCreateContext(context_properties
, 1, &device_id
, NULL
, NULL
, &cl_error
);
612 if (CL_SUCCESS
!= cl_error
)
614 gmx_fatal(FARGS
, "On rank %d failed to create context for GPU #%s:\n OpenCL error %d: %s",
616 devInfo
->device_name
,
617 cl_error
, ocl_get_error_string(cl_error
).c_str());
621 runtimeData
->context
= context
;
624 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
625 static cl_kernel
nbnxn_gpu_create_kernel(gmx_nbnxn_ocl_t
*nb
,
626 const char *kernel_name
)
631 kernel
= clCreateKernel(nb
->dev_rundata
->program
, kernel_name
, &cl_error
);
632 if (CL_SUCCESS
!= cl_error
)
634 gmx_fatal(FARGS
, "Failed to create kernel '%s' for GPU #%s: OpenCL error %d",
636 nb
->dev_info
->device_name
,
643 /*! \brief Clears nonbonded shift force output array and energy outputs on the GPU.
646 nbnxn_ocl_clear_e_fshift(gmx_nbnxn_ocl_t
*nb
)
650 cl_atomdata_t
* adat
= nb
->atdat
;
651 cl_command_queue ls
= nb
->stream
[eintLocal
];
653 size_t local_work_size
[3] = {1, 1, 1};
654 size_t global_work_size
[3] = {1, 1, 1};
656 cl_int shifts
= SHIFTS
*3;
660 cl_kernel zero_e_fshift
= nb
->kernel_zero_e_fshift
;
662 local_work_size
[0] = 64;
663 // Round the total number of threads up from the array size
664 global_work_size
[0] = ((shifts
+ local_work_size
[0] - 1)/local_work_size
[0])*local_work_size
[0];
667 cl_error
= clSetKernelArg(zero_e_fshift
, arg_no
++, sizeof(cl_mem
), &(adat
->fshift
));
668 cl_error
|= clSetKernelArg(zero_e_fshift
, arg_no
++, sizeof(cl_mem
), &(adat
->e_lj
));
669 cl_error
|= clSetKernelArg(zero_e_fshift
, arg_no
++, sizeof(cl_mem
), &(adat
->e_el
));
670 cl_error
|= clSetKernelArg(zero_e_fshift
, arg_no
++, sizeof(cl_uint
), &shifts
);
671 GMX_ASSERT(cl_error
== CL_SUCCESS
, ocl_get_error_string(cl_error
).c_str());
673 cl_error
= clEnqueueNDRangeKernel(ls
, zero_e_fshift
, 3, NULL
, global_work_size
, local_work_size
, 0, NULL
, NULL
);
674 GMX_ASSERT(cl_error
== CL_SUCCESS
, ocl_get_error_string(cl_error
).c_str());
677 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
678 static void nbnxn_gpu_init_kernels(gmx_nbnxn_ocl_t
*nb
)
680 /* Init to 0 main kernel arrays */
681 /* They will be later on initialized in select_nbnxn_kernel */
682 // TODO: consider always creating all variants of the kernels here so that there is no
683 // need for late call to clCreateKernel -- if that gives any advantage?
684 memset(nb
->kernel_ener_noprune_ptr
, 0, sizeof(nb
->kernel_ener_noprune_ptr
));
685 memset(nb
->kernel_ener_prune_ptr
, 0, sizeof(nb
->kernel_ener_prune_ptr
));
686 memset(nb
->kernel_noener_noprune_ptr
, 0, sizeof(nb
->kernel_noener_noprune_ptr
));
687 memset(nb
->kernel_noener_prune_ptr
, 0, sizeof(nb
->kernel_noener_prune_ptr
));
689 /* Init pruning kernels
691 * TODO: we could avoid creating kernels if dynamic pruning is turned off,
692 * but ATM that depends on force flags not passed into the initialization.
694 nb
->kernel_pruneonly
[epruneFirst
] = nbnxn_gpu_create_kernel(nb
, "nbnxn_kernel_prune_opencl");
695 nb
->kernel_pruneonly
[epruneRolling
] = nbnxn_gpu_create_kernel(nb
, "nbnxn_kernel_prune_rolling_opencl");
697 /* Init auxiliary kernels */
698 nb
->kernel_memset_f
= nbnxn_gpu_create_kernel(nb
, "memset_f");
699 nb
->kernel_memset_f2
= nbnxn_gpu_create_kernel(nb
, "memset_f2");
700 nb
->kernel_memset_f3
= nbnxn_gpu_create_kernel(nb
, "memset_f3");
701 nb
->kernel_zero_e_fshift
= nbnxn_gpu_create_kernel(nb
, "zero_e_fshift");
704 /*! \brief Initializes simulation constant data.
706 * Initializes members of the atomdata and nbparam structs and
707 * clears e/fshift output buffers.
709 static void nbnxn_ocl_init_const(gmx_nbnxn_ocl_t
*nb
,
710 const interaction_const_t
*ic
,
711 const NbnxnListParameters
*listParams
,
712 const nbnxn_atomdata_t
*nbat
)
714 init_atomdata_first(nb
->atdat
, nbat
->ntype
, nb
->dev_rundata
);
715 init_nbparam(nb
->nbparam
, ic
, listParams
, nbat
, nb
->dev_rundata
);
719 //! This function is documented in the header file
720 void nbnxn_gpu_init(gmx_nbnxn_ocl_t
**p_nb
,
721 const gmx_device_info_t
*deviceInfo
,
722 const interaction_const_t
*ic
,
723 const NbnxnListParameters
*listParams
,
724 const nbnxn_atomdata_t
*nbat
,
726 gmx_bool bLocalAndNonlocal
)
730 cl_command_queue_properties queue_properties
;
741 snew(nb
->nbparam
, 1);
742 snew(nb
->plist
[eintLocal
], 1);
743 if (bLocalAndNonlocal
)
745 snew(nb
->plist
[eintNonlocal
], 1);
748 nb
->bUseTwoStreams
= bLocalAndNonlocal
;
750 nb
->timers
= new cl_timers_t();
751 snew(nb
->timings
, 1);
753 /* set device info, just point it to the right GPU among the detected ones */
754 nb
->dev_info
= deviceInfo
;
755 snew(nb
->dev_rundata
, 1);
757 /* init to NULL the debug buffer */
758 nb
->debug_buffer
= NULL
;
761 ocl_pmalloc((void**)&nb
->nbst
.e_lj
, sizeof(*nb
->nbst
.e_lj
));
762 ocl_pmalloc((void**)&nb
->nbst
.e_el
, sizeof(*nb
->nbst
.e_el
));
763 ocl_pmalloc((void**)&nb
->nbst
.fshift
, SHIFTS
* sizeof(*nb
->nbst
.fshift
));
765 init_plist(nb
->plist
[eintLocal
]);
767 /* OpenCL timing disabled if GMX_DISABLE_OCL_TIMING is defined. */
768 /* TODO deprecate the first env var in the 2017 release. */
769 nb
->bDoTime
= (getenv("GMX_DISABLE_OCL_TIMING") == NULL
&&
770 getenv("GMX_DISABLE_GPU_TIMING") == NULL
);
772 /* Create queues only after bDoTime has been initialized */
775 queue_properties
= CL_QUEUE_PROFILING_ENABLE
;
779 queue_properties
= 0;
782 nbnxn_gpu_create_context(nb
->dev_rundata
, nb
->dev_info
, rank
);
784 /* local/non-local GPU streams */
785 nb
->stream
[eintLocal
] = clCreateCommandQueue(nb
->dev_rundata
->context
, nb
->dev_info
->ocl_gpu_id
.ocl_device_id
, queue_properties
, &cl_error
);
786 if (CL_SUCCESS
!= cl_error
)
788 gmx_fatal(FARGS
, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
790 nb
->dev_info
->device_name
,
795 if (nb
->bUseTwoStreams
)
797 init_plist(nb
->plist
[eintNonlocal
]);
799 nb
->stream
[eintNonlocal
] = clCreateCommandQueue(nb
->dev_rundata
->context
, nb
->dev_info
->ocl_gpu_id
.ocl_device_id
, queue_properties
, &cl_error
);
800 if (CL_SUCCESS
!= cl_error
)
802 gmx_fatal(FARGS
, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
804 nb
->dev_info
->device_name
,
812 init_timers(nb
->timers
, nb
->bUseTwoStreams
);
813 init_timings(nb
->timings
);
816 nbnxn_ocl_init_const(nb
, ic
, listParams
, nbat
);
818 /* Enable LJ param manual prefetch for AMD or if we request through env. var.
819 * TODO: decide about NVIDIA
821 nb
->bPrefetchLjParam
=
822 (getenv("GMX_OCL_DISABLE_I_PREFETCH") == NULL
) &&
823 ((nb
->dev_info
->vendor_e
== OCL_VENDOR_AMD
) || (getenv("GMX_OCL_ENABLE_I_PREFETCH") != NULL
));
825 /* NOTE: in CUDA we pick L1 cache configuration for the nbnxn kernels here,
826 * but sadly this is not supported in OpenCL (yet?). Consider adding it if
827 * it becomes supported.
829 nbnxn_gpu_compile_kernels(nb
);
830 nbnxn_gpu_init_kernels(nb
);
832 /* clear energy and shift force outputs */
833 nbnxn_ocl_clear_e_fshift(nb
);
839 fprintf(debug
, "Initialized OpenCL data structures.\n");
843 /*! \brief Clears the first natoms_clear elements of the GPU nonbonded force output array.
845 static void nbnxn_ocl_clear_f(gmx_nbnxn_ocl_t
*nb
, int natoms_clear
)
847 if (natoms_clear
== 0)
853 cl_atomdata_t
* adat
= nb
->atdat
;
854 cl_command_queue ls
= nb
->stream
[eintLocal
];
855 cl_float value
= 0.0f
;
857 size_t local_work_size
[3] = {1, 1, 1};
858 size_t global_work_size
[3] = {1, 1, 1};
862 cl_kernel memset_f
= nb
->kernel_memset_f
;
864 cl_uint natoms_flat
= natoms_clear
* (sizeof(rvec
)/sizeof(real
));
866 local_work_size
[0] = 64;
867 // Round the total number of threads up from the array size
868 global_work_size
[0] = ((natoms_flat
+ local_work_size
[0] - 1)/local_work_size
[0])*local_work_size
[0];
872 cl_error
= clSetKernelArg(memset_f
, arg_no
++, sizeof(cl_mem
), &(adat
->f
));
873 cl_error
|= clSetKernelArg(memset_f
, arg_no
++, sizeof(cl_float
), &value
);
874 cl_error
|= clSetKernelArg(memset_f
, arg_no
++, sizeof(cl_uint
), &natoms_flat
);
875 assert(cl_error
== CL_SUCCESS
);
877 cl_error
= clEnqueueNDRangeKernel(ls
, memset_f
, 3, NULL
, global_work_size
, local_work_size
, 0, NULL
, NULL
);
878 assert(cl_error
== CL_SUCCESS
);
881 //! This function is documented in the header file
883 nbnxn_gpu_clear_outputs(gmx_nbnxn_ocl_t
*nb
,
886 nbnxn_ocl_clear_f(nb
, nb
->atdat
->natoms
);
887 /* clear shift force array and energies if the outputs were
888 used in the current step */
889 if (flags
& GMX_FORCE_VIRIAL
)
891 nbnxn_ocl_clear_e_fshift(nb
);
894 /* kick off buffer clearing kernel to ensure concurrency with constraints/update */
895 cl_int gmx_unused cl_error
;
896 cl_error
= clFlush(nb
->stream
[eintLocal
]);
897 assert(CL_SUCCESS
== cl_error
);
900 //! This function is documented in the header file
901 void nbnxn_gpu_init_pairlist(gmx_nbnxn_ocl_t
*nb
,
902 const nbnxn_pairlist_t
*h_plist
,
905 if (canSkipWork(nb
, iloc
))
911 bool bDoTime
= nb
->bDoTime
;
912 cl_command_queue stream
= nb
->stream
[iloc
];
913 cl_plist_t
*d_plist
= nb
->plist
[iloc
];
915 if (d_plist
->na_c
< 0)
917 d_plist
->na_c
= h_plist
->na_ci
;
921 if (d_plist
->na_c
!= h_plist
->na_ci
)
923 sprintf(sbuf
, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
924 d_plist
->na_c
, h_plist
->na_ci
);
931 nb
->timers
->pl_h2d
[iloc
].openTimingRegion(stream
);
932 nb
->timers
->didPairlistH2D
[iloc
] = true;
935 ocl_realloc_buffered(&d_plist
->sci
, h_plist
->sci
, sizeof(nbnxn_sci_t
),
936 &d_plist
->nsci
, &d_plist
->sci_nalloc
,
938 nb
->dev_rundata
->context
,
939 stream
, true, bDoTime
? nb
->timers
->pl_h2d
[iloc
].fetchNextEvent() : nullptr);
941 ocl_realloc_buffered(&d_plist
->cj4
, h_plist
->cj4
, sizeof(nbnxn_cj4_t
),
942 &d_plist
->ncj4
, &d_plist
->cj4_nalloc
,
944 nb
->dev_rundata
->context
,
945 stream
, true, bDoTime
? nb
->timers
->pl_h2d
[iloc
].fetchNextEvent() : nullptr);
947 /* this call only allocates space on the device (no data is transferred) - no timing as well! */
948 ocl_realloc_buffered(&d_plist
->imask
, NULL
, sizeof(unsigned int),
949 &d_plist
->nimask
, &d_plist
->imask_nalloc
,
950 h_plist
->ncj4
*c_nbnxnGpuClusterpairSplit
,
951 nb
->dev_rundata
->context
,
954 ocl_realloc_buffered(&d_plist
->excl
, h_plist
->excl
, sizeof(nbnxn_excl_t
),
955 &d_plist
->nexcl
, &d_plist
->excl_nalloc
,
957 nb
->dev_rundata
->context
,
958 stream
, true, bDoTime
? nb
->timers
->pl_h2d
[iloc
].fetchNextEvent() : nullptr);
962 nb
->timers
->pl_h2d
[iloc
].closeTimingRegion(stream
);
965 /* need to prune the pair list during the next step */
966 d_plist
->haveFreshList
= true;
969 //! This function is documented in the header file
970 void nbnxn_gpu_upload_shiftvec(gmx_nbnxn_ocl_t
*nb
,
971 const nbnxn_atomdata_t
*nbatom
)
973 cl_atomdata_t
*adat
= nb
->atdat
;
974 cl_command_queue ls
= nb
->stream
[eintLocal
];
976 /* only if we have a dynamic box */
977 if (nbatom
->bDynamicBox
|| !adat
->bShiftVecUploaded
)
979 ocl_copy_H2D_async(adat
->shift_vec
, nbatom
->shift_vec
, 0,
980 SHIFTS
* adat
->shift_vec_elem_size
, ls
, NULL
);
981 adat
->bShiftVecUploaded
= true;
985 //! This function is documented in the header file
986 void nbnxn_gpu_init_atomdata(gmx_nbnxn_ocl_t
*nb
,
987 const nbnxn_atomdata_t
*nbat
)
992 bool bDoTime
= nb
->bDoTime
;
993 cl_timers_t
*timers
= nb
->timers
;
994 cl_atomdata_t
*d_atdat
= nb
->atdat
;
995 cl_command_queue ls
= nb
->stream
[eintLocal
];
997 natoms
= nbat
->natoms
;
1002 /* time async copy */
1003 timers
->atdat
.openTimingRegion(ls
);
1006 /* need to reallocate if we have to copy more atoms than the amount of space
1007 available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
1008 if (natoms
> d_atdat
->nalloc
)
1010 nalloc
= over_alloc_small(natoms
);
1012 /* free up first if the arrays have already been initialized */
1013 if (d_atdat
->nalloc
!= -1)
1015 ocl_free_buffered(d_atdat
->f
, &d_atdat
->natoms
, &d_atdat
->nalloc
);
1016 ocl_free_buffered(d_atdat
->xq
, NULL
, NULL
);
1017 ocl_free_buffered(d_atdat
->lj_comb
, NULL
, NULL
);
1018 ocl_free_buffered(d_atdat
->atom_types
, NULL
, NULL
);
1021 d_atdat
->f_elem_size
= sizeof(rvec
);
1023 // TODO: handle errors, check clCreateBuffer flags
1024 d_atdat
->f
= clCreateBuffer(nb
->dev_rundata
->context
, CL_MEM_READ_WRITE
, nalloc
* d_atdat
->f_elem_size
, NULL
, &cl_error
);
1025 assert(CL_SUCCESS
== cl_error
);
1027 // TODO: change the flag to read-only
1028 d_atdat
->xq
= clCreateBuffer(nb
->dev_rundata
->context
, CL_MEM_READ_WRITE
, nalloc
* sizeof(cl_float4
), NULL
, &cl_error
);
1029 assert(CL_SUCCESS
== cl_error
);
1030 // TODO: handle errors, check clCreateBuffer flags
1032 if (useLjCombRule(nb
->nbparam
->vdwtype
))
1034 // TODO: change the flag to read-only
1035 d_atdat
->lj_comb
= clCreateBuffer(nb
->dev_rundata
->context
, CL_MEM_READ_WRITE
, nalloc
* sizeof(cl_float2
), NULL
, &cl_error
);
1036 assert(CL_SUCCESS
== cl_error
);
1037 // TODO: handle errors, check clCreateBuffer flags
1041 // TODO: change the flag to read-only
1042 d_atdat
->atom_types
= clCreateBuffer(nb
->dev_rundata
->context
, CL_MEM_READ_WRITE
, nalloc
* sizeof(int), NULL
, &cl_error
);
1043 assert(CL_SUCCESS
== cl_error
);
1044 // TODO: handle errors, check clCreateBuffer flags
1047 d_atdat
->nalloc
= nalloc
;
1051 d_atdat
->natoms
= natoms
;
1052 d_atdat
->natoms_local
= nbat
->natoms_local
;
1054 /* need to clear GPU f output if realloc happened */
1057 nbnxn_ocl_clear_f(nb
, nalloc
);
1060 if (useLjCombRule(nb
->nbparam
->vdwtype
))
1062 ocl_copy_H2D_async(d_atdat
->lj_comb
, nbat
->lj_comb
, 0,
1063 natoms
*sizeof(cl_float2
), ls
, bDoTime
? timers
->atdat
.fetchNextEvent() : nullptr);
1067 ocl_copy_H2D_async(d_atdat
->atom_types
, nbat
->type
, 0,
1068 natoms
*sizeof(int), ls
, bDoTime
? timers
->atdat
.fetchNextEvent() : nullptr);
1074 timers
->atdat
.closeTimingRegion(ls
);
1077 /* kick off the tasks enqueued above to ensure concurrency with the search */
1078 cl_error
= clFlush(ls
);
1079 assert(CL_SUCCESS
== cl_error
);
1082 /*! \brief Releases an OpenCL kernel pointer */
1083 static void free_kernel(cl_kernel
*kernel_ptr
)
1085 cl_int gmx_unused cl_error
;
1087 assert(NULL
!= kernel_ptr
);
1091 cl_error
= clReleaseKernel(*kernel_ptr
);
1092 assert(cl_error
== CL_SUCCESS
);
1098 /*! \brief Releases a list of OpenCL kernel pointers */
1099 static void free_kernels(cl_kernel
*kernels
, int count
)
1103 for (i
= 0; i
< count
; i
++)
1105 free_kernel(kernels
+ i
);
1109 /*! \brief Free the OpenCL runtime data (context and program).
1111 * The function releases the OpenCL context and program assuciated with the
1112 * device that the calling PP rank is running on.
1114 * \param runData [in] porinter to the structure with runtime data.
1116 static void free_gpu_device_runtime_data(gmx_device_runtime_data_t
*runData
)
1118 if (runData
== NULL
)
1123 cl_int gmx_unused cl_error
;
1125 if (runData
->context
)
1127 cl_error
= clReleaseContext(runData
->context
);
1128 runData
->context
= NULL
;
1129 assert(CL_SUCCESS
== cl_error
);
1132 if (runData
->program
)
1134 cl_error
= clReleaseProgram(runData
->program
);
1135 runData
->program
= NULL
;
1136 assert(CL_SUCCESS
== cl_error
);
1141 //! This function is documented in the header file
1142 void nbnxn_gpu_free(gmx_nbnxn_ocl_t
*nb
)
1150 int kernel_count
= sizeof(nb
->kernel_ener_noprune_ptr
) / sizeof(nb
->kernel_ener_noprune_ptr
[0][0]);
1151 free_kernels((cl_kernel
*)nb
->kernel_ener_noprune_ptr
, kernel_count
);
1153 kernel_count
= sizeof(nb
->kernel_ener_prune_ptr
) / sizeof(nb
->kernel_ener_prune_ptr
[0][0]);
1154 free_kernels((cl_kernel
*)nb
->kernel_ener_prune_ptr
, kernel_count
);
1156 kernel_count
= sizeof(nb
->kernel_noener_noprune_ptr
) / sizeof(nb
->kernel_noener_noprune_ptr
[0][0]);
1157 free_kernels((cl_kernel
*)nb
->kernel_noener_noprune_ptr
, kernel_count
);
1159 kernel_count
= sizeof(nb
->kernel_noener_prune_ptr
) / sizeof(nb
->kernel_noener_prune_ptr
[0][0]);
1160 free_kernels((cl_kernel
*)nb
->kernel_noener_prune_ptr
, kernel_count
);
1162 free_kernel(&(nb
->kernel_memset_f
));
1163 free_kernel(&(nb
->kernel_memset_f2
));
1164 free_kernel(&(nb
->kernel_memset_f3
));
1165 free_kernel(&(nb
->kernel_zero_e_fshift
));
1168 free_ocl_buffer(&(nb
->atdat
->xq
));
1169 free_ocl_buffer(&(nb
->atdat
->f
));
1170 free_ocl_buffer(&(nb
->atdat
->e_lj
));
1171 free_ocl_buffer(&(nb
->atdat
->e_el
));
1172 free_ocl_buffer(&(nb
->atdat
->fshift
));
1173 free_ocl_buffer(&(nb
->atdat
->lj_comb
));
1174 free_ocl_buffer(&(nb
->atdat
->atom_types
));
1175 free_ocl_buffer(&(nb
->atdat
->shift_vec
));
1179 free_ocl_buffer(&(nb
->nbparam
->nbfp_climg2d
));
1180 free_ocl_buffer(&(nb
->nbparam
->nbfp_comb_climg2d
));
1181 free_ocl_buffer(&(nb
->nbparam
->coulomb_tab_climg2d
));
1185 free_ocl_buffer(&(nb
->plist
[eintLocal
]->sci
));
1186 free_ocl_buffer(&(nb
->plist
[eintLocal
]->cj4
));
1187 free_ocl_buffer(&(nb
->plist
[eintLocal
]->imask
));
1188 free_ocl_buffer(&(nb
->plist
[eintLocal
]->excl
));
1189 sfree(nb
->plist
[eintLocal
]);
1190 if (nb
->bUseTwoStreams
)
1192 free_ocl_buffer(&(nb
->plist
[eintNonlocal
]->sci
));
1193 free_ocl_buffer(&(nb
->plist
[eintNonlocal
]->cj4
));
1194 free_ocl_buffer(&(nb
->plist
[eintNonlocal
]->imask
));
1195 free_ocl_buffer(&(nb
->plist
[eintNonlocal
]->excl
));
1196 sfree(nb
->plist
[eintNonlocal
]);
1200 ocl_pfree(nb
->nbst
.e_lj
);
1201 nb
->nbst
.e_lj
= NULL
;
1203 ocl_pfree(nb
->nbst
.e_el
);
1204 nb
->nbst
.e_el
= NULL
;
1206 ocl_pfree(nb
->nbst
.fshift
);
1207 nb
->nbst
.fshift
= NULL
;
1209 /* Free debug buffer */
1210 free_ocl_buffer(&nb
->debug_buffer
);
1212 /* Free command queues */
1213 clReleaseCommandQueue(nb
->stream
[eintLocal
]);
1214 nb
->stream
[eintLocal
] = NULL
;
1215 if (nb
->bUseTwoStreams
)
1217 clReleaseCommandQueue(nb
->stream
[eintNonlocal
]);
1218 nb
->stream
[eintNonlocal
] = NULL
;
1220 /* Free other events */
1221 if (nb
->nonlocal_done
)
1223 clReleaseEvent(nb
->nonlocal_done
);
1224 nb
->nonlocal_done
= NULL
;
1226 if (nb
->misc_ops_and_local_H2D_done
)
1228 clReleaseEvent(nb
->misc_ops_and_local_H2D_done
);
1229 nb
->misc_ops_and_local_H2D_done
= NULL
;
1232 free_gpu_device_runtime_data(nb
->dev_rundata
);
1233 sfree(nb
->dev_rundata
);
1235 /* Free timers and timings */
1242 fprintf(debug
, "Cleaned up OpenCL data structures.\n");
1246 //! This function is documented in the header file
1247 gmx_wallclock_gpu_nbnxn_t
*nbnxn_gpu_get_timings(gmx_nbnxn_ocl_t
*nb
)
1249 return (nb
!= nullptr && nb
->bDoTime
) ? nb
->timings
: nullptr;
1252 //! This function is documented in the header file
1253 void nbnxn_gpu_reset_timings(nonbonded_verlet_t
* nbv
)
1255 if (nbv
->gpu_nbv
&& nbv
->gpu_nbv
->bDoTime
)
1257 init_timings(nbv
->gpu_nbv
->timings
);
1261 //! This function is documented in the header file
1262 int nbnxn_gpu_min_ci_balanced(gmx_nbnxn_ocl_t
*nb
)
1265 gpu_min_ci_balanced_factor
* nb
->dev_info
->compute_units
: 0;
1268 //! This function is documented in the header file
1269 gmx_bool
nbnxn_gpu_is_kernel_ewald_analytical(const gmx_nbnxn_ocl_t
*nb
)
1271 return ((nb
->nbparam
->eeltype
== eelOclEWALD_ANA
) ||
1272 (nb
->nbparam
->eeltype
== eelOclEWALD_ANA_TWIN
));