endif()
if (GMX_COMPILER_WARNINGS)
# 177: function/variable ".." was declared but never referenced
+# 411: class defines no constructor to initialize the following (incorrect for struct, initializer list works)
# 593: variable ".." was set but never used
# 981: operands are evaluated in unspecified order
#1418: external function definition with no prior declaration
#3280: declaration hides member ".."
#11074: Inlining inhibited by limit max-size(/max-total-size)
#11076: To get full report use -opt-report=3 -opt-report-phase ipo (shown for previous remark)
- GMX_TEST_CFLAG(CFLAGS_WARN "-w3 -wd177 -wd593 -wd981 -wd1418 -wd1419 -wd1572 -wd1599 -wd2259 -wd2415 -wd2547 -wd2557 -wd3280 -wd11074 -wd11076" GMXC_CFLAGS)
+ GMX_TEST_CFLAG(CFLAGS_WARN "-w3 -wd177 -wd411 -wd593 -wd981 -wd1418 -wd1419 -wd1572 -wd1599 -wd2259 -wd2415 -wd2547 -wd2557 -wd3280 -wd11074 -wd11076" GMXC_CFLAGS)
endif()
GMX_TEST_CFLAG(CFLAGS_STDGNU "-std=gnu99" GMXC_CFLAGS)
GMX_TEST_CFLAG(CFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -static-intel" GMXC_CFLAGS_RELEASE)
GMX_TEST_CFLAG(CFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- GMX_TEST_CFLAG(CFLAGS_WARN "/W3 /wd177 /wd593 /wd981 /wd1418 /wd1419 /wd1572 /wd1599 /wd2259 /wd2415 /wd2547 /wd2557 /wd3280" GMXC_CFLAGS)
+GMX_TEST_CFLAG(CFLAGS_WARN "/W3 /wd177 /wd411 /wd593 /wd981 /wd1418 /wd1419 /wd1572 /wd1599 /wd2259 /wd2415 /wd2547 /wd2557 /wd3280" GMXC_CFLAGS)
endif()
GMX_TEST_CFLAG(CFLAGS_OPT "/Qip" GMXC_CFLAGS_RELEASE)
endif()
# 383: value copied to temporary, reference to temporary used
# 444: destructor for base class ".." is not virtual
#2282: unrecognized GCC pragma
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "-w3 -wd177 -wd383 -wd444 -wd981 -wd1418 -wd1572 -wd1599 -wd2259 -wd3280 -wd11074 -wd11076 -wd2282" GMXC_CXXFLAGS)
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN "-w3 -wd177 -wd383 -wd411 -wd444 -wd981 -wd1418 -wd1572 -wd1599 -wd2259 -wd3280 -wd11074 -wd11076 -wd2282" GMXC_CXXFLAGS)
endif()
GMX_TEST_CXXFLAG(CXXFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -static-intel" GMXC_CXXFLAGS_RELEASE)
GMX_TEST_CXXFLAG(CXXFLAGS_DEBUG "-O0" GMXC_CXXFLAGS_DEBUG)
GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/W3 /wd177 /wd383 /wd444 /wd981 /wd1418 /wd1572 /wd1599 /wd2259 /wd3280 /wd11074 /wd11076 /wd2282" GMXC_CXXFLAGS)
+GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/W3 /wd177 /wd383 /wd411 /wd444 /wd981 /wd1418 /wd1572 /wd1599 /wd2259 /wd3280 /wd11074 /wd11076 /wd2282" GMXC_CXXFLAGS)
endif()
GMX_TEST_CXXFLAG(CXXFLAGS_OPT "/Qip" GMXC_CXXFLAGS_RELEASE)
endif()
src/gromacs/simd/tests/scalar.cpp: warning: includes "simd.h" unnecessarily
src/gromacs/simd/tests/scalar_math.cpp: warning: includes "simd.h" unnecessarily
src/gromacs/simd/tests/scalar_util.cpp: warning: includes "simd.h" unnecessarily
+src/gromacs/tables/cubicsplinetable.h: warning: includes "simd.h" unnecessarily
+src/gromacs/tables/quadraticsplinetable.h: warning: includes "simd.h" unnecessarily
# These are specific to Folding@Home, and easiest to suppress here
*: warning: includes non-local file as "corewrap.h"
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2015, by the GROMACS development team, led by
+# Copyright (c) 2015,2016, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${TABLE_SOURCES} PARENT_SCOPE)
if (BUILD_TESTING)
-# add_subdirectory(tests)
+ add_subdirectory(tests)
endif()
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements classes for cubic spline table functions
+ *
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#include "gmxpre.h"
+
+#include "cubicsplinetable.h"
+
+#include <cmath>
+
+#include <algorithm>
+#include <functional>
+#include <initializer_list>
+#include <utility>
+#include <vector>
+
+#include "gromacs/tables/tableinput.h"
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/real.h"
+
+#include "splineutil.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+/*! \brief Calculate table elements from function/derivative data
+ *
+ * \param functionValue0 Function value for the present table index
+ * \param functionValue1 Function value for the next table index
+ * \param derivativeValue0 Derivative value for the present table index
+ * \param derivativeValue1 Derivative value for the next table index
+ * \param spacing Distance between table points
+ * \param Y Function value for table index
+ * \param F Component to multiply with offset eps
+ * \param G Component to multiply with eps^2
+ * \param H Component to multiply with eps^3
+ */
+void
+calculateCubicSplineCoefficients(double functionValue0,
+ double functionValue1,
+ double derivativeValue0,
+ double derivativeValue1,
+ double spacing,
+ double *Y,
+ double *F,
+ double *G,
+ double *H)
+{
+ *Y = functionValue0;
+ *F = spacing * derivativeValue0;
+ *G = 3.0*( functionValue1 - functionValue0) - spacing * (derivativeValue1 + 2.0 * derivativeValue0);
+ *H = -2.0*( functionValue1 - functionValue0) + spacing * (derivativeValue1 + derivativeValue0);
+}
+
+/*! \brief Perform cubic spline interpolation in interval from function/derivative
+ *
+ * \param functionValue0 Function value for the present table index
+ * \param functionValue1 Function value for the next table index
+ * \param derivativeValue0 Derivative value for the present table index
+ * \param derivativeValue1 Derivative value for the next table index
+ * \param spacing Distance between table points
+ * \param eps Offset from lower table point for evaluation
+ * \param[out] interpolatedFunctionValue Output function value
+ * \param[out] interpolatedDerivativeValue Output derivative value
+ */
+void
+cubicSplineInterpolationFromFunctionAndDerivative(double functionValue0,
+ double functionValue1,
+ double derivativeValue0,
+ double derivativeValue1,
+ double spacing,
+ double eps,
+ double *interpolatedFunctionValue,
+ double *interpolatedDerivativeValue)
+{
+ double Y, F, G, H;
+
+ calculateCubicSplineCoefficients(functionValue0, functionValue1,
+ derivativeValue0, derivativeValue1,
+ spacing,
+ &Y, &F, &G, &H);
+
+ double Fp = fma(fma(H, eps, G), eps, F);
+
+ *interpolatedFunctionValue = fma(Fp, eps, Y);
+ *interpolatedDerivativeValue = fma(eps, fma(2.0*eps, H, G), Fp)/spacing;
+}
+
+
+
+/*! \brief Construct the data for a single cubic table from analytical functions
+ *
+ * \param[in] function Analytical functiojn
+ * \param[in] derivative Analytical derivative
+ * \param[in] range Upper/lower limit of region to tabulate
+ * \param[in] spacing Distance between table points
+ * \param[out] yfghTableData Output cubic spline table with Y,F,G,H entries
+ */
+void
+fillSingleCubicSplineTableData(const std::function<double(double)> &function,
+ const std::function<double(double)> &derivative,
+ const std::pair<real, real> &range,
+ double spacing,
+ std::vector<real> *yfghTableData)
+{
+ int endIndex = range.second / spacing + 2;
+
+ yfghTableData->resize(4*endIndex);
+
+ double maxMagnitude = 0.0001*GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
+
+ for (int i = endIndex - 1; i >= 0; i--)
+ {
+ double x = i * spacing;
+ double tmpFunctionValue;
+ double tmpDerivativeValue;
+ double nextHigherFunction;
+ double nextHigherDerivative;
+ double Y, F, G, H;
+
+ if (range.first > 0 && i == 0)
+ {
+ // Avoid x==0 if it is not in the range, since it can lead to
+ // singularities even if the value for i==1 was within or required magnitude
+ functionIsInRange = false;
+ }
+
+ if (functionIsInRange)
+ {
+ tmpFunctionValue = function(x);
+ tmpDerivativeValue = derivative(x);
+ nextHigherFunction = ((i+1) < endIndex) ? function(x+spacing) : 0.0;
+ nextHigherDerivative = ((i+1) < endIndex) ? derivative(x+spacing) : 0.0;
+
+ if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
+ {
+ functionIsInRange = false; // Once this happens, it never resets to true again
+ }
+ }
+
+ if (functionIsInRange)
+ {
+ calculateCubicSplineCoefficients(tmpFunctionValue, nextHigherFunction,
+ tmpDerivativeValue, nextHigherDerivative,
+ spacing,
+ &Y, &F, &G, &H);
+ lastIndexInRange--;
+ }
+ else
+ {
+ double lastIndexY = (*yfghTableData)[4*lastIndexInRange];
+ double lastIndexF = (*yfghTableData)[4*lastIndexInRange + 1];
+
+ Y = lastIndexY + lastIndexF * (i - lastIndexInRange);
+ F = lastIndexF;
+ G = 0.0;
+ H = 0.0;
+ }
+
+ (*yfghTableData)[4*i ] = Y;
+ (*yfghTableData)[4*i+1] = F;
+ (*yfghTableData)[4*i+2] = G;
+ (*yfghTableData)[4*i+3] = H;
+ }
+}
+
+
+/*! \brief Construct the data for a single cubic table from vector data
+ *
+ * \param[in] function Input vector with function data
+ * \param[in] derivative Input vector with derivative data
+ * \param[in] inputSpacing Distance between points in input vectors
+ * \param[in] range Upper/lower limit of region to tabulate
+ * \param[in] spacing Distance between table points
+ * \param[out] yfghTableData Output cubic spline table with Y,F,G,H entries
+ */
+void
+fillSingleCubicSplineTableData(ConstArrayRef<double> function,
+ ConstArrayRef<double> derivative,
+ double inputSpacing,
+ const std::pair<real, real> &range,
+ double spacing,
+ std::vector<real> *yfghTableData)
+{
+ int endIndex = range.second / spacing + 2;
+
+ std::vector<double> tmpFunction(endIndex);
+ std::vector<double> tmpDerivative(endIndex);
+
+ double maxMagnitude = 0.0001*GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
+
+ // Interpolate function and derivative values in positions needed for output
+ for (int i = endIndex - 1; i >= 0; i--)
+ {
+ double x = i * spacing;
+ double xtab = x / inputSpacing;
+ int index = xtab;
+ double eps = xtab - index;
+
+ if (range.first > 0 && i == 0)
+ {
+ // Avoid x==0 if it is not in the range, since it can lead to
+ // singularities even if the value for i==1 was within or required magnitude
+ functionIsInRange = false;
+ }
+
+ if (functionIsInRange && (std::abs(function[index]) > maxMagnitude || std::abs(derivative[index]) > maxMagnitude))
+ {
+ functionIsInRange = false; // Once this happens, it never resets to true again
+ }
+
+ if (functionIsInRange)
+ {
+ cubicSplineInterpolationFromFunctionAndDerivative(function[index],
+ function[index+1],
+ derivative[index],
+ derivative[index+1],
+ inputSpacing,
+ eps,
+ &(tmpFunction[i]),
+ &(tmpDerivative[i]));
+ lastIndexInRange--;
+ }
+ else
+ {
+ double lastIndexFunction = tmpFunction[lastIndexInRange];
+ double lastIndexDerivative = tmpDerivative[lastIndexInRange];
+ tmpFunction[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
+ tmpDerivative[i] = lastIndexDerivative;
+ }
+ }
+
+ yfghTableData->resize(4*endIndex);
+
+ for (int i = 0; i < endIndex; i++)
+ {
+ double Y, F, G, H;
+
+ double nextFunction = ((i+1) < endIndex) ? tmpFunction[i+1] : 0.0;
+ double nextDerivative = ((i+1) < endIndex) ? tmpDerivative[i+1] : 0.0;
+
+ calculateCubicSplineCoefficients(tmpFunction[i], nextFunction,
+ tmpDerivative[i], nextDerivative,
+ spacing,
+ &Y, &F, &G, &H);
+ (*yfghTableData)[4*i ] = Y;
+ (*yfghTableData)[4*i+1] = F;
+ (*yfghTableData)[4*i+2] = G;
+ (*yfghTableData)[4*i+3] = H;
+ }
+
+}
+
+} // namespace anonymous
+
+
+
+#if GMX_DOUBLE
+const real
+CubicSplineTable::defaultTolerance = 1e-10;
+#else
+const real
+CubicSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
+#endif
+
+CubicSplineTable::CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
+ const std::pair<real, real> &range,
+ real tolerance)
+ : numFuncInTable_(analyticalInputList.size()), range_(range)
+{
+ // Sanity check on input values
+ if (range_.first < 0.0 || (range_.second-range_.first) < 0.001)
+ {
+ GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ }
+
+ if (tolerance < GMX_REAL_EPS)
+ {
+ GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
+ }
+
+ double minQuotient = GMX_REAL_MAX;
+
+ // loop over all functions to find smallest spacing
+ for (auto thisFuncInput : analyticalInputList)
+ {
+ try
+ {
+ internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, range_);
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating cubic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ // Calculate the required table spacing h. The error we make with linear interpolation
+ // of the derivative will be described by the third-derivative correction term.
+ // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
+ // where f'/f''' is the first and third derivative of the function, respectively.
+
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(thisFuncInput.derivative, range_);
+
+ minQuotient = std::min(minQuotient, thisMinQuotient);
+ }
+
+ double spacing = 0.5 * std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
+
+ tableScale_ = 1.0 / spacing;
+
+ if (range_.second * tableScale_ > 2e6)
+ {
+ GMX_THROW(ToleranceError("Over a million points would be required for table; decrease range or increase tolerance"));
+ }
+
+ // Loop over all tables again.
+ // Here we create the actual table for each function, and then
+ // combine them into a multiplexed table function.
+ std::size_t funcIndex = 0;
+
+ for (auto thisFuncInput : analyticalInputList)
+ {
+ try
+ {
+ std::vector<real> tmpYfghTableData;
+
+ fillSingleCubicSplineTableData(thisFuncInput.function,
+ thisFuncInput.derivative,
+ range_,
+ spacing,
+ &tmpYfghTableData);
+
+ internal::fillMultiplexedTableData(tmpYfghTableData, &yfghMultiTableData_,
+ 4, numFuncInTable_, funcIndex);
+
+ funcIndex++;
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating cubic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ }
+}
+
+
+CubicSplineTable::CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
+ const std::pair<real, real> &range,
+ real tolerance)
+ : numFuncInTable_(numericalInputList.size()), range_(range)
+{
+ // Sanity check on input values
+ if (range.first < 0.0 || (range.second-range.first) < 0.001)
+ {
+ GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ }
+
+ if (tolerance < GMX_REAL_EPS)
+ {
+ GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
+ }
+
+ double minQuotient = GMX_REAL_MAX;
+
+ // loop over all functions to find smallest spacing
+ for (auto thisFuncInput : numericalInputList)
+ {
+ try
+ {
+ // We do not yet know what the margin is, but we need to test that we at least cover
+ // the requested range before starting to calculate derivatives
+ if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
+ {
+ GMX_THROW(InconsistentInputError("Table input vectors must cover requested range, and a margin beyond the upper endpoint"));
+ }
+
+ if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
+ {
+ GMX_THROW(InconsistentInputError("Function and derivative vectors have different lengths"));
+ }
+
+ internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating cubic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ // Calculate the required table spacing h. The error we make with linear interpolation
+ // of the derivative will be described by the third-derivative correction term.
+ // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
+ // where f'/f''' is the first and third derivative of the function, respectively.
+
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(thisFuncInput.derivative, thisFuncInput.spacing, range_);
+
+ minQuotient = std::min(minQuotient, thisMinQuotient);
+ }
+
+ double spacing = std::cbrt(72.0 * std::sqrt(3.0) * tolerance * minQuotient);
+
+ tableScale_ = 1.0 / spacing;
+
+ if (range_.second * tableScale_ > 1e6)
+ {
+ GMX_THROW(ToleranceError("Requested tolerance would require over a million points in table"));
+ }
+
+ // Loop over all tables again.
+ // Here we create the actual table for each function, and then
+ // combine them into a multiplexed table function.
+ std::size_t funcIndex = 0;
+
+ for (auto thisFuncInput : numericalInputList)
+ {
+ try
+ {
+ if (spacing < thisFuncInput.spacing)
+ {
+ GMX_THROW(ToleranceError("Input vector spacing cannot achieve tolerance requested"));
+ }
+
+ std::vector<real> tmpYfghTableData;
+
+ fillSingleCubicSplineTableData(thisFuncInput.function,
+ thisFuncInput.derivative,
+ thisFuncInput.spacing,
+ range,
+ spacing,
+ &tmpYfghTableData);
+
+ internal::fillMultiplexedTableData(tmpYfghTableData, &yfghMultiTableData_,
+ 4, numFuncInTable_, funcIndex);
+
+ funcIndex++;
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating cubic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ }
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ * \brief
+ * Declares classes for cubic spline table
+ *
+ * \inlibraryapi
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#ifndef GMX_TABLES_CUBICSPLINETABLE_H
+#define GMX_TABLES_CUBICSPLINETABLE_H
+
+#include <initializer_list>
+#include <vector>
+
+#include "gromacs/simd/simd.h"
+#include "gromacs/tables/tableinput.h"
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/classhelpers.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
+
+namespace gmx
+{
+
+/*! \libinternal \brief Cubic spline interpolation table.
+ *
+ * This class interpolates a function specified either as an analytical
+ * expression or from user-provided table data.
+ *
+ * At initialization, you provide the reference function of vectors
+ * as a list of tuples that contain a brief name, the function, and
+ * derivative for each function to tabulate. To create a table with
+ * two functions this initializer list can for instance look like
+ *
+ * { {"LJ6", lj6Func, lj6Der}, {"LJ12", lj12Func, lj12Der} }
+ *
+ * The names are only used so exceptions during initialization can
+ * be traced to a specific table.
+ *
+ * When interpolating, there are methods to interpolate either 1, 2, or 3
+ * functions in one go. By default these interpolation routines will
+ * operate on tables with the same number of functions as specified in
+ * the interpolation method (debug builds check that this is consistent with
+ * the table). However, it is also possible to use optional template
+ * parameters that specify the total number of functions in a table, and
+ * what function index to interpolate. For instance, to interpolate the
+ * derivative of the second function (i.e., index 1) in a
+ * multi-function-table with three functions in total, you can write
+ *
+ * table.evaluateDerivative<3,1>(x,&der);
+ *
+ * Here too, debug builds will check that the template parameters are
+ * consistent with the table.
+ *
+ * This class interpolates a function specified either as an analytical
+ * expression or from user-provided table data. The coefficients for each
+ * table point are precalculated such that we simply evaluate
+ *
+ * \f{eqnarray*}{
+ * V(x) = Y + F \epsilon + G \epsilon^2 + H \epsilon^3
+ * V'(x) = (F + 2 G \epsilon + 3 H \epsilon^2)/h
+ * \f}
+ *
+ * Where h is the spacing and epsilon the fractional offset from table point.
+ *
+ * While it is possible to create tables only from function values
+ * (i.e., no derivatives), it is recommended to provide derivatives for higher
+ * accuracy and to avoid issues with numerical differentiation. Note that the
+ * table input should be smooth, i.e. it should not contain noise e.g. from an
+ * (iterative) Boltzmann inversion procedure - you have been warned.
+ *
+ * \note This class is responsible for fundamental interpolation of any function,
+ * which might or might not correspond to a potential. For this reason
+ * both input and output derivatives are proper function derivatives, and
+ * we do not swap any signs to get forces directly from the table.
+ *
+ * \note There will be a small additional accuracy loss from the internal
+ * operation where we calculate the epsilon offset from the nearest table
+ * point, since the integer part we subtract can get large in those cases.
+ *
+ * While this is technically possible to solve with extended precision
+ * arithmetics, that would introduce extra instructions in some highly
+ * performance-sensitive code parts. For typical GROMACS interaction
+ * functions the derivatives will decay faster than the potential, which
+ * means it will never play any role. For other functions it will only
+ * cause a small increase in the relative error for arguments where the
+ * magnitude of the function or derivative is very small.
+ * Since we typically sum several results in GROMACS, this should never
+ * show up in any real cases, and for this reason we choose not to do
+ * the extended precision arithmetics.
+ *
+ * \note These routines are not suitable for table ranges starting far away
+ * from zero, since we allocate memory and calculate indices starting from
+ * range zero for efficiency reasons.
+ */
+class CubicSplineTable
+{
+ private:
+ /*! \brief Change that function value falls inside range when debugging
+ *
+ * \tparam T Lookup argument floating-point type, typically SimdReal or real.
+ * \param r Lookup argument to test
+ *
+ * \throws Debug builds will throw gmx::RangeError for values that are
+ * larger than the upper limit of the range, or smaller than 0.
+ * We allow the table to be called with arguments between 0 and
+ * the lower limit of the range, since this might in theory occur
+ * once-in-a-blue-moon with some algorithms.
+ */
+ template <typename T>
+ void
+ rangeCheck(T gmx_unused r) const
+ {
+#ifndef NDEBUG
+ // Check that all values fall in range when debugging
+ if (anyTrue( r < T(0.0) || T(range_.second) <= r ) )
+ {
+ GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
+ }
+#endif
+ }
+
+ public:
+
+ /*! \brief Default tolerance for cubic spline tables
+ *
+ * This is 10*GMX_FLOAT_EPS in single precision, and
+ * 1e-10 for double precision. It might not be worth setting
+ * this tolerance lower than 1e-10 in double precision, both because
+ * you will end up with very large tables, and because
+ * functions like r^-12 become so large for small values of r the
+ * table generation code will lead to some precision loss even
+ * in double precision.
+ */
+ static const real defaultTolerance;
+
+ /*! \brief Initialize table data from function
+ *
+ * \param analyticalInputList Initializer list with one or more functions to tabulate,
+ * specified as elements with a string description and
+ * the function as well as derivative. The function will also
+ * be called for values smaller than the lower limit of the
+ * range, but we avoid calling it for 0.0 if that value
+ * is not included in the range.
+ * Constructor will throw gmx::APIError for negative values.
+ * Due to the way the numerical derivative evaluation depends
+ * on machine precision internally, this range must be
+ * at least 0.001, or the constructor throws gmx::APIError.
+ * \param range Range over which the function will be tabulated.
+ * Constructor will throw gmx::APIError for negative values,
+ * or if the value/derivative vector does not cover the
+ * range.
+ * \param tolerance Requested accuracy of the table. This will be used to
+ * calculate the required internal spacing. If this cannot
+ * be achieved (for instance because the table would require
+ * too much memory) the constructor will throw gmx::ToleranceError.
+ *
+ * \note The functions are always defined in double precision to avoid
+ * losing accuracy when constructing tables.
+ *
+ * \note Since we fill the table for values below range.first, you can achieve
+ * a smaller table by using a smaller range where the tolerance has to be
+ * met, and accept that a few function calls below range.first do not
+ * quite reach the tolerance.
+ *
+ * \warning For efficiency reasons (since this code is used in some inner
+ * (kernels), we always allocate memory and calculate table indices
+ * for the complete interval [0,range.second], although the data will
+ * not be valid outside the definition range to avoid calling the
+ * function there. This means you should \a not use this class
+ * to tabulate functions for small ranges very far away from zero,
+ * since you would both waste a huge amount of memory and incur
+ * truncation errors when calculating the index.
+ *
+ * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
+ * and gmx::APIError for other incorrect input.
+ */
+ CubicSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
+ const std::pair<real, real> &range,
+ real tolerance = defaultTolerance);
+
+ /*! \brief Initialize table data from tabulated values and derivatives
+ *
+ * \param numericalInputList Initializer list with one or more functions to tabulate,
+ * specified as a string description, vectors with function and
+ * derivative values, and the input spacing. Data points are
+ * separated by the spacing parameter, starting from 0.
+ * Values below the lower limit of the range will be used to
+ * attempt defining the table, but we avoid using index 0
+ * unless 0.0 is included in the range. Some extra points beyond
+ * range.second are required to re-interpolate values, so add
+ * some margin. The constructor will throw gmx::APIError if the
+ * input vectors are too short to cover the requested range
+ * (and they must always be at least five points).
+ * \param range Range over which the function will be tabulated.
+ * Constructor will throw gmx::APIError for negative values,
+ * or if the value/derivative vector does not cover the
+ * range.
+ * \param tolerance Requested accuracy of the table. This will be used to
+ * calculate the required internal spacing and possibly
+ * re-interpolate. The constructor will throw
+ * gmx::ToleranceError if the input spacing is too coarse
+ * to achieve this accuracy.
+ *
+ * \note The input data vectors are always double precision to avoid
+ * losing accuracy when constructing tables.
+ *
+ * \note Since we fill the table for values below range.first, you can achieve
+ * a smaller table by using a smaller range where the tolerance has to be
+ * met, and accept that a few function calls below range.first do not
+ * quite reach the tolerance.
+ *
+ * \warning For efficiency reasons (since this code is used in some inner
+ * (kernels), we always allocate memory and calculate table indices
+ * for the complete interval [0,range.second], although the data will
+ * not be valid outside the definition range to avoid calling the
+ * function there. This means you should \a not use this class
+ * to tabulate functions for small ranges very far away from zero,
+ * since you would both waste a huge amount of memory and incur
+ * truncation errors when calculating the index.
+ */
+ CubicSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
+ const std::pair<real, real> &range,
+ real tolerance = defaultTolerance);
+
+
+ /************************************************************
+ * Evaluation methods for single functions *
+ ************************************************************/
+
+ /*! \brief Evaluate both function and derivative, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue Function value
+ * \param[out] derivativeValue Function derivative
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue,
+ T * derivativeValue) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex, tabIndex, &Y, &F, &G, &H);
+ *functionValue = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /*! \brief Evaluate function value only, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue Function value
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue) const
+ {
+ T der gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
+ }
+
+ /*! \brief Evaluate function derivative only, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue Function derivative
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /************************************************************
+ * Evaluation methods for two functions *
+ ************************************************************/
+
+ /*! \brief Evaluate both function and derivative, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue0 Interpolated value for first function
+ * \param[out] derivativeValue0 Interpolated derivative for first function
+ * \param[out] functionValue1 Interpolated value for second function
+ * \param[out] derivativeValue1 Interpolated derivative for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue0,
+ T * derivativeValue0,
+ T * functionValue1,
+ T * derivativeValue1) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex0, tabIndex, &Y, &F, &G, &H);
+ *functionValue0 = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex1, tabIndex, &Y, &F, &G, &H);
+ *functionValue1 = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /*! \brief Evaluate function value only, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue0 Interpolated value for first function
+ * \param[out] functionValue1 Interpolated value for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue0,
+ T * functionValue1) const
+ {
+ T der0 gmx_unused;
+ T der1 gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(r, functionValue0, &der0, functionValue1, &der1);
+ }
+
+ /*! \brief Evaluate function derivative only, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue0 Interpolated derivative for first function
+ * \param[out] derivativeValue1 Interpolated derivative for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue0,
+ T * derivativeValue1) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /************************************************************
+ * Evaluation methods for three functions *
+ ************************************************************/
+
+
+ /*! \brief Evaluate both function and derivative, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue0 Interpolated value for first function
+ * \param[out] derivativeValue0 Interpolated derivative for first function
+ * \param[out] functionValue1 Interpolated value for second function
+ * \param[out] derivativeValue1 Interpolated derivative for second function
+ * \param[out] functionValue2 Interpolated value for third function
+ * \param[out] derivativeValue2 Interpolated derivative for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue0,
+ T * derivativeValue0,
+ T * functionValue1,
+ T * derivativeValue1,
+ T * functionValue2,
+ T * derivativeValue2) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex0, tabIndex, &Y, &F, &G, &H);
+ *functionValue0 = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex1, tabIndex, &Y, &F, &G, &H);
+ *functionValue1 = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4*funcIndex2, tabIndex, &Y, &F, &G, &H);
+ *functionValue2 = fma(fma(fma(H, eps, G), eps, F), eps, Y);
+ *derivativeValue2 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /*! \brief Evaluate function value only, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue0 Interpolated value for first function
+ * \param[out] functionValue1 Interpolated value for second function
+ * \param[out] functionValue2 Interpolated value for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue0,
+ T * functionValue1,
+ T * functionValue2) const
+ {
+ T der0 gmx_unused;
+ T der1 gmx_unused;
+ T der2 gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(r, functionValue0, &der0, functionValue1, &der1, functionValue2, &der2);
+ }
+
+ /*! \brief Evaluate function derivative only, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue0 Interpolated derivative for first function
+ * \param[out] derivativeValue1 Interpolated derivative for second function
+ * \param[out] derivativeValue2 Interpolated derivative for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue0,
+ T * derivativeValue1,
+ T * derivativeValue2) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T Y, F, G, H;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex0, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue0 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex1, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue1 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(yfghMultiTableData_.data() + 4 * funcIndex2, tabIndex, &Y, &F, &G, &H);
+ *derivativeValue2 = tableScale_ * fma(fma(T(3.0)*H, eps, T(2.0)*G), eps, F);
+ }
+
+ /*! \brief Return the table spacing (distance between points)
+ *
+ * You should never have to use this for normal code, but due to the
+ * way tables are constructed internally we need this in the unit tests
+ * to check relative tolerances over each interval.
+ *
+ * \return table spacing.
+ */
+ real
+ tableSpacing() const { return 1.0 / tableScale_; }
+
+ private:
+
+ std::size_t numFuncInTable_; //!< Number of separate tabluated functions
+ std::pair<real, real> range_; //!< Range for which table evaluation is allowed
+ real tableScale_; //!< Table scale (inverse of spacing between points)
+
+ /*! \brief Vector with combined table data to save calculations after lookup.
+ *
+ * For table point i, this vector contains the four coefficients
+ * Y,F,G,H that we use to express the function value as
+ * V(x) = Y + F e + G e^2 + H e^3, where e is the epsilon offset from
+ * the nearest table point.
+ *
+ * To allow aligned SIMD loads we need to use an aligned allocator for
+ * this container.
+ */
+ std::vector<real, AlignedAllocator<real> > yfghMultiTableData_;
+
+ // There should never be any reason to copy the table since it is read-only
+ GMX_DISALLOW_COPY_AND_ASSIGN(CubicSplineTable);
+};
+
+
+} // namespace gmx
+
+#endif // GMX_TABLES_CUBICSPLINETABLE_H
#ifndef GMX_TABLES_FORCETABLE_H
#define GMX_TABLES_FORCETABLE_H
+/*! \libinternal \file
+ * \brief
+ * Old routines for table generation (will eventually be replaced)
+ *
+ * \inlibraryapi
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+
#include <cstdio>
#include "gromacs/mdtypes/fcdata.h"
#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/utility/real.h"
+/*! \brief Flag to select user tables for make_tables */
#define GMX_MAKETABLES_FORCEUSER (1<<0)
+/*! \brief Flag to only make 1,4 pair tables for make_tables */
#define GMX_MAKETABLES_14ONLY (1<<1)
-/* Index in the tables that says which function to use */
+/*! \brief Enumerated type to describe the interaction types in a table */
enum {
- etiCOUL, etiLJ6, etiLJ12, etiNR
+ etiCOUL, //!< Coulomb
+ etiLJ6, //!< Dispersion
+ etiLJ12, //!< Repulsion
+ etiNR //!< Total number of interaction types
};
-typedef double (*real_space_grid_contribution_computer)(double, double);
-/* Function pointer used to tell table_spline3_fill_ewald_lr whether it
+/*! \brief Function pointer to calculate the grid contribution for coulomb/LJ
+ *
+ * Used to tell table_spline3_fill_ewald_lr whether it
* should calculate the grid contribution for electrostatics or LJ.
*/
+typedef double (*real_space_grid_contribution_computer)(double, double);
+
+/*! \brief Fill tables with the Ewald long-range force interaction
+ *
+ * Fill tables of ntab points with spacing dr with the ewald long-range
+ * (mesh) force.
+ * There are three separate tables with format FDV0, F, and V.
+ * This function interpolates the Ewald mesh potential contribution
+ * with coefficient beta using a quadratic spline.
+ * The force can then be interpolated linearly.
+ *
+ * \param table_F Force table
+ * \param table_V Potential table
+ * \param table_FDV0 Combined table optimized for SIMD loads
+ * \param ntab Number of points in tables
+ * \param dx Spacing
+ * \param beta Ewald splitting paramter
+ * \param v_lr Pointer to function calculating real-space grid contribution
+ */
void table_spline3_fill_ewald_lr(real *table_F,
real *table_V,
real *table_FDV0,
double dx,
real beta,
real_space_grid_contribution_computer v_lr);
-/* Fill tables of ntab points with spacing dr with the ewald long-range
- * (mesh) force.
- * There are three separate tables with format FDV0, F, and V.
- * This function interpolates the Ewald mesh potential contribution
- * with coefficient beta using a quadratic spline.
- * The force can then be interpolated linearly.
- */
+/*! \brief Compute scaling for the Ewald quadratic spline tables.
+ *
+ * \param ic Pointer to interaction constant structure
+ * \return The scaling factor
+ */
real ewald_spline3_table_scale(const interaction_const_t *ic);
-/* Return the scaling for the Ewald quadratic spline tables. */
+/*! \brief Return the real space grid contribution for Ewald
+ *
+ * \param beta Ewald splitting parameter
+ * \param r Distance for which to calculate the real-space contrib
+ * \return Real space grid contribution for Ewald electrostatics
+ */
double v_q_ewald_lr(double beta, double r);
-/* Return the real space grid contribution for Ewald*/
+/*! \brief Return the real space grid contribution for LJ-Ewald
+ *
+ * \param beta Ewald splitting parameter
+ * \param r Distance for which to calculate the real-space contrib
+ * \return Real space grid contribution for Ewald Lennard-Jones interaction
+ */
double v_lj_ewald_lr(double beta, double r);
-/* Return the real space grid contribution for LJ-Ewald*/
+/*! \brief Return tables for inner loops.
+ *
+ * \param fp Log file pointer
+ * \param fr Force record
+ * \param fn File name from which to read user tables
+ * \param rtab Largest interaction distance to tabulate
+ * \param flags Flags to select table settings
+ *
+ * \return Pointer to inner loop table structure
+ */
t_forcetable *make_tables(FILE *fp,
const t_forcerec *fr,
const char *fn, real rtab, int flags);
-/* Return tables for inner loops. */
-bondedtable_t make_bonded_table(FILE *fplog, const char *fn, int angle);
-/* Return a table for bonded interactions,
- * angle should be: bonds 0, angles 1, dihedrals 2
+/*! \brief Return a table for bonded interactions,
+ *
+ * \param fplog Pointer to log file
+ * \param fn File name
+ * \param angle Type of angle: bonds 0, angles 1, dihedrals 2
+ * \return New bonded table datatype
*/
+bondedtable_t make_bonded_table(FILE *fplog, const char *fn, int angle);
-/* Return a table for GB calculations */
+/*! \brief Return a table for GB calculations
+ *
+ * \param fr Force record
+ * \return Pointer to new gb table structure
+ */
t_forcetable *make_gb_table(const t_forcerec *fr);
/*! \brief Construct and return tabulated dispersion and repulsion interactions
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements classes for quadratic spline table functions
+ *
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#include "gmxpre.h"
+
+#include "quadraticsplinetable.h"
+
+#include <cmath>
+
+#include <algorithm>
+#include <functional>
+#include <initializer_list>
+#include <utility>
+#include <vector>
+
+#include "gromacs/tables/tableinput.h"
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/real.h"
+
+#include "splineutil.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+/*! \brief Construct the data for a single quadratic table from analytical functions
+ *
+ * \param[in] function Analytical functiojn
+ * \param[in] derivative Analytical derivative
+ * \param[in] range Upper/lower limit of region to tabulate
+ * \param[in] spacing Distance between table points
+ * \param[out] functionTableData Output table with function data
+ * \param[out] derivativeTableData OUtput table with (adjusted) derivative data
+ */
+void
+fillSingleQuadraticSplineTableData(const std::function<double(double)> &function,
+ const std::function<double(double)> &derivative,
+ const std::pair<real, real> &range,
+ double spacing,
+ std::vector<real> *functionTableData,
+ std::vector<real> *derivativeTableData)
+{
+ std::size_t endIndex = range.second / spacing + 2;
+
+ functionTableData->resize(endIndex);
+ derivativeTableData->resize(endIndex);
+
+ double maxMagnitude = 0.0001*GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
+
+ for (int i = endIndex - 1; i >= 0; i--)
+ {
+ double x = i * spacing;
+ double tmpFunctionValue;
+ double tmpDerivativeValue;
+
+ if (range.first > 0 && i == 0)
+ {
+ // Avoid x==0 if it is not in the range, since it can lead to
+ // singularities even if the value for i==1 was within or required magnitude
+ functionIsInRange = false;
+ }
+
+ if (functionIsInRange)
+ {
+ tmpFunctionValue = function(x);
+
+ // Calculate third derivative term (2nd derivative of the derivative)
+ // Make sure we stay in range. In practice this means we use one-sided
+ // interpolation at the interval endpoints (indentical to an offset for 3-point formula)
+ const double h = std::pow( GMX_DOUBLE_EPS, 0.25 );
+ double y = std::min( std::max(x, range.first + h), range.second - h);
+ double thirdDerivativeValue = ( derivative(y+h) - 2.0 * derivative(y) + derivative(y-h) ) / ( h * h );
+
+ tmpDerivativeValue = derivative(x) - spacing * spacing * thirdDerivativeValue / 12.0;
+
+ if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
+ {
+ functionIsInRange = false; // Once this happens, it never resets to true again
+ }
+ }
+
+ if (functionIsInRange)
+ {
+ (*functionTableData)[i] = tmpFunctionValue;
+ (*derivativeTableData)[i] = tmpDerivativeValue;
+ lastIndexInRange--;
+ }
+ else
+ {
+ // Once the function or derivative (more likely) has reached very large values,
+ // we simply make a linear function from the last in-range value of the derivative.
+ double lastIndexFunction = (*functionTableData)[lastIndexInRange];
+ double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
+ (*functionTableData)[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
+ (*derivativeTableData)[i] = lastIndexDerivative;
+ }
+ }
+}
+
+
+/*! \brief Construct the data for a single quadratic table from vector data
+ *
+ * \param[in] function Input vector with function data
+ * \param[in] derivative Input vector with derivative data
+ * \param[in] inputSpacing Distance between points in input vectors
+ * \param[in] range Upper/lower limit of region to tabulate
+ * \param[in] spacing Distance between table points
+ * \param[out] functionTableData Output table with function data
+ * \param[out] derivativeTableData OUtput table with (adjusted) derivative data
+ */
+void
+fillSingleQuadraticSplineTableData(ConstArrayRef<double> function,
+ ConstArrayRef<double> derivative,
+ double inputSpacing,
+ const std::pair<real, real> &range,
+ double spacing,
+ std::vector<real> *functionTableData,
+ std::vector<real> *derivativeTableData)
+{
+ std::size_t endIndex = range.second / spacing + 2;
+
+ functionTableData->resize(endIndex);
+ derivativeTableData->resize(endIndex);
+
+ std::vector<double> thirdDerivative(internal::vectorSecondDerivative(derivative, inputSpacing));
+
+ double maxMagnitude = 0.0001*GMX_REAL_MAX;
+ bool functionIsInRange = true;
+ std::size_t lastIndexInRange = endIndex - 1;
+
+ for (int i = endIndex - 1; i >= 0; i--)
+ {
+ double x = i * spacing;
+ double tmpFunctionValue;
+ double tmpDerivativeValue;
+
+ if (range.first > 0 && i == 0)
+ {
+ // Avoid x==0 if it is not in the range, since it can lead to
+ // singularities even if the value for i==1 was within or required magnitude
+ functionIsInRange = false;
+ }
+
+ if (functionIsInRange)
+ {
+ // Step 1: Interpolate the function value at x from input table.
+ double inputXTab = x / inputSpacing;
+ int inputIndex = inputXTab;
+ double inputEps = inputXTab - inputIndex;
+
+ // Linear interpolation of input derivative and third derivative
+ double thirdDerivativeValue = (1.0 - inputEps) * thirdDerivative[inputIndex] + inputEps * thirdDerivative[inputIndex+1];
+ double derivativeValue = (1.0 - inputEps) * derivative[inputIndex] + inputEps * derivative[inputIndex+1];
+
+ // Quadratic interpolation for function value
+ tmpFunctionValue = function[inputIndex] + 0.5 * (derivative[inputIndex] + derivativeValue) * inputEps * inputSpacing;
+ tmpDerivativeValue = derivativeValue - spacing * spacing * thirdDerivativeValue / 12.0;
+
+ if (std::abs(tmpFunctionValue) > maxMagnitude || std::abs(tmpDerivativeValue) > maxMagnitude)
+ {
+ functionIsInRange = false; // Once this happens, it never resets to true again
+ }
+ }
+
+ if (functionIsInRange)
+ {
+ (*functionTableData)[i] = tmpFunctionValue;
+ (*derivativeTableData)[i] = tmpDerivativeValue;
+ lastIndexInRange--;
+ }
+ else
+ {
+ // Once the function or derivative (more likely) has reached very large values,
+ // we simply make a linear function from the last in-range value of the derivative.
+ double lastIndexFunction = (*functionTableData)[lastIndexInRange];
+ double lastIndexDerivative = (*derivativeTableData)[lastIndexInRange];
+ (*functionTableData)[i] = lastIndexFunction + lastIndexDerivative * (i - lastIndexInRange) * spacing;
+ (*derivativeTableData)[i] = lastIndexDerivative;
+ }
+ }
+}
+
+/*! \brief Create merged DDFZ vector from function & derivative data
+ *
+ * \param functionTableData Function values
+ * \param derivativeTableData Derivative values. We have already subtracted the
+ * small third derivative component when calling this
+ * function, but in practice it is just an arbitrary
+ * vector here.
+ * \param ddfzTableData Vector four times longer, filled with
+ * the derivative, the difference to the next derivative
+ * point, the function value, and zero.
+ *
+ * \throws If the vector lengths do not match.
+ */
+void
+fillDdfzTableData(const std::vector<real> &functionTableData,
+ const std::vector<real> &derivativeTableData,
+ std::vector<real> *ddfzTableData)
+{
+ GMX_ASSERT(functionTableData.size() == derivativeTableData.size(), "Mismatching vector lengths");
+
+ std::size_t points = functionTableData.size();
+
+ ddfzTableData->resize(4 * points);
+
+ for (std::size_t i = 0; i < points; i++)
+ {
+ (*ddfzTableData)[4*i] = derivativeTableData[i];
+
+ double nextDerivative = ( i < functionTableData.size() - 1 ) ? derivativeTableData[i+1] : 0.0;
+
+ (*ddfzTableData)[4*i + 1] = nextDerivative - derivativeTableData[i];
+ (*ddfzTableData)[4*i + 2] = functionTableData[i];
+ (*ddfzTableData)[4*i + 3] = 0.0;
+ }
+}
+
+} // namespace anonymous
+
+
+
+const real
+QuadraticSplineTable::defaultTolerance = 10.0 * GMX_FLOAT_EPS;
+
+
+QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
+ const std::pair<real, real> &range,
+ real tolerance)
+ : numFuncInTable_(analyticalInputList.size()), range_(range)
+{
+ // Sanity check on input values
+ if (range_.first < 0.0 || (range_.second-range_.first) < 0.001)
+ {
+ GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ }
+
+ if (tolerance < GMX_REAL_EPS)
+ {
+ GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
+ }
+
+ double minQuotient = GMX_REAL_MAX;
+
+ // loop over all functions to find smallest spacing
+ for (auto thisFuncInput : analyticalInputList)
+ {
+ try
+ {
+ internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, range_);
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ // Calculate the required table spacing h. The error we make with linear interpolation
+ // of the derivative will be described by the third-derivative correction term.
+ // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
+ // where f'/f''' is the first and third derivative of the function, respectively.
+
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(thisFuncInput.derivative, range_);
+
+ minQuotient = std::min(minQuotient, thisMinQuotient);
+ }
+
+ double spacing = std::sqrt(12.0 * tolerance * minQuotient);
+
+ halfSpacing_ = 0.5 * spacing;
+ tableScale_ = 1.0 / spacing;
+
+ if (range_.second * tableScale_ > 1e6)
+ {
+ GMX_THROW(ToleranceError("Over a million points would be required for table; decrease range or increase tolerance"));
+ }
+
+ // Loop over all tables again.
+ // Here we create the actual table for each function, and then
+ // combine them into a multiplexed table function.
+ std::size_t funcIndex = 0;
+
+ for (auto thisFuncInput : analyticalInputList)
+ {
+ try
+ {
+ std::vector<real> tmpFuncTableData;
+ std::vector<real> tmpDerTableData;
+ std::vector<real> tmpDdfzTableData;
+
+ fillSingleQuadraticSplineTableData(thisFuncInput.function,
+ thisFuncInput.derivative,
+ range_,
+ spacing,
+ &tmpFuncTableData,
+ &tmpDerTableData);
+
+ fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
+
+ internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_,
+ 1, numFuncInTable_, funcIndex);
+
+ internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_,
+ 4, numFuncInTable_, funcIndex);
+
+ funcIndex++;
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ }
+}
+
+
+QuadraticSplineTable::QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
+ const std::pair<real, real> &range,
+ real tolerance)
+ : numFuncInTable_(numericalInputList.size()), range_(range)
+{
+ // Sanity check on input values
+ if (range.first < 0.0 || (range.second-range.first) < 0.001)
+ {
+ GMX_THROW(InvalidInputError("Range to tabulate cannot include negative values and must span at least 0.001"));
+ }
+
+ if (tolerance < GMX_REAL_EPS)
+ {
+ GMX_THROW(ToleranceError("Table tolerance cannot be smaller than GMX_REAL_EPS"));
+ }
+
+ double minQuotient = GMX_REAL_MAX;
+
+ // loop over all functions to find smallest spacing
+ for (auto thisFuncInput : numericalInputList)
+ {
+ try
+ {
+ // We do not yet know what the margin is, but we need to test that we at least cover
+ // the requested range before starting to calculate derivatives
+ if (thisFuncInput.function.size() < range_.second / thisFuncInput.spacing + 1)
+ {
+ GMX_THROW(InconsistentInputError("Table input vectors must cover requested range, and a margin beyond the upper endpoint"));
+ }
+
+ if (thisFuncInput.function.size() != thisFuncInput.derivative.size())
+ {
+ GMX_THROW(InconsistentInputError("Function and derivative vectors have different lengths"));
+ }
+
+ internal::throwUnlessDerivativeIsConsistentWithFunction(thisFuncInput.function, thisFuncInput.derivative, thisFuncInput.spacing, range_);
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ // Calculate the required table spacing h. The error we make with linear interpolation
+ // of the derivative will be described by the third-derivative correction term.
+ // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
+ // where f'/f''' is the first and third derivative of the function, respectively.
+
+ double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndSecondDerivative(thisFuncInput.derivative, thisFuncInput.spacing, range_);
+
+ minQuotient = std::min(minQuotient, thisMinQuotient);
+ }
+
+ double spacing = std::sqrt(12.0 * tolerance * minQuotient);
+
+ halfSpacing_ = 0.5 * spacing;
+ tableScale_ = 1.0 / spacing;
+
+ if (range_.second * tableScale_ > 1e6)
+ {
+ GMX_THROW(ToleranceError("Requested tolerance would require over a million points in table"));
+ }
+
+ // Loop over all tables again.
+ // Here we create the actual table for each function, and then
+ // combine them into a multiplexed table function.
+ std::size_t funcIndex = 0;
+
+ for (auto thisFuncInput : numericalInputList)
+ {
+ try
+ {
+ if (spacing < thisFuncInput.spacing)
+ {
+ GMX_THROW(ToleranceError("Input vector spacing cannot achieve tolerance requested"));
+ }
+
+ std::vector<real> tmpFuncTableData;
+ std::vector<real> tmpDerTableData;
+ std::vector<real> tmpDdfzTableData;
+
+ fillSingleQuadraticSplineTableData(thisFuncInput.function,
+ thisFuncInput.derivative,
+ thisFuncInput.spacing,
+ range,
+ spacing,
+ &tmpFuncTableData,
+ &tmpDerTableData);
+
+ fillDdfzTableData(tmpFuncTableData, tmpDerTableData, &tmpDdfzTableData);
+
+ internal::fillMultiplexedTableData(tmpDerTableData, &derivativeMultiTableData_,
+ 1, numFuncInTable_, funcIndex);
+
+ internal::fillMultiplexedTableData(tmpDdfzTableData, &ddfzMultiTableData_,
+ 4, numFuncInTable_, funcIndex);
+
+ funcIndex++;
+ }
+ catch (gmx::GromacsException &ex)
+ {
+ ex.prependContext("Error generating quadratic spline table for function '" + thisFuncInput.desc + "'");
+ throw;
+ }
+ }
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal
+ * \defgroup module_tables Classes for table interpolation
+ * \ingroup group_utilitymodules
+ *
+ * \brief Table interpolation from analytical or numerical input
+ *
+ * This module provides quadratic spline interpolation tables used
+ * both for the nonbonded kernels and listed interactions.
+ *
+ * \author Erik Lindahl <erik.lindahl@scilifelab.se>
+ */
+
+
+/*! \libinternal \file
+ * \brief
+ * Declares classes for quadratic spline table
+ *
+ * \inlibraryapi
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#ifndef GMX_TABLES_QUADRATICSPLINETABLE_H
+#define GMX_TABLES_QUADRATICSPLINETABLE_H
+
+#include <functional>
+#include <initializer_list>
+#include <vector>
+
+#include "gromacs/simd/simd.h"
+#include "gromacs/tables/tableinput.h"
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/classhelpers.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
+
+namespace gmx
+{
+
+/*! \libinternal \brief Quadratic spline interpolation table.
+ *
+ * This class interpolates a function specified either as an analytical
+ * expression or from user-provided table data.
+ *
+ * At initialization, you provide the reference function of vectors
+ * as a list of tuples that contain a brief name, the function, and
+ * derivative for each function to tabulate. To create a table with
+ * two functions this initializer list can for instance look like
+ *
+ * { {"LJ6", lj6Func, lj6Der}, {"LJ12", lj12Func, lj12Der} }
+ *
+ * The names are only used so exceptions during initialization can
+ * be traced to a specific table.
+ *
+ * When interpolating, there are methods to interpolate either 1, 2, or 3
+ * functions in one go. By default these interpolation routines will
+ * operate on tables with the same number of functions as specified in
+ * the interpolation method (debug builds check that this is consistent with
+ * the table). However, it is also possible to use optional template
+ * parameters that specify the total number of functions in a table, and
+ * what function index to interpolate. For instance, to interpolate the
+ * derivative of the second function (i.e., index 1) in a
+ * multi-function-table with three functions in total, you can write
+ *
+ * table.evaluateDerivative<3,1>(x,&der);
+ *
+ * Here too, debug builds will check that the template parameters are
+ * consistent with the table.
+ *
+ * The table data is internally adjusted to guarantee that the interpolated
+ * derivative is the true derivative of the interpolated potential, which is
+ * important to avoid systematic errors for the common case when the derivative
+ * is concave/convex in the entire interval.
+ * We do this by expressing the difference in the function value
+ * at a small offset h relative to a reference value in position 0 with a forward
+ * Taylor series expanded around 0, and then doing the opposite of expressing
+ * difference in the function at position 0 relative to a reference value in
+ * position h when using a backward Taylor expansion:
+ *
+ * \f{eqnarray*}{
+ * \Delta V & = & hV'(0) + \frac{1}{2} h^2 V''(0) + \frac{1}{6} h^3 V'''(0) + O(h^4) \\
+ * \Delta V & = & hV'(h) - \frac{1}{2} h^2 V''(h) + \frac{1}{6} h^3 V'''(h) + O(h^4)
+ * \f}
+ *
+ * Summing the equations leads to
+ *
+ * \f[
+ * 2 \Delta V = h(V'(0) + V'(h)) + \frac{1}{2} h^2 (V''(0)-V''(h)) + \frac{1}{6}h^3(V'''(0)+V'''(h)) + O(h^4)
+ * \f]
+ *
+ * To make the second term symmetric too, we can replace it with the average of
+ * the Taylor expansion at 0 and h (i.e., using the third derivative). This gives
+ *
+ * \f[
+ * 2 \Delta V = h(V'(0) + V'(h)) - \frac{1}{12} h^3 (V'''(0)+V'''(h)) + O(h^4)
+ * \f]
+ *
+ * Thus, if we replace the derivative in the internal quadratic table data with
+ *
+ * \f[
+ * V' - \frac{1}{12}h^2 V'''
+ * \f]
+ *
+ * we will cancel the h^3 term in the error. This will make the integral of the
+ * forces match the potential much better (The h^4 term actually disappears, so
+ * when summing over 1/h points the remaining error will be O(h^4).
+ *
+ * While it is possible to create tables only from function values
+ * (i.e., no derivatives), it is recommended to provide derivatives for higher
+ * accuracy and to avoid issues with numerical differentiation. Note that the
+ * table input should be smooth, i.e. it should not contain noise e.g. from an
+ * (iterative) Boltzmann inversion procedure - you have been warned.
+ *
+ * \note This class is responsible for fundamental interpolation of any function,
+ * which might or might not correspond to a potential. For this reason
+ * both input and output derivatives are proper function derivatives, and
+ * we do not swap any signs to get forces directly from the table.
+ *
+ * \note There will be a small additional accuracy loss from the internal
+ * operation where we calculate the epsilon offset from the nearest table
+ * point, since the integer part we subtract can get large in those cases.
+ * The absolute such error both in the function and derivative value will
+ * be roughly f''*x*GMX_REAL_EPS, where x is the argument and f'' the
+ * second derivative.
+ * While this is technically possible to solve with extended precision
+ * arithmetics, that would introduce extra instructions in some highly
+ * performance-sensitive code parts. For typical GROMACS interaction
+ * functions the derivatives will decay faster than the potential, which
+ * means it will never play any role. For other functions it will only
+ * cause a small increase in the relative error for arguments where the
+ * magnitude of the function or derivative is very small.
+ * Since we typically sum several results in GROMACS, this should never
+ * show up in any real cases, and for this reason we choose not to do
+ * the extended precision arithmetics.
+ *
+ * \note These routines are not suitable for table ranges starting far away
+ * from zero, since we allocate memory and calculate indices starting from
+ * range zero for efficiency reasons.
+ */
+class QuadraticSplineTable
+{
+ private:
+ /*! \brief Change that function value falls inside range when debugging
+ *
+ * \tparam T Lookup argument floating-point type, typically SimdReal or real.
+ * \param r Lookup argument to test
+ *
+ * \throws Debug builds will throw gmx::RangeError for values that are
+ * larger than the upper limit of the range, or smaller than 0.
+ * We allow the table to be called with arguments between 0 and
+ * the lower limit of the range, since this might in theory occur
+ * once-in-a-blue-moon with some algorithms.
+ */
+ template <typename T>
+ void
+ rangeCheck(T gmx_unused r) const
+ {
+#ifndef NDEBUG
+ // Check that all values fall in range when debugging
+ if (anyTrue( r < T(0.0) || T(range_.second) <= r ) )
+ {
+ GMX_THROW(RangeError("Interpolation input value falls outside table definition range"));
+ }
+#endif
+ }
+
+ public:
+
+ /*! \brief Default tolerance for tables is 10*GMX_FLOAT_EPS
+ *
+ * \note Even for double precision builds we set the tolerance to
+ * one order of magnitude above the single precision epsilon.
+ */
+ static const real defaultTolerance;
+
+ /*! \brief Initialize table data from function
+ *
+ * \param analyticalInputList Initializer list with one or more functions to tabulate,
+ * specified as pairs containing analytical
+ * functions and their derivatives. The function will also
+ * be called for values smaller than the lower limit of the
+ * range, but we avoid calling it for 0.0 if that value
+ * is not included in the range.
+ * \param range Range over which the function will be tabulated.
+ * Constructor will throw gmx::APIError for negative values.
+ * Due to the way the numerical derivative evaluation depends
+ * on machine precision internally, this range must be
+ * at least 0.001, or the constructor throws gmx::APIError.
+ * \param tolerance Requested accuracy of the table. This will be used to
+ * calculate the required internal spacing. If this cannot
+ * be achieved (for instance because the table would require
+ * too much memory) the constructor will throw gmx::ToleranceError.
+ *
+ * \note The functions are always defined in double precision to avoid
+ * losing accuracy when constructing tables.
+ *
+ * \note Since we fill the table for values below range.first, you can achieve
+ * a smaller table by using a smaller range where the tolerance has to be
+ * met, and accept that a few function calls below range.first do not
+ * quite reach the tolerance.
+ *
+ * \warning For efficiency reasons (since this code is used in some inner
+ * (kernels), we always allocate memory and calculate table indices
+ * for the complete interval [0,range.second], although the data will
+ * not be valid outside the definition range to avoid calling the
+ * function there. This means you should \a not use this class
+ * to tabulate functions for small ranges very far away from zero,
+ * since you would both waste a huge amount of memory and incur
+ * truncation errors when calculating the index.
+ *
+ * \throws gmx::ToleranceError if the requested tolerance cannot be achieved,
+ * and gmx::APIError for other incorrect input.
+ */
+ QuadraticSplineTable(std::initializer_list<AnalyticalSplineTableInput> analyticalInputList,
+ const std::pair<real, real> &range,
+ real tolerance = defaultTolerance);
+
+ /*! \brief Initialize table data from tabulated values and derivatives
+ *
+ * \param numericalInputList Initializer list with one or more functions to tabulate,
+ * specified as pairs containing containing vectors for the
+ * function values and their derivatives. Data points are
+ * separated by the spacing parameter, starting from 0.
+ * Values below the lower limit of the range will be used to
+ * attempt defining the table, but we avoid using index 0
+ * unless 0.0 is included in the range. Some extra points beyond
+ * range.second are required to re-interpolate values, so add
+ * some margin. The constructor will throw gmx::APIError if the
+ * input vectors are too short to cover the requested range
+ * (and they must always be at least five points).
+ * \param range Range over which the function will be tabulated.
+ * Constructor will throw gmx::APIError for negative values,
+ * or if the value/derivative vector does not cover the
+ * range.
+ * \param tolerance Requested accuracy of the table in the range. This will be
+ * used to calculate the required internal spacing and possibly
+ * re-interpolate. The constructor will throw
+ * gmx::ToleranceError if the input spacing is too coarse
+ * to achieve this accuracy.
+ *
+ * \note The input data vectors are always double precision to avoid
+ * losing accuracy when constructing tables.
+ *
+ * \note Since we fill the table for values below range.first, you can achieve
+ * a smaller table by using a smaller range where the tolerance has to be
+ * met, and accept that a few function calls below range.first do not
+ * quite reach the tolerance.
+ *
+ * \warning For efficiency reasons (since this code is used in some inner
+ * (kernels), we always allocate memory and calculate table indices
+ * for the complete interval [0,range.second], although the data will
+ * not be valid outside the definition range to avoid calling the
+ * function there. This means you should \a not use this class
+ * to tabulate functions for small ranges very far away from zero,
+ * since you would both waste a huge amount of memory and incur
+ * truncation errors when calculating the index.
+ */
+ QuadraticSplineTable(std::initializer_list<NumericalSplineTableInput> numericalInputList,
+ const std::pair<real, real> &range,
+ real tolerance = defaultTolerance);
+
+
+ /************************************************************
+ * Evaluation methods for single functions *
+ ************************************************************/
+
+ /*! \brief Evaluate both function and derivative, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue Function value
+ * \param[out] derivativeValue Function derivative
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue,
+ T * derivativeValue) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T t0;
+ T t1;
+ T t2;
+ T t3 gmx_unused;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4 * funcIndex, tabIndex, &t0, &t1, &t2, &t3);
+
+ t1 = t0 + eps * t1;
+ *functionValue = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue = t1;
+ }
+
+ /*! \brief Evaluate function value only, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue Function value
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue) const
+ {
+ T der gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(r, functionValue, &der);
+ }
+
+ /*! \brief Evaluate function derivative only, single table function
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 1
+ * \tparam funcIndex Index of function to evaluate in table, default is 0
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue Function derivative
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 1, int funcIndex = 0, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T t0;
+ T t1;
+ T t2 gmx_unused;
+
+ if (numFuncInTable == 1)
+ {
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex, tabIndex, &t0, &t1); // works for scalar T too
+ }
+ else
+ {
+ // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
+ // only loads a single value from memory to implement it better (will be written)
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ }
+
+ // (1-eps)*t0 + eps*t1
+ *derivativeValue = fma(t1-t0, eps, t0);
+ }
+
+ /************************************************************
+ * Evaluation methods for two functions *
+ ************************************************************/
+
+ /*! \brief Evaluate both function and derivative, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue1 Interpolated value for first function
+ * \param[out] derivativeValue1 Interpolated derivative for first function
+ * \param[out] functionValue2 Interpolated value for second function
+ * \param[out] derivativeValue2 Interpolated derivative for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue1,
+ T * derivativeValue1,
+ T * functionValue2,
+ T * derivativeValue2) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T t0;
+ T t1;
+ T t2;
+ T t3 gmx_unused;
+
+ // Load Derivative, Delta, Function, and Zero values for each table point.
+ // The 4 refers to these four values - not any SIMD width.
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4*funcIndex0, tabIndex, &t0, &t1, &t2, &t3);
+ t1 = t0 + eps * t1;
+ *functionValue1 = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue1 = t1;
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4*funcIndex1, tabIndex, &t0, &t1, &t2, &t3);
+ t1 = t0 + eps * t1;
+ *functionValue2 = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue2 = t1;
+ }
+
+ /*! \brief Evaluate function value only, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue1 Interpolated value for first function
+ * \param[out] functionValue2 Interpolated value for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue1,
+ T * functionValue2) const
+ {
+ T der1 gmx_unused;
+ T der2 gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1>(r, functionValue1, &der1, functionValue2, &der2);
+ }
+
+ /*! \brief Evaluate function derivative only, two table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 2
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue1 Interpolated derivative for first function
+ * \param[out] derivativeValue2 Interpolated derivative for second function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 2, int funcIndex0 = 0, int funcIndex1 = 1, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue1,
+ T * derivativeValue2) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable, "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+
+ if (numFuncInTable == 2 && funcIndex0 == 0 && funcIndex1 == 1)
+ {
+ T t0A, t0B, t1A, t1B;
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data(), tabIndex, &t0A, &t0B); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + 2, tabIndex, &t1A, &t1B); // works for scalar T too
+ *derivativeValue1 = fma(t1A-t0A, eps, t0A);
+ *derivativeValue2 = fma(t1B-t0B, eps, t0B);
+ }
+ else
+ {
+ T t0, t1, t2;
+ // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
+ // only loads a single value from memory to implement it better (will be written)
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ *derivativeValue1 = fma(t1-t0, eps, t0);
+
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ *derivativeValue2 = fma(t1-t0, eps, t0);
+ }
+ }
+
+ /************************************************************
+ * Evaluation methods for three functions *
+ ************************************************************/
+
+
+ /*! \brief Evaluate both function and derivative, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function and derivative
+ * \param[out] functionValue1 Interpolated value for first function
+ * \param[out] derivativeValue1 Interpolated derivative for first function
+ * \param[out] functionValue2 Interpolated value for second function
+ * \param[out] derivativeValue2 Interpolated derivative for second function
+ * \param[out] functionValue3 Interpolated value for third function
+ * \param[out] derivativeValue3 Interpolated derivative for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateFunctionAndDerivative(T r,
+ T * functionValue1,
+ T * derivativeValue1,
+ T * functionValue2,
+ T * derivativeValue2,
+ T * functionValue3,
+ T * derivativeValue3) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
+ "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+ T t0;
+ T t1;
+ T t2;
+ T t3 gmx_unused;
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4*funcIndex0, tabIndex, &t0, &t1, &t2, &t3);
+ t1 = t0 + eps * t1;
+ *functionValue1 = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue1 = t1;
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4*funcIndex1, tabIndex, &t0, &t1, &t2, &t3);
+ t1 = t0 + eps * t1;
+ *functionValue2 = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue2 = t1;
+
+ gatherLoadBySimdIntTranspose<4*numFuncInTable>(ddfzMultiTableData_.data() + 4*funcIndex2, tabIndex, &t0, &t1, &t2, &t3);
+ t1 = t0 + eps * t1;
+ *functionValue3 = fma( eps * T(halfSpacing_), t0 + t1, t2);
+ *derivativeValue3 = t1;
+ }
+
+ /*! \brief Evaluate function value only, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function value
+ * \param[out] functionValue1 Interpolated value for first function
+ * \param[out] functionValue2 Interpolated value for second function
+ * \param[out] functionValue3 Interpolated value for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateFunction(T r,
+ T * functionValue1,
+ T * functionValue2,
+ T * functionValue3) const
+ {
+ T der1 gmx_unused;
+ T der2 gmx_unused;
+ T der3 gmx_unused;
+
+ evaluateFunctionAndDerivative<numFuncInTable, funcIndex0, funcIndex1, funcIndex2>(r, functionValue1, &der1, functionValue2, &der2, functionValue3, &der3);
+ }
+
+ /*! \brief Evaluate function derivative only, three table functions
+ *
+ * This is a templated method where the template can be either real or SimdReal.
+ *
+ * \tparam numFuncInTable Number of separate functions in table, default is 3
+ * \tparam funcIndex0 Index of 1st function to evaluate in table, default is 0
+ * \tparam funcIndex1 Index of 2nd function to evaluate in table, default is 1
+ * \tparam funcIndex2 Index of 3rd function to evaluate in table, default is 2
+ * \tparam T Type (SimdReal or real) of lookup and result
+ * \param r Points for which to evaluate function derivative
+ * \param[out] derivativeValue1 Interpolated derivative for first function
+ * \param[out] derivativeValue2 Interpolated derivative for second function
+ * \param[out] derivativeValue3 Interpolated derivative for third function
+ *
+ * For debug builds we assert that the input values fall in the range
+ * specified when constructing the table.
+ */
+ template <int numFuncInTable = 3, int funcIndex0 = 0, int funcIndex1 = 1, int funcIndex2 = 2, typename T>
+ void
+ evaluateDerivative(T r,
+ T * derivativeValue1,
+ T * derivativeValue2,
+ T * derivativeValue3) const
+ {
+ rangeCheck(r);
+ GMX_ASSERT(numFuncInTable == numFuncInTable_, "Evaluation method not matching number of functions in table");
+ GMX_ASSERT(funcIndex0 < numFuncInTable && funcIndex1 < numFuncInTable && funcIndex2 < numFuncInTable,
+ "Function index not in range of the number of tables");
+
+ T rTable = r * T(tableScale_);
+ auto tabIndex = cvttR2I(rTable); // type is either std::int32_t or SimdInt32
+ T eps = rTable - trunc(rTable);
+
+ if (numFuncInTable == 3 && funcIndex0 == 0 && funcIndex1 == 1 && funcIndex2 == 2)
+ {
+ T t0A, t0B, t0C, t1A, t1B, t1C;
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data(), tabIndex, &t0A, &t0B);
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + 2, tabIndex, &t0C, &t1A);
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + 4, tabIndex, &t1B, &t1C);
+ *derivativeValue1 = fma(t1A-t0A, eps, t0A);
+ *derivativeValue2 = fma(t1B-t0B, eps, t0B);
+ *derivativeValue3 = fma(t1C-t0C, eps, t0C);
+ }
+ else
+ {
+ T t0, t1, t2;
+ // This is not ideal, but we need a version of gatherLoadUBySimdIntTranspose that
+ // only loads a single value from memory to implement it better (will be written)
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex0, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ *derivativeValue1 = fma(t1-t0, eps, t0);
+
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex1, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ *derivativeValue2 = fma(t1-t0, eps, t0);
+
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex2, tabIndex, &t0, &t2); // works for scalar T too
+ gatherLoadUBySimdIntTranspose<numFuncInTable>(derivativeMultiTableData_.data() + funcIndex2, tabIndex + T(1), &t1, &t2); // works for scalar T too
+ *derivativeValue3 = fma(t1-t0, eps, t0);
+ }
+ }
+
+ /*! \brief Return the table spacing (distance between points)
+ *
+ * You should never have to use this for normal code, but due to the
+ * way tables are constructed internally we need this in the unit tests
+ * to check relative tolerances over each interval.
+ *
+ * \return table spacing.
+ */
+ real
+ tableSpacing() const { return 1.0 / tableScale_; }
+
+ private:
+
+ std::size_t numFuncInTable_; //!< Number of separate tabluated functions
+ std::pair<real, real> range_; //!< Range for which table evaluation is allowed
+ real tableScale_; //!< Table scale (inverse of spacing between points)
+ real halfSpacing_; //!< 0.5*spacing (used for DDFZ table data)
+
+ //!< Derivative values only, with the third-derivative subtraction described in the class documentation.
+ std::vector<real> derivativeMultiTableData_;
+
+ /*! \brief Combined derivative, difference to next derivative, value, and zero.
+ *
+ * For table point i, this vector contains the four values:
+ * - derivative[i]
+ * - (derivative[i+1]-derivative[i])
+ * - value[i]
+ * - 0.0
+ *
+ * For the derivative terms we have subtracted the third-derivative term described
+ * in the main class documentation.
+ *
+ * This is typically more efficient than the individual tables, in particular
+ * when using SIMD. The result should be identical outside the class, so this
+ * is merely an internal implementation optimization. However, to allow
+ * aligned SIMD loads we need to use an aligned allocator for this container.
+ * We occasionally abbreviate this data as DDFZ.
+ */
+ std::vector<real, AlignedAllocator<real> > ddfzMultiTableData_;
+
+ // There should never be any reason to copy the table since it is read-only
+ GMX_DISALLOW_COPY_AND_ASSIGN(QuadraticSplineTable);
+};
+
+
+} // namespace gmx
+
+#endif // GMX_TABLES_QUADRATICSPLINETABLE_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements internal utility functions for spline tables
+ *
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#include "gmxpre.h"
+
+#include "splineutil.h"
+
+#include <cmath>
+
+#include <algorithm>
+#include <functional>
+#include <utility>
+#include <vector>
+
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/real.h"
+#include "gromacs/utility/stringutil.h"
+
+namespace gmx
+{
+
+namespace internal
+{
+
+void
+throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)> &function,
+ const std::function<double(double)> &derivative,
+ const std::pair<real, real> &range)
+{
+ // Since the numerical derivative will evaluate extra points
+ // we shrink the interval slightly to avoid calling the function with values
+ // outside the range specified.
+ double h = std::cbrt(GMX_DOUBLE_EPS); // ideal spacing
+ std::pair<double, double> newRange(range.first + h, range.second - h);
+ const int points = 1000; // arbitrary
+ double dx = (newRange.second - newRange.first) / points;
+ bool isConsistent = true;
+ double minFail = newRange.second;
+ double maxFail = newRange.first;
+
+ for (double x = newRange.first; x <= newRange.second; x += dx)
+ {
+ double analyticalDerivative = derivative(x);
+ double numericalDerivative = (function(x+h)-function(x-h))/(2*h);
+ double thirdDerivative = (derivative(x+h)-2.0*derivative(x)+derivative(x-h))/(h*h);
+
+ // We make two types of errors in numerical calculation of the derivative:
+ // - The truncation error: eps * |f| / h
+ // - The rounding error: h * h * |f'''| / 6.0
+ double expectedErr = GMX_DOUBLE_EPS*std::abs(function(x))/h + h*h*std::abs(thirdDerivative)/6.0;
+
+ // To avoid triggering false errors because of compiler optimization or numerical issues
+ // in the function evalulation we allow an extra factor of 10 in the expected error
+ if (std::abs(analyticalDerivative-numericalDerivative) > 10.0*expectedErr)
+ {
+ isConsistent = false;
+ minFail = std::min(minFail, x);
+ maxFail = std::max(maxFail, x);
+ }
+ }
+
+ if (!isConsistent)
+ {
+ GMX_THROW(InconsistentInputError(formatString("Derivative inconsistent with analytical function in range [%d,%d]", minFail, maxFail)));
+ }
+}
+
+
+void
+throwUnlessDerivativeIsConsistentWithFunction(ConstArrayRef<double> function,
+ ConstArrayRef<double> derivative,
+ double inputSpacing,
+ const std::pair<real, real> &range)
+{
+ std::size_t firstIndex = range.first / inputSpacing;
+ std::size_t lastIndex = range.second / inputSpacing;
+ bool isConsistent = true;
+ std::size_t minFail = lastIndex;
+ std::size_t maxFail = firstIndex;
+
+ // The derivative will access one extra point before/after each point, so reduce interval
+ for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
+ {
+ double inputDerivative = derivative[i];
+ double numericalDerivative = (function[i+1] - function[i-1]) / (2.0 * inputSpacing);
+ double thirdDerivative = (derivative[i+1] - 2.0*derivative[i] + derivative[i-1])/(inputSpacing * inputSpacing);
+
+ // We make two types of errors in numerical calculation of the derivative:
+ // - The truncation error: eps * |f| / h
+ // - The rounding error: h * h * |f'''| / 6.0
+ double expectedErr = GMX_DOUBLE_EPS*std::abs(function[i])/inputSpacing + inputSpacing*inputSpacing*std::abs(thirdDerivative)/6.0;
+
+ // To avoid triggering false errors because of compiler optimization or numerical issues
+ // in the function evalulation we allow an extra factor of 10 in the expected error
+ if (std::abs(inputDerivative-numericalDerivative) > 10.0*expectedErr)
+ {
+ isConsistent = false;
+ minFail = std::min(minFail, i);
+ maxFail = std::max(maxFail, i);
+ }
+ }
+ if (!isConsistent)
+ {
+ GMX_THROW(InconsistentInputError(formatString("Derivative inconsistent with numerical vector for elements %d-%d", minFail+1, maxFail+1)));
+ }
+}
+
+
+/*! \brief Update minQuotient if the ratio of this function value and its second derivative is smaller
+ *
+ * This is a utility function used in the functions to find the smallest quotient
+ * in a range.
+ *
+ * \param[in] previousPoint Value of function at x-h.
+ * \param[in] thisPoint Value of function at x.
+ * \param[in] nextPoint Value of function at x+h.
+ * \param[in] spacing Value of h.
+ * \param[inout] minQuotient Current minimum of such quotients, updated if this quotient is smaller.
+ */
+static void
+updateMinQuotientOfFunctionAndSecondDerivative(double previousPoint,
+ double thisPoint,
+ double nextPoint,
+ double spacing,
+ double * minQuotient)
+{
+ double value = std::abs( thisPoint );
+ double secondDerivative = std::abs( (previousPoint - 2.0 * thisPoint + nextPoint) / (spacing * spacing ) );
+
+ // Make sure we do not divide by zero. This limit is arbitrary,
+ // but it doesnt matter since this point will have a very large value,
+ // and the whole routine is searching for the smallest value.
+ secondDerivative = std::max(secondDerivative, static_cast<double>(std::sqrt(GMX_REAL_MIN)));
+
+ *minQuotient = std::min(*minQuotient, value / secondDerivative);
+}
+
+
+real
+findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)> &f,
+ const std::pair<real, real> &range)
+{
+ // Since the numerical second derivative will evaluate extra points
+ // we shrink the interval slightly to avoid calling the function with values
+ // outside the range specified.
+ double h = std::pow( GMX_DOUBLE_EPS, 0.25 );
+ std::pair<double, double> newRange(range.first + h, range.second - h);
+ const int points = 1000; // arbitrary
+ double dx = (newRange.second - newRange.first) / points;
+ double minQuotient = GMX_REAL_MAX;
+
+ for (double x = newRange.first; x <= newRange.second; x += dx)
+ {
+ updateMinQuotientOfFunctionAndSecondDerivative(f(x-h), f(x), f(x+h), h, &minQuotient);
+ }
+ return static_cast<real>(minQuotient);
+}
+
+
+
+
+
+
+real
+findSmallestQuotientOfFunctionAndSecondDerivative(ConstArrayRef<double> function,
+ double inputSpacing,
+ const std::pair<real, real> &range)
+{
+
+ std::size_t firstIndex = range.first / inputSpacing;
+ std::size_t lastIndex = range.second / inputSpacing;
+ double minQuotient = GMX_REAL_MAX;
+
+ for (std::size_t i = firstIndex + 1; (i + 1) < lastIndex; i++)
+ {
+ updateMinQuotientOfFunctionAndSecondDerivative(function[i-1], function[i], function[i+1], inputSpacing, &minQuotient);
+ }
+ return static_cast<real>(minQuotient);
+}
+
+
+
+/*! \brief Update minQuotient if the ratio of this function value and its third derivative is smaller
+ *
+ * This is a utility function used in the functions to find the smallest quotient
+ * in a range.
+ *
+ * \param[in] previousPreviousPoint Value of function at x-2h.
+ * \param[in] previousPoint Value of function at x-h.
+ * \param[in] thisPoint Value of function at x.
+ * \param[in] nextPoint Value of function at x+h.
+ * \param[in] nextNextPoint Value of function at x+2h.
+ * \param[in] spacing Value of h.
+ * \param[inout] minQuotient Current minimum of such quotients, updated if this quotient is smaller.
+ */
+static void
+updateMinQuotientOfFunctionAndThirdDerivative(double previousPreviousPoint,
+ double previousPoint,
+ double thisPoint,
+ double nextPoint,
+ double nextNextPoint,
+ double spacing,
+ double * minQuotient)
+{
+ double value = std::abs( thisPoint );
+ double thirdDerivative = std::abs((nextNextPoint - 2 * nextPoint + 2 * previousPoint - previousPreviousPoint) / (2 * spacing * spacing * spacing));
+
+ // Make sure we do not divide by zero. This limit is arbitrary,
+ // but it doesnt matter since this point will have a very large value,
+ // and the whole routine is searching for the smallest value.
+ thirdDerivative = std::max(thirdDerivative, static_cast<double>(std::sqrt(GMX_REAL_MIN)));
+
+ *minQuotient = std::min(*minQuotient, value / thirdDerivative);
+}
+
+
+real
+findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)> &f,
+ const std::pair<real, real> &range)
+{
+ // Since the numerical second derivative will evaluate extra points
+ // we shrink the interval slightly to avoid calling the function with values
+ // outside the range specified.
+ double h = std::pow( GMX_DOUBLE_EPS, 0.2 ); // optimal spacing for 3rd derivative
+ std::pair<double, double> newRange(range.first + 2*h, range.second - 2*h);
+ const int points = 1000; // arbitrary
+ double dx = (newRange.second - newRange.first) / points;
+ double minQuotient = GMX_REAL_MAX;
+
+ for (double x = newRange.first; x <= newRange.second; x += dx)
+ {
+ updateMinQuotientOfFunctionAndThirdDerivative(f(x-2*h), f(x-h), f(x), f(x+h), f(x+2*h), h, &minQuotient);
+ }
+ return static_cast<real>(minQuotient);
+}
+
+
+real
+findSmallestQuotientOfFunctionAndThirdDerivative(ConstArrayRef<double> function,
+ double inputSpacing,
+ const std::pair<real, real> &range)
+{
+
+ std::size_t firstIndex = range.first / inputSpacing;
+ std::size_t lastIndex = range.second / inputSpacing;
+ double minQuotient = GMX_REAL_MAX;
+
+ for (std::size_t i = firstIndex + 2; (i + 2) < lastIndex; i++)
+ {
+ updateMinQuotientOfFunctionAndThirdDerivative(function[i-2], function[i-1], function[i], function[i+1], function[i+2], inputSpacing, &minQuotient);
+ }
+ return static_cast<real>(minQuotient);
+}
+
+
+
+std::vector<double>
+vectorSecondDerivative(ConstArrayRef<double> f, double spacing)
+{
+ if (f.size() < 5)
+ {
+ GMX_THROW(APIError("Too few points in vector for 5-point differentiation"));
+ }
+
+ std::vector<double> d(f.size());
+ std::size_t i;
+
+ // 5-point formula evaluated for points 0,1
+ i = 0;
+ d[i] = (11 * f[i+4] - 56 * f[i+3] + 114 * f[i+2] - 104 * f[i+1] + 35 * f[i]) / ( 12 * spacing * spacing);
+ i = 1;
+ d[i] = (-f[i+3] + 4 * f[i+2] + 6 * f[i+1] - 20 * f[i] + 11 * f[i-1]) / ( 12 * spacing * spacing);
+
+ for (std::size_t i = 2; i < d.size() - 2; i++)
+ {
+ // 5-point formula evaluated for central point (2)
+ d[i] = (-f[i+2] + 16 * f[i+1] - 30 * f[i] + 16 * f[i-1] - f[i-2]) / (12 * spacing * spacing);
+ }
+
+ // 5-point formula evaluated for points 3,4
+ i = d.size() - 2;
+ d[i] = (11 * f[i+1] - 20 * f[i] + 6 * f[i-1] + 4 * f[i-2] - f[i-3]) / ( 12 * spacing * spacing);
+ i = d.size() - 1;
+ d[i] = (35 * f[i] - 104 * f[i-1] + 114 * f[i-2] - 56 * f[i-3] + 11 * f[i-4]) / ( 12 * spacing * spacing);
+
+ return d;
+}
+
+
+
+} // namespace internal
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Declares internal utility functions for spline tables
+ *
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#ifndef GMX_TABLES_SPLINEUTIL_H
+#define GMX_TABLES_SPLINEUTIL_H
+
+#include <functional>
+#include <utility>
+#include <vector>
+
+#include "gromacs/utility/alignedallocator.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
+
+namespace gmx
+{
+
+namespace internal
+{
+
+/*! \brief Ensure analytical derivative is the derivative of analytical function.
+ *
+ * This routine evaluates the numerical derivative of the function for
+ * a few (1000) points in the interval and checks that the relative difference
+ * between numerical and analytical derivative is within the expected error
+ * for the numerical derivative approximation we use.
+ *
+ * The main point of this routine is to make sure the user has not made a
+ * mistake or sign error when defining the functions.
+ *
+ * \param function Analytical function to differentiate
+ * \param derivative Analytical derivative to compare with
+ * \param range Range to test
+ *
+ * \throws If the provided derivative does not seem to match the function.
+ *
+ * \note The function/derivative are always double-valued to avoid accuracy loss.
+ */
+void
+throwUnlessDerivativeIsConsistentWithFunction(const std::function<double(double)> &function,
+ const std::function<double(double)> &derivative,
+ const std::pair<real, real> &range);
+
+/*! \brief Ensure vector of derivative values is the derivative of function vector.
+ *
+ * This routine differentiates a vector of numerical values and checks
+ * that the relative difference to a provided vector of numerical derivatives
+ * is smaller than the expected error from the numerical differentiation.
+ *
+ * The main point of this routine is to make sure the user has not made a
+ * mistake or sign error when defining the functions.
+ *
+ * To avoid problems if the vectors change from zero to finite values at the
+ * start/end of the interval, we only check inside the range requested.
+ *
+ * \param function Numerical function value vector to differentiate
+ * \param derivative Numerical derivative vector to compare with
+ * \param inputSpacing Distance between input points
+ * \param range Range to test
+ *
+ * \throws If the provided derivative does not seem to match the function.
+ *
+ * \note The function/derivative vectors and spacing are always double-valued
+ * to avoid accuracy loss.
+ */
+void
+throwUnlessDerivativeIsConsistentWithFunction(ConstArrayRef<double> function,
+ ConstArrayRef<double> derivative,
+ double inputSpacing,
+ const std::pair<real, real> &range);
+
+
+
+
+/*! \brief Find smallest quotient between analytical function and its 2nd derivative
+ *
+ * Used to calculate spacing for quadratic spline tables. This function divides the
+ * function value by the second derivative (or a very small number when that is zero),
+ * and returns the smallest such quotient found in the range.
+ *
+ * Our quadratic tables corresponds to linear interpolation of the derivative,
+ * which means the derivative will typically have larger error than the value
+ * when interpolating. The spacing required to reach a particular relative
+ * tolerance in the derivative depends on the quotient between the first
+ * derivative and the third derivative of the function itself.
+ *
+ * You should call this routine with the analytical derivative as the "function"
+ * parameter, and the quotient between "function and second derivative" will
+ * then correspond to the quotient bewteen the derivative and the third derivative
+ * of the actual function we want to tabulate.
+ *
+ * Since all functions that can be tabulated efficiently are reasonably smooth,
+ * we simply check 1,000 points in the interval rather than bother about
+ * implementing any complicated global optimization scheme.
+ *
+ * \param f Analytical function
+ * \param range Interval
+ *
+ * \return Smallest quotient found in range.
+ *
+ * \note The function is always double-valued to avoid accuracy loss.
+ */
+real
+findSmallestQuotientOfFunctionAndSecondDerivative(const std::function<double(double)> &f,
+ const std::pair<real, real> &range);
+
+
+
+
+
+/*! \brief Find smallest quotient between vector of values and its 2nd derivative
+ *
+ * Used to calculate spacing for quadratic spline tables. This function divides the
+ * function value by the second derivative (or a very small number when that is zero),
+ * and returns the smallest such quotient found in the range.
+ *
+ * Our quadratic tables corresponds to linear interpolation of the derivative,
+ * which means the derivative will typically have larger error than the value
+ * when interpolating. The spacing required to reach a particular relative
+ * tolerance in the derivative depends on the quotient between the first
+ * derivative and the third derivative of the function itself.
+ *
+ * You should call this routine with the analytical derivative as the "function"
+ * parameter, and the quotient between "function and second derivative" will
+ * then correspond to the quotient bewteen the derivative and the third derivative
+ * of the actual function we want to tabulate.
+ *
+ * \param function Vector with function values
+ * \param inputSpacing Spacing between function values
+ * \param range Interval to check
+ *
+ * \return Smallest quotient found in range.
+ *
+ * \note The function vector and input spacing are always double-valued to
+ * avoid accuracy loss.
+ */
+real
+findSmallestQuotientOfFunctionAndSecondDerivative(ConstArrayRef<double> function,
+ double inputSpacing,
+ const std::pair<real, real> &range);
+
+
+
+
+
+/*! \brief Find smallest quotient between analytical function and its 3rd derivative
+ *
+ * Used to calculate table spacing. This function divides the function value
+ * by the second derivative (or a very small number when that is zero), and
+ * returns the smallest such quotient found in the range.
+ *
+ * Our quadratic tables corresponds to linear interpolation of the derivative,
+ * which means the derivative will typically have larger error than the value
+ * when interpolating. The spacing required to reach a particular relative
+ * tolerance in the derivative depends on the quotient between the first
+ * derivative and the third derivative of the function itself.
+ *
+ * You should call this routine with the analytical derivative as the "function"
+ * parameter, and the quotient between "function and second derivative" will
+ * then correspond to the quotient bewteen the derivative and the third derivative
+ * of the actual function we want to tabulate.
+ *
+ * Since all functions that can be tabulated efficiently are reasonably smooth,
+ * we simply check 1,000 points in the interval rather than bother about
+ * implementing any complicated global optimization scheme.
+ *
+ * \param f Analytical function
+ * \param range Interval
+ *
+ * \return Smallest quotient found in range.
+ *
+ * \note The function is always double-valued to avoid accuracy loss.
+ */
+real
+findSmallestQuotientOfFunctionAndThirdDerivative(const std::function<double(double)> &f,
+ const std::pair<real, real> &range);
+
+
+
+
+/*! \brief Find smallest quotient between function and 2nd derivative (vectors)
+ *
+ * Used to calculate table spacing. This function divides the function value
+ * by the second derivative (or a very small number when that is zero), and
+ * returns the smallest such quotient found in the range.
+ *
+ * Our quadratic tables corresponds to linear interpolation of the derivative,
+ * which means the derivative will typically have larger error than the value
+ * when interpolating. The spacing required to reach a particular relative
+ * tolerance in the derivative depends on the quotient between the first
+ * derivative and the third derivative of the function itself.
+ *
+ * You should call this routine with the analytical derivative as the "function"
+ * parameter, and the quotient between "function and second derivative" will
+ * then correspond to the quotient bewteen the derivative and the third derivative
+ * of the actual function we want to tabulate.
+ *
+ * \param function Vector with function values
+ * \param inputSpacing Spacing between function values
+ * \param range Interval to check
+ *
+ * \return Smallest quotient found in range.
+ *
+ * \note The function vector and input spacing are always double-valued to
+ * avoid accuracy loss.
+ */
+real
+findSmallestQuotientOfFunctionAndThirdDerivative(ConstArrayRef<double> function,
+ double inputSpacing,
+ const std::pair<real, real> &range);
+
+
+/*! \brief Calculate second derivative of vector and return vector of same length
+ *
+ * 5-point approximations are used, with endpoints using non-center interpolation.
+ *
+ * \param f Vector (function) for which to calculate second derivative
+ * \param spacing Spacing of input data.
+ *
+ * \throws If the input vector has fewer than five data points.
+ *
+ * \note This function always uses double precision arguments since it is meant
+ * to be used on raw user input data for tables, where we want to avoid
+ * accuracy loss (since differentiation can be numerically fragile).
+ */
+std::vector<double>
+vectorSecondDerivative(ConstArrayRef<double> f,
+ double spacing);
+
+
+/*! \brief Copy (temporary) table data into aligned multiplexed vector
+ *
+ * This routine takes the temporary data generated for a single table
+ * and writes multiplexed output into a multiple-table-data vector.
+ * If the output vector is empty we will resize it to fit the data, and
+ * otherwise we assert the size is correct to add out input data.
+ *
+ * \tparam T Type of container for input data
+ * \tparam U Type of container for output data
+ *
+ * \param[in] inputData Input data for single table
+ * \param[inout] multiplexedOutputData Multiplexed output vector, many tables.
+ * \param[in] valuesPerTablePoint Number of real values for each table point,
+ * for instance 4 in DDFZ tables.
+ * \param[in] numTables Number of tables mixed into multiplexed output
+ * \param[in] thisTableIndex Index of this table in multiplexed output
+ *
+ * \note The output container type can be different from the input since the latter
+ * sometimes uses an aligned allocator so the data can be loaded efficiently
+ * in the GROMACS nonbonded kernels.
+ */
+template<class T, class U>
+void
+fillMultiplexedTableData(const T inputData,
+ U * multiplexedOutputData,
+ std::size_t valuesPerTablePoint,
+ std::size_t numTables,
+ std::size_t thisTableIndex)
+{
+ if (multiplexedOutputData->size() == 0)
+ {
+ multiplexedOutputData->resize( inputData.size() * numTables );
+ }
+ else
+ {
+ GMX_ASSERT(inputData.size() * numTables == multiplexedOutputData->size(),
+ "Size mismatch when multiplexing table data");
+ }
+
+ GMX_ASSERT(inputData.size() % valuesPerTablePoint == 0,
+ "Input table size must be a multiple of valuesPerTablePoint");
+
+ std::size_t points = inputData.size() / valuesPerTablePoint;
+
+ for (std::size_t i = 0; i < points; i++)
+ {
+ std::size_t inputOffset = valuesPerTablePoint * i;
+ std::size_t outputOffset = valuesPerTablePoint * ( numTables * i + thisTableIndex );
+
+ for (std::size_t j = 0; j < valuesPerTablePoint; j++)
+ {
+ (*multiplexedOutputData)[outputOffset + j] = inputData[inputOffset + j];
+ }
+ }
+}
+
+
+} // namespace internal
+
+} // namespace gmx
+
+#endif // GMX_TABLES_SPLINEUTIL_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ * \brief
+ * Declares structures for analytical or numerical input data to construct tables
+ *
+ * \inlibraryapi
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+
+#ifndef GMX_TABLES_TABLEINPUT_H
+#define GMX_TABLES_TABLEINPUT_H
+
+#include <functional>
+#include <vector>
+
+#include "gromacs/utility/arrayref.h"
+
+namespace gmx
+{
+
+/*! \libinternal \brief Specification for analytical table function (name, function, derivative)
+ */
+struct
+AnalyticalSplineTableInput
+{
+ const std::string &desc; //!< \libinternal Brief description of function
+ std::function<double(double)> function; //!< \libinternal Analytical form of function
+ std::function<double(double)> derivative; //!< \libinternal Analytical derivative
+};
+
+/*! \libinternal \brief Specification for vector table function (name, function, derivative, spacing)
+ */
+struct
+NumericalSplineTableInput
+{
+ const std::string &desc; //!< \libinternal Brief description of function
+ ConstArrayRef<double> function; //!< \libinternal Vector with function values
+ ConstArrayRef<double> derivative; //!< \libinternal Vector with derivative values
+ double spacing; //!< \libinternal Distance between data points
+};
+
+
+} // namespace gmx
+
+
+#endif // GMX_TABLES_TABLEINPUT_H
--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2016, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+gmx_add_unit_test(TableUnitTests table-test
+ splinetable.cpp
+ )
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for simple math functions.eval
+ *
+ * \author Erik Lindahl <erik.lindahl@gmail.com>
+ * \ingroup module_tables
+ */
+#include "gmxpre.h"
+
+#include <cmath>
+
+#include <algorithm>
+#include <functional>
+#include <utility>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/utilities.h"
+#include "gromacs/options/basicoptions.h"
+#include "gromacs/options/ioptionscontainer.h"
+#include "gromacs/simd/simd.h"
+#include "gromacs/tables/cubicsplinetable.h"
+#include "gromacs/tables/quadraticsplinetable.h"
+
+#include "testutils/testasserts.h"
+#include "testutils/testoptions.h"
+
+
+namespace gmx
+{
+
+namespace test
+{
+
+namespace
+{
+
+class SplineTableTestBase : public ::testing::Test
+{
+ public:
+ static int s_testPoints_; //!< Number of points to use. Public so we can set it as option
+};
+
+int
+SplineTableTestBase::s_testPoints_ = 100;
+
+/*! \cond */
+/*! \brief Command-line option to adjust the number of points used to test SIMD math functions. */
+GMX_TEST_OPTIONS(SplineTableTestOptions, options)
+{
+ options->addOption(::gmx::IntegerOption("npoints")
+ .store(&SplineTableTestBase::s_testPoints_)
+ .description("Number of points to test for spline table functions"));
+}
+/*! \endcond */
+
+
+
+
+/*! \brief Test fixture for table comparision with analytical/numerical functions */
+template <typename T>
+class SplineTableTest : public SplineTableTestBase
+{
+ public:
+ SplineTableTest() : tolerance_(T::defaultTolerance) {}
+
+ /*! \brief Set a new tolerance to be used in table function comparison
+ *
+ * \param tol New tolerance to use
+ */
+ void
+ setTolerance(real tol) { tolerance_ = tol; }
+
+ //! \cond internal
+ /*! \internal \brief
+ * Assertion predicate formatter for comparing table with function/derivative
+ */
+ template<int numFuncInTable = 1, int funcIndex = 0>
+ void
+ testSplineTableAgainstFunctions(const std::string &desc,
+ const std::function<double(double)> &refFunc,
+ const std::function<double(double)> &refDer,
+ const T &table,
+ const std::pair<real, real> &testRange);
+ //! \endcond
+
+ private:
+ real tolerance_; //!< Tolerance to use
+};
+
+template <class T>
+template<int numFuncInTable, int funcIndex>
+void
+SplineTableTest<T>::testSplineTableAgainstFunctions(const std::string &desc,
+ const std::function<double(double)> &refFunc,
+ const std::function<double(double)> &refDer,
+ const T &table,
+ const std::pair<real, real> &testRange)
+{
+ real dx = (testRange.second - testRange.first) / s_testPoints_;
+
+ FloatingPointTolerance funcTolerance(relativeToleranceAsFloatingPoint(0.0, tolerance_));
+
+ for (real x = testRange.first; x < testRange.second; x += dx)
+ {
+ real h = std::sqrt(GMX_REAL_EPS);
+ real secondDerivative = (refDer(x+h)-refDer(x))/h;
+
+ real testFuncValue;
+ real testDerValue;
+
+ table.template evaluateFunctionAndDerivative<numFuncInTable, funcIndex>(x, &testFuncValue, &testDerValue);
+
+ // Check that we get the same values from function/derivative-only methods
+ real tmpFunc, tmpDer;
+
+ table.template evaluateFunction<numFuncInTable, funcIndex>(x, &tmpFunc);
+
+ table.template evaluateDerivative<numFuncInTable, funcIndex>(x, &tmpDer);
+
+ if (testFuncValue != tmpFunc)
+ {
+ ADD_FAILURE()
+ << "Interpolation inconsistency for table " << desc << std::endl
+ << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
+ << "First failure at x = " << x << std::endl
+ << "Function value when evaluating function & derivative: " << testFuncValue << std::endl
+ << "Function value when evaluating only function: " << tmpFunc << std::endl;
+ return;
+ }
+ if (testDerValue != tmpDer)
+ {
+ ADD_FAILURE()
+ << "Interpolation inconsistency for table " << desc << std::endl
+ << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
+ << "First failure at x = " << x << std::endl
+ << "Derivative value when evaluating function & derivative: " << testDerValue << std::endl
+ << "Derivative value when evaluating only derivative: " << tmpDer << std::endl;
+ return;
+ }
+
+ // There are two sources of errors that we need to account for when checking the values,
+ // and we only fail the test if both of these tolerances are violated:
+ //
+ // 1) First, we have the normal relative error of the test vs. reference value. For this
+ // we use the normal testutils relative tolerance checking.
+ // However, there is an additional source of error: When we calculate the forces we
+ // use average higher derivatives over the interval to improve accuracy, but this
+ // also means we won't reproduce values at table points exactly. This is usually not
+ // an issue since the tolerances we have are much larger, but when the reference derivative
+ // value is exactly zero the relative error will be infinite. To account for this, we
+ // extract the spacing from the table and evaluate the reference derivative at a point
+ // this much larger too, and use the largest of the two values as the reference
+ // magnitude for the derivative when setting the relative tolerance.
+ // Note that according to the table function definitions, we should be allowed to evaluate
+ // it one table point beyond the range (this is done already for construction).
+ //
+ // 2) Second, due to the loss-of-accuracy when calculating the index through rtable
+ // there is an internal absolute tolerance that we can calculate.
+ // The source of this error is the subtraction eps=rtab-[rtab], which leaves an
+ // error proportional to eps_machine*rtab=eps_machine*x*tableScale.
+ // To lowest order, the term in the function and derivative values (respectively) that
+ // are proportional to eps will be the next-higher derivative multiplied by the spacing.
+ // This means the truncation error in the value is derivative*x*eps_machine, and in the
+ // derivative the error is 2nd_derivative*x*eps_machine.
+
+ real refFuncValue = refFunc(x);
+ real refDerValue = refDer(x);
+ real nextRefDerValue = refDer(x + table.tableSpacing());
+
+ real derMagnitude = std::max( std::abs(refDerValue), std::abs(nextRefDerValue));
+
+ // Since the reference magnitude will change over each interval we need to re-evaluate
+ // the derivative tolerance inside the loop.
+ FloatingPointTolerance derTolerance(relativeToleranceAsFloatingPoint(derMagnitude, tolerance_));
+
+ FloatingPointDifference funcDiff(refFuncValue, testFuncValue);
+ FloatingPointDifference derDiff(refDerValue, testDerValue);
+
+ real allowedAbsFuncErr = std::abs(refDerValue) * x * GMX_REAL_EPS;
+ real allowedAbsDerErr = std::abs(secondDerivative) * x * GMX_REAL_EPS;
+
+ if ((!funcTolerance.isWithin(funcDiff) && funcDiff.asAbsolute() > allowedAbsFuncErr) ||
+ (!derTolerance.isWithin(derDiff) && derDiff.asAbsolute() > allowedAbsDerErr))
+ {
+ ADD_FAILURE()
+ << "Failing comparison with function for table " << desc << std::endl
+ << numFuncInTable << " function(s) in table, testing index " << funcIndex << std::endl
+ << "Test range is ( " << testRange.first << " , " << testRange.second << " ) " << std::endl
+ << "Tolerance = " << tolerance_ << std::endl
+ << "First failure at x = " << x << std::endl
+ << "Reference function = " << refFuncValue << std::endl
+ << "Test table function = " << testFuncValue << std::endl
+ << "Reference derivative = " << refDerValue << std::endl
+ << "Test table derivative = " << testDerValue << std::endl;
+ return;
+ }
+ }
+}
+
+
+/*! \brief Function similar to coulomb electrostatics
+ *
+ * \param r argument
+ * \return r^-1
+ */
+double
+coulombFunction(double r)
+{
+ return 1.0/r;
+}
+
+/*! \brief Derivative (not force) of coulomb electrostatics
+ *
+ * \param r argument
+ * \return -r^-2
+ */
+double
+coulombDerivative(double r)
+{
+ return -1.0/(r*r);
+}
+
+/*! \brief Function similar to power-6 Lennard-Jones dispersion
+ *
+ * \param r argument
+ * \return r^-6
+ */
+double
+lj6Function(double r)
+{
+ return std::pow(r, -6.0);
+}
+
+/*! \brief Derivative (not force) of the power-6 Lennard-Jones dispersion
+ *
+ * \param r argument
+ * \return -6.0*r^-7
+ */
+double
+lj6Derivative(double r)
+{
+ return -6.0*std::pow(r, -7.0);
+}
+
+/*! \brief Function similar to power-12 Lennard-Jones repulsion
+ *
+ * \param r argument
+ * \return r^-12
+ */
+double
+lj12Function(double r)
+{
+ return std::pow(r, -12.0);
+}
+
+/*! \brief Derivative (not force) of the power-12 Lennard-Jones repulsion
+ *
+ * \param r argument
+ * \return -12.0*r^-13
+ */
+double
+lj12Derivative(double r)
+{
+ return -12.0*std::pow(r, -13.0);
+}
+
+/*! \brief The sinc function, sin(r)/r
+ *
+ * \param r argument
+ * \return sin(r)/r
+ */
+double
+sincFunction(double r)
+{
+ return std::sin(r)/r;
+}
+
+/*! \brief Derivative of the sinc function
+ *
+ * \param r argument
+ * \return derivative of sinc, (r*cos(r)-sin(r))/r^2
+ */
+double
+sincDerivative(double r)
+{
+ return (r*std::cos(r)-std::sin(r))/(r*r);
+}
+
+/*! \brief Function for the direct-space PME correction to 1/r
+ *
+ * \param r argument
+ * \return PME correction function, erf(r)/r
+ */
+double
+pmeCorrFunction(double r)
+{
+ if (r == 0)
+ {
+ return 2.0/std::sqrt(M_PI);
+ }
+ else
+ {
+ return std::erf(r)/r;
+ }
+}
+
+/*! \brief Derivative of the direct-space PME correction to 1/r
+ *
+ * \param r argument
+ * \return Derivative of the PME correction function.
+ */
+double
+pmeCorrDerivative(double r)
+{
+ if (r == 0)
+ {
+ return 0;
+ }
+ else
+ {
+ return (2.0*std::exp(-r*r)/std::sqrt(3.14159265358979323846)*r-erf(r))/(r*r);
+ }
+}
+
+/*! \brief Typed-test list. We test QuadraticSplineTable and CubicSplineTable
+ */
+typedef ::testing::Types<QuadraticSplineTable, CubicSplineTable> SplineTableTypes;
+TYPED_TEST_CASE(SplineTableTest, SplineTableTypes);
+
+
+TYPED_TEST(SplineTableTest, HandlesIncorrectInput)
+{
+ // negative range
+ EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {-1.0, 0.0}), gmx::InvalidInputError);
+
+ // Too small range
+ EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1.0, 1.00001}), gmx::InvalidInputError);
+
+ // bad tolerance
+ EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1.0, 2.0}, 1e-20), gmx::ToleranceError);
+
+ // Range is so close to 0.0 that table would require >1e6 points
+ EXPECT_THROW_GMX(TypeParam( {{"LJ12", lj12Function, lj12Derivative}}, {1e-4, 2.0}), gmx::ToleranceError);
+
+ // mismatching function/derivative
+ EXPECT_THROW_GMX(TypeParam( { {"BadLJ12", lj12Derivative, lj12Function}}, {1.0, 2.0}), gmx::InconsistentInputError);
+}
+
+
+#ifndef NDEBUG
+TYPED_TEST(SplineTableTest, CatchesOutOfRangeValues)
+{
+ TypeParam table( {{"LJ12", lj12Function, lj12Derivative}}, {0.2, 1.0});
+ real x, func, der;
+
+ x = -GMX_REAL_EPS;
+ EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
+
+ x = 1.0;
+ EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
+}
+#endif
+
+
+TYPED_TEST(SplineTableTest, Sinc)
+{
+ std::pair<real, real> range(0.1, 10);
+
+ TypeParam sincTable( {{"Sinc", sincFunction, sincDerivative}}, range);
+
+ TestFixture::testSplineTableAgainstFunctions("Sinc", sincFunction, sincDerivative, sincTable, range);
+}
+
+
+TYPED_TEST(SplineTableTest, LJ12)
+{
+ std::pair<real, real> range(0.2, 2.0);
+
+ TypeParam lj12Table( {{"LJ12", lj12Function, lj12Derivative}}, range);
+
+ TestFixture::testSplineTableAgainstFunctions("LJ12", lj12Function, lj12Derivative, lj12Table, range);
+}
+
+
+TYPED_TEST(SplineTableTest, PmeCorrection)
+{
+ std::pair<real, real> range(0.0, 4.0);
+ real tolerance = 1e-5;
+
+ TypeParam pmeCorrTable( {{"PMECorr", pmeCorrFunction, pmeCorrDerivative}}, range, tolerance);
+
+ TestFixture::setTolerance(tolerance);
+ TestFixture::testSplineTableAgainstFunctions("PMECorr", pmeCorrFunction, pmeCorrDerivative, pmeCorrTable, range);
+}
+
+
+
+TYPED_TEST(SplineTableTest, HandlesIncorrectNumericalInput)
+{
+ // Lengths do not match
+ std::vector<double> functionValues(10);
+ std::vector<double> derivativeValues(20);
+ EXPECT_THROW_GMX(TypeParam( {{"EmptyVectors", functionValues, derivativeValues, 0.001}},
+ {1.0, 2.0}), gmx::InconsistentInputError);
+
+ // Upper range is 2.0, spacing 0.1. This requires at least 21 points. Make sure we get an error for 20.
+ functionValues.resize(20);
+ derivativeValues.resize(20);
+ EXPECT_THROW_GMX(TypeParam( {{"EmptyVectors", functionValues, derivativeValues, 0.1}},
+ {1.0, 2.0}), gmx::InconsistentInputError);
+
+ // Create some test data
+ functionValues.clear();
+ derivativeValues.clear();
+
+ std::vector<double> badDerivativeValues;
+ double spacing = 1e-3;
+
+ for (std::size_t i = 0; i < 1001; i++)
+ {
+ double x = i * spacing;
+ double func = (x >= 0.1) ? lj12Function(x) : 0.0;
+ double der = (x >= 0.1) ? lj12Derivative(x) : 0.0;
+
+ functionValues.push_back(func);
+ derivativeValues.push_back(der);
+ badDerivativeValues.push_back(-der);
+ }
+
+ // Derivatives not consistent with function
+ EXPECT_THROW_GMX(TypeParam( {{"NumericalBadLJ12", functionValues, badDerivativeValues, spacing}},
+ {0.2, 1.0}), gmx::InconsistentInputError);
+
+ // Spacing 1e-3 is not sufficient for r^-12 in range [0.1,1.0]
+ // Make sure we get a tolerance error
+ EXPECT_THROW_GMX(TypeParam( {{"NumericalLJ12", functionValues, derivativeValues, spacing}},
+ {0.2, 1.0}), gmx::ToleranceError);
+}
+
+
+TYPED_TEST(SplineTableTest, NumericalInputPmeCorr)
+{
+ std::pair<real, real> range(0.0, 4.0);
+ std::vector<double> functionValues;
+ std::vector<double> derivativeValues;
+
+ double inputSpacing = 1e-3;
+ real tolerance = 1e-5;
+
+ // We only need data up to the argument 4.0, but add 1% margin
+ for (std::size_t i = 0; i < range.second*1.01/inputSpacing; i++)
+ {
+ double x = i * inputSpacing;
+
+ functionValues.push_back(pmeCorrFunction(x));
+ derivativeValues.push_back(pmeCorrDerivative(x));
+ }
+
+ TypeParam pmeCorrTable( {{"NumericalPMECorr", functionValues, derivativeValues, inputSpacing}},
+ range, tolerance);
+
+ TestFixture::setTolerance(tolerance);
+ TestFixture::testSplineTableAgainstFunctions("NumericalPMECorr", pmeCorrFunction, pmeCorrDerivative, pmeCorrTable, range);
+}
+
+TYPED_TEST(SplineTableTest, TwoFunctions)
+{
+ std::pair<real, real> range(0.2, 2.0);
+
+ TypeParam table( {{"LJ6", lj6Function, lj6Derivative}, {"LJ12", lj12Function, lj12Derivative}}, range);
+
+ // Test entire range for each function. This will use the method that interpolates a single function
+ TestFixture::template testSplineTableAgainstFunctions<2, 0>("LJ6", lj6Function, lj6Derivative, table, range);
+ TestFixture::template testSplineTableAgainstFunctions<2, 1>("LJ12", lj12Function, lj12Derivative, table, range);
+
+ // Test the methods that evaluated both functions for one value
+ real x = 0.5 * (range.first + range.second);
+ real refFunc0 = lj6Function(x);
+ real refDer0 = lj6Derivative(x);
+ real refFunc1 = lj12Function(x);
+ real refDer1 = lj12Derivative(x);
+
+ real tstFunc0, tstDer0, tstFunc1, tstDer1;
+ real tmpFunc0, tmpFunc1, tmpDer0, tmpDer1;
+
+ // test that we reproduce the reference functions
+ table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
+
+ real funcErr0 = std::abs(tstFunc0-refFunc0) / std::abs(refFunc0);
+ real funcErr1 = std::abs(tstFunc1-refFunc1) / std::abs(refFunc1);
+ real derErr0 = std::abs(tstDer0-refDer0) / std::abs(refDer0);
+ real derErr1 = std::abs(tstDer1-refDer1) / std::abs(refDer1);
+
+ // Use asserts, since the following ones compare to these values.
+ ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
+ ASSERT_LT(derErr0, TypeParam::defaultTolerance);
+ ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
+ ASSERT_LT(derErr1, TypeParam::defaultTolerance);
+
+ // Test that function/derivative-only interpolation methods work
+ table.evaluateFunction(x, &tmpFunc0, &tmpFunc1);
+ table.evaluateDerivative(x, &tmpDer0, &tmpDer1);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstDer0, tmpDer0);
+
+ // Test that scrambled order interpolation methods work
+ table.template evaluateFunctionAndDerivative<2, 1, 0>(x, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+
+ // Test scrambled order for function/derivative-only methods
+ table.template evaluateFunction<2, 1, 0>(x, &tmpFunc1, &tmpFunc0);
+ table.template evaluateDerivative<2, 1, 0>(x, &tmpDer1, &tmpDer0);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+}
+
+TYPED_TEST(SplineTableTest, ThreeFunctions)
+{
+ std::pair<real, real> range(0.2, 2.0);
+
+ TypeParam table( {{"Coulomb", coulombFunction, coulombDerivative}, {"LJ6", lj6Function, lj6Derivative}, {"LJ12", lj12Function, lj12Derivative}}, range);
+
+ // Test entire range for each function
+ TestFixture::template testSplineTableAgainstFunctions<3, 0>("Coulomb", coulombFunction, coulombDerivative, table, range);
+ TestFixture::template testSplineTableAgainstFunctions<3, 1>("LJ6", lj6Function, lj6Derivative, table, range);
+ TestFixture::template testSplineTableAgainstFunctions<3, 2>("LJ12", lj12Function, lj12Derivative, table, range);
+
+ // Test the methods that evaluated both functions for one value
+ real x = 0.5 * (range.first + range.second);
+ real refFunc0 = coulombFunction(x);
+ real refDer0 = coulombDerivative(x);
+ real refFunc1 = lj6Function(x);
+ real refDer1 = lj6Derivative(x);
+ real refFunc2 = lj12Function(x);
+ real refDer2 = lj12Derivative(x);
+
+ real tstFunc0, tstDer0, tstFunc1, tstDer1, tstFunc2, tstDer2;
+ real tmpFunc0, tmpFunc1, tmpFunc2, tmpDer0, tmpDer1, tmpDer2;
+
+ // test that we reproduce the reference functions
+ table.evaluateFunctionAndDerivative(x, &tstFunc0, &tstDer0, &tstFunc1, &tstDer1, &tstFunc2, &tstDer2);
+
+ real funcErr0 = std::abs(tstFunc0-refFunc0) / std::abs(refFunc0);
+ real derErr0 = std::abs(tstDer0-refDer0) / std::abs(refDer0);
+ real funcErr1 = std::abs(tstFunc1-refFunc1) / std::abs(refFunc1);
+ real derErr1 = std::abs(tstDer1-refDer1) / std::abs(refDer1);
+ real funcErr2 = std::abs(tstFunc2-refFunc2) / std::abs(refFunc2);
+ real derErr2 = std::abs(tstDer2-refDer2) / std::abs(refDer2);
+
+ // Use asserts, since the following ones compare to these values.
+ ASSERT_LT(funcErr0, TypeParam::defaultTolerance);
+ ASSERT_LT(derErr0, TypeParam::defaultTolerance);
+ ASSERT_LT(funcErr1, TypeParam::defaultTolerance);
+ ASSERT_LT(derErr1, TypeParam::defaultTolerance);
+ ASSERT_LT(funcErr2, TypeParam::defaultTolerance);
+ ASSERT_LT(derErr2, TypeParam::defaultTolerance);
+
+ // Test that function/derivative-only interpolation methods work
+ table.evaluateFunction(x, &tmpFunc0, &tmpFunc1, &tmpFunc2);
+ table.evaluateDerivative(x, &tmpDer0, &tmpDer1, &tmpDer2);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstFunc2, tmpFunc2);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+ EXPECT_EQ(tstDer2, tmpDer2);
+
+ // Test two functions out of three
+ table.template evaluateFunctionAndDerivative<3, 0, 1>(x, &tmpFunc0, &tmpDer0, &tmpFunc1, &tmpDer1);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+
+ // two out of three, function/derivative-only
+ table.template evaluateFunction<3, 0, 1>(x, &tmpFunc0, &tmpFunc1);
+ table.template evaluateDerivative<3, 0, 1>(x, &tmpDer0, &tmpDer1);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+
+ // Test that scrambled order interpolation methods work
+ table.template evaluateFunctionAndDerivative<3, 2, 1, 0>(x, &tstFunc2, &tstDer2, &tstFunc1, &tstDer1, &tstFunc0, &tstDer0);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstFunc2, tmpFunc2);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+ EXPECT_EQ(tstDer2, tmpDer2);
+
+ // Test scrambled order for function/derivative-only methods
+ table.template evaluateFunction<3, 2, 1, 0>(x, &tmpFunc2, &tmpFunc1, &tmpFunc0);
+ table.template evaluateDerivative<3, 2, 1, 0>(x, &tmpDer2, &tmpDer1, &tmpDer0);
+ EXPECT_EQ(tstFunc0, tmpFunc0);
+ EXPECT_EQ(tstFunc1, tmpFunc1);
+ EXPECT_EQ(tstFunc2, tmpFunc2);
+ EXPECT_EQ(tstDer0, tmpDer0);
+ EXPECT_EQ(tstDer1, tmpDer1);
+ EXPECT_EQ(tstDer2, tmpDer2);
+}
+
+#if GMX_SIMD_HAVE_REAL
+TYPED_TEST(SplineTableTest, Simd)
+{
+ std::pair<real, real> range(0.2, 1.0);
+ TypeParam table( {{"LJ12", lj12Function, lj12Derivative}}, range);
+
+ // We already test that the SIMD operations handle the different elements
+ // correctly in the SIMD module, so here we only test that interpolation
+ // works for a single value in the middle of the interval
+
+ real x = 0.5 * (range.first + range.second);
+ real refFunc = lj12Function(x);
+ real refDer = lj12Derivative(x);
+ SimdReal tstFunc, tstDer;
+ real funcErr, derErr;
+ GMX_ALIGNED(real, GMX_SIMD_REAL_WIDTH) alignedMem[GMX_SIMD_REAL_WIDTH];
+
+ table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc, &tstDer);
+
+ store(alignedMem, tstFunc);
+ funcErr = std::abs(alignedMem[0]-refFunc) / std::abs(refFunc);
+
+ store(alignedMem, tstDer);
+ derErr = std::abs(alignedMem[0]-refDer ) / std::abs(refDer );
+
+ EXPECT_LT(funcErr, TypeParam::defaultTolerance);
+ EXPECT_LT(derErr, TypeParam::defaultTolerance);
+}
+
+TYPED_TEST(SplineTableTest, SimdTwoFunctions)
+{
+ std::pair<real, real> range(0.2, 2.0);
+
+ TypeParam table( {{"LJ6", lj6Function, lj6Derivative}, {"LJ12", lj12Function, lj12Derivative}}, range);
+
+ // We already test that the SIMD operations handle the different elements
+ // correctly in the SIMD module, so here we only test that interpolation
+ // works for a single value in the middle of the interval
+
+ real x = 0.5 * (range.first + range.second);
+ real refFunc0 = lj6Function(x);
+ real refDer0 = lj6Derivative(x);
+ real refFunc1 = lj12Function(x);
+ real refDer1 = lj12Derivative(x);
+ SimdReal tstFunc0, tstDer0;
+ SimdReal tstFunc1, tstDer1;
+ real funcErr0, derErr0;
+ real funcErr1, derErr1;
+ GMX_ALIGNED(real, GMX_SIMD_REAL_WIDTH) alignedMem[GMX_SIMD_REAL_WIDTH];
+
+ table.evaluateFunctionAndDerivative(SimdReal(x), &tstFunc0, &tstDer0, &tstFunc1, &tstDer1);
+
+ store(alignedMem, tstFunc0);
+ funcErr0 = std::abs(alignedMem[0]-refFunc0) / std::abs(refFunc0);
+
+ store(alignedMem, tstDer0);
+ derErr0 = std::abs(alignedMem[0]-refDer0 ) / std::abs(refDer0 );
+
+ store(alignedMem, tstFunc1);
+ funcErr1 = std::abs(alignedMem[0]-refFunc1) / std::abs(refFunc1);
+
+ store(alignedMem, tstDer1);
+ derErr1 = std::abs(alignedMem[0]-refDer1 ) / std::abs(refDer1 );
+
+ EXPECT_LT(funcErr0, TypeParam::defaultTolerance);
+ EXPECT_LT(derErr0, TypeParam::defaultTolerance);
+ EXPECT_LT(funcErr1, TypeParam::defaultTolerance);
+ EXPECT_LT(derErr1, TypeParam::defaultTolerance);
+}
+#endif
+
+#if GMX_SIMD_HAVE_REAL && !defined NDEBUG
+TYPED_TEST(SplineTableTest, CatchesOutOfRangeValuesSimd)
+{
+ std::pair<real, real> range(0.2, 1.0);
+ TypeParam table( {{"LJ12", lj12Function, lj12Derivative}}, range);
+ SimdReal x, func, der;
+
+ GMX_ALIGNED(real, GMX_SIMD_REAL_WIDTH) alignedMem[GMX_SIMD_REAL_WIDTH];
+
+ // Make position 1 incorrect if width>=2, otherwise position 0
+ alignedMem[ (GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = -GMX_REAL_EPS;
+ x = load(alignedMem);
+
+ EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
+
+ for (std::size_t i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
+ {
+ alignedMem[i] = range.second*(1.0-GMX_REAL_EPS);
+ }
+ // Make position 1 incorrect if width>=2, otherwise position 0
+ alignedMem[ (GMX_SIMD_REAL_WIDTH >= 2) ? 1 : 0] = range.second;
+ x = load(alignedMem);
+
+ EXPECT_THROW_GMX(table.evaluateFunctionAndDerivative(x, &func, &der), gmx::RangeError);
+}
+#endif
+
+} // namespace
+
+} // namespace test
+
+} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
"System I/O error",
"Error in user input",
"Inconsistency in user input",
+ "Requested tolerance cannot be achieved",
"Simulation instability detected",
"Feature not implemented",
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013,2015, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2015,2016, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
eeInvalidInput,
//! Invalid user input (conflicting or unsupported settings).
eeInconsistentInput,
+ //! Requested tolerance cannot be achieved.
+ eeTolerance,
//! Simulation instability detected.
eeInstability,
return eeInconsistentInput;
}
+int ToleranceError::errorCode() const
+{
+ return eeTolerance;
+}
+
int SimulationInstabilityError::errorCode() const
{
return eeInstability;
return eeAPIError;
}
+int RangeError::errorCode() const
+{
+ return eeRange;
+}
+
int NotImplementedError::errorCode() const
{
return eeNotImplemented;
virtual int errorCode() const;
};
+/*! \brief
+ * Exception class when a specified tolerance cannot be achieved.
+ *
+ * \inpublicapi
+ */
+class ToleranceError : public GromacsException
+{
+ public:
+ /*! \brief
+ * Creates an exception object with the provided initializer/reason.
+ *
+ * \param[in] details Initializer for the exception.
+ * \throws std::bad_alloc if out of memory.
+ *
+ * It is possible to call this constructor either with an explicit
+ * ExceptionInitializer object (useful for more complex cases), or
+ * a simple string if only a reason string needs to be provided.
+ */
+ explicit ToleranceError(const ExceptionInitializer &details)
+ : GromacsException(details) {}
+
+ virtual int errorCode() const;
+};
+
/*! \brief
* Exception class for simulation instabilities.
*
virtual int errorCode() const;
};
+/*! \brief
+ * Exception class for out-of-range values or indices
+ *
+ * \inpublicapi
+ */
+class RangeError : public GromacsException
+{
+ public:
+ //! \copydoc FileIOError::FileIOError()
+ explicit RangeError(const ExceptionInitializer &details)
+ : GromacsException(details) {}
+
+ virtual int errorCode() const;
+};
+
/*! \brief
* Exception class for use of an unimplemented feature.
*
virtual int errorCode() const;
};
-
/*! \brief
* Macro for throwing an exception.
*
// a reasonable thing to do...
const double relativeTolerance
= difference.isDouble() ? doubleRelativeTolerance_ : singleRelativeTolerance_;
+
if (difference.asAbsolute() <= relativeTolerance * difference.termMagnitude())
{
return true;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff