api/nblib/listed_forces/helpers.hpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, 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
  37  * Helper data structures and utility functions for the nblib force calculator.
  38  * Intended for internal use.
  39  *
  40  * \author Victor Holanda <victor.holanda@cscs.ch>
  41  * \author Joe Jordan <ejjordan@kth.se>
  42  * \author Prashanth Kanduri <kanduri@cscs.ch>
  43  * \author Sebastian Keller <keller@cscs.ch>
  44  * \author Artem Zhmurov <zhmurov@gmail.com>
  45  */
  46
  47 #ifndef NBLIB_LISTEDFORCSES_HELPERS_HPP
  48 #define NBLIB_LISTEDFORCSES_HELPERS_HPP
  49
  50 #include <unordered_map>
  51
  52 #include "nblib/pbc.hpp"
  53 #include "definitions.h"
  54 #include "nblib/util/util.hpp"
  55
  56 #define NBLIB_ALWAYS_INLINE __attribute((always_inline))
  57
  58 namespace nblib
  59 {
  60
  61 namespace detail
  62 {
  63 template<class T>
  64 inline void gmxRVecZeroWorkaround([[maybe_unused]] T& value)
  65 {
  66 }
  67
  68 template<>
  69 inline void gmxRVecZeroWorkaround<gmx::RVec>(gmx::RVec& value)
  70 {
  71     for (int i = 0; i < dimSize; ++i)
  72     {
  73         value[i] = 0;
  74     }
  75 }
  76 } // namespace detail
  77
  78 /*! \internal \brief object to store forces for multithreaded listed forces computation
  79  *
  80  */
  81 template<class T>
  82 class ForceBuffer
  83 {
  84     using HashMap = std::unordered_map<int, T>;
  85
  86 public:
  87     ForceBuffer() : rangeStart(0), rangeEnd(0) { }
  88
  89     ForceBuffer(T* mbuf, int rs, int re) :
  90         masterForceBuffer(mbuf),
  91         rangeStart(rs),
  92         rangeEnd(re)
  93     {
  94     }
  95
  96     void clear() { outliers.clear(); }
  97
  98     inline NBLIB_ALWAYS_INLINE T& operator[](int i)
  99     {
 100         if (i >= rangeStart && i < rangeEnd)
 101         {
 102             return masterForceBuffer[i];
 103         }
 104         else
 105         {
 106             if (outliers.count(i) == 0)
 107             {
 108                 T zero = T();
 109                 // if T = gmx::RVec, need to explicitly initialize it to zeros
 110                 detail::gmxRVecZeroWorkaround(zero);
 111                 outliers[i] = zero;
 112             }
 113             return outliers[i];
 114         }
 115     }
 116
 117     typename HashMap::const_iterator begin() { return outliers.begin(); }
 118     typename HashMap::const_iterator end() { return outliers.end(); }
 119
 120     [[nodiscard]] bool inRange(int index) const { return (index >= rangeStart && index < rangeEnd); }
 121
 122 private:
 123     T*  masterForceBuffer;
 124     int rangeStart;
 125     int rangeEnd;
 126
 127     HashMap outliers;
 128 };
 129
 130 namespace detail
 131 {
 132
 133 static int computeChunkIndex(int index, int totalRange, int nSplits)
 134 {
 135     if (totalRange < nSplits)
 136     {
 137         // if there's more threads than particles
 138         return index;
 139     }
 140
 141     int splitLength = totalRange / nSplits;
 142     return std::min(index / splitLength, nSplits - 1);
 143 }
 144
 145 } // namespace detail
 146
 147
 148 /*! \internal \brief splits an interaction tuple into nSplits interaction tuples
 149  *
 150  * \param interactions
 151  * \param totalRange the number of particle sequence coordinates
 152  * \param nSplits number to divide the total work by
 153  * \return
 154  */
 155 inline
 156 std::vector<ListedInteractionData> splitListedWork(const ListedInteractionData& interactions,
 157                                                    int                          totalRange,
 158                                                    int                          nSplits)
 159 {
 160     std::vector<ListedInteractionData> workDivision(nSplits);
 161
 162     auto splitOneElement = [totalRange, nSplits, &workDivision](const auto& inputElement) {
 163         // the index of inputElement in the ListedInteractionsTuple
 164         constexpr int elementIndex =
 165                 FindIndex<std::decay_t<decltype(inputElement)>, ListedInteractionData>{};
 166
 167         // for now, copy all parameters to each split
 168         // Todo: extract only the parameters needed for this split
 169         for (auto& workDivisionSplit : workDivision)
 170         {
 171             std::get<elementIndex>(workDivisionSplit).parameters = inputElement.parameters;
 172         }
 173
 174         // loop over all interactions in inputElement
 175         for (const auto& interactionIndex : inputElement.indices)
 176         {
 177             // each interaction has multiple coordinate indices
 178             // we must pick one of them to assign this interaction to one of the output index ranges
 179             // Todo: count indices outside the current split range in order to minimize the buffer size
 180             int representativeIndex =
 181                     *std::min_element(begin(interactionIndex), end(interactionIndex) - 1);
 182             int splitIndex = detail::computeChunkIndex(representativeIndex, totalRange, nSplits);
 183
 184             std::get<elementIndex>(workDivision[splitIndex]).indices.push_back(interactionIndex);
 185         }
 186     };
 187
 188     // split each interaction type in the input interaction tuple
 189     for_each_tuple(splitOneElement, interactions);
 190
 191     return workDivision;
 192 }
 193
 194 } // namespace nblib
 195
 196 #undef NBLIB_ALWAYS_INLINE
 197
 198 #endif // NBLIB_LISTEDFORCSES_HELPERS_HPP