diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index aa55025..d42ae8a 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -46,8 +46,12 @@ class Rule { // TODO: Define x and w as Eigen::VectorXd std::vector x, w, x_forward, x_backward; T a, b; + // Default constructor + Rule() {}; + // Constructor with x, w, x_forward, x_backward, a, b Rule(const std::vector& x, const std::vector& w, std::vector x_forward, std::vector x_backward, T a = -1, T b = 1) : x(x), w(w), x_forward(x_forward), x_backward(x_backward), a(a), b(b) {} + // Constructor with x, w, a, b Rule(const std::vector& x, const std::vector& w, T a = -1, T b = 1) : x(x), w(w), a(a), b(b) { this->x_forward = x_forward.empty() ? std::vector(x.size(), 0) : x_forward; diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp index 2dd7469..3ac162c 100644 --- a/include/sparseir/sve.hpp +++ b/include/sparseir/sve.hpp @@ -130,7 +130,7 @@ class AbstractSVE { template class SamplingSVE : public AbstractSVE { public: - std::shared_ptr kernel; + K kernel; double epsilon; int n_gauss; int nsvals_hint; @@ -142,25 +142,25 @@ class SamplingSVE : public AbstractSVE { Rule gauss_y; // Constructor - SamplingSVE(std::shared_ptr kernel_, double epsilon_, int n_gauss_ = -1) + SamplingSVE(K kernel_, double epsilon_, int n_gauss_ = -1) : kernel(kernel_), - epsilon(epsilon_), - n_gauss((n_gauss_ > 0) ? n_gauss_ : sve_hints(*kernel, epsilon).ngauss()), - rule(convert(legendre(n_gauss))), - gauss_x(rule.piecewise(segs_x)), // Initialize gauss_x - gauss_y(rule.piecewise(segs_y)) // Initialize gauss_y + epsilon(epsilon_) { - auto hints = sve_hints(*kernel, epsilon); - // TODO: This is not correct, we need to use the segments from the hints + n_gauss = (n_gauss_ > 0) ? n_gauss_ : sve_hints(kernel, epsilon).ngauss(); + // TODO: Implement Rule(n_gauss) + rule = convert(legendre(n_gauss)); + auto hints = sve_hints(kernel, epsilon); nsvals_hint = hints.nsvals(); segs_x = hints.template segments_x(); segs_y = hints.template segments_y(); + gauss_x = rule.piecewise(segs_x); + gauss_y = rule.piecewise(segs_y); } // Compute matrices for SVD std::vector> matrices() const override { std::vector> mats; - Eigen::MatrixX A = matrix_from_gauss(*kernel, gauss_x, gauss_y); + Eigen::MatrixX A = matrix_from_gauss(kernel, gauss_x, gauss_y); // Element-wise multiplication with square roots of weights for (int i = 0; i < gauss_x.w.size(); ++i) { @@ -304,7 +304,7 @@ class SamplingSVE : public AbstractSVE { PiecewiseLegendrePolyVector ulx(polyvec_u); PiecewiseLegendrePolyVector vly(polyvec_v); canonicalize(ulx, vly); - return SVEResult(ulx, s, vly, *kernel, epsilon); + return SVEResult(ulx, s, vly, kernel, epsilon); } }; diff --git a/test/sve.cxx b/test/sve.cxx index 12893f6..6f83c10 100644 --- a/test/sve.cxx +++ b/test/sve.cxx @@ -13,8 +13,6 @@ using std::invalid_argument; -using xprec::DDouble; - // Function to check smoothness void check_smooth(const std::function& u, const std::vector& s, double uscale, double fudge_factor) { /* @@ -58,15 +56,31 @@ void check_smooth(const std::function& u, const std::vector(std::make_shared(sparseir::LogisticKernel(10.0)), 1e-6); - // REQUIRE(sve.nsvals_hint == 10); - // REQUIRE(sve.n_gauss == 10); - // REQUIRE(sve.segs_x.size() == 10); - // REQUIRE(sve.segs_y.size() == 10); - // REQUIRE(sve.gauss_x.size() == 10); - // REQUIRE(sve.gauss_y.size() == 10); + auto hints = sparseir::sve_hints(lk, 1e-6); + int nsvals_hint = hints.nsvals(); + int n_gauss = hints.ngauss(); + std::vector segs_x = hints.template segments_x(); + std::vector segs_y = hints.template segments_y(); + + // Ensure `convert` is declared before this line + sparseir::Rule rule = sparseir::convert(sparseir::legendre(n_gauss)); + + sparseir::Rule gauss_x = rule.piecewise(segs_x); + sparseir::Rule gauss_y = rule.piecewise(segs_y); + auto ssve1 = sparseir::SamplingSVE(lk, 1e-6); + REQUIRE(ssve1.n_gauss == n_gauss); + auto ssve2 = sparseir::SamplingSVE(lk, 1e-6, 12); + REQUIRE(ssve2.n_gauss == 12); + + auto ssve1_double = sparseir::SamplingSVE(lk, 1e-6); + REQUIRE(ssve1_double.n_gauss == n_gauss); + auto ssve2_double = sparseir::SamplingSVE(lk, 1e-6, 12); + REQUIRE(ssve2_double.n_gauss == 12); + + //auto ssve1_ddouble = sparseir::SamplingSVE(lk, 1e-6); + //REQUIRE(ssve1_ddouble.n_gauss == n_gauss); + //auto ssve2_ddouble = sparseir::SamplingSVE(lk, 1e-6, 12); + //REQUIRE(ssve2_ddouble.n_gauss == 12); } TEST_CASE("sve.cpp", "[sve]")