diff --git a/src/interpolation/detail/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp b/src/interpolation/detail/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp index 7596c50d1..cb2c541df 100644 --- a/src/interpolation/detail/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp +++ b/src/interpolation/detail/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp @@ -17,8 +17,8 @@ #include #include #include -#include #include +#include #include #include @@ -110,7 +110,8 @@ class MovingLeastSquaresCoefficientsKernel // We need the inverse of P^T.PHI.P, and because it is symmetric, we can use // the symmetric SVD algorithm to get it. - symmetricPseudoInverseSVDKernel(moment, svd_diag, svd_unit); + ::ArborX::Details::symmetricPseudoInverseSVDKernel(moment, svd_diag, + svd_unit); // Now, the moment has [P^T.PHI.P]^-1 // Finally, the result is produced by computing p(0).[P^T.PHI.P]^-1.P^T.PHI diff --git a/src/interpolation/detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp b/src/misc/ArborX_SymmetricSVD.hpp similarity index 76% rename from src/interpolation/detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp rename to src/misc/ArborX_SymmetricSVD.hpp index e277ee095..4c7184ce2 100644 --- a/src/interpolation/detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp +++ b/src/misc/ArborX_SymmetricSVD.hpp @@ -9,8 +9,8 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ -#ifndef ARBORX_INTERP_SYMMETRIC_PSEUDO_INVERSE_SVD_HPP -#define ARBORX_INTERP_SYMMETRIC_PSEUDO_INVERSE_SVD_HPP +#ifndef ARBORX_SYMMETRIC_SVD_HPP +#define ARBORX_SYMMETRIC_SVD_HPP #include #include @@ -18,7 +18,7 @@ #include #include -namespace ArborX::Interpolation::Details +namespace ArborX::Details { template @@ -77,18 +77,17 @@ KOKKOS_FUNCTION auto argmaxUpperTriangle(Matrix const &mat) return result; } -// Pseudo-inverse of symmetric matrices using SVD +// SVD of a symmetric matrix // We must find U, E (diagonal and positive) and V such that A = U.E.V^T // We also suppose, as the input, that A is symmetric, so U = SV where S is // a sign matrix (only 1 or -1 on the diagonal, 0 elsewhere). -// Thus A = U.ES.U^T and A^-1 = U.[ ES^-1 ].U^T +// Thus A = U.ES.U^T. // // mat <=> initial ES // diag <=> final ES // unit <=> U template -KOKKOS_FUNCTION void symmetricPseudoInverseSVDKernel(Matrix &mat, Diag &diag, - Unit &unit) +KOKKOS_FUNCTION void symmetricSVDKernel(Matrix &mat, Diag &diag, Unit &unit) { ensureIsSquareSymmetricMatrix(mat); static_assert(!std::is_const_v, @@ -203,13 +202,34 @@ KOKKOS_FUNCTION void symmetricPseudoInverseSVDKernel(Matrix &mat, Diag &diag, } } + for (int i = 0; i < size; i++) + diag(i) = mat(i, i); +} + +// Pseudo-inverse of symmetric matrices using SVD +// We must find U, E (diagonal and positive) and V such that A = U.E.V^T +// We also suppose, as the input, that A is symmetric, so U = SV where S is +// a sign matrix (only 1 or -1 on the diagonal, 0 elsewhere). +// Thus A = U.ES.U^T and A^-1 = U.[ ES^-1 ].U^T +// +// mat <=> initial ES +// diag <=> final ES +// unit <=> U +template +KOKKOS_FUNCTION void symmetricPseudoInverseSVDKernel(Matrix &mat, Diag &diag, + Unit &unit) +{ + symmetricSVDKernel(mat, diag, unit); + + int const size = mat.extent(0); + + using Value = typename Matrix::non_const_value_type; + constexpr Value epsilon = Kokkos::Experimental::epsilon_v; + // We compute the max to get a range of the invertible eigenvalues auto max_eigen = epsilon; for (int i = 0; i < size; i++) - { - diag(i) = mat(i, i); max_eigen = Kokkos::max(Kokkos::abs(diag(i)), max_eigen); - } auto const threshold = max_eigen * epsilon; // We invert the diagonal of 'mat', except if "0" is found @@ -226,49 +246,6 @@ KOKKOS_FUNCTION void symmetricPseudoInverseSVDKernel(Matrix &mat, Diag &diag, } } -template -void symmetricPseudoInverseSVD(ExecutionSpace const &space, - InOutMatrices &matrices) -{ - auto guard = - Kokkos::Profiling::ScopedRegion("ArborX::SymmetricPseudoInverseSVD"); - - // InOutMatrices is a list of square symmetric matrices (3D view) - static_assert(Kokkos::is_view_v, "matrices must be a view"); - static_assert(!std::is_const_v, - "matrices must be writable"); - static_assert(InOutMatrices::rank() == 3, - "matrices must be a list of square matrices"); - static_assert( - ArborX::Details::KokkosExt::is_accessible_from< - typename InOutMatrices::memory_space, ExecutionSpace>::value, - "matrices must be accessible from the execution space"); - - KOKKOS_ASSERT(matrices.extent(1) == matrices.extent(2)); // Must be square - - using Value = typename InOutMatrices::non_const_value_type; - using MemorySpace = typename InOutMatrices::memory_space; - - Kokkos::View diags( - Kokkos::view_alloc(space, Kokkos::WithoutInitializing, - "ArborX::SymmetricPseudoInverseSVD::diags"), - matrices.extent(0), matrices.extent(1)); - Kokkos::View units( - Kokkos::view_alloc(space, Kokkos::WithoutInitializing, - "ArborX::SymmetricPseudoInverseSVD::units"), - matrices.extent(0), matrices.extent(1), matrices.extent(2)); - - Kokkos::parallel_for( - "ArborX::SymmetricPseudoInverseSVD::computations", - Kokkos::RangePolicy(space, 0, matrices.extent(0)), - KOKKOS_LAMBDA(int const i) { - auto mat = Kokkos::subview(matrices, i, Kokkos::ALL, Kokkos::ALL); - auto diag = Kokkos::subview(diags, i, Kokkos::ALL); - auto unit = Kokkos::subview(units, i, Kokkos::ALL, Kokkos::ALL); - symmetricPseudoInverseSVDKernel(mat, diag, unit); - }); -} - -} // namespace ArborX::Interpolation::Details +} // namespace ArborX::Details #endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6fdde05a8..ab79b388d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -49,6 +49,7 @@ add_executable(ArborX_Test_DetailsUtils.exe tstAttachIndices.cpp tstDetailsVector.cpp tstDetailsUtils.cpp + tstDetailsSVD.cpp tstDetailsGeometryReducer.cpp utf_main.cpp ) @@ -255,7 +256,6 @@ target_link_libraries(ArborX_Test_BoostAdapters.exe PRIVATE ArborX Boost::unit_t add_test(NAME ArborX_Test_BoostAdapters COMMAND ArborX_Test_BoostAdapters.exe) add_executable(ArborX_Test_InterpMovingLeastSquares.exe - tstInterpDetailsSVD.cpp tstInterpDetailsCompactRadialBasisFunction.cpp tstInterpDetailsPolyBasis.cpp tstInterpDetailsMLSCoefficients.cpp diff --git a/test/tstInterpDetailsSVD.cpp b/test/tstDetailsSVD.cpp similarity index 73% rename from test/tstInterpDetailsSVD.cpp rename to test/tstDetailsSVD.cpp index 1883e4ac8..d95400baf 100644 --- a/test/tstInterpDetailsSVD.cpp +++ b/test/tstDetailsSVD.cpp @@ -11,33 +11,60 @@ #include "ArborX_EnableDeviceTypes.hpp" #include "ArborX_EnableViewComparison.hpp" -#include +#include #include "BoostTest_CUDA_clang_workarounds.hpp" #include -template -void makeCase(ES const &es, V const (&src_arr)[M][N][N], - V const (&ref_arr)[M][N][N]) +template +void makeCase(ExecutionSpace const &exec, Value const (&src_arr)[M][N][N], + Value const (&ref_arr)[M][N][N]) { - using DeviceView = Kokkos::View; + using DeviceView = Kokkos::View; using HostView = typename DeviceView::HostMirror; HostView src("Testing::src"); HostView ref("Testing::ref"); + DeviceView inv("Testing::inv"); + DeviceView unit("Testing::unit"); + Kokkos::View diag("Testing::diag"); for (int i = 0; i < M; i++) + { for (int j = 0; j < N; j++) for (int k = 0; k < N; k++) { - src(i, j, k) = src_arr[i][j][k]; - ref(i, j, k) = ref_arr[i][j][k]; + src(j, k) = src_arr[i][j][k]; + ref(j, k) = ref_arr[i][j][k]; } - Kokkos::deep_copy(es, inv, src); - ArborX::Interpolation::Details::symmetricPseudoInverseSVD(es, inv); - ARBORX_MDVIEW_TEST_TOL(ref, inv, Kokkos::Experimental::epsilon_v); + // Test SVD + Kokkos::deep_copy(exec, inv, src); + Kokkos::parallel_for( + "Testing::run_svd", Kokkos::RangePolicy(exec, 0, 1), + KOKKOS_LAMBDA(int) { + ArborX::Details::symmetricSVDKernel(inv, diag, unit); + for (int p = 0; p < N; ++p) + for (int q = 0; q < N; ++q) + { + inv(p, q) = 0; + for (int s = 0; s < N; ++s) + inv(p, q) += unit(p, s) * diag(s) * unit(q, s); + } + }); + ARBORX_MDVIEW_TEST_TOL(src, inv, Kokkos::Experimental::epsilon_v); + + // Test pseudo-inverse through SVD + Kokkos::deep_copy(exec, inv, src); + Kokkos::parallel_for( + "Testing::run_inverse", Kokkos::RangePolicy(exec, 0, 1), + KOKKOS_LAMBDA(int) { + ArborX::Details::symmetricPseudoInverseSVDKernel(inv, diag, unit); + }); + ARBORX_MDVIEW_TEST_TOL(ref, inv, Kokkos::Experimental::epsilon_v); + } } // Pseudo-inverses were computed using numpy's "linalg.pinv" solver and @@ -115,14 +142,3 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(pseudo_inv_scalar_like, DeviceType, double inv[2][1][1] = {{{1 / 2.}}, {{0}}}; makeCase(space, mat, inv); } - -BOOST_AUTO_TEST_CASE_TEMPLATE(pseudo_inv_empty, DeviceType, ARBORX_DEVICE_TYPES) -{ - using ExecutionSpace = typename DeviceType::execution_space; - using MemorySpace = typename DeviceType::memory_space; - ExecutionSpace space{}; - - Kokkos::View mat("mat", 0, 0, 0); - ArborX::Interpolation::Details::symmetricPseudoInverseSVD(space, mat); - BOOST_TEST(mat.size() == 0); -}