Commit b00f5b44 authored by Nick Thompson's avatar Nick Thompson Committed by Kitware Robot
Browse files

Merge topic 'diff_of_products'

068399b8 Fix build errors on Windows and compiler warnings on Ubuntu.
8d54138d

 Kahan's difference of products algorithm
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland's avatarKenneth Moreland <kmorel@acm.org>
Merge-request: !2460
parents fdb7abe7 068399b8
...@@ -2696,6 +2696,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 ...@@ -2696,6 +2696,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32
return xi - yi; return xi - yi;
} }
// Computes ab - cd.
// See: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html
template<typename T>
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 /// Bitwise operations
/// ///
......
...@@ -1298,6 +1298,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 ...@@ -1298,6 +1298,17 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32
return xi - yi; return xi - yi;
} }
// Computes ab - cd.
// See: https://pharr.org/matt/blog/2019/11/03/difference-of-floats.html
template<typename T>
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 /// Bitwise operations
/// ///
......
...@@ -519,7 +519,7 @@ VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 1, 1>& A) ...@@ -519,7 +519,7 @@ VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 1, 1>& A)
template <typename T> template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 2, 2>& A) VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 2, 2>& 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 <typename T> template <typename T>
......
...@@ -179,7 +179,9 @@ VTKM_EXEC_CONT vtkm::Vec<typename detail::FloatingPointReturnType<T>::Type, 3> C ...@@ -179,7 +179,9 @@ VTKM_EXEC_CONT vtkm::Vec<typename detail::FloatingPointReturnType<T>::Type, 3> C
const vtkm::Vec<T, 3>& y) const vtkm::Vec<T, 3>& y)
{ {
return vtkm::Vec<typename detail::FloatingPointReturnType<T>::Type, 3>( return vtkm::Vec<typename detail::FloatingPointReturnType<T>::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]));
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
......
...@@ -901,6 +901,29 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase ...@@ -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 VTKM_EXEC
void operator()(vtkm::Id) const void operator()(vtkm::Id) const
{ {
......
...@@ -60,31 +60,24 @@ void TestVector(const VectorType& vector) ...@@ -60,31 +60,24 @@ void TestVector(const VectorType& vector)
{ {
using ComponentType = typename vtkm::VecTraits<VectorType>::ComponentType; using ComponentType = typename vtkm::VecTraits<VectorType>::ComponentType;
std::cout << "Testing " << vector << std::endl;
//to do have to implement a norm and normalized call to verify the math ones //to do have to implement a norm and normalized call to verify the math ones
//against //against
std::cout << " Magnitude" << std::endl;
ComponentType magnitude = vtkm::Magnitude(vector); ComponentType magnitude = vtkm::Magnitude(vector);
ComponentType magnitudeCompare = internal::MyMag(vector); ComponentType magnitudeCompare = internal::MyMag(vector);
VTKM_TEST_ASSERT(test_equal(magnitude, magnitudeCompare), "Magnitude failed test."); VTKM_TEST_ASSERT(test_equal(magnitude, magnitudeCompare), "Magnitude failed test.");
std::cout << " Magnitude squared" << std::endl;
ComponentType magnitudeSquared = vtkm::MagnitudeSquared(vector); ComponentType magnitudeSquared = vtkm::MagnitudeSquared(vector);
VTKM_TEST_ASSERT(test_equal(magnitude * magnitude, magnitudeSquared), VTKM_TEST_ASSERT(test_equal(magnitude * magnitude, magnitudeSquared),
"Magnitude squared test failed."); "Magnitude squared test failed.");
if (magnitudeSquared > 0) if (magnitudeSquared > 0)
{ {
std::cout << " Reciprocal magnitude" << std::endl;
ComponentType rmagnitude = vtkm::RMagnitude(vector); ComponentType rmagnitude = vtkm::RMagnitude(vector);
VTKM_TEST_ASSERT(test_equal(1 / magnitude, rmagnitude), "Reciprical magnitude failed."); 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)), VTKM_TEST_ASSERT(test_equal(vtkm::Normal(vector), internal::MyNormal(vector)),
"Normalized vector failed test."); "Normalized vector failed test.");
std::cout << " Normalize" << std::endl;
VectorType normalizedVector = vector; VectorType normalizedVector = vector;
vtkm::Normalize(normalizedVector); vtkm::Normalize(normalizedVector);
VTKM_TEST_ASSERT(test_equal(normalizedVector, internal::MyNormal(vector)), VTKM_TEST_ASSERT(test_equal(normalizedVector, internal::MyNormal(vector)),
...@@ -98,13 +91,11 @@ void TestLerp(const VectorType& a, ...@@ -98,13 +91,11 @@ void TestLerp(const VectorType& a,
const VectorType& w, const VectorType& w,
const typename vtkm::VecTraits<VectorType>::ComponentType& wS) const typename vtkm::VecTraits<VectorType>::ComponentType& wS)
{ {
std::cout << "Linear interpolation: " << a << "-" << b << ": " << w << std::endl;
VectorType vtkmLerp = vtkm::Lerp(a, b, w); VectorType vtkmLerp = vtkm::Lerp(a, b, w);
VectorType otherLerp = internal::MyLerp(a, b, w); VectorType otherLerp = internal::MyLerp(a, b, w);
VTKM_TEST_ASSERT(test_equal(vtkmLerp, otherLerp), VTKM_TEST_ASSERT(test_equal(vtkmLerp, otherLerp),
"Vectors with Vector weight do not lerp() correctly"); "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 lhsS = internal::MyLerp(a, b, wS);
VectorType rhsS = vtkm::Lerp(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"); VTKM_TEST_ASSERT(test_equal(lhsS, rhsS), "Vectors with Scalar weight do not lerp() correctly");
...@@ -113,19 +104,16 @@ void TestLerp(const VectorType& a, ...@@ -113,19 +104,16 @@ void TestLerp(const VectorType& a,
template <typename T> template <typename T>
void TestCross(const vtkm::Vec<T, 3>& x, const vtkm::Vec<T, 3>& y) void TestCross(const vtkm::Vec<T, 3>& x, const vtkm::Vec<T, 3>& y)
{ {
std::cout << "Testing " << x << " x " << y << std::endl;
using Vec3 = vtkm::Vec<T, 3>; using Vec3 = vtkm::Vec<T, 3>;
Vec3 cross = vtkm::Cross(x, y); 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. // 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(vtkm::Abs(vtkm::Dot(cross, x)) <
VTKM_TEST_ASSERT(test_equal(vtkm::Dot(cross, y), T(0.0)), "Cross product not perpendicular."); std::numeric_limits<T>::epsilon() * vtkm::MagnitudeSquared(x),
"Cross product not perpendicular.");
std::cout << " Length" << std::endl; VTKM_TEST_ASSERT(vtkm::Abs(vtkm::Dot(cross, y)) <
std::numeric_limits<T>::epsilon() * vtkm::MagnitudeSquared(y),
"Cross product not perpendicular.");
// The length of cross product should be the lengths of the input vectors // The length of cross product should be the lengths of the input vectors
// times the sin of the angle between them. // times the sin of the angle between them.
T sinAngle = vtkm::Magnitude(cross) * vtkm::RMagnitude(x) * vtkm::RMagnitude(y); T sinAngle = vtkm::Magnitude(cross) * vtkm::RMagnitude(x) * vtkm::RMagnitude(y);
...@@ -139,10 +127,10 @@ void TestCross(const vtkm::Vec<T, 3>& x, const vtkm::Vec<T, 3>& y) ...@@ -139,10 +127,10 @@ void TestCross(const vtkm::Vec<T, 3>& x, const vtkm::Vec<T, 3>& y)
VTKM_TEST_ASSERT(test_equal(sinAngle * sinAngle + cosAngle * cosAngle, T(1.0)), VTKM_TEST_ASSERT(test_equal(sinAngle * sinAngle + cosAngle * cosAngle, T(1.0)),
"Bad cross product length."); "Bad cross product length.");
std::cout << " Triangle normal" << std::endl;
// Test finding the normal to a triangle (similar to cross product). // Test finding the normal to a triangle (similar to cross product).
Vec3 normal = vtkm::TriangleNormal(x, y, Vec3(0, 0, 0)); 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(vtkm::Abs(vtkm::Dot(normal, x - y)) <
std::numeric_limits<T>::epsilon() * vtkm::MagnitudeSquared(x),
"Triangle normal is not really normal."); "Triangle normal is not really normal.");
} }
...@@ -151,13 +139,6 @@ void TestOrthonormalize(const VectorBasisType& inputs, int expectedRank) ...@@ -151,13 +139,6 @@ void TestOrthonormalize(const VectorBasisType& inputs, int expectedRank)
{ {
VectorBasisType outputs; VectorBasisType outputs;
int actualRank = vtkm::Orthonormalize(inputs, 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."); VTKM_TEST_ASSERT(test_equal(actualRank, expectedRank), "Orthonormalized rank is unexpected.");
} }
...@@ -208,6 +189,9 @@ struct TestCrossFunctor ...@@ -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, 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(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)); 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.969f, -30438.8f, -5643.727f));
} }
}; };
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment