src/gromacs/tables/cubicsplinetable.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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
  36 /*! \libinternal \file
  37  * \brief
  38  * Declares classes for cubic spline table
  39  *
  40  * \inlibraryapi
  41  * \author Erik Lindahl <erik.lindahl@gmail.com>
  42  * \ingroup module_tables
  43  */
  44 #ifndef GMX_TABLES_CUBICSPLINETABLE_H
  45 #define GMX_TABLES_CUBICSPLINETABLE_H
  46
  47 #include <initializer_list>
  48 #include <vector>
  49
  50 #include "gromacs/simd/simd.h"
  51 #include "gromacs/tables/tableinput.h"
  52 #include "gromacs/utility/alignedallocator.h"
  53 #include "gromacs/utility/classhelpers.h"
  54 #include "gromacs/utility/exceptions.h"
  55 #include "gromacs/utility/gmxassert.h"
  56 #include "gromacs/utility/real.h"
  57
  58 namespace gmx
  59 {
  60
  61 /*! \libinternal \brief Cubic spline interpolation table.
  62  *
  63  * This class interpolates a function specified either as an analytical
  64  * expression or from user-provided table data.
  65  *
  66  * At initialization, you provide the reference function of vectors
  67  * as a list of tuples that contain a brief name, the function, and
  68  * derivative for each function to tabulate. To create a table with
  69  * two functions this initializer list can for instance look like
  70  *
  71  *     { {"LJ6", lj6Func, lj6Der}, {"LJ12", lj12Func, lj12Der} }
  72  *
  73  * The names are only used so exceptions during initialization can
  74  * be traced to a specific table.
  75  *
  76  * When interpolating, there are methods to interpolate either 1, 2, or 3
  77  * functions in one go. By default these interpolation routines will
  78  * operate on tables with the same number of functions as specified in
  79  * the interpolation method (debug builds check that this is consistent with
  80  * the table). However, it is also possible to use optional template
  81  * parameters that specify the total number of functions in a table, and
  82  * what function index to interpolate. For instance, to interpolate the
  83  * derivative of the second function (i.e., index 1) in a
  84  * multi-function-table with three functions in total, you can write
  85  *
  86  *     table.evaluateDerivative<3,1>(x,&der);
  87  *
  88  * Here too, debug builds will check that the template parameters are
  89  * consistent with the table.
  90  *
  91  * This class interpolates a function specified either as an analytical
  92  * expression or from user-provided table data. The coefficients for each
  93  * table point are precalculated such that we simply evaluate
  94  *
  95  * \f{eqnarray*}{
  96  * V(x)  = Y + F \epsilon + G \epsilon^2 + H \epsilon^3
  97  * V'(x) = (F + 2 G \epsilon + 3 H \epsilon^2)/h
  98  * \f}
  99  *
 100  * Where h is the spacing and epsilon the fractional offset from table point.
 101  *
 102  * While it is possible to create tables only from function values
 103  * (i.e., no derivatives), it is recommended to provide derivatives for higher
 104  * accuracy and to avoid issues with numerical differentiation. Note that the
 105  * table input should be smooth, i.e. it should not contain noise e.g. from an
 106  * (iterative) Boltzmann inversion procedure - you have been warned.
 107  *
 108  * \note This class is responsible for fundamental interpolation of any function,
 109  *       which might or might not correspond to a potential. For this reason
 110  *       both input and output derivatives are proper function derivatives, and
 111  *       we do not swap any signs to get forces directly from the table.
 112  *
 113  * \note There will be a small additional accuracy loss from the internal
 114  *       operation where we calculate the epsilon offset from the nearest table
 115  *       point, since the integer part we subtract can get large in those cases.
 116  *       While this is technically possible to solve with extended precision
 117  *       arithmetics, that would introduce extra instructions in some highly
 118  *       performance-sensitive code parts. For typical GROMACS interaction
 119  *       functions the derivatives will decay faster than the potential, which
 120  *       means it will never play any role. For other functions it will only
 121  *       cause a small increase in the relative error for arguments where the
 122  *       magnitude of the function or derivative is very small.
 123  *       Since we typically sum several results in GROMACS, this should never
 124  *       show up in any real cases, and for this reason we choose not to do
 125  *       the extended precision arithmetics.
 126  *
 127  * \note These routines are not suitable for table ranges starting far away
 128  *       from zero, since we allocate memory and calculate indices starting from
 129  *       range zero for efficiency reasons.
 130  */
 131 class CubicSplineTable
 132 {
 133     private:
 134         /*! \brief Change that function value falls inside range when debugging
 135          *
 136          *  \tparam T   Lookup argument floating-point type, typically SimdReal or real.
 137          *  \param  r   Lookup argument to test
 138          *
 139          *  \throws Debug builds will throw gmx::RangeError for values that are
 140          *          larger than the upper limit of the range, or smaller than 0.
 141          *          We allow the table to be called with arguments between 0 and
 142          *          the lower limit of the range, since this might in theory occur
 143          *          once-in-a-blue-moon with some algorithms.
 144          */
 145         template <typename T>
 146         void
 147         rangeCheck(T gmx_unused r) const
 148         {
 149 #ifndef NDEBUG
 150             // Check that all values fall in range when debugging
 151             if (anyTrue( r < T(0.0) || T(range_.second) <= r ) )
 152             {
 153                 GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
 154             }
 155 #endif
 156         }
 157
 158     public:
 159
 160         /*! \brief Default tolerance for cubic spline tables
 161          *
 162          * This is 10*GMX_FLOAT_EPS in single precision, and
 163          * 1e-10 for double precision. It might not be worth setting
 164          * this tolerance lower than 1e-10 in double precision, both because
 165          * you will end up with very large tables, and because
 166          * functions like r^-12 become so large for small values of r the
 167          * table generation code will lead to some precision loss even
 168          * in double precision.
 169          */
 170         static const real defaultTolerance;
 171
 172         /*! \brief Initialize table data from function
 173          *
 174          * \param analyticalInputList Initializer list with one or more functions to tabulate,
 175          *                            specified as elements with a string description and
 176          *                            the function as well as derivative. The function will also
 177          *                            be called for values smaller than the lower limit of the
 178          *                            range, but we avoid calling it for 0.0 if that value
 179          *                            is not included in the range.
 180          *                            Constructor will throw gmx::APIError for negative values.
 181          *                            Due to the way the numerical derivative evaluation depends
 182          *                            on machine precision internally, this range must be
 183          *                            at least 0.001, or the constructor throws gmx::APIError.
 184          * \param range                Range over which the function will be tabulated.
 185          *                             Constructor will throw gmx::APIError for negative values,
 186          *                             or if the value/derivative vector does not cover the
 187          *                             range.
 188          * \param tolerance           Requested accuracy of the table. This will be used to
 189          *                            calculate the required internal spacing. If this cannot
 190          *                            be achieved (for instance because the table would require
 191          *                            too much memory) the constructor will throw gmx::ToleranceError.
 192          *
 193          * \note The functions are always defined in double precision to avoid
 194          *       losing accuracy when constructing tables.
 195          *
 196          * \note Since we fill the table for values below range.first, you can achieve
 197          *       a smaller table by using a smaller range where the tolerance has to be
 198          *       met, and accept that a few function calls below range.first do not
 199          *       quite reach the tolerance.
 200          *
 201          * \warning For efficiency reasons (since this code is used in some inner
 202          *       (kernels), we always allocate memory and calculate table indices
 203          *       for the complete interval [0,range.second], although the data will
 204          *       not be valid outside the definition range to avoid calling the
 205          *       function there. This means you should \a not use this class
 206          *       to tabulate functions for small ranges very far away from zero,
 207          *       since you would both waste a huge amount of memory and incur
 208          *       truncation errors when calculating the index.
 209          *
 210          * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
 211          *         and gmx::APIError for other incorrect input.
 212          */
 213         CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput>  analyticalInputList,
 214                          const std::pair<real, real>                       &range,
 215                          real                                               tolerance = defaultTolerance);
 216
 217         /*! \brief Initialize table data from tabulated values and derivatives
 218          *
 219          * \param numericalInputList  Initializer list with one or more functions to tabulate,
 220          *                            specified as a string description, vectors with function and
 221          *                            derivative values, and the input spacing. Data points are
 222          *                            separated by the spacing parameter, starting from 0.
 223          *                            Values below the lower limit of the range will be used to
 224          *                            attempt defining the table, but we avoid using index 0
 225          *                            unless 0.0 is included in the range. Some extra points beyond
 226          *                            range.second are required to re-interpolate values, so add
 227          *                            some margin. The constructor will throw gmx::APIError if the
 228          *                            input vectors are too short to cover the requested range
 229          *                            (and they must always be at least five points).
 230          * \param range               Range over which the function will be tabulated.
 231          *                            Constructor will throw gmx::APIError for negative values,
 232          *                            or if the value/derivative vector does not cover the
 233          *                            range.
 234          * \param tolerance           Requested accuracy of the table. This will be used to
 235          *                            calculate the required internal spacing and possibly
 236          *                            re-interpolate. The constructor will throw
 237          *                            gmx::ToleranceError if the input spacing is too coarse
 238          *                            to achieve this accuracy.
 239          *
 240          * \note The input data vectors are always double precision to avoid
 241          *       losing accuracy when constructing tables.
 242          *
 243          * \note Since we fill the table for values below range.first, you can achieve
 244          *       a smaller table by using a smaller range where the tolerance has to be
 245          *       met, and accept that a few function calls below range.first do not
 246          *       quite reach the tolerance.
 247          *
 248          * \warning For efficiency reasons (since this code is used in some inner
 249          *       (kernels), we always allocate memory and calculate table indices
 250          *       for the complete interval [0,range.second], although the data will
 251          *       not be valid outside the definition range to avoid calling the
 252          *       function there. This means you should \a not use this class
 253          *       to tabulate functions for small ranges very far away from zero,
 254          *       since you would both waste a huge amount of memory and incur
 255          *       truncation errors when calculating the index.
 256          */
 257         CubicSplineTable(std::initializer_list<NumericalSplineTableInput>  numericalInputList,
 258                          const std::pair<real, real>                      &range,
 259                          real                                              tolerance = defaultTolerance);
 260
 261
 262         /************************************************************
 263          *           Evaluation methods for single functions        *
 264          ************************************************************/
 265
 266         /*! \brief Evaluate both function and derivative, single table function
 267          *
 268          *  This is a templated method where the template can be either real or SimdReal.
 269          *
 270          *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 271          *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 272          *  \tparam     T               Type (SimdReal or real) of lookup and result
 273          *  \param      r               Points for which to evaluate function and derivative
 274          *  \param[out] functionValue   Function value
 275          *  \param[out] derivativeValue Function derivative
 276          *
 277          *  For debug builds we assert that the input values fall in the range
 278          *  specified when constructing the table.
 279          */
 280         template <int numFuncInTable = 1, int funcIndex = 0, typename T>
 281         void
 282         evaluateFunctionAndDerivative(T     r,
 283                                       T *   functionValue,
 284                                       T *   derivativeValue) const
 285         {
 286             rangeCheck(r);
 287             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 288             GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
 289
 290             T     rTable   = r * T(tableScale_);
 291             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 292             T     eps      = rTable - trunc(rTable);
 293             T     Y, F, G, H;
 294
 295             // Load Derivative, Delta, Function, and Zero values for each table point.
 296             // The 4 refers to these four values - not any SIMD width.
 297             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex, tabIndex, &Y, &F, &G, &H);
 298             *functionValue   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 299             *derivativeValue = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 300         }
 301
 302         /*! \brief Evaluate function value only, single table function
 303          *
 304          *  This is a templated method where the template can be either real or SimdReal.
 305          *
 306          *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 307          *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 308          *  \tparam     T               Type (SimdReal or real) of lookup and result
 309          *  \param      r               Points for which to evaluate function value
 310          *  \param[out] functionValue   Function value
 311          *
 312          *  For debug builds we assert that the input values fall in the range
 313          *  specified when constructing the table.
 314          */
 315         template <int numFuncInTable = 1, int funcIndex = 0, typename T>
 316         void
 317         evaluateFunction(T     r,
 318                          T *   functionValue) const
 319         {
 320             T     der gmx_unused;
 321
 322             evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
 323         }
 324
 325         /*! \brief Evaluate function derivative only, single table function
 326          *
 327          *  This is a templated method where the template can be either real or SimdReal.
 328          *
 329          *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 330          *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 331          *  \tparam     T               Type (SimdReal or real) of lookup and result
 332          *  \param      r               Points for which to evaluate function derivative
 333          *  \param[out] derivativeValue Function derivative
 334          *
 335          *  For debug builds we assert that the input values fall in the range
 336          *  specified when constructing the table.
 337          */
 338         template <int numFuncInTable = 1, int funcIndex = 0, typename T>
 339         void
 340         evaluateDerivative(T     r,
 341                            T *   derivativeValue) const
 342         {
 343             rangeCheck(r);
 344             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 345             GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
 346
 347             T     rTable   = r * T(tableScale_);
 348             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 349             T     eps      = rTable - trunc(rTable);
 350             T     Y, F, G, H;
 351
 352             // Load Derivative, Delta, Function, and Zero values for each table point.
 353             // The 4 refers to these four values - not any SIMD width.
 354             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex, tabIndex, &Y, &F, &G, &H);
 355             *derivativeValue = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 356         }
 357
 358         /************************************************************
 359          *             Evaluation methods for two functions         *
 360          ************************************************************/
 361
 362         /*! \brief Evaluate both function and derivative, two table functions
 363          *
 364          *  This is a templated method where the template can be either real or SimdReal.
 365          *
 366          *  \tparam     numFuncInTable  Number of separate functions in table, default is 2
 367          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 368          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 369          *  \tparam     T                Type (SimdReal or real) of lookup and result
 370          *  \param      r                Points for which to evaluate function and derivative
 371          *  \param[out] functionValue0   Interpolated value for first function
 372          *  \param[out] derivativeValue0 Interpolated derivative for first function
 373          *  \param[out] functionValue1   Interpolated value for second function
 374          *  \param[out] derivativeValue1 Interpolated derivative for second function
 375          *
 376          *  For debug builds we assert that the input values fall in the range
 377          *  specified when constructing the table.
 378          */
 379         template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 380         void
 381         evaluateFunctionAndDerivative(T     r,
 382                                       T *   functionValue0,
 383                                       T *   derivativeValue0,
 384                                       T *   functionValue1,
 385                                       T *   derivativeValue1) const
 386         {
 387             rangeCheck(r);
 388             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 389             GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
 390
 391             T     rTable   = r * T(tableScale_);
 392             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 393             T     eps      = rTable - trunc(rTable);
 394             T     Y, F, G, H;
 395
 396             // Load Derivative, Delta, Function, and Zero values for each table point.
 397             // The 4 refers to these four values - not any SIMD width.
 398             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex0, tabIndex, &Y, &F, &G, &H);
 399             *functionValue0   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 400             *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 401
 402             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex1, tabIndex, &Y, &F, &G, &H);
 403             *functionValue1   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 404             *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 405         }
 406
 407         /*! \brief Evaluate function value only, two table functions
 408          *
 409          *  This is a templated method where the template can be either real or SimdReal.
 410          *
 411          *  \tparam     numFuncInTable   Number of separate functions in table, default is 2
 412          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 413          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 414          *  \tparam     T                Type (SimdReal or real) of lookup and result
 415          *  \param      r                Points for which to evaluate function value
 416          *  \param[out] functionValue0   Interpolated value for first function
 417          *  \param[out] functionValue1   Interpolated value for second function
 418          *
 419          *  For debug builds we assert that the input values fall in the range
 420          *  specified when constructing the table.
 421          */
 422         template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 423         void
 424         evaluateFunction(T     r,
 425                          T *   functionValue0,
 426                          T *   functionValue1) const
 427         {
 428             T     der0 gmx_unused;
 429             T     der1 gmx_unused;
 430
 431             evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(r, functionValue0, &der0, functionValue1, &der1);
 432         }
 433
 434         /*! \brief Evaluate function derivative only, two table functions
 435          *
 436          *  This is a templated method where the template can be either real or SimdReal.
 437          *
 438          *  \tparam     numFuncInTable   Number of separate functions in table, default is 2
 439          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 440          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 441          *  \tparam     T                Type (SimdReal or real) of lookup and result
 442          *  \param      r                Points for which to evaluate function derivative
 443          *  \param[out] derivativeValue0 Interpolated derivative for first function
 444          *  \param[out] derivativeValue1 Interpolated derivative for second function
 445          *
 446          *  For debug builds we assert that the input values fall in the range
 447          *  specified when constructing the table.
 448          */
 449         template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 450         void
 451         evaluateDerivative(T     r,
 452                            T *   derivativeValue0,
 453                            T *   derivativeValue1) const
 454         {
 455             rangeCheck(r);
 456             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 457             GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
 458
 459             T     rTable   = r * T(tableScale_);
 460             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 461             T     eps      = rTable - trunc(rTable);
 462             T     Y, F, G, H;
 463
 464             // Load Derivative, Delta, Function, and Zero values for each table point.
 465             // The 4 refers to these four values - not any SIMD width.
 466             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 467             *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 468
 469             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 470             *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 471         }
 472
 473         /************************************************************
 474          *            Evaluation methods for three functions        *
 475          ************************************************************/
 476
 477
 478         /*! \brief Evaluate both function and derivative, three table functions
 479          *
 480          *  This is a templated method where the template can be either real or SimdReal.
 481          *
 482          *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 483          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 484          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 485          *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 486          *  \tparam     T                Type (SimdReal or real) of lookup and result
 487          *  \param      r                Points for which to evaluate function and derivative
 488          *  \param[out] functionValue0   Interpolated value for first function
 489          *  \param[out] derivativeValue0 Interpolated derivative for first function
 490          *  \param[out] functionValue1   Interpolated value for second function
 491          *  \param[out] derivativeValue1 Interpolated derivative for second function
 492          *  \param[out] functionValue2   Interpolated value for third function
 493          *  \param[out] derivativeValue2 Interpolated derivative for third function
 494          *
 495          *  For debug builds we assert that the input values fall in the range
 496          *  specified when constructing the table.
 497          */
 498         template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 499         void
 500         evaluateFunctionAndDerivative(T     r,
 501                                       T *   functionValue0,
 502                                       T *   derivativeValue0,
 503                                       T *   functionValue1,
 504                                       T *   derivativeValue1,
 505                                       T *   functionValue2,
 506                                       T *   derivativeValue2) const
 507         {
 508             rangeCheck(r);
 509             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 510             GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable, "Function index not in range of the number of tables");
 511
 512             T     rTable   = r * T(tableScale_);
 513             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 514             T     eps      = rTable - trunc(rTable);
 515             T     Y, F, G, H;
 516
 517             // Load Derivative, Delta, Function, and Zero values for each table point.
 518             // The 4 refers to these four values - not any SIMD width.
 519             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex0, tabIndex, &Y, &F, &G, &H);
 520             *functionValue0   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 521             *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 522
 523             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex1, tabIndex, &Y, &F, &G, &H);
 524             *functionValue1   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 525             *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 526
 527             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex2, tabIndex, &Y, &F, &G, &H);
 528             *functionValue2   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 529             *derivativeValue2 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 530         }
 531
 532         /*! \brief Evaluate function value only, three table functions
 533          *
 534          *  This is a templated method where the template can be either real or SimdReal.
 535          *
 536          *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 537          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 538          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 539          *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 540          *  \tparam     T                Type (SimdReal or real) of lookup and result
 541          *  \param      r                Points for which to evaluate function value
 542          *  \param[out] functionValue0   Interpolated value for first function
 543          *  \param[out] functionValue1   Interpolated value for second function
 544          *  \param[out] functionValue2   Interpolated value for third function
 545          *
 546          *  For debug builds we assert that the input values fall in the range
 547          *  specified when constructing the table.
 548          */
 549         template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 550         void
 551         evaluateFunction(T     r,
 552                          T *   functionValue0,
 553                          T *   functionValue1,
 554                          T *   functionValue2) const
 555         {
 556             T     der0 gmx_unused;
 557             T     der1 gmx_unused;
 558             T     der2 gmx_unused;
 559
 560             evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(r, functionValue0, &der0, functionValue1, &der1, functionValue2, &der2);
 561         }
 562
 563         /*! \brief Evaluate function derivative only, three table functions
 564          *
 565          *  This is a templated method where the template can be either real or SimdReal.
 566          *
 567          *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 568          *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 569          *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 570          *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 571          *  \tparam     T                Type (SimdReal or real) of lookup and result
 572          *  \param      r                Points for which to evaluate function derivative
 573          *  \param[out] derivativeValue0 Interpolated derivative for first function
 574          *  \param[out] derivativeValue1 Interpolated derivative for second function
 575          *  \param[out] derivativeValue2 Interpolated derivative for third function
 576          *
 577          *  For debug builds we assert that the input values fall in the range
 578          *  specified when constructing the table.
 579          */
 580         template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 581         void
 582         evaluateDerivative(T     r,
 583                            T *   derivativeValue0,
 584                            T *   derivativeValue1,
 585                            T *   derivativeValue2) const
 586         {
 587             rangeCheck(r);
 588             GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
 589             GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable, "Function index not in range of the number of tables");
 590
 591             T     rTable   = r * T(tableScale_);
 592             auto  tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 593             T     eps      = rTable - trunc(rTable);
 594             T     Y, F, G, H;
 595
 596             // Load Derivative, Delta, Function, and Zero values for each table point.
 597             // The 4 refers to these four values - not any SIMD width.
 598             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 599             *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 600
 601             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 602             *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 603
 604             gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex2, tabIndex, &Y, &F, &G, &H);
 605             *derivativeValue2 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
 606         }
 607
 608         /*! \brief Return the table spacing (distance between points)
 609          *
 610          *  You should never have to use this for normal code, but due to the
 611          *  way tables are constructed internally we need this in the unit tests
 612          *  to check relative tolerances over each interval.
 613          *
 614          *  \return table spacing.
 615          */
 616         real
 617         tableSpacing() const { return 1.0 / tableScale_; }
 618
 619     private:
 620
 621         std::size_t             numFuncInTable_; //!< Number of separate tabluated functions
 622         std::pair<real, real>   range_;          //!< Range for which table evaluation is allowed
 623         real                    tableScale_;     //!< Table scale (inverse of spacing between points)
 624
 625         /*! \brief Vector with combined table data to save calculations after lookup.
 626          *
 627          *  For table point i, this vector contains the four coefficients
 628          *  Y,F,G,H that we use to express the function value as
 629          *  V(x)  = Y + F e + G e^2 + H e^3, where e is the epsilon offset from
 630          *  the nearest table point.
 631          *
 632          *  To allow aligned SIMD loads we need to use an aligned allocator for
 633          *  this container.
 634          */
 635         std::vector<real, AlignedAllocator<real> >  yfghMultiTableData_;
 636
 637         // There should never be any reason to copy the table since it is read-only
 638         GMX_DISALLOW_COPY_AND_ASSIGN(CubicSplineTable);
 639 };
 640
 641
 642 }      // namespace gmx
 643
 644 #endif // GMX_TABLES_CUBICSPLINETABLE_H