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.
40 #include "gromacs/math/utilities.h"
41 #include "gromacs/simd/simd.h"
42 #include "gromacs/simd/simd_math.h"
43 #include "gromacs/options/basicoptions.h"
53 /*! \addtogroup module_simd */
56 #ifdef GMX_SIMD_HAVE_REAL
58 class SimdMathTest : public SimdTest
61 ::testing::AssertionResult
62 compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
63 real refFunc(real x), gmx_simd_real_t simdFunc(gmx_simd_real_t x));
66 /*! \brief Test approximate equality of SIMD vs reference version of a function.
68 * This macro takes vanilla C and SIMD flavors of a function and tests it with
69 * the number of points, range, and tolerances specified by the test fixture class.
71 #define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc) \
72 EXPECT_PRED_FORMAT2(compareSimdMathFunction, refFunc, tstFunc)
74 /*! \brief Implementation routine to compare SIMD vs reference functions.
76 * \param refFuncExpr Description of reference function expression
77 * \param simdFuncExpr Description of SIMD function expression
78 * \param refFunc Reference math function pointer
79 * \param simdFunc SIMD math function pointer
81 * The function will be tested with the range and tolerances specified in
82 * the SimdBaseTest class. You should not never call this function directly,
83 * but use the macro GMX_EXPECT_SIMD_FUNC_NEAR(refFunc,tstFunc) instead.
85 ::testing::AssertionResult
86 SimdMathTest::compareSimdMathFunction(const char * refFuncExpr, const char *simdFuncExpr,
87 real refFunc(real x), gmx_simd_real_t simdFunc(gmx_simd_real_t x))
89 std::vector<real> vx(GMX_SIMD_REAL_WIDTH);
90 std::vector<real> vref(GMX_SIMD_REAL_WIDTH);
91 std::vector<real> vtst(GMX_SIMD_REAL_WIDTH);
93 gmx_int64_t ulpDiff, maxUlpDiff;
95 real refValMaxUlpDiff, simdValMaxUlpDiff;
98 int niter = s_nPoints/GMX_SIMD_REAL_WIDTH;
99 int npoints = niter*GMX_SIMD_REAL_WIDTH;
102 double r; gmx_int64_t i;
106 float r; gmx_int32_t i;
111 dx = (range_.second-range_.first)/npoints;
113 for (iter = 0; iter < niter; iter++)
115 for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
117 vx[i] = range_.first+dx*(iter*GMX_SIMD_REAL_WIDTH+i);
118 vref[i] = refFunc(vx[i]);
120 vtst = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
122 for (i = 0, signOk = true, absOk = true; i < GMX_SIMD_REAL_WIDTH; i++)
124 absDiff = fabs(vref[i]-vtst[i]);
125 absOk = absOk && ( absDiff < absTol_ );
126 signOk = signOk && ( vref[i]*vtst[i] >= 0 );
128 if (absDiff >= absTol_)
130 /* We replicate the trivial ulp differences comparison here rather than
131 * calling the lower-level routine for comparing them, since this enables
132 * us to run through the entire test range and report the largest deviation
133 * without lots of extra glue routines.
137 ulpDiff = llabs(conv0.i-conv1.i);
138 if (ulpDiff > maxUlpDiff)
140 maxUlpDiff = ulpDiff;
141 maxUlpDiffPos = vx[i];
142 refValMaxUlpDiff = vref[i];
143 simdValMaxUlpDiff = vtst[i];
147 if ( (absOk == false) && (signOk == false) )
149 return ::testing::AssertionFailure()
150 << "Failing SIMD math function comparison due to sign differences." << std::endl
151 << "Reference function: " << refFuncExpr << std::endl
152 << "Simd function: " << simdFuncExpr << std::endl
153 << "Test range is ( " << range_.first << " , " << range_.second << " ) " << std::endl
154 << "First sign difference around x=" << std::setprecision(20) << ::testing::PrintToString(vx) << std::endl
155 << "Ref values: " << std::setprecision(20) << ::testing::PrintToString(vref) << std::endl
156 << "SIMD values: " << std::setprecision(20) << ::testing::PrintToString(vtst) << std::endl;
160 if (maxUlpDiff <= ulpTol_)
162 return ::testing::AssertionSuccess();
166 return ::testing::AssertionFailure()
167 << "Failing SIMD math function ulp comparison between " << refFuncExpr << " and " << simdFuncExpr << std::endl
168 << "Requested ulp tolerance: " << ulpTol_ << std::endl
169 << "Requested abs tolerance: " << absTol_ << std::endl
170 << "Largest Ulp difference occurs for x=" << std::setprecision(20) << maxUlpDiffPos << std::endl
171 << "Ref values: " << std::setprecision(20) << refValMaxUlpDiff << std::endl
172 << "SIMD values: " << std::setprecision(20) << simdValMaxUlpDiff << std::endl
173 << "Ulp diff.: " << std::setprecision(20) << maxUlpDiff << std::endl;
181 // Actual math function tests below
187 /*! \cond internal */
188 /*! \addtogroup module_simd */
191 TEST_F(SimdMathTest, gmxSimdXorSignR)
193 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-4, 5, 6), gmx_simd_xor_sign_r(setSimdRealFrom3R(4, 5, 6), setSimdRealFrom3R(-5, 2, 0)));
194 GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(4, -5, -6), gmx_simd_xor_sign_r(setSimdRealFrom3R(-4, -5, -6), setSimdRealFrom3R(-5, 2, 0)));
197 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
204 TEST_F(SimdMathTest, gmxSimdInvsqrtR)
206 setRange(1e-10, 1e10);
207 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, gmx_simd_invsqrt_r);
210 /*! \brief Function wrapper to return first result when testing \ref gmx_simd_invsqrt_pair_r */
212 tst_invsqrt_pair0(gmx_simd_real_t x)
214 gmx_simd_real_t r0, r1;
215 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
219 /*! \brief Function wrapper to return second result when testing \ref gmx_simd_invsqrt_pair_r */
221 tst_invsqrt_pair1(gmx_simd_real_t x)
223 gmx_simd_real_t r0, r1;
224 gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
228 TEST_F(SimdMathTest, gmxSimdInvsqrtPairR)
230 setRange(1e-10, 1e10);
231 // The accuracy conversions lose a bit of extra accuracy compared to
232 // doing the iterations in all-double.
233 setUlpTol(4*ulpTol_);
235 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair0);
236 GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair1);
239 TEST_F(SimdMathTest, gmxSimdSqrtR)
241 // Just make sure sqrt(0)=0 works and isn't evaluated as 0*1/sqrt(0)=NaN
242 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, 2, 3), gmx_simd_sqrt_r(setSimdRealFrom3R(0, 4, 9)));
245 /*! \brief Function wrapper to evaluate reference 1/x */
251 TEST_F(SimdMathTest, gmxSimdInvR)
254 setRange(-1e10, -1e-10);
255 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
256 setRange(1e-10, 1e10);
257 GMX_EXPECT_SIMD_FUNC_NEAR(ref_inv, gmx_simd_inv_r);
260 /*! \brief Function wrapper for log(x), with argument/return in default Gromacs precision */
266 TEST_F(SimdMathTest, gmxSimdLogR)
268 setRange(1e-30, 1e30);
269 GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_r);
272 // MSVC does not support exp2(), so we have no reference to test against
274 /*! \brief Function wrapper for exp2(x), with argument/return in default Gromacs precision */
275 real ref_exp2(real x)
280 TEST_F(SimdMathTest, gmxSimdExp2R)
283 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_r);
287 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
293 TEST_F(SimdMathTest, gmxSimdExpR)
296 GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp, gmx_simd_exp_r);
299 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
301 * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
302 * than the SIMD flavor, so we use double for reference.
309 TEST_F(SimdMathTest, gmxSimdErfR)
312 setAbsTol(GMX_REAL_MIN);
313 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erf, gmx_simd_erf_r);
316 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
318 * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
319 * than the SIMD flavor, so we use double for reference.
321 real ref_erfc(real x)
326 TEST_F(SimdMathTest, gmxSimdErfcR)
329 setAbsTol(GMX_REAL_MIN);
330 // Our erfc algorithm has 4 ulp accuracy, so relax defaultTol a bit
331 setUlpTol(4*ulpTol_);
332 GMX_EXPECT_SIMD_FUNC_NEAR(ref_erfc, gmx_simd_erfc_r);
335 /*! \brief Function wrapper for sin(x), with argument/return in default Gromacs precision */
341 TEST_F(SimdMathTest, gmxSimdSinR)
343 setRange(-8*M_PI, 8*M_PI);
344 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
345 // Range reduction leads to accuracy loss, so we might want higher tolerance here
346 setRange(-10000, 10000);
347 setUlpTol(2*ulpTol_);
348 GMX_EXPECT_SIMD_FUNC_NEAR(ref_sin, gmx_simd_sin_r);
351 /*! \brief Function wrapper for cos(x), with argument/return in default Gromacs precision */
357 TEST_F(SimdMathTest, gmxSimdCosR)
359 setRange(-8*M_PI, 8*M_PI);
360 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
361 // Range reduction leads to accuracy loss, so we might want higher tolerance here
362 setRange(-10000, 10000);
363 setUlpTol(2*ulpTol_);
364 GMX_EXPECT_SIMD_FUNC_NEAR(ref_cos, gmx_simd_cos_r);
367 /*! \brief Function wrapper for tan(x), with argument/return in default Gromacs precision */
373 TEST_F(SimdMathTest, gmxSimdTanR)
375 // Tan(x) is a little sensitive due to the division in the algorithm.
376 // Rather than using lots of extra FP operations, we accept the algorithm
377 // presently only achieves a ~3 ulp error and use the medium tolerance.
378 setRange(-8*M_PI, 8*M_PI);
379 GMX_EXPECT_SIMD_FUNC_NEAR(ref_tan, gmx_simd_tan_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_tan, gmx_simd_tan_r);
386 /*! \brief Function wrapper for asin(x), with argument/return in default Gromacs precision */
387 real ref_asin(real x)
392 TEST_F(SimdMathTest, gmxSimdAsinR)
394 // Our present asin(x) algorithm achieves 2-3 ulp accuracy
396 GMX_EXPECT_SIMD_FUNC_NEAR(ref_asin, gmx_simd_asin_r);
399 /*! \brief Function wrapper for acos(x), with argument/return in default Gromacs precision */
400 real ref_acos(real x)
405 TEST_F(SimdMathTest, gmxSimdAcosR)
407 // Our present acos(x) algorithm achieves 2-3 ulp accuracy
409 GMX_EXPECT_SIMD_FUNC_NEAR(ref_acos, gmx_simd_acos_r);
412 /*! \brief Function wrapper for atan(x), with argument/return in default Gromacs precision */
413 real ref_atan(real x)
418 TEST_F(SimdMathTest, gmxSimdAtanR)
420 // Our present atan(x) algorithm achieves 1 ulp accuracy
421 setRange(-10000, 10000);
422 GMX_EXPECT_SIMD_FUNC_NEAR(ref_atan, gmx_simd_atan_r);
425 TEST_F(SimdMathTest, gmxSimdAtan2R)
427 // test each quadrant
428 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_1_2_3));
429 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_1_2_3));
430 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, -1.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, rSimd_m1_m2_m3));
431 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, -1.0)), gmx_simd_atan2_r(rSimd_1_2_3, rSimd_m1_m2_m3));
432 // cases important for calculating angles
433 // values on coordinate axes
434 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, 1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_1_2_3));
435 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(1.0, 0.0)), gmx_simd_atan2_r(rSimd_1_2_3, gmx_simd_setzero_r()));
436 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(0.0, -1.0)), gmx_simd_atan2_r(gmx_simd_setzero_r(), rSimd_m1_m2_m3));
437 GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(atan2(-1.0, 0.0)), gmx_simd_atan2_r(rSimd_m1_m2_m3, gmx_simd_setzero_r()));
438 // degenerate value (origin) should return 0.0
439 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()));
442 /*! \brief Evaluate reference version of PME force correction. */
443 real ref_pmecorrF(real x)
446 return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
449 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
450 TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
452 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
454 setUlpTol((gmx_int64_t)(1e-10/GMX_REAL_EPS));
456 setUlpTol((gmx_int64_t)(1e-6/GMX_REAL_EPS));
460 setAbsTol(GMX_REAL_EPS);
461 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrF, gmx_simd_pmecorrF_r);
464 /*! \brief Evaluate reference version of PME potential correction. */
465 real ref_pmecorrV(real x)
468 return gmx_erfd(y)/y;
471 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
472 TEST_F(SimdMathTest, gmxSimdPmecorrPotentialR)
474 // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
476 setUlpTol((gmx_int64_t)(1e-10/GMX_REAL_EPS));
478 setUlpTol((gmx_int64_t)(1e-6/GMX_REAL_EPS));
481 setAbsTol(GMX_REAL_EPS);
482 GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrV, gmx_simd_pmecorrV_r);
487 #endif // GMX_SIMD_HAVE_REAL