Provide SVD kernel independently of pseudo-inverse

src/interpolation/details/ArborX_InterpDetailsSymmetricPseudoInverseSVD.hpp

@@ -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

test/tstInterpDetailsSVD.cpp

@@ -39,10 +39,27 @@ void makeCase(ExecutionSpace const &exec, Value const (&src_arr)[M][N][N],
         src(j, k) = src_arr[i][j][k];
         ref(j, k) = ref_arr[i][j][k];
       }
+
+    // Test SVD
     Kokkos::deep_copy(exec, inv, src);
+    Kokkos::parallel_for(
+        "Testing::run_svd", Kokkos::RangePolicy<ExecutionSpace>(exec, 0, 1),
+        KOKKOS_LAMBDA(int) {
+          ArborX::Interpolation::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_case", Kokkos::RangePolicy<ExecutionSpace>(exec, 0, 1),
+        "Testing::run_inverse", Kokkos::RangePolicy<ExecutionSpace>(exec, 0, 1),
         KOKKOS_LAMBDA(int) {
           ArborX::Interpolation::Details::symmetricPseudoInverseSVDKernel(
               inv, diag, unit);