src/gromacs/tables/cubicsplinetable.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2017,2018,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
  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 rangeCheck(T gmx_unused r) const
 147     {
 148 #ifndef NDEBUG
 149         // Check that all values fall in range when debugging
 150         if (anyTrue(r < T(0.0) || T(range_.second) <= r))
 151         {
 152             GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
 153         }
 154 #endif
 155     }
 156
 157 public:
 158     /*! \brief Default tolerance for cubic spline tables
 159      *
 160      * This is 10*GMX_FLOAT_EPS in single precision, and
 161      * 1e-10 for double precision. It might not be worth setting
 162      * this tolerance lower than 1e-10 in double precision, both because
 163      * you will end up with very large tables, and because
 164      * functions like r^-12 become so large for small values of r the
 165      * table generation code will lead to some precision loss even
 166      * in double precision.
 167      */
 168     static const real defaultTolerance;
 169
 170     /*! \brief Initialize table data from function
 171      *
 172      * \param analyticalInputList Initializer list with one or more functions to tabulate,
 173      *                            specified as elements with a string description and
 174      *                            the function as well as derivative. The function will also
 175      *                            be called for values smaller than the lower limit of the
 176      *                            range, but we avoid calling it for 0.0 if that value
 177      *                            is not included in the range.
 178      *                            Constructor will throw gmx::APIError for negative values.
 179      *                            Due to the way the numerical derivative evaluation depends
 180      *                            on machine precision internally, this range must be
 181      *                            at least 0.001, or the constructor throws gmx::APIError.
 182      * \param range                Range over which the function will be tabulated.
 183      *                             Constructor will throw gmx::APIError for negative values,
 184      *                             or if the value/derivative vector does not cover the
 185      *                             range.
 186      * \param tolerance           Requested accuracy of the table. This will be used to
 187      *                            calculate the required internal spacing. If this cannot
 188      *                            be achieved (for instance because the table would require
 189      *                            too much memory) the constructor will throw gmx::ToleranceError.
 190      *
 191      * \note The functions are always defined in double precision to avoid
 192      *       losing accuracy when constructing tables.
 193      *
 194      * \note Since we fill the table for values below range.first, you can achieve
 195      *       a smaller table by using a smaller range where the tolerance has to be
 196      *       met, and accept that a few function calls below range.first do not
 197      *       quite reach the tolerance.
 198      *
 199      * \warning For efficiency reasons (since this code is used in some inner
 200      *       (kernels), we always allocate memory and calculate table indices
 201      *       for the complete interval [0,range.second], although the data will
 202      *       not be valid outside the definition range to avoid calling the
 203      *       function there. This means you should \a not use this class
 204      *       to tabulate functions for small ranges very far away from zero,
 205      *       since you would both waste a huge amount of memory and incur
 206      *       truncation errors when calculating the index.
 207      *
 208      * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
 209      *         and gmx::APIError for other incorrect input.
 210      */
 211     CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
 212                      const std::pair<real, real>&                      range,
 213                      real tolerance = defaultTolerance);
 214
 215     /*! \brief Initialize table data from tabulated values and derivatives
 216      *
 217      * \param numericalInputList  Initializer list with one or more functions to tabulate,
 218      *                            specified as a string description, vectors with function and
 219      *                            derivative values, and the input spacing. Data points are
 220      *                            separated by the spacing parameter, starting from 0.
 221      *                            Values below the lower limit of the range will be used to
 222      *                            attempt defining the table, but we avoid using index 0
 223      *                            unless 0.0 is included in the range. Some extra points beyond
 224      *                            range.second are required to re-interpolate values, so add
 225      *                            some margin. The constructor will throw gmx::APIError if the
 226      *                            input vectors are too short to cover the requested range
 227      *                            (and they must always be at least five points).
 228      * \param range               Range over which the function will be tabulated.
 229      *                            Constructor will throw gmx::APIError for negative values,
 230      *                            or if the value/derivative vector does not cover the
 231      *                            range.
 232      * \param tolerance           Requested accuracy of the table. This will be used to
 233      *                            calculate the required internal spacing and possibly
 234      *                            re-interpolate. The constructor will throw
 235      *                            gmx::ToleranceError if the input spacing is too coarse
 236      *                            to achieve this accuracy.
 237      *
 238      * \note The input data vectors are always double precision to avoid
 239      *       losing accuracy when constructing tables.
 240      *
 241      * \note Since we fill the table for values below range.first, you can achieve
 242      *       a smaller table by using a smaller range where the tolerance has to be
 243      *       met, and accept that a few function calls below range.first do not
 244      *       quite reach the tolerance.
 245      *
 246      * \warning For efficiency reasons (since this code is used in some inner
 247      *       (kernels), we always allocate memory and calculate table indices
 248      *       for the complete interval [0,range.second], although the data will
 249      *       not be valid outside the definition range to avoid calling the
 250      *       function there. This means you should \a not use this class
 251      *       to tabulate functions for small ranges very far away from zero,
 252      *       since you would both waste a huge amount of memory and incur
 253      *       truncation errors when calculating the index.
 254      */
 255     CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
 256                      const std::pair<real, real>&                     range,
 257                      real                                             tolerance = defaultTolerance);
 258
 259
 260     /************************************************************
 261      *           Evaluation methods for single functions        *
 262      ************************************************************/
 263
 264     /*! \brief Evaluate both function and derivative, single table function
 265      *
 266      *  This is a templated method where the template can be either real or SimdReal.
 267      *
 268      *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 269      *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 270      *  \tparam     T               Type (SimdReal or real) of lookup and result
 271      *  \param      r               Points for which to evaluate function and derivative
 272      *  \param[out] functionValue   Function value
 273      *  \param[out] derivativeValue Function derivative
 274      *
 275      *  For debug builds we assert that the input values fall in the range
 276      *  specified when constructing the table.
 277      */
 278     template<int numFuncInTable = 1, int funcIndex = 0, typename T>
 279     void evaluateFunctionAndDerivative(T r, T* functionValue, T* derivativeValue) const
 280     {
 281         rangeCheck(r);
 282         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 283                    "Evaluation method not matching number of functions in table");
 284         GMX_ASSERT(funcIndex < numFuncInTable,
 285                    "Function index not in range of the number of tables");
 286
 287         T    rTable   = r * T(tableScale_);
 288         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 289         T    eps      = rTable - trunc(rTable);
 290         T    Y, F, G, H;
 291
 292         // Load Derivative, Delta, Function, and Zero values for each table point.
 293         // The 4 refers to these four values - not any SIMD width.
 294         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex,
 295                                                          tabIndex, &Y, &F, &G, &H);
 296         *functionValue   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 297         *derivativeValue = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 298     }
 299
 300     /*! \brief Evaluate function value only, single table function
 301      *
 302      *  This is a templated method where the template can be either real or SimdReal.
 303      *
 304      *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 305      *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 306      *  \tparam     T               Type (SimdReal or real) of lookup and result
 307      *  \param      r               Points for which to evaluate function value
 308      *  \param[out] functionValue   Function value
 309      *
 310      *  For debug builds we assert that the input values fall in the range
 311      *  specified when constructing the table.
 312      */
 313     template<int numFuncInTable = 1, int funcIndex = 0, typename T>
 314     void evaluateFunction(T r, T* functionValue) const
 315     {
 316         T der gmx_unused;
 317
 318         evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
 319     }
 320
 321     /*! \brief Evaluate function derivative only, single table function
 322      *
 323      *  This is a templated method where the template can be either real or SimdReal.
 324      *
 325      *  \tparam     numFuncInTable  Number of separate functions in table, default is 1
 326      *  \tparam     funcIndex       Index of function to evaluate in table, default is 0
 327      *  \tparam     T               Type (SimdReal or real) of lookup and result
 328      *  \param      r               Points for which to evaluate function derivative
 329      *  \param[out] derivativeValue Function derivative
 330      *
 331      *  For debug builds we assert that the input values fall in the range
 332      *  specified when constructing the table.
 333      */
 334     template<int numFuncInTable = 1, int funcIndex = 0, typename T>
 335     void evaluateDerivative(T r, T* derivativeValue) const
 336     {
 337         rangeCheck(r);
 338         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 339                    "Evaluation method not matching number of functions in table");
 340         GMX_ASSERT(funcIndex < numFuncInTable,
 341                    "Function index not in range of the number of tables");
 342
 343         T    rTable   = r * T(tableScale_);
 344         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 345         T    eps      = rTable - trunc(rTable);
 346         T    Y, F, G, H;
 347
 348         // Load Derivative, Delta, Function, and Zero values for each table point.
 349         // The 4 refers to these four values - not any SIMD width.
 350         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex,
 351                                                          tabIndex, &Y, &F, &G, &H);
 352         *derivativeValue = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 353     }
 354
 355     /************************************************************
 356      *             Evaluation methods for two functions         *
 357      ************************************************************/
 358
 359     /*! \brief Evaluate both function and derivative, two table functions
 360      *
 361      *  This is a templated method where the template can be either real or SimdReal.
 362      *
 363      *  \tparam     numFuncInTable  Number of separate functions in table, default is 2
 364      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 365      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 366      *  \tparam     T                Type (SimdReal or real) of lookup and result
 367      *  \param      r                Points for which to evaluate function and derivative
 368      *  \param[out] functionValue0   Interpolated value for first function
 369      *  \param[out] derivativeValue0 Interpolated derivative for first function
 370      *  \param[out] functionValue1   Interpolated value for second function
 371      *  \param[out] derivativeValue1 Interpolated derivative for second function
 372      *
 373      *  For debug builds we assert that the input values fall in the range
 374      *  specified when constructing the table.
 375      */
 376     template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 377     void evaluateFunctionAndDerivative(T r, T* functionValue0, T* derivativeValue0, T* functionValue1, T* derivativeValue1) const
 378     {
 379         rangeCheck(r);
 380         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 381                    "Evaluation method not matching number of functions in table");
 382         GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable,
 383                    "Function index not in range of the number of tables");
 384
 385         T    rTable   = r * T(tableScale_);
 386         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 387         T    eps      = rTable - trunc(rTable);
 388         T    Y, F, G, H;
 389
 390         // Load Derivative, Delta, Function, and Zero values for each table point.
 391         // The 4 refers to these four values - not any SIMD width.
 392         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 393                 yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 394         *functionValue0   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 395         *derivativeValue0 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 396
 397         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 398                 yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 399         *functionValue1   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 400         *derivativeValue1 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 401     }
 402
 403     /*! \brief Evaluate function value only, two table functions
 404      *
 405      *  This is a templated method where the template can be either real or SimdReal.
 406      *
 407      *  \tparam     numFuncInTable   Number of separate functions in table, default is 2
 408      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 409      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 410      *  \tparam     T                Type (SimdReal or real) of lookup and result
 411      *  \param      r                Points for which to evaluate function value
 412      *  \param[out] functionValue0   Interpolated value for first function
 413      *  \param[out] functionValue1   Interpolated value for second function
 414      *
 415      *  For debug builds we assert that the input values fall in the range
 416      *  specified when constructing the table.
 417      */
 418     template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 419     void evaluateFunction(T r, T* functionValue0, T* functionValue1) const
 420     {
 421         T der0 gmx_unused;
 422         T der1 gmx_unused;
 423
 424         evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(
 425                 r, functionValue0, &der0, functionValue1, &der1);
 426     }
 427
 428     /*! \brief Evaluate function derivative only, two table functions
 429      *
 430      *  This is a templated method where the template can be either real or SimdReal.
 431      *
 432      *  \tparam     numFuncInTable   Number of separate functions in table, default is 2
 433      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 434      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 435      *  \tparam     T                Type (SimdReal or real) of lookup and result
 436      *  \param      r                Points for which to evaluate function derivative
 437      *  \param[out] derivativeValue0 Interpolated derivative for first function
 438      *  \param[out] derivativeValue1 Interpolated derivative for second function
 439      *
 440      *  For debug builds we assert that the input values fall in the range
 441      *  specified when constructing the table.
 442      */
 443     template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
 444     void evaluateDerivative(T r, T* derivativeValue0, T* derivativeValue1) const
 445     {
 446         rangeCheck(r);
 447         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 448                    "Evaluation method not matching number of functions in table");
 449         GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable,
 450                    "Function index not in range of the number of tables");
 451
 452         T    rTable   = r * T(tableScale_);
 453         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 454         T    eps      = rTable - trunc(rTable);
 455         T    Y, F, G, H;
 456
 457         // Load Derivative, Delta, Function, and Zero values for each table point.
 458         // The 4 refers to these four values - not any SIMD width.
 459         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 460                 yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 461         *derivativeValue0 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 462
 463         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 464                 yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 465         *derivativeValue1 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 466     }
 467
 468     /************************************************************
 469      *            Evaluation methods for three functions        *
 470      ************************************************************/
 471
 472
 473     /*! \brief Evaluate both function and derivative, three table functions
 474      *
 475      *  This is a templated method where the template can be either real or SimdReal.
 476      *
 477      *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 478      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 479      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 480      *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 481      *  \tparam     T                Type (SimdReal or real) of lookup and result
 482      *  \param      r                Points for which to evaluate function and derivative
 483      *  \param[out] functionValue0   Interpolated value for first function
 484      *  \param[out] derivativeValue0 Interpolated derivative for first function
 485      *  \param[out] functionValue1   Interpolated value for second function
 486      *  \param[out] derivativeValue1 Interpolated derivative for second function
 487      *  \param[out] functionValue2   Interpolated value for third function
 488      *  \param[out] derivativeValue2 Interpolated derivative for third function
 489      *
 490      *  For debug builds we assert that the input values fall in the range
 491      *  specified when constructing the table.
 492      */
 493     template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 494     void evaluateFunctionAndDerivative(T  r,
 495                                        T* functionValue0,
 496                                        T* derivativeValue0,
 497                                        T* functionValue1,
 498                                        T* derivativeValue1,
 499                                        T* functionValue2,
 500                                        T* derivativeValue2) const
 501     {
 502         rangeCheck(r);
 503         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 504                    "Evaluation method not matching number of functions in table");
 505         GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
 506                    "Function index not in range of the number of tables");
 507
 508         T    rTable   = r * T(tableScale_);
 509         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 510         T    eps      = rTable - trunc(rTable);
 511         T    Y, F, G, H;
 512
 513         // Load Derivative, Delta, Function, and Zero values for each table point.
 514         // The 4 refers to these four values - not any SIMD width.
 515         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 516                 yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 517         *functionValue0   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 518         *derivativeValue0 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 519
 520         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 521                 yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 522         *functionValue1   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 523         *derivativeValue1 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 524
 525         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 526                 yfghMultiTableData_.data() + 4 * funcIndex2, tabIndex, &Y, &F, &G, &H);
 527         *functionValue2   = fma(fma(fma(H, eps, G), eps, F), eps, Y);
 528         *derivativeValue2 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 529     }
 530
 531     /*! \brief Evaluate function value only, three table functions
 532      *
 533      *  This is a templated method where the template can be either real or SimdReal.
 534      *
 535      *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 536      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 537      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 538      *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 539      *  \tparam     T                Type (SimdReal or real) of lookup and result
 540      *  \param      r                Points for which to evaluate function value
 541      *  \param[out] functionValue0   Interpolated value for first function
 542      *  \param[out] functionValue1   Interpolated value for second function
 543      *  \param[out] functionValue2   Interpolated value for third function
 544      *
 545      *  For debug builds we assert that the input values fall in the range
 546      *  specified when constructing the table.
 547      */
 548     template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 549     void evaluateFunction(T r, T* functionValue0, T* functionValue1, T* functionValue2) const
 550     {
 551         T der0 gmx_unused;
 552         T der1 gmx_unused;
 553         T der2 gmx_unused;
 554
 555         evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(
 556                 r, functionValue0, &der0, functionValue1, &der1, functionValue2, &der2);
 557     }
 558
 559     /*! \brief Evaluate function derivative only, three table functions
 560      *
 561      *  This is a templated method where the template can be either real or SimdReal.
 562      *
 563      *  \tparam     numFuncInTable   Number of separate functions in table, default is 3
 564      *  \tparam     funcIndex0       Index of 1st function to evaluate in table, default is 0
 565      *  \tparam     funcIndex1       Index of 2nd function to evaluate in table, default is 1
 566      *  \tparam     funcIndex2       Index of 3rd function to evaluate in table, default is 2
 567      *  \tparam     T                Type (SimdReal or real) of lookup and result
 568      *  \param      r                Points for which to evaluate function derivative
 569      *  \param[out] derivativeValue0 Interpolated derivative for first function
 570      *  \param[out] derivativeValue1 Interpolated derivative for second function
 571      *  \param[out] derivativeValue2 Interpolated derivative for third function
 572      *
 573      *  For debug builds we assert that the input values fall in the range
 574      *  specified when constructing the table.
 575      */
 576     template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
 577     void evaluateDerivative(T r, T* derivativeValue0, T* derivativeValue1, T* derivativeValue2) const
 578     {
 579         rangeCheck(r);
 580         GMX_ASSERT(numFuncInTable == numFuncInTable_,
 581                    "Evaluation method not matching number of functions in table");
 582         GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
 583                    "Function index not in range of the number of tables");
 584
 585         T    rTable   = r * T(tableScale_);
 586         auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
 587         T    eps      = rTable - trunc(rTable);
 588         T    Y, F, G, H;
 589
 590         // Load Derivative, Delta, Function, and Zero values for each table point.
 591         // The 4 refers to these four values - not any SIMD width.
 592         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 593                 yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
 594         *derivativeValue0 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 595
 596         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 597                 yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
 598         *derivativeValue1 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 599
 600         gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
 601                 yfghMultiTableData_.data() + 4 * funcIndex2, tabIndex, &Y, &F, &G, &H);
 602         *derivativeValue2 = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
 603     }
 604
 605     /*! \brief Return the table spacing (distance between points)
 606      *
 607      *  You should never have to use this for normal code, but due to the
 608      *  way tables are constructed internally we need this in the unit tests
 609      *  to check relative tolerances over each interval.
 610      *
 611      *  \return table spacing.
 612      */
 613     real tableSpacing() const { return 1.0 / tableScale_; }
 614
 615 private:
 616     std::size_t           numFuncInTable_; //!< Number of separate tabluated functions
 617     std::pair<real, real> range_;          //!< Range for which table evaluation is allowed
 618     real                  tableScale_;     //!< Table scale (inverse of spacing between points)
 619
 620     /*! \brief Vector with combined table data to save calculations after lookup.
 621      *
 622      *  For table point i, this vector contains the four coefficients
 623      *  Y,F,G,H that we use to express the function value as
 624      *  V(x)  = Y + F e + G e^2 + H e^3, where e is the epsilon offset from
 625      *  the nearest table point.
 626      *
 627      *  To allow aligned SIMD loads we need to use an aligned allocator for
 628      *  this container.
 629      */
 630     std::vector<real, AlignedAllocator<real>> yfghMultiTableData_;
 631
 632     // There should never be any reason to copy the table since it is read-only
 633     GMX_DISALLOW_COPY_AND_ASSIGN(CubicSplineTable);
 634 };
 635
 636
 637 } // namespace gmx
 638
 639 #endif // GMX_TABLES_CUBICSPLINETABLE_H