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