src/gromacs/simd/tests/simd_math.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 #include "config.h"
  36
  37 #include <vector>
  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"
  42
  43 #include "simd.h"
  44
  45 namespace gmx
  46 {
  47 namespace test
  48 {
  49
  50 /*! \cond internal */
  51 /*! \addtogroup module_simd */
  52 /*! \{ */
  53
  54 #ifdef GMX_SIMD_HAVE_REAL
  55
  56 class SimdMathTest : public SimdTest
  57 {
  58     public:
  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));
  62 };
  63
  64 /*! \brief Test approximate equality of SIMD vs reference version of a function.
  65  *
  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.
  68  */
  69 #define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc) \
  70     EXPECT_PRED_FORMAT2(compareSimdMathFunction, refFunc, tstFunc)
  71
  72 /*! \brief Implementation routine to compare SIMD vs reference functions.
  73  *
  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
  78  *
  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.
  82  */
  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))
  86 {
  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);
  90     real                         dx, absDiff;
  91     gmx_int64_t                  ulpDiff, maxUlpDiff;
  92     real                         maxUlpDiffPos;
  93     real                         refValMaxUlpDiff, simdValMaxUlpDiff;
  94     bool                         absOk, signOk;
  95     int                          i, iter;
  96     int                          niter   = s_nPoints/GMX_SIMD_REAL_WIDTH;
  97     int                          npoints = niter*GMX_SIMD_REAL_WIDTH;
  98 #    ifdef GMX_DOUBLE
  99     union {
 100         double r; gmx_int64_t i;
 101     } conv0, conv1;
 102 #    else
 103     union {
 104         float  r; gmx_int32_t i;
 105     } conv0, conv1;
 106 #    endif
 107
 108     maxUlpDiff = 0;
 109     dx         = (range_.second-range_.first)/npoints;
 110
 111     for (iter = 0; iter < niter; iter++)
 112     {
 113         for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 114         {
 115             vx[i]   = range_.first+dx*(iter*GMX_SIMD_REAL_WIDTH+i);
 116             vref[i] = refFunc(vx[i]);
 117         }
 118         vtst  = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
 119
 120         for (i = 0, signOk = true, absOk = true; i < GMX_SIMD_REAL_WIDTH; i++)
 121         {
 122             absDiff = fabs(vref[i]-vtst[i]);
 123             absOk   = absOk  && ( absDiff < absTol_ );
 124             signOk  = signOk && ( vref[i]*vtst[i] >= 0 );
 125
 126             if (absDiff >= absTol_)
 127             {
 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.
 132                  */
 133                 conv0.r           = vref[i];
 134                 conv1.r           = vtst[i];
 135                 ulpDiff           = llabs(conv0.i-conv1.i);
 136                 if (ulpDiff > maxUlpDiff)
 137                 {
 138                     maxUlpDiff        = ulpDiff;
 139                     maxUlpDiffPos     = vx[i];
 140                     refValMaxUlpDiff  = vref[i];
 141                     simdValMaxUlpDiff = vtst[i];
 142                 }
 143             }
 144         }
 145         if ( (absOk == false) && (signOk == false) )
 146         {
 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;
 155         }
 156     }
 157
 158     if (maxUlpDiff <= ulpTol_)
 159     {
 160         return ::testing::AssertionSuccess();
 161     }
 162     else
 163     {
 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;
 172     }
 173 }
 174
 175 /*! \} */
 176 /*! \endcond */
 177
 178
 179 // Actual math function tests below
 180
 181
 182 namespace
 183 {
 184
 185 /*! \cond internal */
 186 /*! \addtogroup module_simd */
 187 /*! \{ */
 188
 189 TEST_F(SimdMathTest, gmxSimdXorSignR)
 190 {
 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)));
 193 }
 194
 195 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
 196 static real
 197 ref_invsqrt(real x)
 198 {
 199     return 1.0/sqrt(x);
 200 }
 201
 202 TEST_F(SimdMathTest, gmxSimdInvsqrtR)
 203 {
 204     setRange(1e-10, 1e10);
 205     GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, gmx_simd_invsqrt_r);
 206 }
 207
 208 /*! \brief Function wrapper to return first result when testing \ref gmx_simd_invsqrt_pair_r */
 209 gmx_simd_real_t
 210 tst_invsqrt_pair0(gmx_simd_real_t x)
 211 {
 212     gmx_simd_real_t r0, r1;
 213     gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
 214     return r0;
 215 }
 216
 217 /*! \brief Function wrapper to return second result when testing \ref gmx_simd_invsqrt_pair_r */
 218 gmx_simd_real_t
 219 tst_invsqrt_pair1(gmx_simd_real_t x)
 220 {
 221     gmx_simd_real_t r0, r1;
 222     gmx_simd_invsqrt_pair_r(x, x, &r0, &r1);
 223     return r1;
 224 }
 225
 226 TEST_F(SimdMathTest, gmxSimdInvsqrtPairR)
 227 {
 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_);
 232
 233     GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair0);
 234     GMX_EXPECT_SIMD_FUNC_NEAR(ref_invsqrt, tst_invsqrt_pair1);
 235 }
 236
 237 TEST_F(SimdMathTest, gmxSimdSqrtR)
 238 {
 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)));
 241 }
 242
 243 /*! \brief Function wrapper to evaluate reference 1/x */
 244 real ref_inv(real x)
 245 {
 246     return 1.0/x;
 247 }
 248
 249 TEST_F(SimdMathTest, gmxSimdInvR)
 250 {
 251     // test <0
 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);
 256 }
 257
 258 /*! \brief Function wrapper for log(x), with argument/return in default Gromacs precision */
 259 real ref_log(real x)
 260 {
 261     return log(x);
 262 }
 263
 264 TEST_F(SimdMathTest, gmxSimdLogR)
 265 {
 266     setRange(1e-30, 1e30);
 267     GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_r);
 268 }
 269
 270 // MSVC does not support exp2(), so we have no reference to test against
 271 #ifndef _MSC_VER
 272 /*! \brief Function wrapper for exp2(x), with argument/return in default Gromacs precision */
 273 real ref_exp2(real x)
 274 {
 275     return exp2(x);
 276 }
 277
 278 TEST_F(SimdMathTest, gmxSimdExp2R)
 279 {
 280     setRange(-100, 100);
 281     GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_r);
 282 }
 283 #endif
 284
 285 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
 286 real ref_exp(real x)
 287 {
 288     return exp(x);
 289 }
 290
 291 TEST_F(SimdMathTest, gmxSimdExpR)
 292 {
 293     setRange(-75, 75);
 294     GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp, gmx_simd_exp_r);
 295 }
 296
 297 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
 298  *
 299  * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
 300  * than the SIMD flavor, so we use double for reference.
 301  */
 302 real ref_erf(real x)
 303 {
 304     return gmx_erfd(x);
 305 }
 306
 307 TEST_F(SimdMathTest, gmxSimdErfR)
 308 {
 309     setRange(-9, 9);
 310     setAbsTol(GMX_REAL_MIN);
 311     GMX_EXPECT_SIMD_FUNC_NEAR(ref_erf, gmx_simd_erf_r);
 312 }
 313
 314 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
 315  *
 316  * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
 317  * than the SIMD flavor, so we use double for reference.
 318  */
 319 real ref_erfc(real x)
 320 {
 321     return gmx_erfcd(x);
 322 }
 323
 324 TEST_F(SimdMathTest, gmxSimdErfcR)
 325 {
 326     setRange(-9, 9);
 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);
 331 }
 332
 333 /*! \brief Function wrapper for sin(x), with argument/return in default Gromacs precision */
 334 real ref_sin(real x)
 335 {
 336     return sin(x);
 337 }
 338
 339 TEST_F(SimdMathTest, gmxSimdSinR)
 340 {
 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);
 347 }
 348
 349 /*! \brief Function wrapper for cos(x), with argument/return in default Gromacs precision */
 350 real ref_cos(real x)
 351 {
 352     return cos(x);
 353 }
 354
 355 TEST_F(SimdMathTest, gmxSimdCosR)
 356 {
 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);
 363 }
 364
 365 /*! \brief Function wrapper for tan(x), with argument/return in default Gromacs precision */
 366 real ref_tan(real x)
 367 {
 368     return tan(x);
 369 }
 370
 371 TEST_F(SimdMathTest, gmxSimdTanR)
 372 {
 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);
 382 }
 383
 384 /*! \brief Function wrapper for asin(x), with argument/return in default Gromacs precision */
 385 real ref_asin(real x)
 386 {
 387     return asin(x);
 388 }
 389
 390 TEST_F(SimdMathTest, gmxSimdAsinR)
 391 {
 392     // Our present asin(x) algorithm achieves 2-3 ulp accuracy
 393     setRange(-1, 1);
 394     GMX_EXPECT_SIMD_FUNC_NEAR(ref_asin, gmx_simd_asin_r);
 395 }
 396
 397 /*! \brief Function wrapper for acos(x), with argument/return in default Gromacs precision */
 398 real ref_acos(real x)
 399 {
 400     return acos(x);
 401 }
 402
 403 TEST_F(SimdMathTest, gmxSimdAcosR)
 404 {
 405     // Our present acos(x) algorithm achieves 2-3 ulp accuracy
 406     setRange(-1, 1);
 407     GMX_EXPECT_SIMD_FUNC_NEAR(ref_acos, gmx_simd_acos_r);
 408 }
 409
 410 /*! \brief Function wrapper for atan(x), with argument/return in default Gromacs precision */
 411 real ref_atan(real x)
 412 {
 413     return atan(x);
 414 }
 415
 416 TEST_F(SimdMathTest, gmxSimdAtanR)
 417 {
 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);
 421 }
 422
 423 TEST_F(SimdMathTest, gmxSimdAtan2R)
 424 {
 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()));
 438 }
 439
 440 /*! \brief Evaluate reference version of PME force correction. */
 441 real ref_pmecorrF(real x)
 442 {
 443     real y = sqrt(x);
 444     return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
 445 }
 446
 447 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
 448 TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
 449 {
 450     // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
 451 #ifdef GMX_DOUBLE
 452     setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
 453 #else
 454     setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
 455 #endif
 456
 457     setRange(0.15, 4);
 458     setAbsTol(GMX_REAL_EPS);
 459     GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrF, gmx_simd_pmecorrF_r);
 460 }
 461
 462 /*! \brief Evaluate reference version of PME potential correction. */
 463 real ref_pmecorrV(real x)
 464 {
 465     real y = sqrt(x);
 466     return gmx_erfd(y)/y;
 467 }
 468
 469 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
 470 TEST_F(SimdMathTest, gmxSimdPmecorrPotentialR)
 471 {
 472     // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
 473 #ifdef GMX_DOUBLE
 474     setUlpTol((gmx_int64_t)(5e-10/GMX_REAL_EPS));
 475 #else
 476     setUlpTol((gmx_int64_t)(5e-6/GMX_REAL_EPS));
 477 #endif
 478     setRange(0.15, 4);
 479     setAbsTol(GMX_REAL_EPS);
 480     GMX_EXPECT_SIMD_FUNC_NEAR(ref_pmecorrV, gmx_simd_pmecorrV_r);
 481 }
 482
 483 }      // namespace
 484
 485 #endif // GMX_SIMD_HAVE_REAL
 486
 487 /*! \} */
 488 /*! \endcond */
 489
 490 }      // namespace
 491 }      // namespace