diff --git a/examples/moving_least_squares/CMakeLists.txt b/examples/moving_least_squares/CMakeLists.txt index 9f3ec8d3f..a22bdbe30 100644 --- a/examples/moving_least_squares/CMakeLists.txt +++ b/examples/moving_least_squares/CMakeLists.txt @@ -1,3 +1,3 @@ add_executable(ArborX_Example_MovingLeastSquares.exe moving_least_squares.cpp) target_link_libraries(ArborX_Example_MovingLeastSquares.exe ArborX::ArborX) -add_test(NAME ArborX_Example_MovingLeastSquares COMMAND ArborX_Example_MovingLeastSquares.exe) \ No newline at end of file +add_test(NAME ArborX_Example_MovingLeastSquares COMMAND ArborX_Example_MovingLeastSquares.exe) diff --git a/examples/moving_least_squares/moving_least_squares.cpp b/examples/moving_least_squares/moving_least_squares.cpp index ef7cf161e..dcc15cc12 100644 --- a/examples/moving_least_squares/moving_least_squares.cpp +++ b/examples/moving_least_squares/moving_least_squares.cpp @@ -110,4 +110,4 @@ int main(int argc, char *argv[]) << diff(2) << '\n'; return 0; -} \ No newline at end of file +} diff --git a/src/geometry/ArborX_GeometryTraits.hpp b/src/geometry/ArborX_GeometryTraits.hpp index 6315e5e3a..7c13b9bb7 100644 --- a/src/geometry/ArborX_GeometryTraits.hpp +++ b/src/geometry/ArborX_GeometryTraits.hpp @@ -90,7 +90,7 @@ struct is_triangle : std::is_same::type, TriangleTag> {}; template -KOKKOS_FUNCTION constexpr void check_valid_geometry_traits(Geometry const &) +void check_valid_geometry_traits(Geometry const &) { static_assert( !Kokkos::is_detected{}, diff --git a/src/interpolation/ArborX_InterpMovingLeastSquares.hpp b/src/interpolation/ArborX_InterpMovingLeastSquares.hpp index d540d675a..c6f418301 100644 --- a/src/interpolation/ArborX_InterpMovingLeastSquares.hpp +++ b/src/interpolation/ArborX_InterpMovingLeastSquares.hpp @@ -31,40 +31,58 @@ namespace ArborX::Interpolation::Details { // This is done to avoid a clash with another predicates access trait -template -struct MLSTargetPointsPredicateWrapper +template +struct MLSPredicateWrapper { - Points target_points; + TargetAccess target_access; int num_neighbors; }; +// Functor used in the tree query to create the 2D source view and indices +template +struct MLSSearchNeighborsCallback +{ + SourceView source_view; + IndicesView indices; + CounterView counter; + + using SourcePoint = typename SourceView::non_const_value_type; + + template + KOKKOS_FUNCTION void + operator()(Predicate const &predicate, + ArborX::PairValueIndex const &primitive) const + { + int const target = getData(predicate); + int const source = primitive.index; + auto count = Kokkos::atomic_fetch_add(&counter(target), 1); + indices(target, count) = source; + source_view(target, count) = primitive.value; + } +}; + } // namespace ArborX::Interpolation::Details namespace ArborX { -template -struct AccessTraits< - Interpolation::Details::MLSTargetPointsPredicateWrapper, - PredicatesTag> +template +struct AccessTraits, + PredicatesTag> { - KOKKOS_INLINE_FUNCTION static auto size( - Interpolation::Details::MLSTargetPointsPredicateWrapper const &tp) + using Self = Interpolation::Details::MLSPredicateWrapper; + + KOKKOS_FUNCTION static auto size(Self const &tp) { - return AccessTraits::size(tp.target_points); + return tp.target_access.size(); } - KOKKOS_INLINE_FUNCTION static auto - get(Interpolation::Details::MLSTargetPointsPredicateWrapper const &tp, - int const i) + KOKKOS_FUNCTION static auto get(Self const &tp, int const i) { - return nearest( - AccessTraits::get(tp.target_points, i), - tp.num_neighbors); + return attach(nearest(tp.target_access(i), tp.num_neighbors), i); } - using memory_space = - typename AccessTraits::memory_space; + using memory_space = typename TargetAccess::memory_space; }; } // namespace ArborX @@ -77,11 +95,11 @@ class MovingLeastSquares { public: template , + typename TargetPoints, typename CRBFunc = CRBF::Wendland<0>, typename PolynomialDegree = PolynomialDegree<2>> MovingLeastSquares(ExecutionSpace const &space, SourcePoints const &source_points, - TargetPoints const &target_points, CRBF = {}, + TargetPoints const &target_points, CRBFunc = {}, PolynomialDegree = {}, std::optional num_neighbors = std::nullopt) { @@ -95,105 +113,53 @@ class MovingLeastSquares // SourcePoints is an access trait of points ArborX::Details::check_valid_access_traits(PrimitivesTag{}, source_points); - using src_acc = AccessTraits; - static_assert(KokkosExt::is_accessible_from::value, - "Source points must be accessible from the execution space"); - using src_point = - typename ArborX::Details::AccessTraitsHelper::type; - GeometryTraits::check_valid_geometry_traits(src_point{}); - static_assert(GeometryTraits::is_point::value, + using SourceAccess = + ArborX::Details::AccessValues; + static_assert( + KokkosExt::is_accessible_from::value, + "Source points must be accessible from the execution space"); + using SourcePoint = typename SourceAccess::value_type; + GeometryTraits::check_valid_geometry_traits(SourcePoint{}); + static_assert(GeometryTraits::is_point::value, "Source points elements must be points"); - static constexpr int dimension = GeometryTraits::dimension_v; + static constexpr int dimension = GeometryTraits::dimension_v; // TargetPoints is an access trait of points ArborX::Details::check_valid_access_traits(PrimitivesTag{}, target_points); - using tgt_acc = AccessTraits; - static_assert(KokkosExt::is_accessible_from::value, - "Target points must be accessible from the execution space"); - using tgt_point = - typename ArborX::Details::AccessTraitsHelper::type; - GeometryTraits::check_valid_geometry_traits(tgt_point{}); - static_assert(GeometryTraits::is_point::value, + using TargetAccess = + ArborX::Details::AccessValues; + static_assert( + KokkosExt::is_accessible_from::value, + "Target points must be accessible from the execution space"); + using TargetPoint = typename TargetAccess::value_type; + GeometryTraits::check_valid_geometry_traits(TargetPoint{}); + static_assert(GeometryTraits::is_point::value, "Target points elements must be points"); - static_assert(dimension == GeometryTraits::dimension_v, + static_assert(dimension == GeometryTraits::dimension_v, "Target and source points must have the same dimension"); - int num_neighbors_val = + _num_neighbors = num_neighbors ? *num_neighbors : Details::polynomialBasisSize(); - int const num_targets = tgt_acc::size(target_points); - _source_size = src_acc::size(source_points); + TargetAccess target_access{target_points}; + SourceAccess source_access{source_points}; + + _num_targets = target_access.size(); + _source_size = source_access.size(); // There must be enough source points - KOKKOS_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size); + KOKKOS_ASSERT(0 < _num_neighbors && _num_neighbors <= _source_size); - // Organize the source points as a tree - BoundingVolumeHierarchy> - source_tree(space, ArborX::Experimental::attach_indices(source_points)); - - // Create the predicates - Details::MLSTargetPointsPredicateWrapper predicates{ - target_points, num_neighbors_val}; - - // Query the source - Kokkos::View indices( - "ArborX::MovingLeastSquares::indices", 0); - Kokkos::View offsets( - "ArborX::MovingLeastSquares::offsets", 0); - source_tree.query(space, predicates, - ArborX::Details::LegacyDefaultCallback{}, indices, - offsets); - - // Fill in the value indices object so values can be transferred from a 1D - // source data to a properly distributed 2D array for each target. - auto const source_view = fillValuesIndicesAndGetSourceView( - space, indices, offsets, num_targets, num_neighbors_val, source_points); + // Search for neighbors and get the arranged source points + auto source_view = searchNeighbors(space, source_access, target_access); // Compute the moving least squares coefficients - _coeffs = Details::movingLeastSquaresCoefficients< - CRBF, PolynomialDegree, FloatingCalculationType, MemorySpace>( - space, source_view, target_points); - } - - template - Kokkos::View>::type **, - MemorySpace> - fillValuesIndicesAndGetSourceView( - ExecutionSpace const &space, - Kokkos::View const &indices, - Kokkos::View const &offsets, int const num_targets, - int const num_neighbors, SourcePoints const &source_points) - { - auto guard = Kokkos::Profiling::ScopedRegion( - "ArborX::MovingLeastSquares::fillValuesIndicesAndGetSourceView"); - - using src_acc = AccessTraits; - using src_point = - typename ArborX::Details::AccessTraitsHelper::type; - - _values_indices = Kokkos::View( - Kokkos::view_alloc(Kokkos::WithoutInitializing, - "ArborX::MovingLeastSquares::values_indices"), - num_targets, num_neighbors); - Kokkos::View source_view( - Kokkos::view_alloc(Kokkos::WithoutInitializing, - "ArborX::MovingLeastSquares::source_view"), - num_targets, num_neighbors); - Kokkos::parallel_for( - "ArborX::MovingLeastSquares::values_indices_and_source_view_fill", - Kokkos::MDRangePolicy>( - space, {0, 0}, {num_targets, num_neighbors}), - KOKKOS_CLASS_LAMBDA(int const i, int const j) { - auto index = indices(offsets(i) + j); - _values_indices(i, j) = index; - source_view(i, j) = src_acc::get(source_points, index); - }); - - return source_view; + _coeffs = Details::movingLeastSquaresCoefficients( + space, source_view, target_access._values); } template (space, 0, num_targets), + Kokkos::RangePolicy(space, 0, _num_targets), KOKKOS_CLASS_LAMBDA(int const i) { - value_t tmp = 0; - for (int j = 0; j < num_neighbors; j++) - tmp += _coeffs(i, j) * source_values(_values_indices(i, j)); + Value tmp = 0; + for (int j = 0; j < _num_neighbors; j++) + tmp += _coeffs(i, j) * source_values(_indices(i, j)); approx_values(i) = tmp; }); } private: + template + auto searchNeighbors(ExecutionSpace const &space, + SourceAccess const &source_access, + TargetAccess const &target_access) + { + auto guard = Kokkos::Profiling::ScopedRegion( + "ArborX::MovingLeastSquares::searchNeighbors"); + + // Organize the source points as a tree + using SourcePoint = typename SourceAccess::value_type; + BoundingVolumeHierarchy> + source_tree(space, ArborX::Experimental::attach_indices(source_access)); + + // Create the predicates + Details::MLSPredicateWrapper predicates{target_access, + _num_neighbors}; + + // Create the callback + Kokkos::View source_view( + Kokkos::view_alloc(space, Kokkos::WithoutInitializing, + "ArborX::MovingLeastSquares::source_view"), + _num_targets, _num_neighbors); + _indices = Kokkos::View( + Kokkos::view_alloc(space, Kokkos::WithoutInitializing, + "ArborX::MovingLeastSquares::indices"), + _num_targets, _num_neighbors); + Kokkos::View counter( + "ArborX::MovingLeastSquares::counter", _num_targets); + Details::MLSSearchNeighborsCallback + callback{source_view, _indices, counter}; + + // Query the source tree + source_tree.query(space, predicates, callback); + + return source_view; + } + Kokkos::View _coeffs; - Kokkos::View _values_indices; + Kokkos::View _indices; + int _num_targets; + int _num_neighbors; int _source_size; }; diff --git a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp index 23d12bd14..019cb29e3 100644 --- a/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsCompactRadialBasisFunction.hpp @@ -12,6 +12,9 @@ #ifndef ARBORX_INTERP_DETAILS_COMPACT_RADIAL_BASIS_FUNCTION_HPP #define ARBORX_INTERP_DETAILS_COMPACT_RADIAL_BASIS_FUNCTION_HPP +#include +#include + #include #include @@ -40,20 +43,22 @@ namespace CRBF { #define CRBF_DECL(NAME) \ - template \ + template \ struct NAME; #define CRBF_DEF(NAME, N, FUNC) \ template <> \ struct NAME \ { \ - template \ + template \ KOKKOS_INLINE_FUNCTION static constexpr T evaluate(T const y) \ { \ /* We force the input to be between 0 and 1. \ Because CRBF(-a) = CRBF(a) = CRBF(|a|), we take the absolute value \ and clamp the range to [0, 1] before entering in the definition of \ - the CRBF. */ \ + the CRBF. \ + We also template the internal function on the dimension as CRBFs \ + depend on the point's dimensionality. */ \ T const x = Kokkos::min(Kokkos::abs(y), T(1)); \ return Kokkos::abs(FUNC); \ } \ @@ -90,8 +95,20 @@ CRBF_DEF(Buhmann, 4, #undef CRBF_DEF #undef CRBF_DECL +template +KOKKOS_INLINE_FUNCTION constexpr auto +evaluate(Point const &point, + typename GeometryTraits::coordinate_type::type const radius) +{ + static_assert(GeometryTraits::is_point::value, + "point must be a point"); + constexpr std::size_t dim = GeometryTraits::dimension_v; + return CRBFunc::template evaluate( + ArborX::Details::distance(point, Point{}) / radius); +} + } // namespace CRBF } // namespace ArborX::Interpolation -#endif \ No newline at end of file +#endif diff --git a/src/interpolation/details/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp b/src/interpolation/details/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp index 35fd048fa..8261520e7 100644 --- a/src/interpolation/details/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -25,19 +26,20 @@ namespace ArborX::Interpolation::Details { -template -Kokkos::View +template +Kokkos::View movingLeastSquaresCoefficients(ExecutionSpace const &space, SourcePoints const &source_points, - TargetPoints const &target_points) + TargetAccess const &target_access) { auto guard = Kokkos::Profiling::ScopedRegion("ArborX::MovingLeastSquaresCoefficients"); namespace KokkosExt = ::ArborX::Details::KokkosExt; + using MemorySpace = typename SourcePoints::memory_space; static_assert( KokkosExt::is_accessible_from::value, "Memory space must be accessible from the execution space"); @@ -49,32 +51,31 @@ movingLeastSquaresCoefficients(ExecutionSpace const &space, KokkosExt::is_accessible_from::value, "source points must be accessible from the execution space"); - using src_point = typename SourcePoints::non_const_value_type; - GeometryTraits::check_valid_geometry_traits(src_point{}); - static_assert(GeometryTraits::is_point::value, + using SourcePoint = typename SourcePoints::non_const_value_type; + GeometryTraits::check_valid_geometry_traits(SourcePoint{}); + static_assert(GeometryTraits::is_point::value, "source points elements must be points"); - static constexpr int dimension = GeometryTraits::dimension_v; - - // TargetPoints is an access trait of points - ArborX::Details::check_valid_access_traits(PrimitivesTag{}, target_points); - using tgt_acc = AccessTraits; - static_assert(KokkosExt::is_accessible_from::value, - "target points must be accessible from the execution space"); - using tgt_point = typename ArborX::Details::AccessTraitsHelper::type; - GeometryTraits::check_valid_geometry_traits(tgt_point{}); - static_assert(GeometryTraits::is_point::value, - "target points elements must be points"); - static_assert(dimension == GeometryTraits::dimension_v, + static constexpr int dimension = GeometryTraits::dimension_v; + + // TargetAccess is an access values of points + static_assert( + KokkosExt::is_accessible_from::value, + "target access must be accessible from the execution space"); + using TargetPoint = typename TargetAccess::value_type; + GeometryTraits::check_valid_geometry_traits(TargetPoint{}); + static_assert(GeometryTraits::is_point::value, + "target access elements must be points"); + static_assert(dimension == GeometryTraits::dimension_v, "target and source points must have the same dimension"); - int const num_targets = tgt_acc::size(target_points); + int const num_targets = target_access.size(); int const num_neighbors = source_points.extent(1); // There must be a set of neighbors for each target KOKKOS_ASSERT(num_targets == source_points.extent_int(0)); - using point_t = ExperimentalHyperGeometry::Point; + using Point = ExperimentalHyperGeometry::Point; static constexpr auto epsilon = Kokkos::Experimental::epsilon_v; static constexpr int degree = PolynomialDegree::value; @@ -93,7 +94,7 @@ movingLeastSquaresCoefficients(ExecutionSpace const &space, // We first change the origin of the evaluation to be at the target point. // This lets us use p(0) which is [1 0 ... 0]. - Kokkos::View source_ref_target( + Kokkos::View source_ref_target( Kokkos::view_alloc( space, Kokkos::WithoutInitializing, "ArborX::MovingLeastSquaresCoefficients::source_ref_target"), @@ -104,8 +105,8 @@ movingLeastSquaresCoefficients(ExecutionSpace const &space, space, {0, 0}, {num_targets, num_neighbors}), KOKKOS_LAMBDA(int const i, int const j) { auto src = source_points(i, j); - auto tgt = tgt_acc::get(target_points, i); - point_t t{}; + auto tgt = target_access(i); + Point t{}; for (int k = 0; k < dimension; k++) t[k] = src[k] - tgt[k]; @@ -132,7 +133,7 @@ movingLeastSquaresCoefficients(ExecutionSpace const &space, for (int j = 0; j < num_neighbors; j++) { CoefficientsType norm = - ArborX::Details::distance(source_ref_target(i, j), point_t{}); + ArborX::Details::distance(source_ref_target(i, j), Point{}); radius = Kokkos::max(radius, norm); } @@ -154,9 +155,7 @@ movingLeastSquaresCoefficients(ExecutionSpace const &space, Kokkos::MDRangePolicy>( space, {0, 0}, {num_targets, num_neighbors}), KOKKOS_LAMBDA(int const i, int const j) { - CoefficientsType norm = - ArborX::Details::distance(source_ref_target(i, j), point_t{}); - phi(i, j) = CRBF::evaluate(norm / radii(i)); + phi(i, j) = CRBF::evaluate(source_ref_target(i, j), radii(i)); }); Kokkos::Profiling::popRegion(); diff --git a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp index 84776973f..2efed72f2 100644 --- a/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp @@ -107,18 +107,17 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) static_assert(GeometryTraits::is_point::value, "point must be a point"); static constexpr std::size_t DIM = GeometryTraits::dimension_v; - using value_t = typename GeometryTraits::coordinate_type::type; + using Value = typename GeometryTraits::coordinate_type::type; static_assert(DIM > 0, "Polynomial basis with no dimension is invalid"); - Kokkos::Array()> arr{}; - arr[0] = value_t(1); + Kokkos::Array()> arr{}; + arr[0] = Value(1); if constexpr (Degree > 0) { // Cannot use structured binding with constexpr - static constexpr auto slice_lengths_struct = + static constexpr auto slice_lengths = polynomialBasisSliceLengths(); - auto &slice_lengths = slice_lengths_struct.arr; std::size_t prev_col = 0; std::size_t curr_col = 1; @@ -129,10 +128,10 @@ KOKKOS_FUNCTION auto evaluatePolynomialBasis(Point const &p) for (std::size_t dim = 0; dim < DIM; dim++) { // copy the previous column and multply by p[dim] - for (std::size_t i = 0; i < slice_lengths[deg][dim]; i++) + for (std::size_t i = 0; i < slice_lengths.arr[deg][dim]; i++) arr[loc_offset + i] = arr[prev_col + i] * p[dim]; - loc_offset += slice_lengths[deg][dim]; + loc_offset += slice_lengths.arr[deg][dim]; } prev_col = curr_col; diff --git a/src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp b/src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp index 47a1b0cd0..029408001 100644 --- a/src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp +++ b/src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp @@ -52,11 +52,11 @@ template KOKKOS_FUNCTION auto argmaxUpperTriangle(Matrix const &mat) { ensureIsSquareMatrix(mat); - using value_t = typename Matrix::non_const_value_type; + using Value = typename Matrix::non_const_value_type; struct { - value_t max = 0; + Value max = 0; int row = 0; int col = 0; } result; @@ -65,7 +65,7 @@ KOKKOS_FUNCTION auto argmaxUpperTriangle(Matrix const &mat) for (int i = 0; i < size; i++) for (int j = i + 1; j < size; j++) { - value_t val = Kokkos::abs(mat(i, j)); + Value val = Kokkos::abs(mat(i, j)); if (result.max < val) { result.max = val; @@ -101,18 +101,18 @@ symmetricPseudoInverseSVDSerialKernel(AMatrix &A, ESMatrix &ES, UMatrix &U) typename UMatrix::value_type>, "All input matrices must have the same value type"); KOKKOS_ASSERT(A.extent(0) == ES.extent(0) && ES.extent(0) == U.extent(0)); - using value_t = typename AMatrix::non_const_value_type; + using Value = typename AMatrix::non_const_value_type; int const size = A.extent(0); // We first initialize U as the identity matrix and copy A to ES for (int i = 0; i < size; i++) for (int j = 0; j < size; j++) { - U(i, j) = value_t(i == j); + U(i, j) = Value(i == j); ES(i, j) = A(i, j); } - static constexpr value_t epsilon = Kokkos::Experimental::epsilon_v; + static constexpr Value epsilon = Kokkos::Experimental::epsilon_v; while (true) { // We have a guarantee that p < q @@ -138,13 +138,13 @@ symmetricPseudoInverseSVDSerialKernel(AMatrix &A, ESMatrix &ES, UMatrix &U) // | b | c | | 0 | y | // +---+---+ +---+---+ - value_t cos_theta; - value_t sin_theta; - value_t x; - value_t y; + Value cos_theta; + Value sin_theta; + Value x; + Value y; if (a == c) { - cos_theta = Kokkos::sqrt(value_t(2)) / 2; + cos_theta = Kokkos::sqrt(Value(2)) / 2; sin_theta = cos_theta; x = a + b; y = a - b; @@ -212,7 +212,7 @@ symmetricPseudoInverseSVDSerialKernel(AMatrix &A, ESMatrix &ES, UMatrix &U) for (int i = 0; i < size; i++) for (int j = 0; j < size; j++) { - value_t tmp = 0; + Value tmp = 0; for (int k = 0; k < size; k++) tmp += ES(k, k) * U(i, k) * U(j, k); A(i, j) = tmp; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4e6cf2016..b1cfa4c18 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -253,7 +253,7 @@ add_test(NAME ArborX_Test_BoostAdapters COMMAND ArborX_Test_BoostAdapters.exe) add_executable(ArborX_Test_InterpMovingLeastSquares.exe tstInterpDetailsSVD.cpp - tstInterpDetailsCRBF.cpp + tstInterpDetailsCompactRadialBasisFunction.cpp tstInterpDetailsPolyBasis.cpp tstInterpDetailsMLSCoefficients.cpp tstInterpMovingLeastSquares.cpp diff --git a/test/tstInterpDetailsCRBF.cpp b/test/tstInterpDetailsCompactRadialBasisFunction.cpp similarity index 90% rename from test/tstInterpDetailsCRBF.cpp rename to test/tstInterpDetailsCompactRadialBasisFunction.cpp index 39d1ed457..d17419bdb 100644 --- a/test/tstInterpDetailsCRBF.cpp +++ b/test/tstInterpDetailsCompactRadialBasisFunction.cpp @@ -11,6 +11,8 @@ #include "ArborX_EnableDeviceTypes.hpp" #include "ArborX_EnableViewComparison.hpp" +#include +#include #include #include "BoostTest_CUDA_clang_workarounds.hpp" @@ -20,6 +22,7 @@ template void makeCase(ES const &es, std::function const &tf, T tol = 1e-5) { + using Point = ArborX::ExperimentalHyperGeometry::Point<1, T>; using View = Kokkos::View; using HostView = typename View::HostMirror; static constexpr int range = 15; @@ -37,7 +40,10 @@ void makeCase(ES const &es, std::function const &tf, T tol = 1e-5) Kokkos::deep_copy(es, eval, input); Kokkos::parallel_for( "Testing::eval_crbf", Kokkos::RangePolicy(es, 0, 4 * range), - KOKKOS_LAMBDA(int const i) { eval(i) = CRBF::evaluate(eval(i)); }); + KOKKOS_LAMBDA(int const i) { + eval(i) = + ArborX::Interpolation::CRBF::evaluate(Point{eval(i)}, 1); + }); if (bool(tf)) { @@ -59,8 +65,8 @@ void makeCase(ES const &es, std::function const &tf, T tol = 1e-5) } #define MAKE_TEST(F, I, TF) \ - BOOST_AUTO_TEST_CASE_TEMPLATE(crbf_##F##_##I, DeviceType, \ - ARBORX_DEVICE_TYPES) \ + BOOST_AUTO_TEST_CASE_TEMPLATE(compact_radial_basis_function_##F##_##I, \ + DeviceType, ARBORX_DEVICE_TYPES) \ { \ makeCase>( \ typename DeviceType::execution_space{}, TF); \ @@ -94,4 +100,4 @@ MAKE_TEST_NONE(Buhmann, 4) #undef MAKE_TEST_NONE #undef MAKE_TEST_POLY -#undef MAKE_TEST \ No newline at end of file +#undef MAKE_TEST diff --git a/test/tstInterpDetailsMLSCoefficients.cpp b/test/tstInterpDetailsMLSCoefficients.cpp index ee2ff5473..c94ec70f6 100644 --- a/test/tstInterpDetailsMLSCoefficients.cpp +++ b/test/tstInterpDetailsMLSCoefficients.cpp @@ -51,9 +51,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) // -------0---------------> // SRC: 0 2 4 6 // TGT: 1 3 5 - using point0 = ArborX::ExperimentalHyperGeometry::Point<1, double>; - Kokkos::View srcp0("Testing::srcp0", 3, 2); - Kokkos::View tgtp0("Testing::tgtp0", 3); + using Point0 = ArborX::ExperimentalHyperGeometry::Point<1, double>; + Kokkos::View srcp0("Testing::srcp0", 3, 2); + Kokkos::View tgtp0("Testing::tgtp0", 3); Kokkos::View srcv0("Testing::srcv0", 3, 2); Kokkos::View tgtv0("Testing::tgtv0", 3); Kokkos::parallel_for( @@ -64,7 +64,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) srcp0(i, 1) = {{2. * i + 2}}; tgtp0(i) = {{2. * i + 1}}; - auto f = [](const point0 &) { return 3.; }; + auto f = [](const Point0 &) { return 3.; }; srcv0(i, 0) = f(srcp0(i, 0)); srcv0(i, 1) = f(srcp0(i, 1)); @@ -72,8 +72,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) }); auto coeffs0 = ArborX::Interpolation::Details::movingLeastSquaresCoefficients< ArborX::Interpolation::CRBF::Wendland<0>, - ArborX::Interpolation::PolynomialDegree<1>, double, MemorySpace>( - space, srcp0, tgtp0); + ArborX::Interpolation::PolynomialDegree<1>, double>(space, srcp0, tgtp0); auto eval0 = interpolate(space, srcv0, coeffs0); ARBORX_MDVIEW_TEST_TOL(eval0, tgtv0, Kokkos::Experimental::epsilon_v); @@ -86,9 +85,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) // T | T // S S S // | - using point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp1("Testing::srcp1", 4, 8); - Kokkos::View tgtp1("Testing::tgtp1", 4); + using Point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp1("Testing::srcp1", 4, 8); + Kokkos::View tgtp1("Testing::tgtp1", 4); Kokkos::View srcv1("Testing::srcv1", 4, 8); Kokkos::View tgtv1("Testing::tgtv1", 4); Kokkos::parallel_for( @@ -109,7 +108,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) } tgtp1(i) = {{double(u), double(v)}}; - auto f = [](const point1 &p) { return p[0] * p[1] + 4 * p[0]; }; + auto f = [](const Point1 &p) { return p[0] * p[1] + 4 * p[0]; }; for (int j = 0; j < 8; j++) srcv1(i, j) = f(srcp1(i, j)); @@ -117,8 +116,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients, DeviceType, ARBORX_DEVICE_TYPES) }); auto coeffs1 = ArborX::Interpolation::Details::movingLeastSquaresCoefficients< ArborX::Interpolation::CRBF::Wendland<2>, - ArborX::Interpolation::PolynomialDegree<2>, double, MemorySpace>( - space, srcp1, tgtp1); + ArborX::Interpolation::PolynomialDegree<2>, double>(space, srcp1, tgtp1); auto eval1 = interpolate(space, srcv1, coeffs1); ARBORX_MDVIEW_TEST_TOL(eval1, tgtv1, Kokkos::Experimental::epsilon_v); } @@ -131,9 +129,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients_edge_cases, DeviceType, ExecutionSpace space{}; // Case 1: Same as previous case 1, but points are 2D and locked on y=0 - using point0 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp0("Testing::srcp0", 3, 2); - Kokkos::View tgtp0("Testing::tgtp0", 3); + using Point0 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp0("Testing::srcp0", 3, 2); + Kokkos::View tgtp0("Testing::tgtp0", 3); Kokkos::View srcv0("Testing::srcv0", 3, 2); Kokkos::View tgtv0("Testing::tgtv0", 3); Kokkos::parallel_for( @@ -144,7 +142,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients_edge_cases, DeviceType, srcp0(i, 1) = {{2. * i + 2, 0.}}; tgtp0(i) = {{2. * i + 1, 0.}}; - auto f = [](const point0 &) { return 3.; }; + auto f = [](const Point0 &) { return 3.; }; srcv0(i, 0) = f(srcp0(i, 0)); srcv0(i, 1) = f(srcp0(i, 1)); @@ -152,15 +150,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients_edge_cases, DeviceType, }); auto coeffs0 = ArborX::Interpolation::Details::movingLeastSquaresCoefficients< ArborX::Interpolation::CRBF::Wendland<0>, - ArborX::Interpolation::PolynomialDegree<1>, double, MemorySpace>( - space, srcp0, tgtp0); + ArborX::Interpolation::PolynomialDegree<1>, double>(space, srcp0, tgtp0); auto eval0 = interpolate(space, srcv0, coeffs0); ARBORX_MDVIEW_TEST_TOL(eval0, tgtv0, Kokkos::Experimental::epsilon_v); // Case 2: Same but corner source points are also targets - using point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp1("Testing::srcp1", 4, 8); - Kokkos::View tgtp1("Testing::tgtp1", 4); + using Point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp1("Testing::srcp1", 4, 8); + Kokkos::View tgtp1("Testing::tgtp1", 4); Kokkos::View srcv1("Testing::srcv1", 4, 8); Kokkos::View tgtv1("Testing::tgtv1", 4); Kokkos::parallel_for( @@ -181,7 +178,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients_edge_cases, DeviceType, } tgtp1(i) = {{u * 2., v * 2.}}; - auto f = [](const point1 &p) { return p[0] * p[1] + 4 * p[0]; }; + auto f = [](const Point1 &p) { return p[0] * p[1] + 4 * p[0]; }; for (int j = 0; j < 8; j++) srcv1(i, j) = f(srcp1(i, j)); @@ -189,8 +186,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(mls_coefficients_edge_cases, DeviceType, }); auto coeffs1 = ArborX::Interpolation::Details::movingLeastSquaresCoefficients< ArborX::Interpolation::CRBF::Wendland<2>, - ArborX::Interpolation::PolynomialDegree<2>, double, MemorySpace>( - space, srcp1, tgtp1); + ArborX::Interpolation::PolynomialDegree<2>, double>(space, srcp1, tgtp1); auto eval1 = interpolate(space, srcv1, coeffs1); ARBORX_MDVIEW_TEST_TOL(eval1, tgtv1, Kokkos::Experimental::epsilon_v); -} \ No newline at end of file +} diff --git a/test/tstInterpDetailsPolyBasis.cpp b/test/tstInterpDetailsPolyBasis.cpp index 4aba902c1..22f9499f6 100644 --- a/test/tstInterpDetailsPolyBasis.cpp +++ b/test/tstInterpDetailsPolyBasis.cpp @@ -17,21 +17,21 @@ BOOST_AUTO_TEST_CASE(polynomial_basis_slice_lengths) { - using view = Kokkos::View; + using View = Kokkos::View; auto [arr0] = ArborX::Interpolation::Details::polynomialBasisSliceLengths<5, 3>(); std::size_t ref0[3][5] = { {1, 1, 1, 1, 1}, {1, 2, 3, 4, 5}, {1, 3, 6, 10, 15}}; - view arr0_view(&arr0[0][0], 3, 5); - view ref0_view(&ref0[0][0], 3, 5); + View arr0_view(&arr0[0][0], 3, 5); + View ref0_view(&ref0[0][0], 3, 5); ARBORX_MDVIEW_TEST(arr0_view, ref0_view); auto [arr1] = ArborX::Interpolation::Details::polynomialBasisSliceLengths<2, 3>(); std::size_t ref1[3][2] = {{1, 1}, {1, 2}, {1, 3}}; - view arr1_view(&arr1[0][0], 3, 2); - view ref1_view(&ref1[0][0], 3, 2); + View arr1_view(&arr1[0][0], 3, 2); + View ref1_view(&ref1[0][0], 3, 2); ARBORX_MDVIEW_TEST(arr1_view, ref1_view); } @@ -47,7 +47,7 @@ BOOST_AUTO_TEST_CASE(polynomial_basis_size) BOOST_AUTO_TEST_CASE(polynomial_basis) { - using view = Kokkos::View; + using View = Kokkos::View; ArborX::ExperimentalHyperGeometry::Point<5, double> point0 = {1, 2, 3, 4, 5}; auto arr0 = @@ -56,15 +56,15 @@ BOOST_AUTO_TEST_CASE(polynomial_basis) 12, 16, 5, 10, 15, 20, 25, 1, 2, 4, 8, 3, 6, 12, 9, 18, 27, 4, 8, 16, 12, 24, 36, 16, 32, 48, 64, 5, 10, 20, 15, 30, 45, 20, 40, 60, 80, 25, 50, 75, 100, 125}; - view arr0_view(arr0.data(), 56); - view ref0_view(&ref0[0], 56); + View arr0_view(arr0.data(), 56); + View ref0_view(&ref0[0], 56); ARBORX_MDVIEW_TEST(arr0_view, ref0_view); ArborX::ExperimentalHyperGeometry::Point<2, double> point1 = {-2, 0}; auto arr1 = ArborX::Interpolation::Details::evaluatePolynomialBasis<3>(point1); double ref1[10] = {1, -2, 0, 4, 0, 0, -8, 0, 0, 0}; - view arr1_view(arr1.data(), 10); - view ref1_view(&ref1[0], 10); + View arr1_view(arr1.data(), 10); + View ref1_view(&ref1[0], 10); ARBORX_MDVIEW_TEST(arr1_view, ref1_view); -} \ No newline at end of file +} diff --git a/test/tstInterpDetailsSVD.cpp b/test/tstInterpDetailsSVD.cpp index f0f56af1f..df8b26033 100644 --- a/test/tstInterpDetailsSVD.cpp +++ b/test/tstInterpDetailsSVD.cpp @@ -20,12 +20,12 @@ template void makeCase(ES const &es, V const (&src_arr)[M][N][N], V const (&ref_arr)[M][N][N]) { - using device_view = Kokkos::View; - using host_view = typename device_view::HostMirror; + using DeviceView = Kokkos::View; + using HostView = typename DeviceView::HostMirror; - host_view src("Testing::src"); - host_view ref("Testing::ref"); - device_view inv("Testing::inv"); + HostView src("Testing::src"); + HostView ref("Testing::ref"); + DeviceView inv("Testing::inv"); for (int i = 0; i < M; i++) for (int j = 0; j < N; j++) @@ -125,4 +125,4 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(pseudo_inv_empty, DeviceType, ARBORX_DEVICE_TYPES) Kokkos::View mat("mat", 0, 0, 0); ArborX::Interpolation::Details::symmetricPseudoInverseSVD(space, mat); BOOST_TEST(mat.size() == 0); -} \ No newline at end of file +} diff --git a/test/tstInterpMovingLeastSquares.cpp b/test/tstInterpMovingLeastSquares.cpp index ff9e4a767..ceee3e7b4 100644 --- a/test/tstInterpMovingLeastSquares.cpp +++ b/test/tstInterpMovingLeastSquares.cpp @@ -29,17 +29,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares, DeviceType, // -------0---------------> // SRC: 0 2 4 6 // TGT: 1 3 5 - using point0 = ArborX::ExperimentalHyperGeometry::Point<1, double>; - Kokkos::View srcp0("Testing::srcp0", 4); - Kokkos::View tgtp0("Testing::tgtp0", 3); + using Point0 = ArborX::ExperimentalHyperGeometry::Point<1, double>; + Kokkos::View srcp0("Testing::srcp0", 4); + Kokkos::View tgtp0("Testing::tgtp0", 3); Kokkos::View srcv0("Testing::srcv0", 4); Kokkos::View tgtv0("Testing::tgtv0", 3); - Kokkos::View eval0("Testing::eval0", 0); + Kokkos::View eval0("Testing::eval0", 3); Kokkos::parallel_for( "Testing::moving_least_squares::for0", Kokkos::RangePolicy(space, 0, 4), KOKKOS_LAMBDA(int const i) { - auto f = [](const point0 &) { return 3.; }; + auto f = [](const Point0 &) { return 3.; }; srcp0(i) = {{2. * i}}; srcv0(i) = f(srcp0(i)); @@ -65,12 +65,12 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares, DeviceType, // T | T // S S S // | - using point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp1("Testing::srcp1", 9); - Kokkos::View tgtp1("Testing::tgtp1", 4); + using Point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp1("Testing::srcp1", 9); + Kokkos::View tgtp1("Testing::tgtp1", 4); Kokkos::View srcv1("Testing::srcv1", 9); Kokkos::View tgtv1("Testing::tgtv1", 4); - Kokkos::View eval1("Testing::eval1", 0); + Kokkos::View eval1("Testing::eval1", 4); Kokkos::parallel_for( "Testing::moving_least_squares::for1", Kokkos::RangePolicy(space, 0, 9), @@ -79,7 +79,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares, DeviceType, int v = (i % 2) * 2 - 1; int x = (i / 3) - 1; int y = (i % 3) - 1; - auto f = [](const point1 &p) { return p[0] * p[1] + 4 * p[0]; }; + auto f = [](const Point1 &p) { return p[0] * p[1] + 4 * p[0]; }; srcp1(i) = {{x * 2., y * 2.}}; srcv1(i) = f(srcp1(i)); @@ -104,17 +104,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares_edge_cases, DeviceType, ExecutionSpace space{}; // Case 1: Same as previous case 1, but points are 2D and locked on y=0 - using point0 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp0("Testing::srcp0", 4); - Kokkos::View tgtp0("Testing::tgtp0", 3); + using Point0 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp0("Testing::srcp0", 4); + Kokkos::View tgtp0("Testing::tgtp0", 3); Kokkos::View srcv0("Testing::srcv0", 4); Kokkos::View tgtv0("Testing::tgtv0", 3); - Kokkos::View eval0("Testing::eval0", 0); + Kokkos::View eval0("Testing::eval0", 3); Kokkos::parallel_for( "Testing::moving_least_squares_edge_cases::for0", Kokkos::RangePolicy(space, 0, 4), KOKKOS_LAMBDA(int const i) { - auto f = [](const point0 &) { return 3.; }; + auto f = [](const Point0 &) { return 3.; }; srcp0(i) = {{2. * i, 0.}}; srcv0(i) = f(srcp0(i)); @@ -132,12 +132,12 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares_edge_cases, DeviceType, ARBORX_MDVIEW_TEST_TOL(eval0, tgtv0, Kokkos::Experimental::epsilon_v); // Case 2: Same but corner source points are also targets - using point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; - Kokkos::View srcp1("Testing::srcp1", 9); - Kokkos::View tgtp1("Testing::tgtp1", 4); + using Point1 = ArborX::ExperimentalHyperGeometry::Point<2, double>; + Kokkos::View srcp1("Testing::srcp1", 9); + Kokkos::View tgtp1("Testing::tgtp1", 4); Kokkos::View srcv1("Testing::srcv1", 9); Kokkos::View tgtv1("Testing::tgtv1", 4); - Kokkos::View eval1("Testing::eval1", 0); + Kokkos::View eval1("Testing::eval1", 4); Kokkos::parallel_for( "Testing::moving_least_squares_edge_cases::for0", Kokkos::RangePolicy(space, 0, 9), @@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares_edge_cases, DeviceType, int v = (i % 2) * 2 - 1; int x = (i / 3) - 1; int y = (i % 3) - 1; - auto f = [](const point1 &p) { return p[0] * p[1] + 4 * p[0]; }; + auto f = [](const Point1 &p) { return p[0] * p[1] + 4 * p[0]; }; srcp1(i) = {{x * 2., y * 2.}}; srcv1(i) = f(srcp1(i)); @@ -161,4 +161,4 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(moving_least_squares_edge_cases, DeviceType, ArborX::Interpolation::PolynomialDegree<2>{}, 8); mls1.interpolate(space, srcv1, eval1); ARBORX_MDVIEW_TEST_TOL(eval1, tgtv1, Kokkos::Experimental::epsilon_v); -} \ No newline at end of file +}