From c358bbc923235c14139fc178f5cbfaafd6c3f59e Mon Sep 17 00:00:00 2001
From: Sebastian Grimberg <sjg@amazon.com>
Date: Thu, 25 Apr 2024 12:00:27 -0700
Subject: [PATCH] Add singular value wrappers to dense matrix functionality

---
 palace/linalg/CMakeLists.txt       |  1 +
 palace/linalg/densematrix.cpp      | 41 ++++++++++++++++++++++++++++++
 palace/linalg/densematrix.hpp      |  4 +++
 palace/models/materialoperator.cpp |  4 +--
 4 files changed, 48 insertions(+), 2 deletions(-)

diff --git a/palace/linalg/CMakeLists.txt b/palace/linalg/CMakeLists.txt
index 421324738..86eae9f82 100644
--- 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
index 5e37f9e5e..2c990bbb7 100644
--- 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
index 080f45168..41eb40fcd 100644
--- 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
index 4907775af..2af6efdd5 100644
--- 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);