src/gromacs/nbnxm/nbnxm_setup.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, 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 Common functions for the different NBNXN GPU implementations.
  37  *
  38  * \author Berk Hess <hess@kth.se>
  39  *
  40  * \ingroup module_nbnxm
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "gromacs/domdec/domdec.h"
  46 #include "gromacs/domdec/domdec_struct.h"
  47 #include "gromacs/hardware/hw_info.h"
  48 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  49 #include "gromacs/mdtypes/commrec.h"
  50 #include "gromacs/mdtypes/forcerec.h"
  51 #include "gromacs/mdtypes/inputrec.h"
  52 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  53 #include "gromacs/nbnxm/nbnxm.h"
  54 #include "gromacs/nbnxm/pairlist_tuning.h"
  55 #include "gromacs/simd/simd.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/logger.h"
  58
  59 #include "atomdata.h"
  60 #include "gpu_types.h"
  61 #include "grid.h"
  62 #include "nbnxm_geometry.h"
  63 #include "nbnxm_simd.h"
  64 #include "pairlist.h"
  65 #include "pairlistset.h"
  66 #include "pairlistsets.h"
  67 #include "pairsearch.h"
  68
  69 namespace Nbnxm
  70 {
  71
  72 /*! \brief Resources that can be used to execute non-bonded kernels on */
  73 enum class NonbondedResource : int
  74 {
  75     Cpu,
  76     Gpu,
  77     EmulateGpu
  78 };
  79
  80 /*! \brief Returns whether CPU SIMD support exists for the given inputrec
  81  *
  82  * If the return value is FALSE and fplog/cr != NULL, prints a fallback
  83  * message to fplog/stderr.
  84  */
  85 static gmx_bool nbnxn_simd_supported(const gmx::MDLogger &mdlog,
  86                                      const t_inputrec    *ir)
  87 {
  88     if (ir->vdwtype == evdwPME && ir->ljpme_combination_rule == eljpmeLB)
  89     {
  90         /* LJ PME with LB combination rule does 7 mesh operations.
  91          * This so slow that we don't compile SIMD non-bonded kernels
  92          * for that. */
  93         GMX_LOG(mdlog.warning).asParagraph().appendText("LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling back to plain C kernels");
  94         return FALSE;
  95     }
  96
  97     return TRUE;
  98 }
  99
 100 /*! \brief Returns the most suitable CPU kernel type and Ewald handling */
 101 static KernelSetup
 102 pick_nbnxn_kernel_cpu(const t_inputrec gmx_unused    *ir,
 103                       const gmx_hw_info_t gmx_unused &hardwareInfo)
 104 {
 105     KernelSetup kernelSetup;
 106
 107     if (!GMX_SIMD)
 108     {
 109         kernelSetup.kernelType         = KernelType::Cpu4x4_PlainC;
 110         kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
 111     }
 112     else
 113     {
 114 #ifdef GMX_NBNXN_SIMD_4XN
 115         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 116 #endif
 117 #ifdef GMX_NBNXN_SIMD_2XNN
 118         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 119 #endif
 120
 121 #if defined GMX_NBNXN_SIMD_2XNN && defined GMX_NBNXN_SIMD_4XN
 122         /* We need to choose if we want 2x(N+N) or 4xN kernels.
 123          * This is based on the SIMD acceleration choice and CPU information
 124          * detected at runtime.
 125          *
 126          * 4xN calculates more (zero) interactions, but has less pair-search
 127          * work and much better kernel instruction scheduling.
 128          *
 129          * Up till now we have only seen that on Intel Sandy/Ivy Bridge,
 130          * which doesn't have FMA, both the analytical and tabulated Ewald
 131          * kernels have similar pair rates for 4x8 and 2x(4+4), so we choose
 132          * 2x(4+4) because it results in significantly fewer pairs.
 133          * For RF, the raw pair rate of the 4x8 kernel is higher than 2x(4+4),
 134          * 10% with HT, 50% without HT. As we currently don't detect the actual
 135          * use of HT, use 4x8 to avoid a potential performance hit.
 136          * On Intel Haswell 4x8 is always faster.
 137          */
 138         kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 139
 140         if (!GMX_SIMD_HAVE_FMA && (EEL_PME_EWALD(ir->coulombtype) ||
 141                                    EVDW_PME(ir->vdwtype)))
 142         {
 143             /* We have Ewald kernels without FMA (Intel Sandy/Ivy Bridge).
 144              * There are enough instructions to make 2x(4+4) efficient.
 145              */
 146             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 147         }
 148
 149         if (hardwareInfo.haveAmdZen1Cpu)
 150         {
 151             /* One 256-bit FMA per cycle makes 2xNN faster */
 152             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 153         }
 154 #endif      /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
 155
 156
 157         if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
 158         {
 159 #ifdef GMX_NBNXN_SIMD_4XN
 160             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_4xN;
 161 #else
 162             gmx_fatal(FARGS, "SIMD 4xN kernels requested, but GROMACS has been compiled without support for these kernels");
 163 #endif
 164         }
 165         if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
 166         {
 167 #ifdef GMX_NBNXN_SIMD_2XNN
 168             kernelSetup.kernelType = KernelType::Cpu4xN_Simd_2xNN;
 169 #else
 170             gmx_fatal(FARGS, "SIMD 2x(N+N) kernels requested, but GROMACS has been compiled without support for these kernels");
 171 #endif
 172         }
 173
 174         /* Analytical Ewald exclusion correction is only an option in
 175          * the SIMD kernel.
 176          * Since table lookup's don't parallelize with SIMD, analytical
 177          * will probably always be faster for a SIMD width of 8 or more.
 178          * With FMA analytical is sometimes faster for a width if 4 as well.
 179          * In single precision, this is faster on Bulldozer.
 180          * On AMD Zen, tabulated Ewald kernels are faster on all 4 combinations
 181          * of single or double precision and 128 or 256-bit AVX2.
 182          */
 183         if (
 184 #if GMX_SIMD
 185             (GMX_SIMD_REAL_WIDTH >= 8 ||
 186              (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)) &&
 187 #endif
 188             !hardwareInfo.haveAmdZen1Cpu)
 189         {
 190             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 191         }
 192         else
 193         {
 194             kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
 195         }
 196         if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
 197         {
 198             kernelSetup.ewaldExclusionType = EwaldExclusionType::Table;
 199         }
 200         if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
 201         {
 202             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 203         }
 204
 205     }
 206
 207     return kernelSetup;
 208 }
 209
 210 const char *lookup_kernel_name(const KernelType kernelType)
 211 {
 212     const char *returnvalue = nullptr;
 213     switch (kernelType)
 214     {
 215         case KernelType::NotSet:
 216             returnvalue = "not set";
 217             break;
 218         case KernelType::Cpu4x4_PlainC:
 219             returnvalue = "plain C";
 220             break;
 221         case KernelType::Cpu4xN_Simd_4xN:
 222         case KernelType::Cpu4xN_Simd_2xNN:
 223 #if GMX_SIMD
 224             returnvalue = "SIMD";
 225 #else  // GMX_SIMD
 226             returnvalue = "not available";
 227 #endif // GMX_SIMD
 228             break;
 229         case KernelType::Gpu8x8x8: returnvalue        = "GPU"; break;
 230         case KernelType::Cpu8x8x8_PlainC: returnvalue = "plain C"; break;
 231
 232         default:
 233             gmx_fatal(FARGS, "Illegal kernel type selected");
 234     }
 235     return returnvalue;
 236 };
 237
 238 /*! \brief Returns the most suitable kernel type and Ewald handling */
 239 static KernelSetup
 240 pick_nbnxn_kernel(const gmx::MDLogger     &mdlog,
 241                   gmx_bool                 use_simd_kernels,
 242                   const gmx_hw_info_t     &hardwareInfo,
 243                   const NonbondedResource &nonbondedResource,
 244                   const t_inputrec        *ir,
 245                   gmx_bool                 bDoNonbonded)
 246 {
 247     KernelSetup kernelSetup;
 248
 249     if (nonbondedResource == NonbondedResource::EmulateGpu)
 250     {
 251         kernelSetup.kernelType         = KernelType::Cpu8x8x8_PlainC;
 252         kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
 253
 254         if (bDoNonbonded)
 255         {
 256             GMX_LOG(mdlog.warning).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
 257         }
 258     }
 259     else if (nonbondedResource == NonbondedResource::Gpu)
 260     {
 261         kernelSetup.kernelType         = KernelType::Gpu8x8x8;
 262         kernelSetup.ewaldExclusionType = EwaldExclusionType::DecidedByGpuModule;
 263     }
 264     else
 265     {
 266         if (use_simd_kernels &&
 267             nbnxn_simd_supported(mdlog, ir))
 268         {
 269             kernelSetup = pick_nbnxn_kernel_cpu(ir, hardwareInfo);
 270         }
 271         else
 272         {
 273             kernelSetup.kernelType         = KernelType::Cpu4x4_PlainC;
 274             kernelSetup.ewaldExclusionType = EwaldExclusionType::Analytical;
 275         }
 276     }
 277
 278     if (bDoNonbonded)
 279     {
 280         GMX_LOG(mdlog.info).asParagraph().appendTextFormatted(
 281                 "Using %s %dx%d nonbonded short-range kernels",
 282                 lookup_kernel_name(kernelSetup.kernelType),
 283                 IClusterSizePerKernelType[kernelSetup.kernelType],
 284                 JClusterSizePerKernelType[kernelSetup.kernelType]);
 285
 286         if (KernelType::Cpu4x4_PlainC == kernelSetup.kernelType ||
 287             KernelType::Cpu8x8x8_PlainC == kernelSetup.kernelType)
 288         {
 289             GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 290                     "WARNING: Using the slow %s kernels. This should\n"
 291                     "not happen during routine usage on supported platforms.",
 292                     lookup_kernel_name(kernelSetup.kernelType));
 293         }
 294     }
 295
 296     GMX_RELEASE_ASSERT(kernelSetup.kernelType != KernelType::NotSet &&
 297                        kernelSetup.ewaldExclusionType != EwaldExclusionType::NotSet,
 298                        "All kernel setup parameters should be set here");
 299
 300     return kernelSetup;
 301 }
 302
 303 } // namespace Nbnxm
 304
 305 PairlistSets::PairlistSets(const PairlistParams &pairlistParams,
 306                            const bool            haveMultipleDomains,
 307                            const int             minimumIlistCountForGpuBalancing) :
 308     params_(pairlistParams),
 309     minimumIlistCountForGpuBalancing_(minimumIlistCountForGpuBalancing)
 310 {
 311     localSet_ = std::make_unique<PairlistSet>(Nbnxm::InteractionLocality::Local,
 312                                               params_);
 313
 314     if (haveMultipleDomains)
 315     {
 316         nonlocalSet_ = std::make_unique<PairlistSet>(Nbnxm::InteractionLocality::NonLocal,
 317                                                      params_);
 318     }
 319 }
 320
 321 namespace Nbnxm
 322 {
 323
 324 /*! \brief Gets and returns the minimum i-list count for balacing based on the GPU used or env.var. when set */
 325 static int getMinimumIlistCountForGpuBalancing(gmx_nbnxn_gpu_t *nbnxmGpu)
 326 {
 327     int minimumIlistCount;
 328
 329     if (const char *env = getenv("GMX_NB_MIN_CI"))
 330     {
 331         char *end;
 332
 333         minimumIlistCount = strtol(env, &end, 10);
 334         if (!end || (*end != 0) || minimumIlistCount < 0)
 335         {
 336             gmx_fatal(FARGS, "Invalid value passed in GMX_NB_MIN_CI=%s, non-negative integer required", env);
 337         }
 338
 339         if (debug)
 340         {
 341             fprintf(debug, "Neighbor-list balancing parameter: %d (passed as env. var.)\n",
 342                     minimumIlistCount);
 343         }
 344     }
 345     else
 346     {
 347         minimumIlistCount = gpu_min_ci_balanced(nbnxmGpu);
 348         if (debug)
 349         {
 350             fprintf(debug, "Neighbor-list balancing parameter: %d (auto-adjusted to the number of GPU multi-processors)\n",
 351                     minimumIlistCount);
 352         }
 353     }
 354
 355     return minimumIlistCount;
 356 }
 357
 358 std::unique_ptr<nonbonded_verlet_t>
 359 init_nb_verlet(const gmx::MDLogger     &mdlog,
 360                gmx_bool                 bFEP_NonBonded,
 361                const t_inputrec        *ir,
 362                const t_forcerec        *fr,
 363                const t_commrec         *cr,
 364                const gmx_hw_info_t     &hardwareInfo,
 365                const gmx_device_info_t *deviceInfo,
 366                const gmx_mtop_t        *mtop,
 367                matrix                   box,
 368                gmx_wallcycle           *wcycle)
 369 {
 370     const bool          emulateGpu = (getenv("GMX_EMULATE_GPU") != nullptr);
 371     const bool          useGpu     = deviceInfo != nullptr;
 372
 373     GMX_RELEASE_ASSERT(!(emulateGpu && useGpu), "When GPU emulation is active, there cannot be a GPU assignment");
 374
 375     NonbondedResource nonbondedResource;
 376     if (useGpu)
 377     {
 378         nonbondedResource = NonbondedResource::Gpu;
 379     }
 380     else if (emulateGpu)
 381     {
 382         nonbondedResource = NonbondedResource::EmulateGpu;
 383     }
 384     else
 385     {
 386         nonbondedResource = NonbondedResource::Cpu;
 387     }
 388
 389     Nbnxm::KernelSetup kernelSetup =
 390         pick_nbnxn_kernel(mdlog, fr->use_simd_kernels, hardwareInfo,
 391                           nonbondedResource, ir,
 392                           fr->bNonbonded);
 393
 394     const bool          haveMultipleDomains = (DOMAINDECOMP(cr) && cr->dd->nnodes > 1);
 395
 396     PairlistParams      pairlistParams(kernelSetup.kernelType,
 397                                        bFEP_NonBonded,
 398                                        ir->rlist,
 399                                        havePPDomainDecomposition(cr));
 400
 401     setupDynamicPairlistPruning(mdlog, ir, mtop, box, fr->ic,
 402                                 &pairlistParams);
 403
 404     int      enbnxninitcombrule;
 405     if (fr->ic->vdwtype == evdwCUT &&
 406         (fr->ic->vdw_modifier == eintmodNONE ||
 407          fr->ic->vdw_modifier == eintmodPOTSHIFT) &&
 408         getenv("GMX_NO_LJ_COMB_RULE") == nullptr)
 409     {
 410         /* Plain LJ cut-off: we can optimize with combination rules */
 411         enbnxninitcombrule = enbnxninitcombruleDETECT;
 412     }
 413     else if (fr->ic->vdwtype == evdwPME)
 414     {
 415         /* LJ-PME: we need to use a combination rule for the grid */
 416         if (fr->ljpme_combination_rule == eljpmeGEOM)
 417         {
 418             enbnxninitcombrule = enbnxninitcombruleGEOM;
 419         }
 420         else
 421         {
 422             enbnxninitcombrule = enbnxninitcombruleLB;
 423         }
 424     }
 425     else
 426     {
 427         /* We use a full combination matrix: no rule required */
 428         enbnxninitcombrule = enbnxninitcombruleNONE;
 429     }
 430
 431     auto pinPolicy = (useGpu ? gmx::PinningPolicy::PinnedIfSupported : gmx::PinningPolicy::CannotBePinned);
 432
 433     auto nbat      = std::make_unique<nbnxn_atomdata_t>(pinPolicy);
 434
 435     int  mimimumNumEnergyGroupNonbonded = ir->opts.ngener;
 436     if (ir->opts.ngener - ir->nwall == 1)
 437     {
 438         /* We have only one non-wall energy group, we do not need energy group
 439          * support in the non-bondeds kernels, since all non-bonded energy
 440          * contributions go to the first element of the energy group matrix.
 441          */
 442         mimimumNumEnergyGroupNonbonded = 1;
 443     }
 444     nbnxn_atomdata_init(mdlog,
 445                         nbat.get(),
 446                         kernelSetup.kernelType,
 447                         enbnxninitcombrule,
 448                         fr->ntype, fr->nbfp,
 449                         mimimumNumEnergyGroupNonbonded,
 450                         (useGpu || emulateGpu) ? 1 : gmx_omp_nthreads_get(emntNonbonded));
 451
 452     gmx_nbnxn_gpu_t *gpu_nbv = nullptr;
 453     int              minimumIlistCountForGpuBalancing = 0;
 454     if (useGpu)
 455     {
 456         /* init the NxN GPU data; the last argument tells whether we'll have
 457          * both local and non-local NB calculation on GPU */
 458         gpu_nbv = gpu_init(deviceInfo,
 459                            fr->ic,
 460                            pairlistParams,
 461                            nbat.get(),
 462                            cr->nodeid,
 463                            haveMultipleDomains);
 464
 465         minimumIlistCountForGpuBalancing = getMinimumIlistCountForGpuBalancing(gpu_nbv);
 466     }
 467
 468     auto pairlistSets =
 469         std::make_unique<PairlistSets>(pairlistParams,
 470                                        haveMultipleDomains,
 471                                        minimumIlistCountForGpuBalancing);
 472
 473     auto pairSearch =
 474         std::make_unique<PairSearch>(ir->ePBC,
 475                                      EI_TPI(ir->eI),
 476                                      DOMAINDECOMP(cr) ? &cr->dd->nc : nullptr,
 477                                      DOMAINDECOMP(cr) ? domdec_zones(cr->dd) : nullptr,
 478                                      pairlistParams.pairlistType,
 479                                      bFEP_NonBonded,
 480                                      gmx_omp_nthreads_get(emntPairsearch),
 481                                      pinPolicy);
 482
 483     return std::make_unique<nonbonded_verlet_t>(std::move(pairlistSets),
 484                                                 std::move(pairSearch),
 485                                                 std::move(nbat),
 486                                                 kernelSetup,
 487                                                 gpu_nbv,
 488                                                 wcycle);
 489 }
 490
 491 } // namespace Nbnxm
 492
 493 nonbonded_verlet_t::nonbonded_verlet_t(std::unique_ptr<PairlistSets>      pairlistSets,
 494                                        std::unique_ptr<PairSearch>        pairSearch,
 495                                        std::unique_ptr<nbnxn_atomdata_t>  nbat_in,
 496                                        const Nbnxm::KernelSetup          &kernelSetup,
 497                                        gmx_nbnxn_gpu_t                   *gpu_nbv_ptr,
 498                                        gmx_wallcycle                     *wcycle) :
 499     pairlistSets_(std::move(pairlistSets)),
 500     pairSearch_(std::move(pairSearch)),
 501     nbat(std::move(nbat_in)),
 502     kernelSetup_(kernelSetup),
 503     wcycle_(wcycle),
 504     gpu_nbv(gpu_nbv_ptr)
 505 {
 506     GMX_RELEASE_ASSERT(pairlistSets_, "Need valid pairlistSets");
 507     GMX_RELEASE_ASSERT(pairSearch_, "Need valid search object");
 508     GMX_RELEASE_ASSERT(nbat, "Need valid atomdata object");
 509 }
 510
 511 nonbonded_verlet_t::~nonbonded_verlet_t()
 512 {
 513     Nbnxm::gpu_free(gpu_nbv);
 514 }