2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2016,2017,2018,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.
36 /*! \libinternal \file
38 * Declares classes for cubic spline table
41 * \author Erik Lindahl <erik.lindahl@gmail.com>
42 * \ingroup module_tables
44 #ifndef GMX_TABLES_CUBICSPLINETABLE_H
45 #define GMX_TABLES_CUBICSPLINETABLE_H
47 #include <initializer_list>
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"
61 /*! \libinternal \brief Cubic spline interpolation table.
63 * This class interpolates a function specified either as an analytical
64 * expression or from user-provided table data.
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
71 * { {"LJ6", lj6Func, lj6Der}, {"LJ12", lj12Func, lj12Der} }
73 * The names are only used so exceptions during initialization can
74 * be traced to a specific table.
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
86 * table.evaluateDerivative<3,1>(x,&der);
88 * Here too, debug builds will check that the template parameters are
89 * consistent with the table.
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
96 * V(x) = Y + F \epsilon + G \epsilon^2 + H \epsilon^3
97 * V'(x) = (F + 2 G \epsilon + 3 H \epsilon^2)/h
100 * Where h is the spacing and epsilon the fractional offset from table point.
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.
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.
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.
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.
131 class CubicSplineTable
134 /*! \brief Change that function value falls inside range when debugging
136 * \tparam T Lookup argument floating-point type, typically SimdReal or real.
137 * \param r Lookup argument to test
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.
146 void rangeCheck(T gmx_unused r) const
149 // Check that all values fall in range when debugging
150 if (anyTrue(r < T(0.0) || T(range_.second) <= r))
152 GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
158 /*! \brief Default tolerance for cubic spline tables
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.
168 static const real defaultTolerance;
170 /*! \brief Initialize table data from function
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
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.
191 * \note The functions are always defined in double precision to avoid
192 * losing accuracy when constructing tables.
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.
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.
208 * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
209 * and gmx::APIError for other incorrect input.
211 CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
212 const std::pair<real, real>& range,
213 real tolerance = defaultTolerance);
215 /*! \brief Initialize table data from tabulated values and derivatives
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
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.
238 * \note The input data vectors are always double precision to avoid
239 * losing accuracy when constructing tables.
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.
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.
255 CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
256 const std::pair<real, real>& range,
257 real tolerance = defaultTolerance);
260 /************************************************************
261 * Evaluation methods for single functions *
262 ************************************************************/
264 /*! \brief Evaluate both function and derivative, single table function
266 * This is a templated method where the template can be either real or SimdReal.
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
275 * For debug builds we assert that the input values fall in the range
276 * specified when constructing the table.
278 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
279 void evaluateFunctionAndDerivative(T r, T* functionValue, T* derivativeValue) const
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");
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);
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>(
295 yfghMultiTableData_.data() + 4 * funcIndex, 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);
300 /*! \brief Evaluate function value only, 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 value
308 * \param[out] functionValue Function value
310 * For debug builds we assert that the input values fall in the range
311 * specified when constructing the table.
313 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
314 void evaluateFunction(T r, T* functionValue) const
318 evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
321 /*! \brief Evaluate function derivative only, single table function
323 * This is a templated method where the template can be either real or SimdReal.
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
331 * For debug builds we assert that the input values fall in the range
332 * specified when constructing the table.
334 template<int numFuncInTable = 1, int funcIndex = 0, typename T>
335 void evaluateDerivative(T r, T* derivativeValue) const
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");
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);
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>(
351 yfghMultiTableData_.data() + 4 * funcIndex, tabIndex, &Y, &F, &G, &H);
352 *derivativeValue = tableScale_ * fma(fma(T(3.0) * H, eps, T(2.0) * G), eps, F);
355 /************************************************************
356 * Evaluation methods for two functions *
357 ************************************************************/
359 /*! \brief Evaluate both function and derivative, two table functions
361 * This is a templated method where the template can be either real or SimdReal.
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
373 * For debug builds we assert that the input values fall in the range
374 * specified when constructing the table.
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
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");
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);
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);
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);
403 /*! \brief Evaluate function value only, two table functions
405 * This is a templated method where the template can be either real or SimdReal.
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
415 * For debug builds we assert that the input values fall in the range
416 * specified when constructing the table.
418 template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
419 void evaluateFunction(T r, T* functionValue0, T* functionValue1) const
424 evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(
425 r, functionValue0, &der0, functionValue1, &der1);
428 /*! \brief Evaluate function derivative only, two table functions
430 * This is a templated method where the template can be either real or SimdReal.
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
440 * For debug builds we assert that the input values fall in the range
441 * specified when constructing the table.
443 template<int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
444 void evaluateDerivative(T r, T* derivativeValue0, T* derivativeValue1) const
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");
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);
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);
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);
468 /************************************************************
469 * Evaluation methods for three functions *
470 ************************************************************/
473 /*! \brief Evaluate both function and derivative, three table functions
475 * This is a templated method where the template can be either real or SimdReal.
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
490 * For debug builds we assert that the input values fall in the range
491 * specified when constructing the table.
493 template<int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
494 void evaluateFunctionAndDerivative(T r,
500 T* derivativeValue2) const
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");
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);
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);
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);
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);
531 /*! \brief Evaluate function value only, three table functions
533 * This is a templated method where the template can be either real or SimdReal.
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
545 * For debug builds we assert that the input values fall in the range
546 * specified when constructing the table.
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
555 evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(
556 r, functionValue0, &der0, functionValue1, &der1, functionValue2, &der2);
559 /*! \brief Evaluate function derivative only, three table functions
561 * This is a templated method where the template can be either real or SimdReal.
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
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 evaluateDerivative(T r, T* derivativeValue0, T* derivativeValue1, T* derivativeValue2) const
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");
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);
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);
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);
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);
605 /*! \brief Return the table spacing (distance between points)
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.
611 * \return table spacing.
613 real tableSpacing() const { return 1.0 / tableScale_; }
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)
620 /*! \brief Vector with combined table data to save calculations after lookup.
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.
627 * To allow aligned SIMD loads we need to use an aligned allocator for
630 std::vector<real, AlignedAllocator<real>> yfghMultiTableData_;
632 // There should never be any reason to copy the table since it is read-only
633 GMX_DISALLOW_COPY_AND_ASSIGN(CubicSplineTable);
639 #endif // GMX_TABLES_CUBICSPLINETABLE_H