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