2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2019,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.
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.
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.
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.
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.
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.
37 * \defgroup module_tables Classes for table interpolation
38 * \ingroup group_utilitymodules
40 * \brief Table interpolation from analytical or numerical input
42 * This module provides quadratic spline interpolation tables used
43 * both for the nonbonded kernels and listed interactions.
45 * \author Erik Lindahl <erik.lindahl@scilifelab.se>
49 /*! \libinternal \file
51 * Declares classes for quadratic spline table
54 * \author Erik Lindahl <erik.lindahl@gmail.com>
55 * \ingroup module_tables
57 #ifndef GMX_TABLES_QUADRATICSPLINETABLE_H
58 #define GMX_TABLES_QUADRATICSPLINETABLE_H
61 #include <initializer_list>
64 #include "gromacs/simd/simd.h"
65 #include "gromacs/tables/tableinput.h"
66 #include "gromacs/utility/alignedallocator.h"
67 #include "gromacs/utility/classhelpers.h"
68 #include "gromacs/utility/exceptions.h"
69 #include "gromacs/utility/gmxassert.h"
70 #include "gromacs/utility/real.h"
75 /*! \libinternal \brief Quadratic spline interpolation table.
77 * This class interpolates a function specified either as an analytical
78 * expression or from user-provided table data.
80 * At initialization, you provide the reference function of vectors
81 * as a list of tuples that contain a brief name, the function, and
82 * derivative for each function to tabulate. To create a table with
83 * two functions this initializer list can for instance look like
85 * { {"LJ6", lj6Func, lj6Der}, {"LJ12", lj12Func, lj12Der} }
87 * The names are only used so exceptions during initialization can
88 * be traced to a specific table.
90 * When interpolating, there are methods to interpolate either 1, 2, or 3
91 * functions in one go. By default these interpolation routines will
92 * operate on tables with the same number of functions as specified in
93 * the interpolation method (debug builds check that this is consistent with
94 * the table). However, it is also possible to use optional template
95 * parameters that specify the total number of functions in a table, and
96 * what function index to interpolate. For instance, to interpolate the
97 * derivative of the second function (i.e., index 1) in a
98 * multi-function-table with three functions in total, you can write
100 * table.evaluateDerivative<3,1>(x,&der);
102 * Here too, debug builds will check that the template parameters are
103 * consistent with the table.
105 * The table data is internally adjusted to guarantee that the interpolated
106 * derivative is the true derivative of the interpolated potential, which is
107 * important to avoid systematic errors for the common case when the derivative
108 * is concave/convex in the entire interval.
109 * We do this by expressing the difference in the function value
110 * at a small offset h relative to a reference value in position 0 with a forward
111 * Taylor series expanded around 0, and then doing the opposite of expressing
112 * difference in the function at position 0 relative to a reference value in
113 * position h when using a backward Taylor expansion:
116 * \Delta V & = & hV'(0) + \frac{1}{2} h^2 V''(0) + \frac{1}{6} h^3 V'''(0) + O(h^4) \\
117 * \Delta V & = & hV'(h) - \frac{1}{2} h^2 V''(h) + \frac{1}{6} h^3 V'''(h) + O(h^4)
120 * Summing the equations leads to
123 * 2 \Delta V = h(V'(0) + V'(h)) + \frac{1}{2} h^2 (V''(0)-V''(h)) +
124 * \frac{1}{6}h^3(V'''(0)+V'''(h)) + O(h^4) \f]
126 * To make the second term symmetric too, we can replace it with the average of
127 * the Taylor expansion at 0 and h (i.e., using the third derivative). This gives
130 * 2 \Delta V = h(V'(0) + V'(h)) - \frac{1}{12} h^3 (V'''(0)+V'''(h)) + O(h^4)
133 * Thus, if we replace the derivative in the internal quadratic table data with
136 * V' - \frac{1}{12}h^2 V'''
139 * we will cancel the h^3 term in the error. This will make the integral of the
140 * forces match the potential much better (The h^4 term actually disappears, so
141 * when summing over 1/h points the remaining error will be O(h^4).
143 * While it is possible to create tables only from function values
144 * (i.e., no derivatives), it is recommended to provide derivatives for higher
145 * accuracy and to avoid issues with numerical differentiation. Note that the
146 * table input should be smooth, i.e. it should not contain noise e.g. from an
147 * (iterative) Boltzmann inversion procedure - you have been warned.
149 * \note This class is responsible for fundamental interpolation of any function,
150 * which might or might not correspond to a potential. For this reason
151 * both input and output derivatives are proper function derivatives, and
152 * we do not swap any signs to get forces directly from the table.
154 * \note There will be a small additional accuracy loss from the internal
155 * operation where we calculate the epsilon offset from the nearest table
156 * point, since the integer part we subtract can get large in those cases.
157 * The absolute such error both in the function and derivative value will
158 * be roughly f''*x*GMX_REAL_EPS, where x is the argument and f'' the
160 * While this is technically possible to solve with extended precision
161 * arithmetics, that would introduce extra instructions in some highly
162 * performance-sensitive code parts. For typical GROMACS interaction
163 * functions the derivatives will decay faster than the potential, which
164 * means it will never play any role. For other functions it will only
165 * cause a small increase in the relative error for arguments where the
166 * magnitude of the function or derivative is very small.
167 * Since we typically sum several results in GROMACS, this should never
168 * show up in any real cases, and for this reason we choose not to do
169 * the extended precision arithmetics.
171 * \note These routines are not suitable for table ranges starting far away
172 * from zero, since we allocate memory and calculate indices starting from
173 * range zero for efficiency reasons.
175 class QuadraticSplineTable
178 /*! \brief Change that function value falls inside range when debugging
180 * \tparam T Lookup argument floating-point type, typically SimdReal or real.
181 * \param r Lookup argument to test
183 * \throws Debug builds will throw gmx::RangeError for values that are
184 * larger than the upper limit of the range, or smaller than 0.
185 * We allow the table to be called with arguments between 0 and
186 * the lower limit of the range, since this might in theory occur
187 * once-in-a-blue-moon with some algorithms.
190 void rangeCheck(T gmx_unused r) const
193 // Check that all values fall in range when debugging
194 if (anyTrue(r < T(0.0) || T(range_.second) <= r))
196 GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
202 /*! \brief Default tolerance for tables is 10*GMX_FLOAT_EPS
204 * \note Even for double precision builds we set the tolerance to
205 * one order of magnitude above the single precision epsilon.
207 static const real defaultTolerance;
209 /*! \brief Initialize table data from function
211 * \param analyticalInputList Initializer list with one or more functions to tabulate,
212 * specified as pairs containing analytical
213 * functions and their derivatives. The function will also
214 * be called for values smaller than the lower limit of the
215 * range, but we avoid calling it for 0.0 if that value
216 * is not included in the range.
217 * \param range Range over which the function will be tabulated.
218 * Constructor will throw gmx::APIError for negative values.
219 * Due to the way the numerical derivative evaluation depends
220 * on machine precision internally, this range must be
221 * at least 0.001, or the constructor throws gmx::APIError.
222 * \param tolerance Requested accuracy of the table. This will be used to
223 * calculate the required internal spacing. If this cannot
224 * be achieved (for instance because the table would require
225 * too much memory) the constructor will throw gmx::ToleranceError.
227 * \note The functions are always defined in double precision to avoid
228 * losing accuracy when constructing tables.
230 * \note Since we fill the table for values below range.first, you can achieve
231 * a smaller table by using a smaller range where the tolerance has to be
232 * met, and accept that a few function calls below range.first do not
233 * quite reach the tolerance.
235 * \warning For efficiency reasons (since this code is used in some inner
236 * (kernels), we always allocate memory and calculate table indices
237 * for the complete interval [0,range.second], although the data will
238 * not be valid outside the definition range to avoid calling the
239 * function there. This means you should \a not use this class
240 * to tabulate functions for small ranges very far away from zero,
241 * since you would both waste a huge amount of memory and incur
242 * truncation errors when calculating the index.
244 * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
245 * and gmx::APIError for other incorrect input.
247 QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
248 const std::pair<real, real>& range,
249 real tolerance = defaultTolerance);
251 /*! \brief Initialize table data from tabulated values and derivatives
253 * \param numericalInputList Initializer list with one or more functions to tabulate,
254 * specified as pairs containing containing vectors for the
255 * function values and their derivatives. Data points are
256 * separated by the spacing parameter, starting from 0.
257 * Values below the lower limit of the range will be used to
258 * attempt defining the table, but we avoid using index 0
259 * unless 0.0 is included in the range. Some extra points beyond
260 * range.second are required to re-interpolate values, so add
261 * some margin. The constructor will throw gmx::APIError if the
262 * input vectors are too short to cover the requested range
263 * (and they must always be at least five points).
264 * \param range Range over which the function will be tabulated.
265 * Constructor will throw gmx::APIError for negative values,
266 * or if the value/derivative vector does not cover the
268 * \param tolerance Requested accuracy of the table in the range. This will be
269 * used to calculate the required internal spacing and possibly
270 * re-interpolate. The constructor will throw
271 * gmx::ToleranceError if the input spacing is too coarse
272 * to achieve this accuracy.
274 * \note The input data vectors are always double precision to avoid
275 * losing accuracy when constructing tables.
277 * \note Since we fill the table for values below range.first, you can achieve
278 * a smaller table by using a smaller range where the tolerance has to be
279 * met, and accept that a few function calls below range.first do not
280 * quite reach the tolerance.
282 * \warning For efficiency reasons (since this code is used in some inner
283 * (kernels), we always allocate memory and calculate table indices
284 * for the complete interval [0,range.second], although the data will
285 * not be valid outside the definition range to avoid calling the
286 * function there. This means you should \a not use this class
287 * to tabulate functions for small ranges very far away from zero,
288 * since you would both waste a huge amount of memory and incur
289 * truncation errors when calculating the index.
291 QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
292 const std::pair<real, real>& range,
293 real tolerance = defaultTolerance);
296 /************************************************************
297 * Evaluation methods for single functions *
298 ************************************************************/
300 /*! \brief Evaluate both function and derivative, single table function
302 * This is a templated method where the template can be either real or SimdReal.
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 and derivative
308 * \param[out] functionValue Function value
309 * \param[out] derivativeValue Function derivative
311 * For debug builds we assert that the input values fall in the range
312 * specified when constructing the table.
314 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
315 void evaluateFunctionAndDerivative(T r, T* functionValue, T* derivativeValue) const
318 GMX_ASSERT(numFuncInTable == numFuncInTable_,
319 "Evaluation method not matching number of functions in table");
320 GMX_ASSERT(funcIndex < numFuncInTable,
321 "Function index not in range of the number of tables");
323 T rTable = r * T(tableScale_);
324 auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
325 T eps = rTable - trunc(rTable);
331 // Load Derivative, Delta, Function, and Zero values for each table point.
332 // The 4 refers to these four values - not any SIMD width.
333 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
334 ddfzMultiTableData_.data() + 4 * funcIndex, tabIndex, &t0, &t1, &t2, &t3);
337 *functionValue = fma(eps * T(halfSpacing_), t0 + t1, t2);
338 *derivativeValue = t1;
341 /*! \brief Evaluate function value only, single table function
343 * This is a templated method where the template can be either real or SimdReal.
345 * \tparam numFuncInTable Number of separate functions in table, default is 1
346 * \tparam funcIndex Index of function to evaluate in table, default is 0
347 * \tparam T Type (SimdReal or real) of lookup and result
348 * \param r Points for which to evaluate function value
349 * \param[out] functionValue Function value
351 * For debug builds we assert that the input values fall in the range
352 * specified when constructing the table.
354 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
355 void evaluateFunction(T r, T* functionValue) const
359 evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
362 /*! \brief Evaluate function derivative only, single table function
364 * This is a templated method where the template can be either real or SimdReal.
366 * \tparam numFuncInTable Number of separate functions in table, default is 1
367 * \tparam funcIndex Index of function to evaluate in table, default is 0
368 * \tparam T Type (SimdReal or real) of lookup and result
369 * \param r Points for which to evaluate function derivative
370 * \param[out] derivativeValue Function derivative
372 * For debug builds we assert that the input values fall in the range
373 * specified when constructing the table.
375 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
376 void evaluateDerivative(T r, T* derivativeValue) const
379 GMX_ASSERT(numFuncInTable == numFuncInTable_,
380 "Evaluation method not matching number of functions in table");
381 GMX_ASSERT(funcIndex < numFuncInTable,
382 "Function index not in range of the number of tables");
384 T rTable = r * T(tableScale_);
385 auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
386 T eps = rTable - trunc(rTable);
391 if (numFuncInTable == 1)
393 gatherLoadUBySimdIntTranspose<numFuncInTable>(
394 derivativeMultiTableData_.data() + funcIndex, tabIndex, &t0, &t1); // works for scalar T too
398 // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
399 // only loads a single value from memory to implement it better (will be written)
400 gatherLoadUBySimdIntTranspose<numFuncInTable>(
401 derivativeMultiTableData_.data() + funcIndex, tabIndex, &t0, &t2); // works for scalar T too
402 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex,
405 &t2); // works for scalar T too
408 // (1-eps)*t0 + eps*t1
409 *derivativeValue = fma(t1 - t0, eps, t0);
412 /************************************************************
413 * Evaluation methods for two functions *
414 ************************************************************/
416 /*! \brief Evaluate both function and derivative, two table functions
418 * This is a templated method where the template can be either real or SimdReal.
420 * \tparam numFuncInTable Number of separate functions in table, default is 2
421 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
422 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
423 * \tparam T Type (SimdReal or real) of lookup and result
424 * \param r Points for which to evaluate function and derivative
425 * \param[out] functionValue1 Interpolated value for first function
426 * \param[out] derivativeValue1 Interpolated derivative for first function
427 * \param[out] functionValue2 Interpolated value for second function
428 * \param[out] derivativeValue2 Interpolated derivative for second function
430 * For debug builds we assert that the input values fall in the range
431 * specified when constructing the table.
433 template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
434 void evaluateFunctionAndDerivative(T r, T* functionValue1, T* derivativeValue1, T* functionValue2, T* derivativeValue2) const
437 GMX_ASSERT(numFuncInTable == numFuncInTable_,
438 "Evaluation method not matching number of functions in table");
439 GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable,
440 "Function index not in range of the number of tables");
442 T rTable = r * T(tableScale_);
443 auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
444 T eps = rTable - trunc(rTable);
450 // Load Derivative, Delta, Function, and Zero values for each table point.
451 // The 4 refers to these four values - not any SIMD width.
452 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
453 ddfzMultiTableData_.data() + 4 * funcIndex0, tabIndex, &t0, &t1, &t2, &t3);
455 *functionValue1 = fma(eps * T(halfSpacing_), t0 + t1, t2);
456 *derivativeValue1 = t1;
458 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
459 ddfzMultiTableData_.data() + 4 * funcIndex1, tabIndex, &t0, &t1, &t2, &t3);
461 *functionValue2 = fma(eps * T(halfSpacing_), t0 + t1, t2);
462 *derivativeValue2 = t1;
465 /*! \brief Evaluate function value only, two table functions
467 * This is a templated method where the template can be either real or SimdReal.
469 * \tparam numFuncInTable Number of separate functions in table, default is 2
470 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
471 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
472 * \tparam T Type (SimdReal or real) of lookup and result
473 * \param r Points for which to evaluate function value
474 * \param[out] functionValue1 Interpolated value for first function
475 * \param[out] functionValue2 Interpolated value for second function
477 * For debug builds we assert that the input values fall in the range
478 * specified when constructing the table.
480 template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
481 void evaluateFunction(T r, T* functionValue1, T* functionValue2) const
486 evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(
487 r, functionValue1, &der1, functionValue2, &der2);
490 /*! \brief Evaluate function derivative only, two table functions
492 * This is a templated method where the template can be either real or SimdReal.
494 * \tparam numFuncInTable Number of separate functions in table, default is 2
495 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
496 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
497 * \tparam T Type (SimdReal or real) of lookup and result
498 * \param r Points for which to evaluate function derivative
499 * \param[out] derivativeValue1 Interpolated derivative for first function
500 * \param[out] derivativeValue2 Interpolated derivative for second function
502 * For debug builds we assert that the input values fall in the range
503 * specified when constructing the table.
505 template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
506 void evaluateDerivative(T r, T* derivativeValue1, T* derivativeValue2) const
509 GMX_ASSERT(numFuncInTable == numFuncInTable_,
510 "Evaluation method not matching number of functions in table");
511 GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable,
512 "Function index not in range of the number of tables");
514 T rTable = r * T(tableScale_);
515 auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
516 T eps = rTable - trunc(rTable);
518 if (numFuncInTable == 2 && funcIndex0 == 0 && funcIndex1 == 1)
520 T t0A, t0B, t1A, t1B;
521 gatherLoadUBySimdIntTranspose<numFuncInTable>(
522 derivativeMultiTableData_.data(), tabIndex, &t0A, &t0B); // works for scalar T too
523 gatherLoadUBySimdIntTranspose<numFuncInTable>(
524 derivativeMultiTableData_.data() + 2, tabIndex, &t1A, &t1B); // works for scalar T too
525 *derivativeValue1 = fma(t1A - t0A, eps, t0A);
526 *derivativeValue2 = fma(t1B - t0B, eps, t0B);
531 // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
532 // only loads a single value from memory to implement it better (will be written)
533 gatherLoadUBySimdIntTranspose<numFuncInTable>(
534 derivativeMultiTableData_.data() + funcIndex0, tabIndex, &t0, &t2); // works for scalar T too
535 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0,
538 &t2); // works for scalar T too
539 *derivativeValue1 = fma(t1 - t0, eps, t0);
541 gatherLoadUBySimdIntTranspose<numFuncInTable>(
542 derivativeMultiTableData_.data() + funcIndex1, tabIndex, &t0, &t2); // works for scalar T too
543 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1,
546 &t2); // works for scalar T too
547 *derivativeValue2 = fma(t1 - t0, eps, t0);
551 /************************************************************
552 * Evaluation methods for three functions *
553 ************************************************************/
556 /*! \brief Evaluate both function and derivative, three table functions
558 * This is a templated method where the template can be either real or SimdReal.
560 * \tparam numFuncInTable Number of separate functions in table, default is 3
561 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
562 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
563 * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
564 * \tparam T Type (SimdReal or real) of lookup and result
565 * \param r Points for which to evaluate function and derivative
566 * \param[out] functionValue1 Interpolated value for first function
567 * \param[out] derivativeValue1 Interpolated derivative for first function
568 * \param[out] functionValue2 Interpolated value for second function
569 * \param[out] derivativeValue2 Interpolated derivative for second function
570 * \param[out] functionValue3 Interpolated value for third function
571 * \param[out] derivativeValue3 Interpolated derivative for third function
573 * For debug builds we assert that the input values fall in the range
574 * specified when constructing the table.
576 template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
577 void evaluateFunctionAndDerivative(T r,
583 T* derivativeValue3) const
586 GMX_ASSERT(numFuncInTable == numFuncInTable,
587 "Evaluation method not matching number of functions in table");
588 GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
589 "Function index not in range of the number of tables");
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);
599 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
600 ddfzMultiTableData_.data() + 4 * funcIndex0, tabIndex, &t0, &t1, &t2, &t3);
602 *functionValue1 = fma(eps * T(halfSpacing_), t0 + t1, t2);
603 *derivativeValue1 = t1;
605 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
606 ddfzMultiTableData_.data() + 4 * funcIndex1, tabIndex, &t0, &t1, &t2, &t3);
608 *functionValue2 = fma(eps * T(halfSpacing_), t0 + t1, t2);
609 *derivativeValue2 = t1;
611 gatherLoadBySimdIntTranspose<4 * numFuncInTable>(
612 ddfzMultiTableData_.data() + 4 * funcIndex2, tabIndex, &t0, &t1, &t2, &t3);
614 *functionValue3 = fma(eps * T(halfSpacing_), t0 + t1, t2);
615 *derivativeValue3 = t1;
618 /*! \brief Evaluate function value only, three table functions
620 * This is a templated method where the template can be either real or SimdReal.
622 * \tparam numFuncInTable Number of separate functions in table, default is 3
623 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
624 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
625 * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
626 * \tparam T Type (SimdReal or real) of lookup and result
627 * \param r Points for which to evaluate function value
628 * \param[out] functionValue1 Interpolated value for first function
629 * \param[out] functionValue2 Interpolated value for second function
630 * \param[out] functionValue3 Interpolated value for third function
632 * For debug builds we assert that the input values fall in the range
633 * specified when constructing the table.
635 template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
636 void evaluateFunction(T r, T* functionValue1, T* functionValue2, T* functionValue3) const
642 evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(
643 r, functionValue1, &der1, functionValue2, &der2, functionValue3, &der3);
646 /*! \brief Evaluate function derivative only, three table functions
648 * This is a templated method where the template can be either real or SimdReal.
650 * \tparam numFuncInTable Number of separate functions in table, default is 3
651 * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
652 * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
653 * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
654 * \tparam T Type (SimdReal or real) of lookup and result
655 * \param r Points for which to evaluate function derivative
656 * \param[out] derivativeValue1 Interpolated derivative for first function
657 * \param[out] derivativeValue2 Interpolated derivative for second function
658 * \param[out] derivativeValue3 Interpolated derivative for third function
660 * For debug builds we assert that the input values fall in the range
661 * specified when constructing the table.
663 template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
664 void evaluateDerivative(T r, T* derivativeValue1, T* derivativeValue2, T* derivativeValue3) const
667 GMX_ASSERT(numFuncInTable == numFuncInTable_,
668 "Evaluation method not matching number of functions in table");
669 GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
670 "Function index not in range of the number of tables");
672 T rTable = r * T(tableScale_);
673 auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
674 T eps = rTable - trunc(rTable);
676 if (numFuncInTable == 3 && funcIndex0 == 0 && funcIndex1 == 1 && funcIndex2 == 2)
678 T t0A, t0B, t0C, t1A, t1B, t1C;
679 gatherLoadUBySimdIntTranspose<numFuncInTable>(
680 derivativeMultiTableData_.data(), tabIndex, &t0A, &t0B);
681 gatherLoadUBySimdIntTranspose<numFuncInTable>(
682 derivativeMultiTableData_.data() + 2, tabIndex, &t0C, &t1A);
683 gatherLoadUBySimdIntTranspose<numFuncInTable>(
684 derivativeMultiTableData_.data() + 4, tabIndex, &t1B, &t1C);
685 *derivativeValue1 = fma(t1A - t0A, eps, t0A);
686 *derivativeValue2 = fma(t1B - t0B, eps, t0B);
687 *derivativeValue3 = fma(t1C - t0C, eps, t0C);
692 // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
693 // only loads a single value from memory to implement it better (will be written)
694 gatherLoadUBySimdIntTranspose<numFuncInTable>(
695 derivativeMultiTableData_.data() + funcIndex0, tabIndex, &t0, &t2); // works for scalar T too
696 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0,
699 &t2); // works for scalar T too
700 *derivativeValue1 = fma(t1 - t0, eps, t0);
702 gatherLoadUBySimdIntTranspose<numFuncInTable>(
703 derivativeMultiTableData_.data() + funcIndex1, tabIndex, &t0, &t2); // works for scalar T too
704 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1,
707 &t2); // works for scalar T too
708 *derivativeValue2 = fma(t1 - t0, eps, t0);
710 gatherLoadUBySimdIntTranspose<numFuncInTable>(
711 derivativeMultiTableData_.data() + funcIndex2, tabIndex, &t0, &t2); // works for scalar T too
712 gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex2,
715 &t2); // works for scalar T too
716 *derivativeValue3 = fma(t1 - t0, eps, t0);
720 /*! \brief Return the table spacing (distance between points)
722 * You should never have to use this for normal code, but due to the
723 * way tables are constructed internally we need this in the unit tests
724 * to check relative tolerances over each interval.
726 * \return table spacing.
728 real tableSpacing() const { return 1.0 / tableScale_; }
731 std::size_t numFuncInTable_; //!< Number of separate tabluated functions
732 std::pair<real, real> range_; //!< Range for which table evaluation is allowed
733 real tableScale_; //!< Table scale (inverse of spacing between points)
734 real halfSpacing_; //!< 0.5*spacing (used for DDFZ table data)
736 //!< Derivative values only, with the third-derivative subtraction described in the class documentation.
737 std::vector<real> derivativeMultiTableData_;
739 /*! \brief Combined derivative, difference to next derivative, value, and zero.
741 * For table point i, this vector contains the four values:
743 * - (derivative[i+1]-derivative[i])
747 * For the derivative terms we have subtracted the third-derivative term described
748 * in the main class documentation.
750 * This is typically more efficient than the individual tables, in particular
751 * when using SIMD. The result should be identical outside the class, so this
752 * is merely an internal implementation optimization. However, to allow
753 * aligned SIMD loads we need to use an aligned allocator for this container.
754 * We occasionally abbreviate this data as DDFZ.
756 std::vector<real, AlignedAllocator<real>> ddfzMultiTableData_;
758 // There should never be any reason to copy the table since it is read-only
759 GMX_DISALLOW_COPY_AND_ASSIGN(QuadraticSplineTable);
765 #endif // GMX_TABLES_QUADRATICSPLINETABLE_H