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.
37 #include "gromacs/simd/simd_math.h"
43 #include "gromacs/math/utilities.h"
44 #include "gromacs/options/basicoptions.h"
45 #include "gromacs/simd/simd.h"
55 /*! \addtogroup module_simd */
58 #ifdef GMX_SIMD_HAVE_REAL
60 class SimdMathTest : public SimdTest
63 ::testing::AssertionResult
64 compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
65 real refFunc(real x), gmx_simd_real_t gmx_simdcall simdFunc(gmx_simd_real_t x));
68 /*! \brief Test approximate equality of SIMD vs reference version of a function.
70 * This macro takes vanilla C and SIMD flavors of a function and tests it with
71 * the number of points, range, and tolerances specified by the test fixture class.
73 #define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc) \
74 EXPECT_PRED_FORMAT2(compareSimdMathFunction, refFunc, tstFunc)
76 /*! \brief Implementation routine to compare SIMD vs reference functions.
78 * \param refFuncExpr Description of reference function expression
79 * \param simdFuncExpr Description of SIMD function expression
80 * \param refFunc Reference math function pointer
81 * \param simdFunc SIMD math function pointer
83 * The function will be tested with the range and tolerances specified in
84 * the SimdBaseTest class. You should not never call this function directly,
85 * but use the macro GMX_EXPECT_SIMD_FUNC_NEAR(refFunc,tstFunc) instead.
87 ::testing::AssertionResult
88 SimdMathTest::compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
89 real refFunc(real x), gmx_simd_real_t gmx_simdcall simdFunc(gmx_simd_real_t x))
91 std::vector<real> vx(GMX_SIMD_REAL_WIDTH);
92 std::vector<real> vref(GMX_SIMD_REAL_WIDTH);
93 std::vector<real> vtst(GMX_SIMD_REAL_WIDTH);
95 gmx_int64_t ulpDiff, maxUlpDiff;
97 real refValMaxUlpDiff, simdValMaxUlpDiff;
100 int niter = s_nPoints/GMX_SIMD_REAL_WIDTH;
101 int npoints = niter*GMX_SIMD_REAL_WIDTH;
104 double r; gmx_int64_t i;
108 float r; gmx_int32_t i;
113 dx = (range_.second-range_.first)/npoints;
115 for (iter = 0; iter < niter; iter++)
117 for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
119 vx[i] = range_.first+dx*(iter*GMX_SIMD_REAL_WIDTH+i);
120 vref[i] = refFunc(vx[i]);
122 vtst = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
124 for (i = 0, signOk = true, absOk = true; i < GMX_SIMD_REAL_WIDTH; i++)
126 absDiff = fabs(vref[i]-vtst[i]);
127 absOk = absOk && ( absDiff < absTol_ );
128 signOk = signOk && ( vref[i]*vtst[i] >= 0 );
130 if (absDiff >= absTol_)
132 /* We replicate the trivial ulp differences comparison here rather than
133 * calling the lower-level routine for comparing them, since this enables
134 * us to run through the entire test range and report the largest deviation
135 * without lots of extra glue routines.
139 ulpDiff = llabs(conv0.i-conv1.i);
140 if (ulpDiff > maxUlpDiff)
142 maxUlpDiff = ulpDiff;
143 maxUlpDiffPos = vx[i];
144 refValMaxUlpDiff = vref[i];
145 simdValMaxUlpDiff = vtst[i];
149 if ( (absOk == false) && (signOk == false) )
151 return ::testing::AssertionFailure()
152 << "Failing SIMD math function comparison due to sign differences." << std::endl
153 << "Reference function: " << refFuncExpr << std::endl
154 << "Simd function: " << simdFuncExpr << std::endl
155 << "Test range is ( " << range_.first << " , " << range_.second << " ) " << std::endl
156 << "First sign difference around x=" << std::setprecision(20) << ::testing::PrintToString(vx) << std::endl
157 << "Ref values: " << std::setprecision(20) << ::testing::PrintToString(vref) << std::endl
158 << "SIMD values: " << std::setprecision(20) << ::testing::PrintToString(vtst) << std::endl;
162 if (maxUlpDiff <= ulpTol_)
164 return ::testing::AssertionSuccess();
168 return ::testing::AssertionFailure()
169 << "Failing SIMD math function ulp comparison between " << refFuncExpr << " and " << simdFuncExpr << std::endl
170 << "Requested ulp tolerance: " << ulpTol_ << std::endl
171 << "Requested abs tolerance: " << absTol_ << std::endl
172 << "Largest Ulp difference occurs for x=" << std::setprecision(20) << maxUlpDiffPos << std::endl
173 << "Ref values: " << std::setprecision(20) << refValMaxUlpDiff << std::endl
174 << "SIMD values: " << std::setprecision(20) << simdValMaxUlpDiff << std::endl
175 << "Ulp diff.: " << std::setprecision(20) << maxUlpDiff << std::endl;
183 // Actual math function tests below
189 /*! \cond internal */
190 /*! \addtogroup module_simd */
193 TEST_F(SimdMathTest, gmxSimdXorSignR)
195 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-4, 5, 6), gmx_simd_xor_sign_r(setSimdRealFrom3R(4, 5, 6), setSimdRealFrom3R(-5, 2, 0)));
196 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(4, -5, -6), gmx_simd_xor_sign_r(setSimdRealFrom3R(-4, -5, -6), setSimdRealFrom3R(-5, 2, 0)));
199 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
206 TEST_F(SimdMathTest, gmxSimdInvsqrtR)
208 setRange(1e-10, 1e10);
209 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, gmx_simd_invsqrt_r);
212 /*! \brief Function wrapper to return first result when testing \ref gmx_simd_invsqrt_pair_r */
213 gmx_simd_real_t gmx_simdcall
214 tst_invsqrt_pair0(gmx_simd_real_t x)
216 gmx_simd_real_t r0, r1;
217 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
221 /*! \brief Function wrapper to return second result when testing \ref gmx_simd_invsqrt_pair_r */
222 gmx_simd_real_t gmx_simdcall
223 tst_invsqrt_pair1(gmx_simd_real_t x)
225 gmx_simd_real_t r0, r1;
226 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
230 TEST_F(SimdMathTest, gmxSimdInvsqrtPairR)
232 setRange(1e-10, 1e10);
233 // The accuracy conversions lose a bit of extra accuracy compared to
234 // doing the iterations in all-double.
235 setUlpTol(4*ulpTol_);
237 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair0);
238 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair1);
241 TEST_F(SimdMathTest, gmxSimdSqrtR)
243 // Just make sure sqrt(0)=0 works and isn't evaluated as 0*1/sqrt(0)=NaN
244 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, 2, 3), gmx_simd_sqrt_r(setSimdRealFrom3R(0, 4, 9)));
247 /*! \brief Function wrapper to evaluate reference 1/x */
253 TEST_F(SimdMathTest, gmxSimdInvR)
256 setRange(-1e10, -1e-10);
257 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
258 setRange(1e-10, 1e10);
259 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
262 /*! \brief Function wrapper for log(x), with argument/return in default Gromacs precision */
268 TEST_F(SimdMathTest, gmxSimdLogR)
270 setRange(1e-30, 1e30);
271 GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_r);
274 // MSVC does not support exp2(), so we have no reference to test against
276 /*! \brief Function wrapper for exp2(x), with argument/return in default Gromacs precision */
277 real ref_exp2(real x)
282 TEST_F(SimdMathTest, gmxSimdExp2R)
285 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_r);
287 // We do not care about the SIMD implementation getting denormal values right,
288 // but they must be clamped to zero rather than producing garbage.
289 // Check by setting the absolute tolerance to machine precision.
290 setAbsTol(GMX_REAL_EPS);
292 // First two values will have denormal results in single, third value in double too.
293 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-150.0), ref_exp2(-300.0), ref_exp2(-1050.0)),
294 gmx_simd_exp2_r(setSimdRealFrom3R(-150.0, -300.0, -1050.0)));
296 // Reset absolute tolerance to enforce ULP checking
299 // Make sure that underflowing values are set to zero.
300 // First two values underflow in single, third value in double too.
301 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-200.0), ref_exp2(-600.0), ref_exp2(-1500.0)),
302 gmx_simd_exp2_r(setSimdRealFrom3R(-200.0, -600.0, -1500.0)));
306 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
312 TEST_F(SimdMathTest, gmxSimdExpR)
315 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp, gmx_simd_exp_r);
317 // We do not care about the SIMD implementation getting denormal values right,
318 // but they must be clamped to zero rather than producing garbage.
319 // Check by setting the absolute tolerance to machine precision.
320 setAbsTol(GMX_REAL_EPS);
321 // First two values will have denormal results in single, third value in double too.
322 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp(-90.0), ref_exp(-100.0), ref_exp(-725.0)),
323 gmx_simd_exp_r(setSimdRealFrom3R(-90.0, -100.0, -725.0)));
325 // Reset absolute tolerance to enforce ULP checking
328 // Make sure that underflowing values are set to zero.
329 // First two values underflow in single, third value in double too.
330 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp(-150.0), ref_exp(-300.0), ref_exp(-800.0)),
331 gmx_simd_exp_r(setSimdRealFrom3R(-150.0, -300.0, -800.0)));
334 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
336 * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
337 * than the SIMD flavor, so we use double for reference.
344 TEST_F(SimdMathTest, gmxSimdErfR)
347 setAbsTol(GMX_REAL_MIN);
348 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erf, gmx_simd_erf_r);
351 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
353 * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
354 * than the SIMD flavor, so we use double for reference.
356 real ref_erfc(real x)
361 TEST_F(SimdMathTest, gmxSimdErfcR)
364 setAbsTol(GMX_REAL_MIN);
365 // Our erfc algorithm has 4 ulp accuracy, so relax defaultTol a bit
366 setUlpTol(4*ulpTol_);
367 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erfc, gmx_simd_erfc_r);
370 /*! \brief Function wrapper for sin(x), with argument/return in default Gromacs precision */
376 TEST_F(SimdMathTest, gmxSimdSinR)
378 setRange(-8*M_PI, 8*M_PI);
379 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
380 // Range reduction leads to accuracy loss, so we might want higher tolerance here
381 setRange(-10000, 10000);
382 setUlpTol(2*ulpTol_);
383 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
386 /*! \brief Function wrapper for cos(x), with argument/return in default Gromacs precision */
392 TEST_F(SimdMathTest, gmxSimdCosR)
394 setRange(-8*M_PI, 8*M_PI);
395 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
396 // Range reduction leads to accuracy loss, so we might want higher tolerance here
397 setRange(-10000, 10000);
398 setUlpTol(2*ulpTol_);
399 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
402 /*! \brief Function wrapper for tan(x), with argument/return in default Gromacs precision */
408 TEST_F(SimdMathTest, gmxSimdTanR)
410 // Tan(x) is a little sensitive due to the division in the algorithm.
411 // Rather than using lots of extra FP operations, we accept the algorithm
412 // presently only achieves a ~3 ulp error and use the medium tolerance.
413 setRange(-8*M_PI, 8*M_PI);
414 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
415 // Range reduction leads to accuracy loss, so we might want higher tolerance here
416 setRange(-10000, 10000);
417 setUlpTol(2*ulpTol_);
418 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_r);
421 /*! \brief Function wrapper for asin(x), with argument/return in default Gromacs precision */
422 real ref_asin(real x)
427 TEST_F(SimdMathTest, gmxSimdAsinR)
429 // Our present asin(x) algorithm achieves 2-3 ulp accuracy
431 GMX_EXPECT_SIMD_FUNC_NEAR(ref_asin, gmx_simd_asin_r);
434 /*! \brief Function wrapper for acos(x), with argument/return in default Gromacs precision */
435 real ref_acos(real x)
440 TEST_F(SimdMathTest, gmxSimdAcosR)
442 // Our present acos(x) algorithm achieves 2-3 ulp accuracy
444 GMX_EXPECT_SIMD_FUNC_NEAR(ref_acos, gmx_simd_acos_r);
447 /*! \brief Function wrapper for atan(x), with argument/return in default Gromacs precision */
448 real ref_atan(real x)
453 TEST_F(SimdMathTest, gmxSimdAtanR)
455 // Our present atan(x) algorithm achieves 1 ulp accuracy
456 setRange(-10000, 10000);
457 GMX_EXPECT_SIMD_FUNC_NEAR(ref_atan, gmx_simd_atan_r);
460 TEST_F(SimdMathTest, gmxSimdAtan2R)
462 // test each quadrant
463 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_1_2_3));
464 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_1_2_3));
465 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, -1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_m1_m2_m3));
466 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, -1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_m1_m2_m3));
467 // cases important for calculating angles
468 // values on coordinate axes
469 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, 1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_1_2_3));
470 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 0.0)), gmx_simd_atan2_r(rSimd_1_2_3, gmx_simd_setzero_r()));
471 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, -1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_m1_m2_m3));
472 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 0.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, gmx_simd_setzero_r()));
473 // degenerate value (origin) should return 0.0
474 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()));
477 /*! \brief Evaluate reference version of PME force correction. */
478 real ref_pmecorrF(real x)
481 return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
484 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
485 TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
487 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
489 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
491 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
495 setAbsTol(GMX_REAL_EPS);
496 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrF, gmx_simd_pmecorrF_r);
499 /*! \brief Evaluate reference version of PME potential correction. */
500 real ref_pmecorrV(real x)
503 return gmx_erfd(y)/y;
506 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
507 TEST_F(SimdMathTest, gmxSimdPmecorrPotentialR)
509 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
511 setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
513 setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
516 setAbsTol(GMX_REAL_EPS);
517 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrV, gmx_simd_pmecorrV_r);
522 #endif // GMX_SIMD_HAVE_REAL