Add singular value wrappers to dense matrix functionality

awslabs · May 1, 2024 · c358bbc · c358bbc
1 parent a7c996d
commit c358bbc
Show file tree

Hide file tree

Showing 4 changed files with 48 additions and 2 deletions.
diff --git a/palace/linalg/CMakeLists.txt b/palace/linalg/CMakeLists.txt
@@ -35,6 +35,7 @@ target_sources(${LIB_TARGET_NAME}
 set(TARGET_SOURCES_DEVICE
   ${TARGET_SOURCES_DEVICE}
   ${CMAKE_CURRENT_SOURCE_DIR}/chebyshev.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/densematrix.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/jacobi.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/operator.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/vector.cpp

diff --git a/palace/linalg/densematrix.cpp b/palace/linalg/densematrix.cpp
@@ -6,6 +6,7 @@
 #include <functional>
 #include <limits>
 #include <mfem.hpp>
+#include <mfem/linalg/kernels.hpp>
 
 namespace palace
 {
@@ -222,6 +223,46 @@ mfem::DenseTensor MatrixPow(const mfem::DenseTensor &T, double p)
   return S;
 }
 
+double SingularValueMax(const mfem::DenseMatrix &M)
+{
+  MFEM_ASSERT(
+      M.Height() == M.Width() && M.Height() > 0 && M.Height() <= 3,
+      "Matrix singular values only available for square matrices of dimension <= 3!");
+  const int N = M.Height();
+  if (N == 1)
+  {
+    return M(0, 0);
+  }
+  else if (N == 2)
+  {
+    return mfem::kernels::CalcSingularvalue<2>(M.Data(), 0);
+  }
+  else
+  {
+    return mfem::kernels::CalcSingularvalue<3>(M.Data(), 0);
+  }
+}
+
+double SingularValueMin(const mfem::DenseMatrix &M)
+{
+  MFEM_ASSERT(
+      M.Height() == M.Width() && M.Height() > 0 && M.Height() <= 3,
+      "Matrix singular values only available for square matrices of dimension <= 3!");
+  const int N = M.Height();
+  if (N == 1)
+  {
+    return M(0, 0);
+  }
+  else if (N == 2)
+  {
+    return mfem::kernels::CalcSingularvalue<2>(M.Data(), 1);
+  }
+  else
+  {
+    return mfem::kernels::CalcSingularvalue<3>(M.Data(), 2);
+  }
+}
+
 mfem::DenseTensor Mult(const mfem::DenseTensor &A, const mfem::DenseTensor &B)
 {
   MFEM_VERIFY(A.SizeK() == B.SizeK(),

diff --git a/palace/linalg/densematrix.hpp b/palace/linalg/densematrix.hpp
@@ -29,6 +29,10 @@ mfem::DenseMatrix MatrixPow(const mfem::DenseMatrix &M, double p);
 
 mfem::DenseTensor MatrixPow(const mfem::DenseTensor &T, double p);
 
+double SingularValueMax(const mfem::DenseMatrix &M);
+
+double SingularValueMin(const mfem::DenseMatrix &M);
+
 mfem::DenseTensor Mult(const mfem::DenseTensor &A, const mfem::DenseTensor &B);
 
 }  // namespace palace::linalg

diff --git a/palace/models/materialoperator.cpp b/palace/models/materialoperator.cpp
@@ -262,8 +262,8 @@ void MaterialOperator::SetUpMaterialProperties(const IoData &iodata,
     // √((μ ε)⁻¹)
     Mult(mat_mu, mat_epsilon(count), T);
     mat_c0(count) = linalg::MatrixPow(T, -0.5);
-    mat_c0_min[count] = mat_c0(count).CalcSingularvalue(sdim - 1);
-    mat_c0_max[count] = mat_c0(count).CalcSingularvalue(0);
+    mat_c0_min[count] = linalg::SingularValueMin(mat_c0(count));
+    mat_c0_max[count] = linalg::SingularValueMax(mat_c0(count));
 
     // Electrical conductivity, σ
     mat_sigma(count) = ToDenseMatrix(data.sigma);