Implement svd.hpp

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);
 }
 

include/sparseir/_root.hpp

@@ -1,3 +1,5 @@
+#pragma once
+
 #include <vector>
 #include <cmath>
 #include <algorithm>

include/sparseir/poly.hpp

@@ -1,3 +1,4 @@
+#pragma once
 // C++ Standard Library headers
 #include <algorithm>
 #include <cassert>

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"
+#include "./kernel.hpp"
+#include "./svd.hpp"

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);
+    }
+}

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)

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);
+    */
+}