Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Return a vec<T,2> rather than a vtkm::Pair<T,T>.

Browse Source
This commit is contained in:
Nick Thompson 2021-04-08 14:45:59 -04:00
parent 91bec19e97
commit dd8863637a
3 changed files with 22 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
vtkm/Math.h
View File

@ -2701,12 +2701,6 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32
template<typename T> template<typename T>
inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d)
{ {
// This is a bit awkward. clang-11 does not define FP_FAST_FMA, but still compiles this to the correct assembly.
// Windows, however, generates truly horrendous assembly from this, with no fmas, to the extent I assume it could
// contort itself into a performance bug.
// That said, MSVC converts even a*b - c*d into horrible assembly, so it may be a wash.
// You'd want to just use #ifdef FP_FAST_FMA, but then you'd lose the (correct) generated assembly on clang.
// See: https://stackoverflow.com/a/40765925/904050
T cd = c * d; T cd = c * d;
T err = std::fma(-c, d, cd); T err = std::fma(-c, d, cd);
T dop = std::fma(a, b, -cd); T dop = std::fma(a, b, -cd);
@ -2719,7 +2713,7 @@ inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d)
// If there are no real roots, both elements are NaNs. // If there are no real roots, both elements are NaNs.
// The error should be at most 3 ulps. // The error should be at most 3 ulps.
template<typename T> template<typename T>
inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c) inline VTKM_EXEC_CONT vtkm::Vec<T, 2> QuadraticRoots(T a, T b, T c)
{ {
if (a == 0) if (a == 0)
{ {
@ -2728,19 +2722,19 @@ inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c)
if (c == 0) if (c == 0)
{ {
// A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use. // A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use.
return vtkm::Pair<T, T>(0,0); return vtkm::Vec<T,2>(0,0);
} }
else else
{ {
return vtkm::Pair<T, T>(vtkm::Nan<T>(), vtkm::Nan<T>()); return vtkm::Vec<T,2>(vtkm::Nan<T>(), vtkm::Nan<T>());
} }
} }
return vtkm::Pair<T, T>(-c/b, -c/b); return vtkm::Vec<T,2>(-c/b, -c/b);
} }
T delta = DifferenceOfProducts(b, b, 4*a, c); T delta = DifferenceOfProducts(b, b, 4*a, c);
if (delta < 0) if (delta < 0)
{ {
return vtkm::Pair<T, T>(vtkm::Nan<T>(), vtkm::Nan<T>()); return vtkm::Vec<T,2>(vtkm::Nan<T>(), vtkm::Nan<T>());
} }
T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2; T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2;
@ -2748,9 +2742,9 @@ inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c)
T r1 = c / q; T r1 = c / q;
if (r0 < r1) if (r0 < r1)
{ {
return vtkm::Pair<T, T>(r0, r1); return vtkm::Vec<T,2>(r0, r1);
} }
return vtkm::Pair<T, T>(r1, r0); return vtkm::Vec<T,2>(r1, r0);
} }
/// Bitwise operations /// Bitwise operations

20
vtkm/Math.h.in
View File

@ -1303,12 +1303,6 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32
template<typename T> template<typename T>
inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d)
{ {
// This is a bit awkward. clang-11 does not define FP_FAST_FMA, but still compiles this to the correct assembly.
// Windows, however, generates truly horrendous assembly from this, with no fmas, to the extent I assume it could
// contort itself into a performance bug.
// That said, MSVC converts even a*b - c*d into horrible assembly, so it may be a wash.
// You'd want to just use #ifdef FP_FAST_FMA, but then you'd lose the (correct) generated assembly on clang.
// See: https://stackoverflow.com/a/40765925/904050
T cd = c * d; T cd = c * d;
T err = std::fma(-c, d, cd); T err = std::fma(-c, d, cd);
T dop = std::fma(a, b, -cd); T dop = std::fma(a, b, -cd);
@ -1321,7 +1315,7 @@ inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d)
// If there are no real roots, both elements are NaNs. // If there are no real roots, both elements are NaNs.
// The error should be at most 3 ulps. // The error should be at most 3 ulps.
template<typename T> template<typename T>
inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c) inline VTKM_EXEC_CONT vtkm::Vec<T, 2> QuadraticRoots(T a, T b, T c)
{ {
if (a == 0) if (a == 0)
{ {
@ -1330,19 +1324,19 @@ inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c)
if (c == 0) if (c == 0)
{ {
// A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use. // A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use.
return vtkm::Pair<T, T>(0,0); return vtkm::Vec<T,2>(0,0);
} }
else else
{ {
return vtkm::Pair<T, T>(vtkm::Nan<T>(), vtkm::Nan<T>()); return vtkm::Vec<T,2>(vtkm::Nan<T>(), vtkm::Nan<T>());
} }
} }
return vtkm::Pair<T, T>(-c/b, -c/b); return vtkm::Vec<T,2>(-c/b, -c/b);
} }
T delta = DifferenceOfProducts(b, b, 4*a, c); T delta = DifferenceOfProducts(b, b, 4*a, c);
if (delta < 0) if (delta < 0)
{ {
return vtkm::Pair<T, T>(vtkm::Nan<T>(), vtkm::Nan<T>()); return vtkm::Vec<T,2>(vtkm::Nan<T>(), vtkm::Nan<T>());
} }
T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2; T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2;
@ -1350,9 +1344,9 @@ inline VTKM_EXEC_CONT vtkm::Pair<T, T> QuadraticRoots(T a, T b, T c)
T r1 = c / q; T r1 = c / q;
if (r0 < r1) if (r0 < r1)
{ {
return vtkm::Pair<T, T>(r0, r1); return vtkm::Vec<T,2>(r0, r1);
} }
return vtkm::Pair<T, T>(r1, r0); return vtkm::Vec<T,2>(r1, r0);
} }
/// Bitwise operations /// Bitwise operations

