2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
38 #include "gromacs/math/utilities.h"
39 #include "gromacs/simd/simd.h"
40 #include "gromacs/simd/simd_math.h"
41 #include "gromacs/options/basicoptions.h"
51 /*! \addtogroup module_simd */
54 #ifdef GMX_SIMD_HAVE_REAL
56 class SimdMathTest : public SimdTest
59 ::testing::AssertionResult
60 compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
61 real refFunc(real x), gmx_simd_real_t simdFunc(gmx_simd_real_t x));
64 /*! \brief Test approximate equality of SIMD vs reference version of a function.
66 * This macro takes vanilla C and SIMD flavors of a function and tests it with
67 * the number of points, range, and tolerances specified by the test fixture class.
69 #define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc) \
70 EXPECT_PRED_FORMAT2(compareSimdMathFunction, refFunc, tstFunc)
72 /*! \brief Implementation routine to compare SIMD vs reference functions.
74 * \param refFuncExpr Description of reference function expression
75 * \param simdFuncExpr Description of SIMD function expression
76 * \param refFunc Reference math function pointer
77 * \param simdFunc SIMD math function pointer
79 * The function will be tested with the range and tolerances specified in
80 * the SimdBaseTest class. You should not never call this function directly,
81 * but use the macro GMX_EXPECT_SIMD_FUNC_NEAR(refFunc,tstFunc) instead.
83 ::testing::AssertionResult
84 SimdMathTest::compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
85 real refFunc(real x), gmx_simd_real_t simdFunc(gmx_simd_real_t x))
87 std::vector<real> vx(GMX_SIMD_REAL_WIDTH);
88 std::vector<real> vref(GMX_SIMD_REAL_WIDTH);
89 std::vector<real> vtst(GMX_SIMD_REAL_WIDTH);
91 gmx_int64_t ulpDiff, maxUlpDiff;
93 real refValMaxUlpDiff, simdValMaxUlpDiff;
96 int niter = s_nPoints/GMX_SIMD_REAL_WIDTH;
97 int npoints = niter*GMX_SIMD_REAL_WIDTH;
100 double r; gmx_int64_t i;
104 float r; gmx_int32_t i;
109 dx = (range_.second-range_.first)/npoints;
111 for (iter = 0; iter < niter; iter++)
113 for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
115 vx[i] = range_.first+dx*(iter*GMX_SIMD_REAL_WIDTH+i);
116 vref[i] = refFunc(vx[i]);
118 vtst = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
120 for (i = 0, signOk = true, absOk = true; i < GMX_SIMD_REAL_WIDTH; i++)
122 absDiff = fabs(vref[i]-vtst[i]);
123 absOk = absOk && ( absDiff < absTol_ );
124 signOk = signOk && ( vref[i]*vtst[i] >= 0 );
126 if (absDiff >= absTol_)
128 /* We replicate the trivial ulp differences comparison here rather than
129 * calling the lower-level routine for comparing them, since this enables
130 * us to run through the entire test range and report the largest deviation
131 * without lots of extra glue routines.
135 ulpDiff = llabs(conv0.i-conv1.i);
136 if (ulpDiff > maxUlpDiff)
138 maxUlpDiff = ulpDiff;
139 maxUlpDiffPos = vx[i];
140 refValMaxUlpDiff = vref[i];
141 simdValMaxUlpDiff = vtst[i];
145 if ( (absOk == false) && (signOk == false) )
147 return ::testing::AssertionFailure()
148 << "Failing SIMD math function comparison due to sign differences." << std::endl
149 << "Reference function: " << refFuncExpr << std::endl
150 << "Simd function: " << simdFuncExpr << std::endl
151 << "Test range is ( " << range_.first << " , " << range_.second << " ) " << std::endl
152 << "First sign difference around x=" << std::setprecision(20) << ::testing::PrintToString(vx) << std::endl
153 << "Ref values: " << std::setprecision(20) << ::testing::PrintToString(vref) << std::endl
154 << "SIMD values: " << std::setprecision(20) << ::testing::PrintToString(vtst) << std::endl;
158 if (maxUlpDiff <= ulpTol_)
160 return ::testing::AssertionSuccess();
164 return ::testing::AssertionFailure()
165 << "Failing SIMD math function ulp comparison between " << refFuncExpr << " and " << simdFuncExpr << std::endl
166 << "Requested ulp tolerance: " << ulpTol_ << std::endl
167 << "Requested abs tolerance: " << absTol_ << std::endl
168 << "Largest Ulp difference occurs for x=" << std::setprecision(20) << maxUlpDiffPos << std::endl
169 << "Ref values: " << std::setprecision(20) << refValMaxUlpDiff << std::endl
170 << "SIMD values: " << std::setprecision(20) << simdValMaxUlpDiff << std::endl
171 << "Ulp diff.: " << std::setprecision(20) << maxUlpDiff << std::endl;
179 // Actual math function tests below
185 /*! \cond internal */
186 /*! \addtogroup module_simd */
189 TEST_F(SimdMathTest, gmxSimdXorSignR)
191 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-4, 5, 6), gmx_simd_xor_sign_r(setSimdRealFrom3R(4, 5, 6), setSimdRealFrom3R(-5, 2, 0)));
192 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(4, -5, -6), gmx_simd_xor_sign_r(setSimdRealFrom3R(-4, -5, -6), setSimdRealFrom3R(-5, 2, 0)));
195 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
202 TEST_F(SimdMathTest, gmxSimdInvsqrtR)
204 setRange(1e-10, 1e10);
205 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, gmx_simd_invsqrt_r);
208 /*! \brief Function wrapper to return first result when testing \ref gmx_simd_invsqrt_pair_r */
210 tst_invsqrt_pair0(gmx_simd_real_t x)
212 gmx_simd_real_t r0, r1;
213 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
217 /*! \brief Function wrapper to return second result when testing \ref gmx_simd_invsqrt_pair_r */
219 tst_invsqrt_pair1(gmx_simd_real_t x)
221 gmx_simd_real_t r0, r1;
222 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
226 TEST_F(SimdMathTest, gmxSimdInvsqrtPairR)
228 setRange(1e-10, 1e10);
229 // The accuracy conversions lose a bit of extra accuracy compared to
230 // doing the iterations in all-double.
231 setUlpTol(4*ulpTol_);
233 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair0);
234 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair1);
237 TEST_F(SimdMathTest, gmxSimdSqrtR)
239 // Just make sure sqrt(0)=0 works and isn't evaluated as 0*1/sqrt(0)=NaN
240 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, 2, 3), gmx_simd_sqrt_r(setSimdRealFrom3R(0, 4, 9)));
243 /*! \brief Function wrapper to evaluate reference 1/x */
249 TEST_F(SimdMathTest, gmxSimdInvR)
252 setRange(-1e10, -1e-10);
253 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
254 setRange(1e-10, 1e10);
255 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
258 /*! \brief Function wrapper for log(x), with argument/return in default Gromacs precision */
264 TEST_F(SimdMathTest, gmxSimdLogR)
266 setRange(1e-30, 1e30);
267 GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_r);
270 // MSVC does not support exp2(), so we have no reference to test against
272 /*! \brief Function wrapper for exp2(x), with argument/return in default Gromacs precision */
273 real ref_exp2(real x)
278 TEST_F(SimdMathTest, gmxSimdExp2R)
281 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_r);
285 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
291 TEST_F(SimdMathTest, gmxSimdExpR)
294 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp, gmx_simd_exp_r);
297 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
299 * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
300 * than the SIMD flavor, so we use double for reference.
307 TEST_F(SimdMathTest, gmxSimdErfR)
310 setAbsTol(GMX_REAL_MIN);
311 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erf, gmx_simd_erf_r);
314 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
316 * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
317 * than the SIMD flavor, so we use double for reference.
319 real ref_erfc(real x)
324 TEST_F(SimdMathTest, gmxSimdErfcR)
327 setAbsTol(GMX_REAL_MIN);
328 // Our erfc algorithm has 4 ulp accuracy, so relax defaultTol a bit
329 setUlpTol(4*ulpTol_);
330 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erfc, gmx_simd_erfc_r);
333 /*! \brief Function wrapper for sin(x), with argument/return in default Gromacs precision */
339 TEST_F(SimdMathTest, gmxSimdSinR)
341 setRange(-8*M_PI, 8*M_PI);
342 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
343 // Range reduction leads to accuracy loss, so we might want higher tolerance here
344 setRange(-10000, 10000);
345 setUlpTol(2*ulpTol_);
346 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
349 /*! \brief Function wrapper for cos(x), with argument/return in default Gromacs precision */
355 TEST_F(SimdMathTest, gmxSimdCosR)
357 setRange(-8*M_PI, 8*M_PI);
358 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
359 // Range reduction leads to accuracy loss, so we might want higher tolerance here
360 setRange(-10000, 10000);
361 setUlpTol(2*ulpTol_);
362 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
365 /*! \brief Function wrapper for tan(x), with argument/return in default Gromacs precision */
371 TEST_F(SimdMathTest, gmxSimdTanR)
373 // Tan(x) is a little sensitive due to the division in the algorithm.
374 // Rather than using lots of extra FP operations, we accept the algorithm
375 // presently only achieves a ~3 ulp error and use the medium tolerance.
376 setRange(-8*M_PI, 8*M_PI);
377 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
378 // Range reduction leads to accuracy loss, so we might want higher tolerance here
379 setRange(-10000, 10000);
380 setUlpTol(2*ulpTol_);
381 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
384 /*! \brief Function wrapper for asin(x), with argument/return in default Gromacs precision */
385 real ref_asin(real x)
390 TEST_F(SimdMathTest, gmxSimdAsinR)
392 // Our present asin(x) algorithm achieves 2-3 ulp accuracy
394 GMX_EXPECT_SIMD_FUNC_NEAR(ref_asin, gmx_simd_asin_r);
397 /*! \brief Function wrapper for acos(x), with argument/return in default Gromacs precision */
398 real ref_acos(real x)
403 TEST_F(SimdMathTest, gmxSimdAcosR)
405 // Our present acos(x) algorithm achieves 2-3 ulp accuracy
407 GMX_EXPECT_SIMD_FUNC_NEAR(ref_acos, gmx_simd_acos_r);
410 /*! \brief Function wrapper for atan(x), with argument/return in default Gromacs precision */
411 real ref_atan(real x)
416 TEST_F(SimdMathTest, gmxSimdAtanR)
418 // Our present atan(x) algorithm achieves 1 ulp accuracy
419 setRange(-10000, 10000);
420 GMX_EXPECT_SIMD_FUNC_NEAR(ref_atan, gmx_simd_atan_r);
423 TEST_F(SimdMathTest, gmxSimdAtan2R)
425 // test each quadrant
426 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_1_2_3));
427 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_1_2_3));
428 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, -1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_m1_m2_m3));
429 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, -1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_m1_m2_m3));
430 // cases important for calculating angles
431 // values on coordinate axes
432 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, 1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_1_2_3));
433 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 0.0)), gmx_simd_atan2_r(rSimd_1_2_3, gmx_simd_setzero_r()));
434 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, -1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_m1_m2_m3));
435 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 0.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, gmx_simd_setzero_r()));
436 // degenerate value (origin) should return 0.0
437 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, 0.0)), gmx_simd_atan2_r(setSimdRealFrom3R(0.0, 0.0, 0.0), gmx_simd_setzero_r()));
440 /*! \brief Evaluate reference version of PME force correction. */
441 real ref_pmecorrF(real x)
444 return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
447 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
448 TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
450 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
452 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
454 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
458 setAbsTol(GMX_REAL_EPS);
459 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrF, gmx_simd_pmecorrF_r);
462 /*! \brief Evaluate reference version of PME potential correction. */
463 real ref_pmecorrV(real x)
466 return gmx_erfd(y)/y;
469 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
470 TEST_F(SimdMathTest, gmxSimdPmecorrPotentialR)
472 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
474 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
476 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
479 setAbsTol(GMX_REAL_EPS);
480 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrV, gmx_simd_pmecorrV_r);
485 #endif // GMX_SIMD_HAVE_REAL