From cb472caff28a4af6732f2d4b20f14da06e962005 Mon Sep 17 00:00:00 2001
From: Hiroshi Shinaoka <h.shinaoka@gmail.com>
Date: Thu, 5 Dec 2024 10:06:41 -0500
Subject: [PATCH] Implement svd.hpp

---
 include/sparseir/_linalg.hpp              |  8 ++--
 include/sparseir/_root.hpp                |  2 +
 include/sparseir/poly.hpp                 |  1 +
 include/sparseir/sparseir-header-only.hpp |  4 +-
 include/sparseir/svd.hpp                  | 37 ++++++++++++++++
 test/CMakeLists.txt                       |  1 +
 test/svd.cxx                              | 54 +++++++++++++++++++++++
 7 files changed, 103 insertions(+), 4 deletions(-)
 create mode 100644 include/sparseir/svd.hpp
 create mode 100644 test/svd.cxx

diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp
index 60d187c..5473c50 100644
--- a/include/sparseir/_linalg.hpp
+++ b/include/sparseir/_linalg.hpp
@@ -1,3 +1,5 @@
+#pragma once
+
 #include <iostream>
 
 #include <Eigen/Dense>
@@ -38,7 +40,7 @@ int argmax(const Eigen::MatrixBase<Derived>& vec) {
 /*
 This function ports Julia's implementation of the `invperm` function to C++.
 */
-Eigen::VectorXi invperm(const Eigen::VectorXi& a) {
+inline Eigen::VectorXi invperm(const Eigen::VectorXi& a) {
     int n = a.size();
     Eigen::VectorXi b(n);
     b.setConstant(-1);
@@ -192,12 +194,12 @@ MatrixX<T> getPropertyR(const QRPivoted<T>& F) {
 }
 
 // General template for _copysign, handles standard floating-point types like double and float
-double _copysign(double x, double y) {
+inline double _copysign(double x, double y) {
     return std::copysign(x, y);
 }
 
 // Specialization for xprec::DDouble type, assuming xprec::copysign is defined
-xprec::DDouble _copysign(xprec::DDouble x, xprec::DDouble y) {
+inline xprec::DDouble _copysign(xprec::DDouble x, xprec::DDouble y) {
     return xprec::copysign(x, y);
 }
 
diff --git a/include/sparseir/_root.hpp b/include/sparseir/_root.hpp
index eb33d13..322f6c9 100644
--- a/include/sparseir/_root.hpp
+++ b/include/sparseir/_root.hpp
@@ -1,3 +1,5 @@
+#pragma once
+
 #include <vector>
 #include <cmath>
 #include <algorithm>
diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp
index 1ee0faf..413c75f 100644
--- a/include/sparseir/poly.hpp
+++ b/include/sparseir/poly.hpp
@@ -1,3 +1,4 @@
+#pragma once
 // C++ Standard Library headers
 #include <algorithm>
 #include <cassert>
diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp
index 86d5095..31e5193 100644
--- a/include/sparseir/sparseir-header-only.hpp
+++ b/include/sparseir/sparseir-header-only.hpp
@@ -3,8 +3,10 @@
 #include "./_specfuncs.hpp"
 
 #include "./_root.hpp"
+#include "./_linalg.hpp"
 
 #include "./gauss.hpp"
 #include "./freq.hpp"
 #include "./poly.hpp"
-#include "./kernel.hpp"
\ No newline at end of file
+#include "./kernel.hpp"
+#include "./svd.hpp"
\ No newline at end of file
diff --git a/include/sparseir/svd.hpp b/include/sparseir/svd.hpp
new file mode 100644
index 0000000..0eda606
--- /dev/null
+++ b/include/sparseir/svd.hpp
@@ -0,0 +1,37 @@
+#pragma once
+
+// C++ Standard Library headers
+#include <stdexcept>
+#include <iostream>
+
+// Eigen headers
+#include <Eigen/Dense>
+
+namespace sparseir
+{
+
+    using namespace Eigen;
+
+    template <typename T>
+    std::tuple<MatrixX<T>, VectorX<T>, MatrixX<T>> compute_svd(const MatrixX<T> &A, int n_sv_hint = 0, std::string strategy = "default")
+    {
+        if (n_sv_hint != 0)
+        {
+            std::cout << "n_sv_hint is set but will not be used in the current implementation!" << std::endl;
+        }
+
+        if (strategy != "default")
+        {
+            std::cout << "strategy is set but will not be used in the current implementation!" << std::endl;
+        }
+
+        MatrixX<T> A_copy = A; // create a copy of A
+        return tsvd<T>(A_copy);
+        // auto svd_result = tsvd<T>(A_copy);
+        // MatrixX<T> u = std::get<0>(svd_result);
+        // VectorX<T> s = std::get<1>(svd_result);
+        // MatrixX<T> v = std::get<2>(svd_result);
+
+        // return std::make_tuple(u, s, v);
+    }
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 1282b4d..5bfe076 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -20,6 +20,7 @@ add_executable(libsparseirtests
     freq.cxx
     poly.cxx
     kernel.cxx
+    svd.cxx
 )
 
 target_link_libraries(libsparseirtests PRIVATE Catch2::Catch2WithMain)
diff --git a/test/svd.cxx b/test/svd.cxx
new file mode 100644
index 0000000..07e490c
--- /dev/null
+++ b/test/svd.cxx
@@ -0,0 +1,54 @@
+#include <Eigen/Dense>
+#include <vector>
+#include <algorithm>
+#include <random>
+#include <cmath>
+#include <stdexcept>
+#include <iostream>
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <xprec/ddouble-header-only.hpp>
+#include <sparseir/sparseir-header-only.hpp>
+
+// test_piecewise_legendre_poly.cpp
+
+#include <catch2/catch_test_macros.hpp>
+#include <Eigen/Dense>
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <random>
+#include <functional>
+
+TEST_CASE("svd.cpp")
+{
+    using namespace sparseir;
+    using namespace Eigen;
+    using namespace xprec;
+
+    // Create a matrix of Float64x2 equivalent (here just Eigen::MatrixXd for simplicity)
+    /**
+    Eigen::MatrixXd mat64x2 = Eigen::MatrixXd::Random(4, 6);
+    REQUIRE_NOTHROW(sparseir::compute_svd(mat64x2, "accurate", 2));
+    std::cout << "n_sv_hint is set but will not be used in the current implementation!" << std::endl;
+
+    REQUIRE_NOTHROW({
+        compute_svd(mat64x2, "accurate");
+        std::cout << "strategy is set but will not be used in the current implementation!" << std::endl;
+    });
+    */
+
+    // Create a standard matrix
+    MatrixX<DDouble> mat = MatrixX<DDouble>::Random(5, 6);
+    auto svd_result = compute_svd(mat, 0, "default");
+    auto U = std::get<0>(svd_result);
+    auto S = std::get<1>(svd_result);
+    auto V = std::get<2>(svd_result);
+    auto diff = (mat - U * S.asDiagonal() * V.transpose()).norm() / mat.norm();
+    REQUIRE(diff < 1e-28);
+
+    /*
+    REQUIRE_THROWS_AS(compute_svd(mat, "fast"), std::domain_error);
+    */
+}
\ No newline at end of file