clang-tidy modernize
[alexxy/gromacs.git] / src / gromacs / taskassignment / resourcedivision.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2015,2016,2017,2018, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*! \internal \file
36  * \brief Defines utility functionality for dividing resources and
37  * checking for consistency and usefulness.
38  *
39  * \author Mark Abraham <mark.j.abraham@gmail.com>
40  * \ingroup module_taskassignment
41  */
42
43 #include "gmxpre.h"
44
45 #include "resourcedivision.h"
46
47 #include "config.h"
48
49 #include <cstdlib>
50 #include <cstring>
51
52 #include <algorithm>
53
54 #include "gromacs/ewald/pme.h"
55 #include "gromacs/hardware/cpuinfo.h"
56 #include "gromacs/hardware/detecthardware.h"
57 #include "gromacs/hardware/hardwaretopology.h"
58 #include "gromacs/hardware/hw_info.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/mdlib/gmx_omp_nthreads.h"
61 #include "gromacs/mdtypes/commrec.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/utility/baseversion.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/logger.h"
70 #include "gromacs/utility/physicalnodecommunicator.h"
71 #include "gromacs/utility/stringutil.h"
72
73
74 /* DISCLAIMER: All the atom count and thread numbers below are heuristic.
75  * The real switching points will depend on the system simulation,
76  * the algorithms used and the hardware it's running on, as well as if there
77  * are other jobs running on the same machine. We try to take into account
78  * factors that have a large influence, such as recent Intel CPUs being
79  * much better at wide multi-threading. The remaining factors should
80  * (hopefully) have a small influence, such that the performance just before
81  * and after a switch point doesn't change too much.
82  */
83
84 //! Constant used to help minimize preprocessed code
85 static const bool bHasOmpSupport = GMX_OPENMP;
86
87 /*! \brief The minimum number of atoms per thread-MPI thread when GPUs
88  * are present. With fewer atoms than this, the number of thread-MPI
89  * ranks will get lowered.
90  */
91 static const int min_atoms_per_mpi_thread =  90;
92 /*! \brief The minimum number of atoms per GPU with thread-MPI
93  * active. With fewer atoms than this, the number of thread-MPI ranks
94  * will get lowered.
95  */
96 static const int min_atoms_per_gpu        = 900;
97
98 /**@{*/
99 /*! \brief Constants for implementing default divisions of threads */
100
101 /* TODO choose nthreads_omp based on hardware topology
102    when we have a hardware topology detection library */
103 /* First we consider the case of no MPI (1 MPI rank).
104  * In general, when running up to 8 threads, OpenMP should be faster.
105  * Note: on AMD Bulldozer we should avoid running OpenMP over two dies.
106  * On Intel>=Nehalem running OpenMP on a single CPU is always faster,
107  * even on two CPUs it's usually faster (but with many OpenMP threads
108  * it could be faster not to use HT, currently we always use HT).
109  * On Nehalem/Westmere we want to avoid running 16 threads over
110  * two CPUs with HT, so we need a limit<16; thus we use 12.
111  * A reasonable limit for Intel Sandy and Ivy bridge,
112  * not knowing the topology, is 16 threads.
113  * Below we check for Intel and AVX, which for now includes
114  * Sandy/Ivy Bridge, Has/Broadwell. By checking for AVX instead of
115  * model numbers we ensure also future Intel CPUs are covered.
116  */
117 const int nthreads_omp_faster_default   =  8;
118 const int nthreads_omp_faster_Nehalem   = 12;
119 const int nthreads_omp_faster_Intel_AVX = 16;
120 const int nthreads_omp_faster_AMD_Ryzen = 16;
121 /* For CPU only runs the fastest options are usually MPI or OpenMP only.
122  * With one GPU, using MPI only is almost never optimal, so we need to
123  * compare running pure OpenMP with combined MPI+OpenMP. This means higher
124  * OpenMP threads counts can still be ok. Multiplying the numbers above
125  * by a factor of 2 seems to be a good estimate.
126  */
127 const int nthreads_omp_faster_gpu_fac   =  2;
128
129 /* This is the case with MPI (2 or more MPI PP ranks).
130  * By default we will terminate with a fatal error when more than 8
131  * OpenMP thread are (indirectly) requested, since using less threads
132  * nearly always results in better performance.
133  * With thread-mpi and multiple GPUs or one GPU and too many threads
134  * we first try 6 OpenMP threads and then less until the number of MPI ranks
135  * is divisible by the number of GPUs.
136  */
137 #if GMX_OPENMP && GMX_MPI
138 const int nthreads_omp_mpi_ok_max              =  8;
139 const int nthreads_omp_mpi_ok_min_cpu          =  1;
140 #endif
141 const int nthreads_omp_mpi_ok_min_gpu          =  2;
142 const int nthreads_omp_mpi_target_max          =  6;
143
144 /**@}*/
145
146 /*! \brief Returns the maximum OpenMP thread count for which using a single MPI rank
147  * should be faster than using multiple ranks with the same total thread count.
148  */
149 static int nthreads_omp_faster(const gmx::CpuInfo &cpuInfo, gmx_bool bUseGPU)
150 {
151     int nth;
152
153     if (cpuInfo.vendor() == gmx::CpuInfo::Vendor::Intel &&
154         cpuInfo.feature(gmx::CpuInfo::Feature::X86_Avx))
155     {
156         nth = nthreads_omp_faster_Intel_AVX;
157     }
158     else if (gmx::cpuIsX86Nehalem(cpuInfo))
159     {
160         // Intel Nehalem
161         nth = nthreads_omp_faster_Nehalem;
162     }
163     else if (cpuInfo.vendor() == gmx::CpuInfo::Vendor::Amd && cpuInfo.family() >= 23)
164     {
165         // AMD Ryzen
166         nth = nthreads_omp_faster_AMD_Ryzen;
167     }
168     else
169     {
170         nth = nthreads_omp_faster_default;
171     }
172
173     if (bUseGPU)
174     {
175         nth *= nthreads_omp_faster_gpu_fac;
176     }
177
178     nth = std::min(nth, GMX_OPENMP_MAX_THREADS);
179
180     return nth;
181 }
182
183 /*! \brief Returns that maximum OpenMP thread count that passes the efficiency check */
184 gmx_unused static int nthreads_omp_efficient_max(int gmx_unused       nrank,
185                                                  const gmx::CpuInfo  &cpuInfo,
186                                                  gmx_bool             bUseGPU)
187 {
188 #if GMX_OPENMP && GMX_MPI
189     if (nrank > 1)
190     {
191         return nthreads_omp_mpi_ok_max;
192     }
193     else
194 #endif
195     {
196         return nthreads_omp_faster(cpuInfo, bUseGPU);
197     }
198 }
199
200 /*! \brief Return the number of thread-MPI ranks to use.
201  * This is chosen such that we can always obey our own efficiency checks.
202  */
203 gmx_unused static int get_tmpi_omp_thread_division(const gmx_hw_info_t *hwinfo,
204                                                    const gmx_hw_opt_t  &hw_opt,
205                                                    int                  nthreads_tot,
206                                                    int                  ngpu)
207 {
208     int                 nrank;
209     const gmx::CpuInfo &cpuInfo = *hwinfo->cpuInfo;
210
211     GMX_RELEASE_ASSERT(nthreads_tot > 0, "There must be at least one thread per rank");
212
213     /* There are no separate PME nodes here, as we ensured in
214      * check_and_update_hw_opt that nthreads_tmpi>0 with PME nodes
215      * and a conditional ensures we would not have ended up here.
216      * Note that separate PME nodes might be switched on later.
217      */
218     if (ngpu > 0)
219     {
220         if (hw_opt.nthreads_omp > 0)
221         {
222             /* In this case it is unclear if we should use 1 rank per GPU
223              * or more or less, so we require also setting the number of ranks.
224              */
225             gmx_fatal(FARGS, "When using GPUs, setting the number of OpenMP threads without specifying the number of ranks can lead to conflicting demands. Please specify the number of thread-MPI ranks as well (option -ntmpi).");
226         }
227
228         nrank = ngpu;
229
230         /* When the user sets nthreads_omp, we can end up oversubscribing CPU cores
231          * if we simply start as many ranks as GPUs. To avoid this, we start as few
232          * tMPI ranks as necessary to avoid oversubscription and instead leave GPUs idle.
233          * If the user does not set the number of OpenMP threads, nthreads_omp==0 and
234          * this code has no effect.
235          */
236         GMX_RELEASE_ASSERT(hw_opt.nthreads_omp >= 0, "nthreads_omp is negative, but previous checks should have prevented this");
237         while (nrank*hw_opt.nthreads_omp > hwinfo->nthreads_hw_avail && nrank > 1)
238         {
239             nrank--;
240         }
241
242         if (nthreads_tot < nrank)
243         {
244             /* #thread < #gpu is very unlikely, but if so: waste gpu(s) */
245             nrank = nthreads_tot;
246         }
247         else if (nthreads_tot > nthreads_omp_faster(cpuInfo, ngpu > 0) ||
248                  (ngpu > 1 && nthreads_tot/ngpu > nthreads_omp_mpi_target_max))
249         {
250             /* The high OpenMP thread count will likely result in sub-optimal
251              * performance. Increase the rank count to reduce the thread count
252              * per rank. This will lead to GPU sharing by MPI ranks/threads.
253              */
254             int nshare;
255
256             /* Increase the rank count as long as have we more than 6 OpenMP
257              * threads per rank or the number of hardware threads is not
258              * divisible by the rank count. Don't go below 2 OpenMP threads.
259              */
260             nshare = 1;
261             do
262             {
263                 nshare++;
264                 nrank = ngpu*nshare;
265             }
266             while (nthreads_tot/nrank > nthreads_omp_mpi_target_max ||
267                    (nthreads_tot/(ngpu*(nshare + 1)) >= nthreads_omp_mpi_ok_min_gpu && nthreads_tot % nrank != 0));
268         }
269     }
270     else if (hw_opt.nthreads_omp > 0)
271     {
272         /* Here we could oversubscribe, when we do, we issue a warning later */
273         nrank = std::max(1, nthreads_tot/hw_opt.nthreads_omp);
274     }
275     else
276     {
277         if (nthreads_tot <= nthreads_omp_faster(cpuInfo, ngpu > 0))
278         {
279             /* Use pure OpenMP parallelization */
280             nrank = 1;
281         }
282         else
283         {
284             /* Don't use OpenMP parallelization */
285             nrank = nthreads_tot;
286         }
287     }
288
289     return nrank;
290 }
291
292 //! Return whether hyper threading is enabled.
293 static bool
294 gmxSmtIsEnabled(const gmx::HardwareTopology &hwTop)
295 {
296     return (hwTop.supportLevel() >= gmx::HardwareTopology::SupportLevel::Basic && hwTop.machine().sockets[0].cores[0].hwThreads.size() > 1);
297 }
298
299 namespace
300 {
301
302 //! Handles checks for algorithms that must use a single rank.
303 class SingleRankChecker
304 {
305     public:
306         //! Constructor
307         SingleRankChecker() : value_(false) {}
308         /*! \brief Call this function for each possible condition
309             under which a single rank is required, along with a string
310             describing the constraint when it is applied. */
311         void applyConstraint(bool condition, const char *description)
312         {
313             if (condition)
314             {
315                 value_ = true;
316                 reasons_.push_back(gmx::formatString("%s only supports a single rank.", description));
317             }
318         }
319         //! After applying any conditions, is a single rank required?
320         bool mustUseOneRank() const
321         {
322             return value_;
323         }
324         /*! \brief Return a formatted string to use when writing a
325             message when a single rank is required, (or empty if no
326             constraint exists.) */
327         std::string getMessage() const
328         {
329             return formatAndJoin(reasons_, "\n", gmx::IdentityFormatter());
330         }
331     private:
332         bool                     value_;
333         std::vector<std::string> reasons_;
334 };
335
336 } // namespace
337
338 /* Get the number of MPI ranks to use for thread-MPI based on how many
339  * were requested, which algorithms we're using,
340  * and how many particles there are.
341  * At the point we have already called check_and_update_hw_opt.
342  * Thus all options should be internally consistent and consistent
343  * with the hardware, except that ntmpi could be larger than #GPU.
344  */
345 int get_nthreads_mpi(const gmx_hw_info_t    *hwinfo,
346                      gmx_hw_opt_t           *hw_opt,
347                      const std::vector<int> &gpuIdsToUse,
348                      bool                    nonbondedOnGpu,
349                      bool                    pmeOnGpu,
350                      const t_inputrec       *inputrec,
351                      const gmx_mtop_t       *mtop,
352                      const gmx::MDLogger    &mdlog,
353                      bool                    doMembed)
354 {
355     int                          nthreads_hw, nthreads_tot_max, nrank, ngpu;
356     int                          min_atoms_per_mpi_rank;
357
358     const gmx::CpuInfo          &cpuInfo = *hwinfo->cpuInfo;
359     const gmx::HardwareTopology &hwTop   = *hwinfo->hardwareTopology;
360
361     if (pmeOnGpu)
362     {
363         GMX_RELEASE_ASSERT((EEL_PME(inputrec->coulombtype) || EVDW_PME(inputrec->vdwtype)) &&
364                            pme_gpu_supports_build(nullptr) && pme_gpu_supports_input(inputrec, nullptr),
365                            "PME can't be on GPUs unless we are using PME");
366
367         // PME on GPUs supports a single PME rank with PP running on the same or few other ranks.
368         // For now, let's treat separate PME GPU rank as opt-in.
369         if (hw_opt->nthreads_tmpi < 1)
370         {
371             return 1;
372         }
373     }
374
375     {
376         /* Check if an algorithm does not support parallel simulation.  */
377         // TODO This might work better if e.g. implemented algorithms
378         // had to define a function that returns such requirements,
379         // and a description string.
380         SingleRankChecker checker;
381         checker.applyConstraint(inputrec->eI == eiLBFGS, "L-BFGS minimization");
382         checker.applyConstraint(inputrec->coulombtype == eelEWALD, "Plain Ewald electrostatics");
383         checker.applyConstraint(doMembed, "Membrane embedding");
384         bool useOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
385         checker.applyConstraint(useOrientationRestraints, "Orientation restraints");
386         if (checker.mustUseOneRank())
387         {
388             std::string message = checker.getMessage();
389             if (hw_opt->nthreads_tmpi > 1)
390             {
391                 gmx_fatal(FARGS, "%s However, you asked for more than 1 thread-MPI rank, so mdrun cannot continue. Choose a single rank, or a different algorithm.", message.c_str());
392             }
393             GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted("%s Choosing to use only a single thread-MPI rank.", message.c_str());
394             return 1;
395         }
396     }
397
398     if (hw_opt->nthreads_tmpi > 0)
399     {
400         /* Trivial, return the user's choice right away */
401         return hw_opt->nthreads_tmpi;
402     }
403
404     // Now implement automatic selection of number of thread-MPI ranks
405     nthreads_hw = hwinfo->nthreads_hw_avail;
406
407     if (nthreads_hw <= 0)
408     {
409         /* This should normally not happen, but if it does, we handle it */
410         gmx_fatal(FARGS, "The number of available hardware threads can not be detected, please specify the number of MPI ranks and the number of OpenMP threads (if supported) manually with options -ntmpi and -ntomp, respectively");
411     }
412
413     /* How many total (#tMPI*#OpenMP) threads can we start? */
414     if (hw_opt->nthreads_tot > 0)
415     {
416         nthreads_tot_max = hw_opt->nthreads_tot;
417     }
418     else
419     {
420         nthreads_tot_max = nthreads_hw;
421     }
422
423     /* nonbondedOnGpu might be false e.g. because this simulation uses
424      * the group scheme, or is a rerun with energy groups. */
425     ngpu = (nonbondedOnGpu ? static_cast<int>(gpuIdsToUse.size()) : 0);
426
427     if (inputrec->cutoff_scheme == ecutsGROUP)
428     {
429         /* We checked this before, but it doesn't hurt to do it once more */
430         GMX_RELEASE_ASSERT(hw_opt->nthreads_omp == 1, "The group scheme only supports one OpenMP thread per rank");
431     }
432
433     nrank =
434         get_tmpi_omp_thread_division(hwinfo, *hw_opt, nthreads_tot_max, ngpu);
435
436     if (inputrec->eI == eiNM || EI_TPI(inputrec->eI))
437     {
438         /* Dims/steps are divided over the nodes iso splitting the atoms.
439          * With NM we can't have more ranks than #atoms*#dim. With TPI it's
440          * unlikely we have fewer atoms than ranks, and if so, communication
441          * would become a bottleneck, so we set the limit to 1 atom/rank.
442          */
443         min_atoms_per_mpi_rank = 1;
444     }
445     else
446     {
447         if (ngpu >= 1)
448         {
449             min_atoms_per_mpi_rank = min_atoms_per_gpu;
450         }
451         else
452         {
453             min_atoms_per_mpi_rank = min_atoms_per_mpi_thread;
454         }
455     }
456
457     if (mtop->natoms/nrank < min_atoms_per_mpi_rank)
458     {
459         int nrank_new;
460
461         /* the rank number was chosen automatically, but there are too few
462            atoms per rank, so we need to reduce the rank count */
463         nrank_new = std::max(1, mtop->natoms/min_atoms_per_mpi_rank);
464
465         /* Avoid partial use of Hyper-Threading */
466         if (gmxSmtIsEnabled(hwTop) &&
467             nrank_new > nthreads_hw/2 && nrank_new < nthreads_hw)
468         {
469             nrank_new = nthreads_hw/2;
470         }
471
472         /* If the user specified the total thread count, ensure this is
473          * divisible by the number of ranks.
474          * It is quite likely that we have too many total threads compared
475          * to the size of the system, but if the user asked for this many
476          * threads we should respect that.
477          */
478         while (hw_opt->nthreads_tot > 0 &&
479                hw_opt->nthreads_tot % nrank_new != 0)
480         {
481             nrank_new--;
482         }
483
484         /* Avoid large prime numbers in the rank count */
485         if (nrank_new >= 6)
486         {
487             /* Use only 6,8,10 with additional factors of 2 */
488             int fac;
489
490             fac = 2;
491             while (3*fac*2 <= nrank_new)
492             {
493                 fac *= 2;
494             }
495
496             nrank_new = (nrank_new/fac)*fac;
497         }
498         else
499         {
500             /* Avoid 5, since small system won't fit 5 domains along
501              * a dimension. This might lead to waisting some cores, but this
502              * will have a small impact in this regime of very small systems.
503              */
504             if (nrank_new == 5)
505             {
506                 nrank_new = 4;
507             }
508         }
509
510         nrank = nrank_new;
511
512         /* We reduced the number of tMPI ranks, which means we might violate
513          * our own efficiency checks if we simply use all hardware threads.
514          */
515         if (bHasOmpSupport && hw_opt->nthreads_omp <= 0 && hw_opt->nthreads_tot <= 0)
516         {
517             /* The user set neither the total nor the OpenMP thread count,
518              * we should use all hardware threads, unless we will violate
519              * our own efficiency limitation on the thread count.
520              */
521             int  nt_omp_max;
522
523             nt_omp_max = nthreads_omp_efficient_max(nrank, cpuInfo, ngpu >= 1);
524
525             if (nrank*nt_omp_max < hwinfo->nthreads_hw_avail)
526             {
527                 /* Limit the number of OpenMP threads to start */
528                 hw_opt->nthreads_omp = nt_omp_max;
529             }
530         }
531
532         fprintf(stderr, "\n");
533         fprintf(stderr, "NOTE: Parallelization is limited by the small number of atoms,\n");
534         fprintf(stderr, "      only starting %d thread-MPI ranks.\n", nrank);
535         fprintf(stderr, "      You can use the -nt and/or -ntmpi option to optimize the number of threads.\n\n");
536     }
537
538     return nrank;
539 }
540
541
542 void check_resource_division_efficiency(const gmx_hw_info_t *hwinfo,
543                                         bool                 willUsePhysicalGpu,
544                                         gmx_bool             bNtOmpOptionSet,
545                                         t_commrec           *cr,
546                                         const gmx::MDLogger &mdlog)
547 {
548     GMX_UNUSED_VALUE(hwinfo);
549 #if GMX_OPENMP && GMX_MPI
550     int         nth_omp_min, nth_omp_max;
551     char        buf[1000];
552 #if GMX_THREAD_MPI
553     const char *mpi_option = " (option -ntmpi)";
554 #else
555     const char *mpi_option = "";
556 #endif
557
558     /* This function should be called after thread-MPI (when configured) and
559      * OpenMP have been initialized. Check that here.
560      */
561 #if GMX_THREAD_MPI
562     GMX_RELEASE_ASSERT(nthreads_omp_faster_default >= nthreads_omp_mpi_ok_max, "Inconsistent OpenMP thread count default values");
563 #endif
564     GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) >= 1, "Must have at least one OpenMP thread");
565
566     nth_omp_min = gmx_omp_nthreads_get(emntDefault);
567     nth_omp_max = gmx_omp_nthreads_get(emntDefault);
568
569     bool anyRankIsUsingGpus = willUsePhysicalGpu;
570     /* Thread-MPI seems to have a bug with reduce on 1 node, so use a cond. */
571     if (cr->nnodes + cr->npmenodes > 1)
572     {
573         int count[3], count_max[3];
574
575         count[0] = -nth_omp_min;
576         count[1] =  nth_omp_max;
577         count[2] =  willUsePhysicalGpu;
578
579         MPI_Allreduce(count, count_max, 3, MPI_INT, MPI_MAX, cr->mpi_comm_mysim);
580
581         /* In case of an inhomogeneous run setup we use the maximum counts */
582         nth_omp_min        = -count_max[0];
583         nth_omp_max        =  count_max[1];
584         anyRankIsUsingGpus = count_max[2] > 0;
585     }
586
587     int nthreads_omp_mpi_ok_min;
588
589     if (!anyRankIsUsingGpus)
590     {
591         nthreads_omp_mpi_ok_min = nthreads_omp_mpi_ok_min_cpu;
592     }
593     else
594     {
595         /* With GPUs we set the minimum number of OpenMP threads to 2 to catch
596          * cases where the user specifies #ranks == #cores.
597          */
598         nthreads_omp_mpi_ok_min = nthreads_omp_mpi_ok_min_gpu;
599     }
600
601     if (DOMAINDECOMP(cr))
602     {
603         if (nth_omp_max < nthreads_omp_mpi_ok_min ||
604             nth_omp_max > nthreads_omp_mpi_ok_max)
605         {
606             /* Note that we print target_max here, not ok_max */
607             sprintf(buf, "Your choice of number of MPI ranks and amount of resources results in using %d OpenMP threads per rank, which is most likely inefficient. The optimum is usually between %d and %d threads per rank.",
608                     nth_omp_max,
609                     nthreads_omp_mpi_ok_min,
610                     nthreads_omp_mpi_target_max);
611
612             if (bNtOmpOptionSet)
613             {
614                 GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted("NOTE: %s", buf);
615             }
616             else
617             {
618                 /* This fatal error, and the one below, is nasty, but it's
619                  * probably the only way to ensure that all users don't waste
620                  * a lot of resources, since many users don't read logs/stderr.
621                  */
622                 gmx_fatal(FARGS, "%s If you want to run with this setup, specify the -ntomp option. But we suggest to change the number of MPI ranks%s.", buf, mpi_option);
623             }
624         }
625     }
626 #else /* GMX_OPENMP && GMX_MPI */
627       /* No OpenMP and/or MPI: it doesn't make much sense to check */
628     GMX_UNUSED_VALUE(bNtOmpOptionSet);
629     GMX_UNUSED_VALUE(willUsePhysicalGpu);
630     GMX_UNUSED_VALUE(cr);
631     /* Check if we have more than 1 physical core, if detected,
632      * or more than 1 hardware thread if physical cores were not detected.
633      */
634     if (!GMX_OPENMP && !GMX_MPI && hwinfo->hardwareTopology->numberOfCores() > 1)
635     {
636         GMX_LOG(mdlog.warning).asParagraph().appendText("NOTE: GROMACS was compiled without OpenMP and (thread-)MPI support, can only use a single CPU core");
637     }
638 #endif /* GMX_OPENMP && GMX_MPI */
639 }
640
641
642 //! Dump a \c hw_opt to \c fp.
643 static void print_hw_opt(FILE *fp, const gmx_hw_opt_t *hw_opt)
644 {
645     fprintf(fp, "hw_opt: nt %d ntmpi %d ntomp %d ntomp_pme %d gpu_id '%s' gputasks '%s'\n",
646             hw_opt->nthreads_tot,
647             hw_opt->nthreads_tmpi,
648             hw_opt->nthreads_omp,
649             hw_opt->nthreads_omp_pme,
650             hw_opt->gpuIdsAvailable.c_str(),
651             hw_opt->userGpuTaskAssignment.c_str());
652 }
653
654 void check_and_update_hw_opt_1(const gmx::MDLogger &mdlog,
655                                gmx_hw_opt_t        *hw_opt,
656                                const t_commrec     *cr,
657                                int                  nPmeRanks)
658 {
659     /* Currently hw_opt only contains default settings or settings supplied
660      * by the user on the command line.
661      */
662     if (hw_opt->nthreads_omp < 0)
663     {
664         gmx_fatal(FARGS, "The number of OpenMP threads supplied on the command line is %d, which is negative and not allowed", hw_opt->nthreads_omp);
665     }
666
667     /* Check for OpenMP settings stored in environment variables, which can
668      * potentially be different on different MPI ranks.
669      */
670     gmx_omp_nthreads_read_env(mdlog, &hw_opt->nthreads_omp);
671
672     /* Check restrictions on the user supplied options before modifying them.
673      * TODO: Put the user values in a const struct and preserve them.
674      */
675 #if !GMX_THREAD_MPI
676     if (hw_opt->nthreads_tot > 0)
677     {
678         gmx_fatal(FARGS, "Setting the total number of threads is only supported with thread-MPI and GROMACS was compiled without thread-MPI");
679     }
680     if (hw_opt->nthreads_tmpi > 0)
681     {
682         gmx_fatal(FARGS, "Setting the number of thread-MPI ranks is only supported with thread-MPI and GROMACS was compiled without thread-MPI");
683     }
684 #endif
685
686     /* With thread-MPI the master thread sets hw_opt->totNumThreadsIsAuto.
687      * The other threads receive a partially processed hw_opt from the master
688      * thread and should not set hw_opt->totNumThreadsIsAuto again.
689      */
690     if (!GMX_THREAD_MPI || SIMMASTER(cr))
691     {
692         /* Check if mdrun is free to choose the total number of threads */
693         hw_opt->totNumThreadsIsAuto = (hw_opt->nthreads_omp == 0 && hw_opt->nthreads_omp_pme == 0 && hw_opt->nthreads_tot == 0);
694     }
695
696     if (bHasOmpSupport)
697     {
698         /* Check restrictions on PME thread related options set by the user */
699
700         if (hw_opt->nthreads_omp_pme > 0 && hw_opt->nthreads_omp <= 0)
701         {
702             gmx_fatal(FARGS, "You need to specify -ntomp in addition to -ntomp_pme");
703         }
704
705         if (hw_opt->nthreads_omp_pme >= 1 &&
706             hw_opt->nthreads_omp_pme != hw_opt->nthreads_omp &&
707             nPmeRanks <= 0)
708         {
709             /* This can result in a fatal error on many MPI ranks,
710              * but since the thread count can differ per rank,
711              * we can't easily avoid this.
712              */
713             gmx_fatal(FARGS, "You need to explicitly specify the number of PME ranks (-npme) when using different number of OpenMP threads for PP and PME ranks");
714         }
715     }
716     else
717     {
718         /* GROMACS was configured without OpenMP support */
719
720         if (hw_opt->nthreads_omp > 1 || hw_opt->nthreads_omp_pme > 1)
721         {
722             gmx_fatal(FARGS, "More than 1 OpenMP thread requested, but GROMACS was compiled without OpenMP support");
723         }
724         hw_opt->nthreads_omp     = 1;
725         hw_opt->nthreads_omp_pme = 1;
726     }
727
728     if (hw_opt->nthreads_tot > 0 && hw_opt->nthreads_omp_pme <= 0)
729     {
730         /* We have the same number of OpenMP threads for PP and PME ranks,
731          * thus we can perform several consistency checks.
732          */
733         if (hw_opt->nthreads_tmpi > 0 &&
734             hw_opt->nthreads_omp > 0 &&
735             hw_opt->nthreads_tot != hw_opt->nthreads_tmpi*hw_opt->nthreads_omp)
736         {
737             gmx_fatal(FARGS, "The total number of threads requested (%d) does not match the thread-MPI ranks (%d) times the OpenMP threads (%d) requested",
738                       hw_opt->nthreads_tot, hw_opt->nthreads_tmpi, hw_opt->nthreads_omp);
739         }
740
741         if (hw_opt->nthreads_tmpi > 0 &&
742             hw_opt->nthreads_tot % hw_opt->nthreads_tmpi != 0)
743         {
744             gmx_fatal(FARGS, "The total number of threads requested (%d) is not divisible by the number of thread-MPI ranks requested (%d)",
745                       hw_opt->nthreads_tot, hw_opt->nthreads_tmpi);
746         }
747
748         if (hw_opt->nthreads_omp > 0 &&
749             hw_opt->nthreads_tot % hw_opt->nthreads_omp != 0)
750         {
751             gmx_fatal(FARGS, "The total number of threads requested (%d) is not divisible by the number of OpenMP threads requested (%d)",
752                       hw_opt->nthreads_tot, hw_opt->nthreads_omp);
753         }
754     }
755
756     if (hw_opt->nthreads_tot > 0)
757     {
758         if (hw_opt->nthreads_omp > hw_opt->nthreads_tot)
759         {
760             gmx_fatal(FARGS, "You requested %d OpenMP threads with %d total threads. Choose a total number of threads that is a multiple of the number of OpenMP threads.",
761                       hw_opt->nthreads_omp, hw_opt->nthreads_tot);
762         }
763
764         if (hw_opt->nthreads_tmpi > hw_opt->nthreads_tot)
765         {
766             gmx_fatal(FARGS, "You requested %d thread-MPI ranks with %d total threads. Choose a total number of threads that is a multiple of the number of thread-MPI ranks.",
767                       hw_opt->nthreads_tmpi, hw_opt->nthreads_tot);
768         }
769     }
770
771     if (debug)
772     {
773         print_hw_opt(debug, hw_opt);
774     }
775
776     /* Asserting this simplifies the hardware resource division later
777      * on. */
778     GMX_RELEASE_ASSERT(!(hw_opt->nthreads_omp_pme >= 1 && hw_opt->nthreads_omp <= 0),
779                        "PME thread count should only be set when the normal thread count is also set");
780 }
781
782 void check_and_update_hw_opt_2(gmx_hw_opt_t *hw_opt,
783                                int           cutoff_scheme)
784 {
785     if (cutoff_scheme == ecutsGROUP)
786     {
787         /* We only have OpenMP support for PME only nodes */
788         if (hw_opt->nthreads_omp > 1)
789         {
790             gmx_fatal(FARGS, "OpenMP threads have been requested with cut-off scheme %s, but these are only supported with cut-off scheme %s",
791                       ecutscheme_names[cutoff_scheme],
792                       ecutscheme_names[ecutsVERLET]);
793         }
794         hw_opt->nthreads_omp = 1;
795     }
796 }
797
798 void checkAndUpdateRequestedNumOpenmpThreads(gmx_hw_opt_t         *hw_opt,
799                                              const gmx_hw_info_t  &hwinfo,
800                                              const t_commrec      *cr,
801                                              const gmx_multisim_t *ms,
802                                              int                   numRanksOnThisNode,
803                                              PmeRunMode            pmeRunMode,
804                                              const gmx_mtop_t     &mtop)
805 {
806 #if GMX_THREAD_MPI
807     GMX_RELEASE_ASSERT(hw_opt->nthreads_tmpi >= 1, "Must have at least one thread-MPI rank");
808
809     /* If the user set the total number of threads on the command line
810      * and did not specify the number of OpenMP threads, set the latter here.
811      */
812     if (hw_opt->nthreads_tot > 0 && hw_opt->nthreads_omp <= 0)
813     {
814         hw_opt->nthreads_omp = hw_opt->nthreads_tot/hw_opt->nthreads_tmpi;
815
816         if (!bHasOmpSupport && hw_opt->nthreads_omp > 1)
817         {
818             gmx_fatal(FARGS, "You (indirectly) asked for OpenMP threads by setting -nt > -ntmpi, but GROMACS was compiled without OpenMP support");
819         }
820     }
821 #endif
822
823     /* With both non-bonded and PME on GPU, the work left on the CPU is often
824      * (much) slower with SMT than without SMT. This is mostly the case with
825      * few atoms per core. Thus, if the number of threads is set to auto,
826      * we turn off SMT in that case. Note that PME on GPU implies that also
827      * the non-bonded are computed on the GPU.
828      * We only need to do this when the number of hardware theads is larger
829      * than the number of cores. Note that a queuing system could limit
830      * the number of hardware threads available, but we are not trying to be
831      * too smart here in that case.
832      */
833     /* The thread reduction and synchronization costs go up roughy quadratically
834      * with the threads count, so we apply a threshold quadratic in #cores.
835      * Also more cores per GPU usually means the CPU gets faster than the GPU.
836      * The number 1000 atoms per core^2 is a reasonable threshold
837      * for Intel x86 and AMD Threadripper.
838      */
839     constexpr int c_numAtomsPerCoreSquaredSmtThreshold = 1000;
840
841     /* Prepare conditions for deciding if we should disable SMT.
842      * We currently only limit SMT for simulations using a single rank.
843      * TODO: Consider limiting also for multi-rank simulations.
844      */
845     bool canChooseNumOpenmpThreads      = (bHasOmpSupport && hw_opt->nthreads_omp <= 0);
846     bool haveSmtSupport                 = (hwinfo.hardwareTopology->supportLevel() >= gmx::HardwareTopology::SupportLevel::Basic &&
847                                            hwinfo.hardwareTopology->machine().logicalProcessorCount > hwinfo.hardwareTopology->numberOfCores());
848     bool simRunsSingleRankNBAndPmeOnGpu = (cr->nnodes == 1 && pmeRunMode == PmeRunMode::GPU);
849
850     if (canChooseNumOpenmpThreads && haveSmtSupport &&
851         simRunsSingleRankNBAndPmeOnGpu)
852     {
853         /* Note that the queing system might have limited us from using
854          * all detected ncore_tot physical cores. We are currently not
855          * checking for that here.
856          */
857         int numRanksTot     = cr->nnodes*(isMultiSim(ms) ? ms->nsim : 1);
858         int numAtomsPerRank = mtop.natoms/cr->nnodes;
859         int numCoresPerRank = hwinfo.ncore_tot/numRanksTot;
860         if (numAtomsPerRank < c_numAtomsPerCoreSquaredSmtThreshold*gmx::square(numCoresPerRank))
861         {
862             /* Choose one OpenMP thread per physical core */
863             hw_opt->nthreads_omp = std::max(1, hwinfo.hardwareTopology->numberOfCores()/numRanksOnThisNode);
864         }
865     }
866
867     GMX_RELEASE_ASSERT(bHasOmpSupport || hw_opt->nthreads_omp == 1, "Without OpenMP support, only one thread per rank can be used");
868
869     /* We are done with updating nthreads_omp, we can set nthreads_omp_pme */
870     if (hw_opt->nthreads_omp_pme <= 0 && hw_opt->nthreads_omp > 0)
871     {
872         hw_opt->nthreads_omp_pme = hw_opt->nthreads_omp;
873     }
874
875     if (debug)
876     {
877         print_hw_opt(debug, hw_opt);
878     }
879 }
880
881 namespace gmx
882 {
883
884 void checkHardwareOversubscription(int                             numThreadsOnThisRank,
885                                    int                             rank,
886                                    const HardwareTopology         &hwTop,
887                                    const PhysicalNodeCommunicator &comm,
888                                    const MDLogger                 &mdlog)
889 {
890     if (hwTop.supportLevel() < HardwareTopology::SupportLevel::LogicalProcessorCount)
891     {
892         /* There is nothing we can check */
893         return;
894     }
895
896     int numRanksOnThisNode   = comm.size_;
897     int numThreadsOnThisNode = numThreadsOnThisRank;
898     /* Avoid MPI calls with uninitialized thread-MPI communicators */
899     if (comm.size_ > 1)
900     {
901 #if GMX_MPI
902         /* Count the threads within this physical node */
903         MPI_Allreduce(&numThreadsOnThisRank, &numThreadsOnThisNode, 1, MPI_INT, MPI_SUM, comm.comm_);
904 #endif
905     }
906
907     if (numThreadsOnThisNode > hwTop.machine().logicalProcessorCount)
908     {
909         std::string mesg = "WARNING: ";
910         if (GMX_LIB_MPI)
911         {
912             mesg += formatString("On rank %d: o", rank);
913         }
914         else
915         {
916             mesg += "O";
917         }
918         mesg     += formatString("versubscribing the available %d logical CPU cores", hwTop.machine().logicalProcessorCount);
919         if (GMX_LIB_MPI)
920         {
921             mesg += " per node";
922         }
923         mesg     += formatString(" with %d ", numThreadsOnThisNode);
924         if (numRanksOnThisNode == numThreadsOnThisNode)
925         {
926             if (GMX_THREAD_MPI)
927             {
928                 mesg += "thread-MPI threads.";
929             }
930             else
931             {
932                 mesg += "MPI processes.";
933             }
934         }
935         else
936         {
937             mesg += "threads.";
938         }
939         mesg     += "\n         This will cause considerable performance loss.";
940         /* Note that only the master rank logs to stderr and only ranks
941          * with an open log file write to log.
942          * TODO: When we have a proper parallel logging framework,
943          *       the framework should add the rank and node numbers.
944          */
945         GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted("%s", mesg.c_str());
946     }
947 }
948
949 }  // namespace gmx