From 8d54138d3e87cfa20d41f6a355879192a016dc13 Mon Sep 17 00:00:00 2001 From: Nick Thompson <4nt@ornl.gov> Date: Mon, 5 Apr 2021 15:20:56 -0400 Subject: [PATCH] Kahan's difference of products algorithm --- vtkm/Math.h | 11 +++++++ vtkm/Math.h.in | 11 +++++++ vtkm/Matrix.h | 2 +- vtkm/VectorAnalysis.h | 4 ++- vtkm/testing/UnitTestMath.cxx | 23 +++++++++++++++ vtkm/testing/UnitTestVectorAnalysis.cxx | 38 +++++++------------------ 6 files changed, 60 insertions(+), 29 deletions(-) diff --git a/vtkm/Math.h b/vtkm/Math.h index 586f543cf..9e7caf23b 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -2696,6 +2696,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 return xi - yi; } +// Computes ab - cd. +// See: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html +template +inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) +{ + T cd = c * d; + T err = std::fma(-c, d, cd); + T dop = std::fma(a, b, -cd); + return dop + err; +} + /// Bitwise operations /// diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index 6e5145fbf..d0f3f68d8 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -1298,6 +1298,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 return xi - yi; } +// Computes ab - cd. +// See: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html +template +inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) +{ + T cd = c * d; + T err = std::fma(-c, d, cd); + T dop = std::fma(a, b, -cd); + return dop + err; +} + /// Bitwise operations /// diff --git a/vtkm/Matrix.h b/vtkm/Matrix.h index 2ed4f102b..33c93bfff 100644 --- a/vtkm/Matrix.h +++ b/vtkm/Matrix.h @@ -519,7 +519,7 @@ VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix& A) template VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix& A) { - return A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1); + return vtkm::DifferenceOfProducts(A(0, 0), A(1, 1), A(1, 0), A(0, 1)); } template diff --git a/vtkm/VectorAnalysis.h b/vtkm/VectorAnalysis.h index 44e78e21b..691eee18d 100644 --- a/vtkm/VectorAnalysis.h +++ b/vtkm/VectorAnalysis.h @@ -179,7 +179,9 @@ VTKM_EXEC_CONT vtkm::Vec::Type, 3> C const vtkm::Vec& y) { return vtkm::Vec::Type, 3>( - x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0]); + DifferenceOfProducts(x[1], y[2], x[2], y[1]), + DifferenceOfProducts(x[2], y[0], x[0], y[2]), + DifferenceOfProducts(x[0], y[1], x[1], y[0])); } //----------------------------------------------------------------------------- diff --git a/vtkm/testing/UnitTestMath.cxx b/vtkm/testing/UnitTestMath.cxx index 85ebd4bc0..5e2e1a52f 100644 --- a/vtkm/testing/UnitTestMath.cxx +++ b/vtkm/testing/UnitTestMath.cxx @@ -901,6 +901,29 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase } } + VTKM_EXEC + void TestDifferenceOfProducts() const + { + { + // Example taken from: + // https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html + vtkm::Float32 a = 33962.035f; + vtkm::Float32 b = -30438.8f; + vtkm::Float32 c = 41563.4f; + vtkm::Float32 d = -24871.969f; + vtkm::Float32 computed = vtkm::DifferenceOfProducts(a, b, c, d); + // Expected result, computed in double precision and cast back to float: + vtkm::Float32 expected = 5.376600027084351f; + + vtkm::UInt64 dist = vtkm::FloatDistance(expected, computed); + VTKM_MATH_ASSERT( + dist < 2, + "Float distance for difference of products exceeds 1.5; this is in violation of a theorem " + "proved by Jeannerod in doi.org/10.1090/S0025-5718-2013-02679-8. Is your build compiled " + "with fma's enabled?"); + } + } + VTKM_EXEC void operator()(vtkm::Id) const { diff --git a/vtkm/testing/UnitTestVectorAnalysis.cxx b/vtkm/testing/UnitTestVectorAnalysis.cxx index 1b8bb3bbf..dd83473d9 100644 --- a/vtkm/testing/UnitTestVectorAnalysis.cxx +++ b/vtkm/testing/UnitTestVectorAnalysis.cxx @@ -60,31 +60,24 @@ void TestVector(const VectorType& vector) { using ComponentType = typename vtkm::VecTraits::ComponentType; - std::cout << "Testing " << vector << std::endl; - //to do have to implement a norm and normalized call to verify the math ones //against - std::cout << " Magnitude" << std::endl; ComponentType magnitude = vtkm::Magnitude(vector); ComponentType magnitudeCompare = internal::MyMag(vector); VTKM_TEST_ASSERT(test_equal(magnitude, magnitudeCompare), "Magnitude failed test."); - std::cout << " Magnitude squared" << std::endl; ComponentType magnitudeSquared = vtkm::MagnitudeSquared(vector); VTKM_TEST_ASSERT(test_equal(magnitude * magnitude, magnitudeSquared), "Magnitude squared test failed."); if (magnitudeSquared > 0) { - std::cout << " Reciprocal magnitude" << std::endl; ComponentType rmagnitude = vtkm::RMagnitude(vector); VTKM_TEST_ASSERT(test_equal(1 / magnitude, rmagnitude), "Reciprical magnitude failed."); - std::cout << " Normal" << std::endl; VTKM_TEST_ASSERT(test_equal(vtkm::Normal(vector), internal::MyNormal(vector)), "Normalized vector failed test."); - std::cout << " Normalize" << std::endl; VectorType normalizedVector = vector; vtkm::Normalize(normalizedVector); VTKM_TEST_ASSERT(test_equal(normalizedVector, internal::MyNormal(vector)), @@ -98,13 +91,11 @@ void TestLerp(const VectorType& a, const VectorType& w, const typename vtkm::VecTraits::ComponentType& wS) { - std::cout << "Linear interpolation: " << a << "-" << b << ": " << w << std::endl; VectorType vtkmLerp = vtkm::Lerp(a, b, w); VectorType otherLerp = internal::MyLerp(a, b, w); VTKM_TEST_ASSERT(test_equal(vtkmLerp, otherLerp), "Vectors with Vector weight do not lerp() correctly"); - std::cout << "Linear interpolation: " << a << "-" << b << ": " << wS << std::endl; VectorType lhsS = internal::MyLerp(a, b, wS); VectorType rhsS = vtkm::Lerp(a, b, wS); VTKM_TEST_ASSERT(test_equal(lhsS, rhsS), "Vectors with Scalar weight do not lerp() correctly"); @@ -113,19 +104,16 @@ void TestLerp(const VectorType& a, template void TestCross(const vtkm::Vec& x, const vtkm::Vec& y) { - std::cout << "Testing " << x << " x " << y << std::endl; - using Vec3 = vtkm::Vec; Vec3 cross = vtkm::Cross(x, y); - std::cout << " = " << cross << std::endl; - - std::cout << " Orthogonality" << std::endl; // The cross product result should be perpendicular to input vectors. - VTKM_TEST_ASSERT(test_equal(vtkm::Dot(cross, x), T(0.0)), "Cross product not perpendicular."); - VTKM_TEST_ASSERT(test_equal(vtkm::Dot(cross, y), T(0.0)), "Cross product not perpendicular."); - - std::cout << " Length" << std::endl; + VTKM_TEST_ASSERT(abs(vtkm::Dot(cross, x)) < + std::numeric_limits::epsilon() * vtkm::MagnitudeSquared(x), + "Cross product not perpendicular."); + VTKM_TEST_ASSERT(abs(vtkm::Dot(cross, y)) < + std::numeric_limits::epsilon() * vtkm::MagnitudeSquared(y), + "Cross product not perpendicular."); // The length of cross product should be the lengths of the input vectors // times the sin of the angle between them. T sinAngle = vtkm::Magnitude(cross) * vtkm::RMagnitude(x) * vtkm::RMagnitude(y); @@ -139,10 +127,10 @@ void TestCross(const vtkm::Vec& x, const vtkm::Vec& y) VTKM_TEST_ASSERT(test_equal(sinAngle * sinAngle + cosAngle * cosAngle, T(1.0)), "Bad cross product length."); - std::cout << " Triangle normal" << std::endl; // Test finding the normal to a triangle (similar to cross product). Vec3 normal = vtkm::TriangleNormal(x, y, Vec3(0, 0, 0)); - VTKM_TEST_ASSERT(test_equal(vtkm::Dot(normal, x - y), T(0.0)), + VTKM_TEST_ASSERT(abs(vtkm::Dot(normal, x - y)) < + std::numeric_limits::epsilon() * vtkm::MagnitudeSquared(x), "Triangle normal is not really normal."); } @@ -151,13 +139,6 @@ void TestOrthonormalize(const VectorBasisType& inputs, int expectedRank) { VectorBasisType outputs; int actualRank = vtkm::Orthonormalize(inputs, outputs); - std::cout << "Testing orthonormalize\n" - << " Rank " << actualRank << " expected " << expectedRank << "\n" - << " Basis vectors:\n"; - for (int i = 0; i < actualRank; ++i) - { - std::cout << " " << i << " " << outputs[i] << "\n"; - } VTKM_TEST_ASSERT(test_equal(actualRank, expectedRank), "Orthonormalized rank is unexpected."); } @@ -208,6 +189,9 @@ struct TestCrossFunctor TestCross(VectorType(1.0f, 0.0f, 0.0f), VectorType(0.0f, 1.0f, 0.0f)); TestCross(VectorType(1.0f, 2.0f, 3.0f), VectorType(-3.0f, -1.0f, 1.0f)); TestCross(VectorType(0.0f, 0.0f, 1.0f), VectorType(0.001f, 0.01f, 2.0f)); + // Example from: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html + TestCross(VectorType(33962.035f, 41563.4f, 7706.415f), + VectorType(-24871.969, -30438.8, -5643.727f)); } };