17
vtkm/testing/UnitTestMath.cxx
View File

@ -916,7 +916,6 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase
vtkm::Float32 expected = 5.376600027084351f; vtkm::Float32 expected = 5.376600027084351f;
vtkm::UInt64 dist = vtkm::FloatDistance(expected, computed); vtkm::UInt64 dist = vtkm::FloatDistance(expected, computed);
std::cout << "Dist = " << dist << "\n";
VTKM_MATH_ASSERT( VTKM_MATH_ASSERT(
dist < 2, dist < 2,
"Float distance for difference of products is which exceeds 1.5; this is in violation of a " "Float distance for difference of products is which exceeds 1.5; this is in violation of a "
@ -934,17 +933,17 @@ void TestQuadraticRoots() const
// (x-1)(x+1) = x² - 1: // (x-1)(x+1) = x² - 1:
auto roots = vtkm::QuadraticRoots(1.0f, 0.0f, -1.0f); auto roots = vtkm::QuadraticRoots(1.0f, 0.0f, -1.0f);
vtkm::UInt64 dist = vtkm::FloatDistance(-1.0f, roots.first); vtkm::UInt64 dist = vtkm::FloatDistance(-1.0f, roots[0]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
dist = vtkm::FloatDistance(1.0f, roots.second); dist = vtkm::FloatDistance(1.0f, roots[1]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
// No real roots: // No real roots:
roots = vtkm::QuadraticRoots(1.0f, 0.0f, 1.0f); roots = vtkm::QuadraticRoots(1.0f, 0.0f, 1.0f);
VTKM_MATH_ASSERT(vtkm::IsNan(roots.first), VTKM_MATH_ASSERT(vtkm::IsNan(roots[0]),
"Roots should be Nan for a quadratic with complex roots."); "Roots should be Nan for a quadratic with complex roots.");
VTKM_MATH_ASSERT(vtkm::IsNan(roots.second), VTKM_MATH_ASSERT(vtkm::IsNan(roots[1]),
"Roots should be Nan for a quadratic with complex roots."); "Roots should be Nan for a quadratic with complex roots.");
#ifdef FP_FAST_FMA #ifdef FP_FAST_FMA
@ -952,18 +951,18 @@ void TestQuadraticRoots() const
// x² + 200x - 0.000015 = 0 has roots // x² + 200x - 0.000015 = 0 has roots
// -200.000000075, 7.5e-8 // -200.000000075, 7.5e-8
roots = vtkm::QuadraticRoots(1.0f, 200.0f, -0.000015f); roots = vtkm::QuadraticRoots(1.0f, 200.0f, -0.000015f);
dist = vtkm::FloatDistance(-200.000000075f, roots.first); dist = vtkm::FloatDistance(-200.000000075f, roots[0]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
dist = vtkm::FloatDistance(7.5e-8f, roots.second); dist = vtkm::FloatDistance(7.5e-8f, roots[1]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
// Kahan's example: // Kahan's example:
auto roots64 = vtkm::QuadraticRoots(94906265.625, 94906267.000, 94906268.375); auto roots64 = vtkm::QuadraticRoots(94906265.625, 94906267.000, 94906268.375);
dist = vtkm::FloatDistance(1.0, roots64.first); dist = vtkm::FloatDistance(1.0, roots64[0]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
dist = vtkm::FloatDistance(1.000000028975958, roots64.second); dist = vtkm::FloatDistance(1.000000028975958, roots64[1]);
VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps.");
#endif #endif
} }