2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,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.
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 * Tests for simple math functions.eval
39 * \author Erik Lindahl <erik.lindahl@gmail.com>
40 * \ingroup module_tables
50 #include <gtest/gtest.h>
52 #include "gromacs/math/utilities.h"
53 #include "gromacs/options/basicoptions.h"
54 #include "gromacs/options/ioptionscontainer.h"
55 #include "gromacs/simd/simd.h"
56 #include "gromacs/tables/cubicsplinetable.h"
57 #include "gromacs/tables/quadraticsplinetable.h"
59 #include "testutils/testasserts.h"
60 #include "testutils/testoptions.h"
72 class SplineTableTestBase : public ::testing::Test
75 static int s_testPoints_; //!< Number of points to use. Public so we can set it as option
78 int SplineTableTestBase::s_testPoints_ = 100;
81 /*! \brief Command-line option to adjust the number of points used to test SIMD math functions. */
82 GMX_TEST_OPTIONS(SplineTableTestOptions, options)
84 options->addOption(::gmx::IntegerOption("npoints")
85 .store(&SplineTableTestBase::s_testPoints_)
86 .description("Number of points to test for spline table functions"));
91 /*! \brief Test fixture for table comparision with analytical/numerical functions */
93 class SplineTableTest : public SplineTableTestBase
96 SplineTableTest() : tolerance_(T::defaultTolerance) {}
98 /*! \brief Set a new tolerance to be used in table function comparison
100 * \param tol New tolerance to use
102 void setTolerance(real tol) { tolerance_ = tol; }
106 * Assertion predicate formatter for comparing table with function/derivative
108 template<int numFuncInTable = 1, int funcIndex = 0>
109 void testSplineTableAgainstFunctions(const std::string& desc,
110 const std::function<double(double)>& refFunc,
111 const std::function<double(double)>& refDer,
113 const std::pair<real, real>& testRange);
117 real tolerance_; //!< Tolerance to use
121 template<int numFuncInTable, int funcIndex>
122 void SplineTableTest<T>::testSplineTableAgainstFunctions(const std::string& desc,
123 const std::function<double(double)>& refFunc,
124 const std::function<double(double)>& refDer,
126 const std::pair<real, real>& testRange)
128 real dx = (testRange.second - testRange.first) / s_testPoints_;
130 FloatingPointTolerance funcTolerance(relativeToleranceAsFloatingPoint(0.0, tolerance_));
132 for (real x = testRange.first; x < testRange.second; x += dx) // NOLINT(clang-analyzer-security.FloatLoopCounter)
134 real h = std::sqrt(GMX_REAL_EPS);
135 real secondDerivative = (refDer(x + h) - refDer(x)) / h;
140 table.template evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(x, &testFuncValue,
143 // Check that we get the same values from function/derivative-only methods
144 real tmpFunc, tmpDer;
146 table.template evaluateFunction<numFuncInTable, funcIndex>(x, &tmpFunc);
148 table.template evaluateDerivative<numFuncInTable, funcIndex>(x, &tmpDer);
150 // Before we even start to think about errors related to the table interpolation
151 // accuracy, we want to test that the interpolations are consistent whether we
152 // call the routine that evaluates both the function and derivative or only one
154 // Note that for these tests the relevant tolerance is NOT the default one
155 // provided based on the requested accuracy of the table, but a tolerance related
156 // to the floating-point precision used. For now we only allow deviations up
157 // to 4 ulp (one for the FMA order, and then some margin).
158 FloatingPointTolerance consistencyTolerance(ulpTolerance(4));
160 FloatingPointDifference evaluateFuncDiff(tmpFunc, testFuncValue);
161 if (!consistencyTolerance.isWithin(evaluateFuncDiff))
163 ADD_FAILURE() << "Interpolation inconsistency for table " << desc << std::endl
164 << numFuncInTable << " function(s) in table, testing index " << funcIndex
166 << "First failure at x = " << x << std::endl
167 << "Function value when evaluating function & derivative: " << testFuncValue
169 << "Function value when evaluating only function: " << tmpFunc
174 FloatingPointDifference evaluateDerDiff(tmpDer, testDerValue);
175 if (!consistencyTolerance.isWithin(evaluateDerDiff))
177 ADD_FAILURE() << "Interpolation inconsistency for table " << desc << std::endl
178 << numFuncInTable << " function(s) in table, testing index " << funcIndex
180 << "First failure at x = " << x << std::endl
181 << "Derivative value when evaluating function & derivative: " << testDerValue
183 << "Derivative value when evaluating only derivative: " << tmpDer
188 // Next, we should examine that the table is exact enough relative
189 // to the requested accuracy in the interpolation.
191 // There are two sources of errors that we need to account for when checking the values,
192 // and we only fail the test if both of these tolerances are violated:
194 // 1) First, we have the normal relative error of the test vs. reference value. For this
195 // we use the normal testutils relative tolerance checking.
196 // However, there is an additional source of error: When we calculate the forces we
197 // use average higher derivatives over the interval to improve accuracy, but this
198 // also means we won't reproduce values at table points exactly. This is usually not
199 // an issue since the tolerances we have are much larger, but when the reference derivative
200 // value is exactly zero the relative error will be infinite. To account for this, we
201 // extract the spacing from the table and evaluate the reference derivative at a point
202 // this much larger too, and use the largest of the two values as the reference
203 // magnitude for the derivative when setting the relative tolerance.
204 // Note that according to the table function definitions, we should be allowed to evaluate
205 // it one table point beyond the range (this is done already for construction).
207 // 2) Second, due to the loss-of-accuracy when calculating the index through rtable
208 // there is an internal absolute tolerance that we can calculate.
209 // The source of this error is the subtraction eps=rtab-[rtab], which leaves an
210 // error proportional to eps_machine*rtab=eps_machine*x*tableScale.
211 // To lowest order, the term in the function and derivative values (respectively) that
212 // are proportional to eps will be the next-higher derivative multiplied by the spacing.
213 // This means the truncation error in the value is derivative*x*eps_machine, and in the
214 // derivative the error is 2nd_derivative*x*eps_machine.
216 real refFuncValue = refFunc(x);
217 real refDerValue = refDer(x);
218 real nextRefDerValue = refDer(x + table.tableSpacing());
220 real derMagnitude = std::max(std::abs(refDerValue), std::abs(nextRefDerValue));
222 // Since the reference magnitude will change over each interval we need to re-evaluate
223 // the derivative tolerance inside the loop.
224 FloatingPointTolerance derTolerance(relativeToleranceAsFloatingPoint(derMagnitude, tolerance_));
226 FloatingPointDifference funcDiff(refFuncValue, testFuncValue);
227 FloatingPointDifference derDiff(refDerValue, testDerValue);
229 real allowedAbsFuncErr = std::abs(refDerValue) * x * GMX_REAL_EPS;
230 real allowedAbsDerErr = std::abs(secondDerivative) * x * GMX_REAL_EPS;
232 if ((!funcTolerance.isWithin(funcDiff) && funcDiff.asAbsolute() > allowedAbsFuncErr)
233 || (!derTolerance.isWithin(derDiff) && derDiff.asAbsolute() > allowedAbsDerErr))
235 ADD_FAILURE() << "Failing comparison with function for table " << desc << std::endl
236 << numFuncInTable << " function(s) in table, testing index " << funcIndex
238 << "Test range is ( " << testRange.first << " , " << testRange.second
239 << " ) " << std::endl
240 << "Tolerance = " << tolerance_ << std::endl
241 << "First failure at x = " << x << std::endl
242 << "Reference function = " << refFuncValue << std::endl
243 << "Test table function = " << testFuncValue << std::endl
244 << "Allowed abs func err. = " << allowedAbsFuncErr << std::endl
245 << "Reference derivative = " << refDerValue << std::endl
246 << "Test table derivative = " << testDerValue << std::endl
247 << "Allowed abs der. err. = " << allowedAbsDerErr << std::endl
248 << "Actual abs der. err. = " << derDiff.asAbsolute() << std::endl;
255 /*! \brief Function similar to coulomb electrostatics
260 double coulombFunction(double r)
265 /*! \brief Derivative (not force) of coulomb electrostatics
270 double coulombDerivative(double r)
272 return -1.0 / (r * r);
275 /*! \brief Function similar to power-6 Lennard-Jones dispersion
280 double lj6Function(double r)
282 return std::pow(r, -6.0);
285 /*! \brief Derivative (not force) of the power-6 Lennard-Jones dispersion
290 double lj6Derivative(double r)
292 return -6.0 * std::pow(r, -7.0);
295 /*! \brief Function similar to power-12 Lennard-Jones repulsion
300 double lj12Function(double r)
302 return std::pow(r, -12.0);
305 /*! \brief Derivative (not force) of the power-12 Lennard-Jones repulsion
308 * \return -12.0*r^-13
310 double lj12Derivative(double r)
312 return -12.0 * std::pow(r, -13.0);
315 /*! \brief The sinc function, sin(r)/r
320 double sincFunction(double r)
322 return std::sin(r) / r;
325 /*! \brief Derivative of the sinc function
328 * \return derivative of sinc, (r*cos(r)-sin(r))/r^2
330 double sincDerivative(double r)
332 return (r * std::cos(r) - std::sin(r)) / (r * r);
335 /*! \brief Function for the direct-space PME correction to 1/r
338 * \return PME correction function, erf(r)/r
340 double pmeCorrFunction(double r)
344 return 2.0 / std::sqrt(M_PI);
348 return std::erf(r) / r;
352 /*! \brief Derivative of the direct-space PME correction to 1/r
355 * \return Derivative of the PME correction function.
357 double pmeCorrDerivative(double r)
365 return (2.0 * std::exp(-r * r) / std::sqrt(3.14159265358979323846) * r - erf(r)) / (r * r);
369 /*! \brief Typed-test list. We test QuadraticSplineTable and CubicSplineTable
371 typedef ::testing::Types<QuadraticSplineTable, CubicSplineTable> SplineTableTypes;
372 TYPED_TEST_CASE(SplineTableTest, SplineTableTypes);
375 TYPED_TEST(SplineTableTest, HandlesIncorrectInput)
378 EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { -1.0, 0.0 }),
379 gmx::InvalidInputError);
382 EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1.0, 1.00001 }),
383 gmx::InvalidInputError);
386 EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1.0, 2.0 }, 1e-20),
387 gmx::ToleranceError);
389 // Range is so close to 0.0 that table would require >1e6 points
390 EXPECT_THROW_GMX(TypeParam({ { "LJ12", lj12Function, lj12Derivative } }, { 1e-4, 2.0 }),
391 gmx::ToleranceError);
393 // mismatching function/derivative
394 EXPECT_THROW_GMX(TypeParam({ { "BadLJ12", lj12Derivative, lj12Function } }, { 1.0, 2.0 }),
395 gmx::InconsistentInputError);
400 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValues)
402 TypeParam table({ { "LJ12", lj12Function, lj12Derivative } }, { 0.2, 1.0 });
406 EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
409 EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
414 TYPED_TEST(SplineTableTest, Sinc)
416 // Sinc hits some sensitive parts of the table construction code which means
417 // we will not have full relative accuracy close to the zeros in the
418 // derivative. Since this is intentially a pathological function we reduce
419 // the interval slightly for now.
420 std::pair<real, real> range(0.1, 3.1);
422 TypeParam sincTable({ { "Sinc", sincFunction, sincDerivative } }, range);
424 TestFixture::testSplineTableAgainstFunctions("Sinc", sincFunction, sincDerivative, sincTable, range);
428 TYPED_TEST(SplineTableTest, LJ12)
430 std::pair<real, real> range(0.2, 2.0);
432 TypeParam lj12Table({ { "LJ12", lj12Function, lj12Derivative } }, range);
434 TestFixture::testSplineTableAgainstFunctions("LJ12", lj12Function, lj12Derivative, lj12Table, range);
438 TYPED_TEST(SplineTableTest, PmeCorrection)
440 std::pair<real, real> range(0.0, 4.0);
441 real tolerance = 1e-5;
443 TypeParam pmeCorrTable({ { "PMECorr", pmeCorrFunction, pmeCorrDerivative } }, range, tolerance);
445 TestFixture::setTolerance(tolerance);
446 TestFixture::testSplineTableAgainstFunctions("PMECorr", pmeCorrFunction, pmeCorrDerivative,
447 pmeCorrTable, range);
451 TYPED_TEST(SplineTableTest, HandlesIncorrectNumericalInput)
453 // Lengths do not match
454 std::vector<double> functionValues(10);
455 std::vector<double> derivativeValues(20);
457 TypeParam({ { "EmptyVectors", functionValues, derivativeValues, 0.001 } }, { 1.0, 2.0 }),
458 gmx::InconsistentInputError);
460 // Upper range is 2.0, spacing 0.1. This requires at least 21 points. Make sure we get an error for 20.
461 functionValues.resize(20);
462 derivativeValues.resize(20);
463 EXPECT_THROW_GMX(TypeParam({ { "EmptyVectors", functionValues, derivativeValues, 0.1 } }, { 1.0, 2.0 }),
464 gmx::InconsistentInputError);
466 // Create some test data
467 functionValues.clear();
468 derivativeValues.clear();
470 std::vector<double> badDerivativeValues;
471 double spacing = 1e-3;
473 for (std::size_t i = 0; i < 1001; i++)
475 double x = i * spacing;
476 double func = (x >= 0.1) ? lj12Function(x) : 0.0;
477 double der = (x >= 0.1) ? lj12Derivative(x) : 0.0;
479 functionValues.push_back(func);
480 derivativeValues.push_back(der);
481 badDerivativeValues.push_back(-der);
484 // Derivatives not consistent with function
485 EXPECT_THROW_GMX(TypeParam({ { "NumericalBadLJ12", functionValues, badDerivativeValues, spacing } },
487 gmx::InconsistentInputError);
489 // Spacing 1e-3 is not sufficient for r^-12 in range [0.1,1.0]
490 // Make sure we get a tolerance error
492 TypeParam({ { "NumericalLJ12", functionValues, derivativeValues, spacing } }, { 0.2, 1.0 }),
493 gmx::ToleranceError);
497 TYPED_TEST(SplineTableTest, NumericalInputPmeCorr)
499 std::pair<real, real> range(0.0, 4.0);
500 std::vector<double> functionValues;
501 std::vector<double> derivativeValues;
503 double inputSpacing = 1e-3;
504 real tolerance = 1e-5;
506 // We only need data up to the argument 4.0, but add 1% margin
507 for (std::size_t i = 0; i < range.second * 1.01 / inputSpacing; i++)
509 double x = i * inputSpacing;
511 functionValues.push_back(pmeCorrFunction(x));
512 derivativeValues.push_back(pmeCorrDerivative(x));
515 TypeParam pmeCorrTable({ { "NumericalPMECorr", functionValues, derivativeValues, inputSpacing } },
518 TestFixture::setTolerance(tolerance);
519 TestFixture::testSplineTableAgainstFunctions("NumericalPMECorr", pmeCorrFunction,
520 pmeCorrDerivative, pmeCorrTable, range);
523 TYPED_TEST(SplineTableTest, TwoFunctions)
525 std::pair<real, real> range(0.2, 2.0);
527 TypeParam table({ { "LJ6", lj6Function, lj6Derivative }, { "LJ12", lj12Function, lj12Derivative } },
530 // Test entire range for each function. This will use the method that interpolates a single function
531 TestFixture::template testSplineTableAgainstFunctions<2, 0>("LJ6", lj6Function, lj6Derivative,
533 TestFixture::template testSplineTableAgainstFunctions<2, 1>("LJ12", lj12Function,
534 lj12Derivative, table, range);
536 // Test the methods that evaluated both functions for one value
537 real x = 0.5 * (range.first + range.second);
538 real refFunc0 = lj6Function(x);
539 real refDer0 = lj6Derivative(x);
540 real refFunc1 = lj12Function(x);
541 real refDer1 = lj12Derivative(x);
543 real tstFunc0, tstDer0, tstFunc1, tstDer1;
544 real tmpFunc0, tmpFunc1, tmpDer0, tmpDer1;
546 // test that we reproduce the reference functions
547 table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
549 real funcErr0 = std::abs(tstFunc0 - refFunc0) / std::abs(refFunc0);
550 real funcErr1 = std::abs(tstFunc1 - refFunc1) / std::abs(refFunc1);
551 real derErr0 = std::abs(tstDer0 - refDer0) / std::abs(refDer0);
552 real derErr1 = std::abs(tstDer1 - refDer1) / std::abs(refDer1);
554 // Use asserts, since the following ones compare to these values.
555 ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
556 ASSERT_LT(derErr0, TypeParam::defaultTolerance);
557 ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
558 ASSERT_LT(derErr1, TypeParam::defaultTolerance);
560 // Test that function/derivative-only interpolation methods work
561 table.evaluateFunction(x, &tmpFunc0, &tmpFunc1);
562 table.evaluateDerivative(x, &tmpDer0, &tmpDer1);
563 EXPECT_EQ(tstFunc0, tmpFunc0);
564 EXPECT_EQ(tstFunc1, tmpFunc1);
565 EXPECT_EQ(tstDer0, tmpDer0);
567 // Test that scrambled order interpolation methods work
568 table.template evaluateFunctionAndDerivative<2, 1, 0>(x, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
569 EXPECT_EQ(tstFunc0, tmpFunc0);
570 EXPECT_EQ(tstFunc1, tmpFunc1);
571 EXPECT_EQ(tstDer0, tmpDer0);
572 EXPECT_EQ(tstDer1, tmpDer1);
574 // Test scrambled order for function/derivative-only methods
575 table.template evaluateFunction<2, 1, 0>(x, &tmpFunc1, &tmpFunc0);
576 table.template evaluateDerivative<2, 1, 0>(x, &tmpDer1, &tmpDer0);
577 EXPECT_EQ(tstFunc0, tmpFunc0);
578 EXPECT_EQ(tstFunc1, tmpFunc1);
579 EXPECT_EQ(tstDer0, tmpDer0);
580 EXPECT_EQ(tstDer1, tmpDer1);
583 TYPED_TEST(SplineTableTest, ThreeFunctions)
585 std::pair<real, real> range(0.2, 2.0);
587 TypeParam table({ { "Coulomb", coulombFunction, coulombDerivative },
588 { "LJ6", lj6Function, lj6Derivative },
589 { "LJ12", lj12Function, lj12Derivative } },
592 // Test entire range for each function
593 TestFixture::template testSplineTableAgainstFunctions<3, 0>("Coulomb", coulombFunction,
594 coulombDerivative, table, range);
595 TestFixture::template testSplineTableAgainstFunctions<3, 1>("LJ6", lj6Function, lj6Derivative,
597 TestFixture::template testSplineTableAgainstFunctions<3, 2>("LJ12", lj12Function,
598 lj12Derivative, table, range);
600 // Test the methods that evaluated both functions for one value
601 real x = 0.5 * (range.first + range.second);
602 real refFunc0 = coulombFunction(x);
603 real refDer0 = coulombDerivative(x);
604 real refFunc1 = lj6Function(x);
605 real refDer1 = lj6Derivative(x);
606 real refFunc2 = lj12Function(x);
607 real refDer2 = lj12Derivative(x);
609 real tstFunc0, tstDer0, tstFunc1, tstDer1, tstFunc2, tstDer2;
610 real tmpFunc0, tmpFunc1, tmpFunc2, tmpDer0, tmpDer1, tmpDer2;
612 // test that we reproduce the reference functions
613 table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1, &tstFunc2, &tstDer2);
615 real funcErr0 = std::abs(tstFunc0 - refFunc0) / std::abs(refFunc0);
616 real derErr0 = std::abs(tstDer0 - refDer0) / std::abs(refDer0);
617 real funcErr1 = std::abs(tstFunc1 - refFunc1) / std::abs(refFunc1);
618 real derErr1 = std::abs(tstDer1 - refDer1) / std::abs(refDer1);
619 real funcErr2 = std::abs(tstFunc2 - refFunc2) / std::abs(refFunc2);
620 real derErr2 = std::abs(tstDer2 - refDer2) / std::abs(refDer2);
622 // Use asserts, since the following ones compare to these values.
623 ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
624 ASSERT_LT(derErr0, TypeParam::defaultTolerance);
625 ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
626 ASSERT_LT(derErr1, TypeParam::defaultTolerance);
627 ASSERT_LT(funcErr2, TypeParam::defaultTolerance);
628 ASSERT_LT(derErr2, TypeParam::defaultTolerance);
630 // Test that function/derivative-only interpolation methods work
631 table.evaluateFunction(x, &tmpFunc0, &tmpFunc1, &tmpFunc2);
632 table.evaluateDerivative(x, &tmpDer0, &tmpDer1, &tmpDer2);
633 EXPECT_EQ(tstFunc0, tmpFunc0);
634 EXPECT_EQ(tstFunc1, tmpFunc1);
635 EXPECT_EQ(tstFunc2, tmpFunc2);
636 EXPECT_EQ(tstDer0, tmpDer0);
637 EXPECT_EQ(tstDer1, tmpDer1);
638 EXPECT_EQ(tstDer2, tmpDer2);
640 // Test two functions out of three
641 table.template evaluateFunctionAndDerivative<3, 0, 1>(x, &tmpFunc0, &tmpDer0, &tmpFunc1, &tmpDer1);
642 EXPECT_EQ(tstFunc0, tmpFunc0);
643 EXPECT_EQ(tstFunc1, tmpFunc1);
644 EXPECT_EQ(tstDer0, tmpDer0);
645 EXPECT_EQ(tstDer1, tmpDer1);
647 // two out of three, function/derivative-only
648 table.template evaluateFunction<3, 0, 1>(x, &tmpFunc0, &tmpFunc1);
649 table.template evaluateDerivative<3, 0, 1>(x, &tmpDer0, &tmpDer1);
650 EXPECT_EQ(tstFunc0, tmpFunc0);
651 EXPECT_EQ(tstFunc1, tmpFunc1);
652 EXPECT_EQ(tstDer0, tmpDer0);
653 EXPECT_EQ(tstDer1, tmpDer1);
655 // Test that scrambled order interpolation methods work
656 table.template evaluateFunctionAndDerivative<3, 2, 1, 0>(x, &tstFunc2, &tstDer2, &tstFunc1,
657 &tstDer1, &tstFunc0, &tstDer0);
658 EXPECT_EQ(tstFunc0, tmpFunc0);
659 EXPECT_EQ(tstFunc1, tmpFunc1);
660 EXPECT_EQ(tstFunc2, tmpFunc2);
661 EXPECT_EQ(tstDer0, tmpDer0);
662 EXPECT_EQ(tstDer1, tmpDer1);
663 EXPECT_EQ(tstDer2, tmpDer2);
665 // Test scrambled order for function/derivative-only methods
666 table.template evaluateFunction<3, 2, 1, 0>(x, &tmpFunc2, &tmpFunc1, &tmpFunc0);
667 table.template evaluateDerivative<3, 2, 1, 0>(x, &tmpDer2, &tmpDer1, &tmpDer0);
668 EXPECT_EQ(tstFunc0, tmpFunc0);
669 EXPECT_EQ(tstFunc1, tmpFunc1);
670 EXPECT_EQ(tstFunc2, tmpFunc2);
671 EXPECT_EQ(tstDer0, tmpDer0);
672 EXPECT_EQ(tstDer1, tmpDer1);
673 EXPECT_EQ(tstDer2, tmpDer2);
676 #if GMX_SIMD_HAVE_REAL
677 TYPED_TEST(SplineTableTest, Simd)
679 std::pair<real, real> range(0.2, 1.0);
680 TypeParam table({ { "LJ12", lj12Function, lj12Derivative } }, range);
682 // We already test that the SIMD operations handle the different elements
683 // correctly in the SIMD module, so here we only test that interpolation
684 // works for a single value in the middle of the interval
686 real x = 0.5 * (range.first + range.second);
687 real refFunc = lj12Function(x);
688 real refDer = lj12Derivative(x);
689 SimdReal tstFunc, tstDer;
690 real funcErr, derErr;
691 alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
693 table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc, &tstDer);
695 store(alignedMem, tstFunc);
696 funcErr = std::abs(alignedMem[0] - refFunc) / std::abs(refFunc);
698 store(alignedMem, tstDer);
699 derErr = std::abs(alignedMem[0] - refDer) / std::abs(refDer);
701 EXPECT_LT(funcErr, TypeParam::defaultTolerance);
702 EXPECT_LT(derErr, TypeParam::defaultTolerance);
705 TYPED_TEST(SplineTableTest, SimdTwoFunctions)
707 std::pair<real, real> range(0.2, 2.0);
709 TypeParam table({ { "LJ6", lj6Function, lj6Derivative }, { "LJ12", lj12Function, lj12Derivative } },
712 // We already test that the SIMD operations handle the different elements
713 // correctly in the SIMD module, so here we only test that interpolation
714 // works for a single value in the middle of the interval
716 real x = 0.5 * (range.first + range.second);
717 real refFunc0 = lj6Function(x);
718 real refDer0 = lj6Derivative(x);
719 real refFunc1 = lj12Function(x);
720 real refDer1 = lj12Derivative(x);
721 SimdReal tstFunc0, tstDer0;
722 SimdReal tstFunc1, tstDer1;
723 real funcErr0, derErr0;
724 real funcErr1, derErr1;
725 alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
727 table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
729 store(alignedMem, tstFunc0);
730 funcErr0 = std::abs(alignedMem[0] - refFunc0) / std::abs(refFunc0);
732 store(alignedMem, tstDer0);
733 derErr0 = std::abs(alignedMem[0] - refDer0) / std::abs(refDer0);
735 store(alignedMem, tstFunc1);
736 funcErr1 = std::abs(alignedMem[0] - refFunc1) / std::abs(refFunc1);
738 store(alignedMem, tstDer1);
739 derErr1 = std::abs(alignedMem[0] - refDer1) / std::abs(refDer1);
741 EXPECT_LT(funcErr0, TypeParam::defaultTolerance);
742 EXPECT_LT(derErr0, TypeParam::defaultTolerance);
743 EXPECT_LT(funcErr1, TypeParam::defaultTolerance);
744 EXPECT_LT(derErr1, TypeParam::defaultTolerance);
748 #if GMX_SIMD_HAVE_REAL && !defined NDEBUG
749 TYPED_TEST(SplineTableTest, CatchesOutOfRangeValuesSimd)
751 std::pair<real, real> range(0.2, 1.0);
752 TypeParam table({ { "LJ12", lj12Function, lj12Derivative } }, range);
753 SimdReal x, func, der;
755 AlignedArray<real, GMX_SIMD_REAL_WIDTH> alignedMem;
757 alignedMem.fill(range.first);
758 // Make position 1 incorrect if width>=2, otherwise position 0
759 // range.first-GMX_REAL_EPS is not invalid. See comment in table.
760 alignedMem[(GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = -GMX_REAL_EPS;
761 x = load<SimdReal>(alignedMem);
763 EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
765 // Make position 1 incorrect if width>=2, otherwise position 0
766 alignedMem[(GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = range.second;
767 x = load<SimdReal>(alignedMem);
769 EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
772 TYPED_TEST(SplineTableTest, AcceptsInRangeValuesSimd)
774 std::pair<real, real> range(0.2, 1.0);
775 TypeParam table({ { "LJ12", lj12Function, lj12Derivative } }, range);
776 SimdReal x, func, der;
778 alignas(GMX_SIMD_ALIGNMENT) real alignedMem[GMX_SIMD_REAL_WIDTH];
780 // Test all values between 0 and range.second
781 for (std::size_t i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
783 alignedMem[i] = range.second * (1.0 - GMX_REAL_EPS) * i / (GMX_SIMD_REAL_WIDTH - 1);
785 x = load<SimdReal>(alignedMem);
787 EXPECT_NO_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der));