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);
283 // We do not care about the SIMD implementation getting denormal values right,
284 // but they must be clamped to zero rather than producing garbage.
285 // Check by setting the absolute tolerance to machine precision.
286 setAbsTol(GMX_REAL_EPS);
288 // First two values will have denormal results in single, third value in double too.
289 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-150.0), ref_exp2(-300.0), ref_exp2(-1050.0)),
290 gmx_simd_exp2_r(setSimdRealFrom3R(-150.0, -300.0, -1050.0)));
292 // Reset absolute tolerance to enforce ULP checking
295 // Make sure that underflowing values are set to zero.
296 // First two values underflow in single, third value in double too.
297 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-200.0), ref_exp2(-600.0), ref_exp2(-1500.0)),
298 gmx_simd_exp2_r(setSimdRealFrom3R(-200.0, -600.0, -1500.0)));
302 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
308 TEST_F(SimdMathTest, gmxSimdExpR)
311 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp, gmx_simd_exp_r);
313 // We do not care about the SIMD implementation getting denormal values right,
314 // but they must be clamped to zero rather than producing garbage.
315 // Check by setting the absolute tolerance to machine precision.
316 setAbsTol(GMX_REAL_EPS);
317 // First two values will have denormal results in single, third value in double too.
318 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp(-90.0), ref_exp(-100.0), ref_exp(-725.0)),
319 gmx_simd_exp_r(setSimdRealFrom3R(-90.0, -100.0, -725.0)));
321 // Reset absolute tolerance to enforce ULP checking
324 // Make sure that underflowing values are set to zero.
325 // First two values underflow in single, third value in double too.
326 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp(-150.0), ref_exp(-300.0), ref_exp(-800.0)),
327 gmx_simd_exp_r(setSimdRealFrom3R(-150.0, -300.0, -800.0)));
330 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
332 * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
333 * than the SIMD flavor, so we use double for reference.
340 TEST_F(SimdMathTest, gmxSimdErfR)
343 setAbsTol(GMX_REAL_MIN);
344 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erf, gmx_simd_erf_r);
347 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
349 * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
350 * than the SIMD flavor, so we use double for reference.
352 real ref_erfc(real x)
357 TEST_F(SimdMathTest, gmxSimdErfcR)
360 setAbsTol(GMX_REAL_MIN);
361 // Our erfc algorithm has 4 ulp accuracy, so relax defaultTol a bit
362 setUlpTol(4*ulpTol_);
363 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erfc, gmx_simd_erfc_r);
366 /*! \brief Function wrapper for sin(x), with argument/return in default Gromacs precision */
372 TEST_F(SimdMathTest, gmxSimdSinR)
374 setRange(-8*M_PI, 8*M_PI);
375 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
376 // Range reduction leads to accuracy loss, so we might want higher tolerance here
377 setRange(-10000, 10000);
378 setUlpTol(2*ulpTol_);
379 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
382 /*! \brief Function wrapper for cos(x), with argument/return in default Gromacs precision */
388 TEST_F(SimdMathTest, gmxSimdCosR)
390 setRange(-8*M_PI, 8*M_PI);
391 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
392 // Range reduction leads to accuracy loss, so we might want higher tolerance here
393 setRange(-10000, 10000);
394 setUlpTol(2*ulpTol_);
395 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
398 /*! \brief Function wrapper for tan(x), with argument/return in default Gromacs precision */
404 TEST_F(SimdMathTest, gmxSimdTanR)
406 // Tan(x) is a little sensitive due to the division in the algorithm.
407 // Rather than using lots of extra FP operations, we accept the algorithm
408 // presently only achieves a ~3 ulp error and use the medium tolerance.
409 setRange(-8*M_PI, 8*M_PI);
410 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
411 // Range reduction leads to accuracy loss, so we might want higher tolerance here
412 setRange(-10000, 10000);
413 setUlpTol(2*ulpTol_);
414 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
417 /*! \brief Function wrapper for asin(x), with argument/return in default Gromacs precision */
418 real ref_asin(real x)
423 TEST_F(SimdMathTest, gmxSimdAsinR)
425 // Our present asin(x) algorithm achieves 2-3 ulp accuracy
427 GMX_EXPECT_SIMD_FUNC_NEAR(ref_asin, gmx_simd_asin_r);
430 /*! \brief Function wrapper for acos(x), with argument/return in default Gromacs precision */
431 real ref_acos(real x)
436 TEST_F(SimdMathTest, gmxSimdAcosR)
438 // Our present acos(x) algorithm achieves 2-3 ulp accuracy
440 GMX_EXPECT_SIMD_FUNC_NEAR(ref_acos, gmx_simd_acos_r);
443 /*! \brief Function wrapper for atan(x), with argument/return in default Gromacs precision */
444 real ref_atan(real x)
449 TEST_F(SimdMathTest, gmxSimdAtanR)
451 // Our present atan(x) algorithm achieves 1 ulp accuracy
452 setRange(-10000, 10000);
453 GMX_EXPECT_SIMD_FUNC_NEAR(ref_atan, gmx_simd_atan_r);
456 TEST_F(SimdMathTest, gmxSimdAtan2R)
458 // test each quadrant
459 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_1_2_3));
460 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_1_2_3));
461 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, -1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_m1_m2_m3));
462 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, -1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_m1_m2_m3));
463 // cases important for calculating angles
464 // values on coordinate axes
465 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, 1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_1_2_3));
466 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 0.0)), gmx_simd_atan2_r(rSimd_1_2_3, gmx_simd_setzero_r()));
467 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, -1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_m1_m2_m3));
468 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 0.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, gmx_simd_setzero_r()));
469 // degenerate value (origin) should return 0.0
470 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()));
473 /*! \brief Evaluate reference version of PME force correction. */
474 real ref_pmecorrF(real x)
477 return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
480 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
481 TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
483 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
485 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
487 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
491 setAbsTol(GMX_REAL_EPS);
492 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrF, gmx_simd_pmecorrF_r);
495 /*! \brief Evaluate reference version of PME potential correction. */
496 real ref_pmecorrV(real x)
499 return gmx_erfd(y)/y;
502 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
503 TEST_F(SimdMathTest, gmxSimdPmecorrPotentialR)
505 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
507 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
509 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
512 setAbsTol(GMX_REAL_EPS);
513 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrV, gmx_simd_pmecorrV_r);
518 #endif // GMX_SIMD_HAVE_REAL