Merge pull request arborx#1128 from aprokop/sym_svd

Author: aprokop

src/interpolation/detail/ArborX_InterpDetailsMovingLeastSquaresCoefficients.hpp

@@ -17,8 +17,8 @@
 #include <detail/ArborX_AccessTraits.hpp>
 #include <detail/ArborX_InterpDetailsCompactRadialBasisFunction.hpp>
 #include <detail/ArborX_InterpDetailsPolynomialBasis.hpp>
-#include <detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp>
 #include <kokkos_ext/ArborX_KokkosExtAccessibilityTraits.hpp>
+#include <misc/ArborX_SymmetricSVD.hpp>
 
 #include <Kokkos_Core.hpp>
 #include <Kokkos_Profiling_ScopedRegion.hpp>
@@ -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

src/interpolation/detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp src/misc/ArborX_SymmetricSVD.hpp

@@ -9,16 +9,16 @@
  * 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 <kokkos_ext/ArborX_KokkosExtAccessibilityTraits.hpp>
 #include <misc/ArborX_Exception.hpp>
 
 #include <Kokkos_Core.hpp>
 #include <Kokkos_Profiling_ScopedRegion.hpp>
 
-namespace ArborX::Interpolation::Details
+namespace ArborX::Details
 {
 
 template <typename Matrix>
@@ -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 <typename Matrix, typename Diag, typename Unit>
-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<typename Matrix::value_type>,
@@ -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 <typename Matrix, typename Diag, typename Unit>
+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<float>;
+
   // 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 <typename ExecutionSpace, typename InOutMatrices>
-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<InOutMatrices>, "matrices must be a view");
-  static_assert(!std::is_const_v<typename InOutMatrices::value_type>,
-                "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<Value **, MemorySpace> diags(
-      Kokkos::view_alloc(space, Kokkos::WithoutInitializing,
-                         "ArborX::SymmetricPseudoInverseSVD::diags"),
-      matrices.extent(0), matrices.extent(1));
-  Kokkos::View<Value ***, MemorySpace> 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

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

test/tstInterpDetailsSVD.cpp test/tstDetailsSVD.cpp

@@ -11,33 +11,60 @@
 
 #include "ArborX_EnableDeviceTypes.hpp"
 #include "ArborX_EnableViewComparison.hpp"
-#include <detail/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp>
+#include <misc/ArborX_SymmetricSVD.hpp>
 
 #include "BoostTest_CUDA_clang_workarounds.hpp"
 #include <boost/test/unit_test.hpp>
 
-template <typename MS, typename ES, typename V, int M, int N>
-void makeCase(ES const &es, V const (&src_arr)[M][N][N],
-              V const (&ref_arr)[M][N][N])
+template <typename MemorySpace, typename ExecutionSpace, typename Value, int M,
+          int N>
+void makeCase(ExecutionSpace const &exec, Value const (&src_arr)[M][N][N],
+              Value const (&ref_arr)[M][N][N])
 {
-  using DeviceView = Kokkos::View<V[M][N][N], MS>;
+  using DeviceView = Kokkos::View<Value[N][N], MemorySpace>;
   using HostView = typename DeviceView::HostMirror;
 
   HostView src("Testing::src");
   HostView ref("Testing::ref");
+
   DeviceView inv("Testing::inv");
+  DeviceView unit("Testing::unit");
+  Kokkos::View<Value[N], MemorySpace> 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<float>);
+    // Test SVD
+    Kokkos::deep_copy(exec, inv, src);
+    Kokkos::parallel_for(
+        "Testing::run_svd", Kokkos::RangePolicy<ExecutionSpace>(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<float>);
+
+    // Test pseudo-inverse through SVD
+    Kokkos::deep_copy(exec, inv, src);
+    Kokkos::parallel_for(
+        "Testing::run_inverse", Kokkos::RangePolicy<ExecutionSpace>(exec, 0, 1),
+        KOKKOS_LAMBDA(int) {
+          ArborX::Details::symmetricPseudoInverseSVDKernel(inv, diag, unit);
+        });
+    ARBORX_MDVIEW_TEST_TOL(ref, inv, Kokkos::Experimental::epsilon_v<float>);
+  }
 }
 
 // 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<MemorySpace>(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<double ***, MemorySpace> mat("mat", 0, 0, 0);
-  ArborX::Interpolation::Details::symmetricPseudoInverseSVD(space, mat);
-  BOOST_TEST(mat.size() == 0);
-}