diff --git a/include/sparseir/_impl.hpp b/include/sparseir/_impl.hpp new file mode 100644 index 0000000..fd201be --- /dev/null +++ b/include/sparseir/_impl.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "xprec/ddouble.hpp" + + +namespace sparseir { +// General case +template +inline T sqrt_impl(const T& x) { + return std::sqrt(x); +} + +// Specialization for DDouble +template<> +inline xprec::DDouble sqrt_impl(const xprec::DDouble& x) { + return xprec::sqrt(x); +} + +template +inline T cosh_impl(const T& x) { + return std::cosh(x); +} + +template<> +inline xprec::DDouble cosh_impl(const xprec::DDouble& x) { + return xprec::cosh(x); +} + +template +inline T sinh_impl(const T& x) { + return std::sinh(x); +} + +template<> +inline xprec::DDouble sinh_impl(const xprec::DDouble& x) { + return xprec::sinh(x); +} + +template +inline T exp_impl(const T& x) { + return std::exp(x); +} + +template<> +inline xprec::DDouble exp_impl(const xprec::DDouble& x) { + return xprec::exp(x); +} + +} \ No newline at end of file diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 1181f75..89105e6 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -357,14 +357,14 @@ namespace sparseir if (v >= 0) { - numerator = std::exp(-u_plus * abs_v); + numerator = exp_impl(-u_plus * abs_v); } else { - numerator = std::exp(-u_minus * abs_v); + numerator = exp_impl(-u_minus * abs_v); } - denominator = 1.0 + std::exp(-abs_v); + denominator = 1.0 + exp_impl(-abs_v); return numerator / denominator; } @@ -736,7 +736,7 @@ namespace sparseir // Calculate diffs using the inverse hyperbolic cosine std::vector diffs(nzeros); for (int i = 0; i < nzeros; ++i) { - diffs[i] = 1.0 / std::cosh(temp[i]); + diffs[i] = 1.0 / cosh_impl(temp[i]); } // Calculate cumulative sum of diffs @@ -781,7 +781,7 @@ namespace sparseir // Calculate trailing differences for (int i = 20; i < nzeros; ++i) { T x = (T)0.141 * i; - diffs.push_back(0.25 * std::exp(-x)); + diffs.push_back(0.25 * exp_impl(-x)); } // Calculate cumulative sum of diffs @@ -945,7 +945,7 @@ namespace sparseir bool sinh_range = 1e-200 < v_half && v_half < 85; if (xy_small && sinh_range) { - return y * std::sinh(xv_half) / std::sinh(v_half); + return y * sinh_impl(xv_half) / sinh_impl(v_half); } else { @@ -1049,8 +1049,14 @@ namespace sparseir{ for (size_t i = 0; i < n; ++i) { threads.emplace_back([&, i]() { for (size_t j = 0; j < m; ++j) { - res(i, j) = kernel(gauss_x.x[i], gauss_y.x[j], - gauss_x.x_forward[i], gauss_x.x_backward[i]); + // TODO: return type should be T + double ko = kernel( + static_cast(gauss_x.x[i]), + static_cast(gauss_y.x[j]), + static_cast(gauss_x.x_forward[i]), + static_cast(gauss_x.x_backward[i]) + ); + res(i, j) = T(ko); } }); } diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp index 46c7615..9672402 100644 --- a/include/sparseir/sparseir-header-only.hpp +++ b/include/sparseir/sparseir-header-only.hpp @@ -1,5 +1,6 @@ #pragma once +#include "./_impl.hpp" #include "./_linalg.hpp" #include "./_root.hpp" #include "./_specfuncs.hpp" diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp index 7f9a4e3..565f46a 100644 --- a/include/sparseir/sve.hpp +++ b/include/sparseir/sve.hpp @@ -164,11 +164,11 @@ class SamplingSVE : public AbstractSVE { // Element-wise multiplication with square roots of weights for (int i = 0; i < gauss_x.w.size(); ++i) { - A.row(i) *= std::sqrt(gauss_x.w[i]); + A.row(i) *= sqrt_impl(gauss_x.w[i]); } for (int j = 0; j < gauss_y.w.size(); ++j) { - A.col(j) *= std::sqrt(gauss_y.w[j]); + A.col(j) *= sqrt_impl(gauss_y.w[j]); } mats.push_back(A); return mats; @@ -253,7 +253,7 @@ class SamplingSVE : public AbstractSVE { { for (int i = 0; i < u_data.dimension(0); ++i) { - u_data(i, j, k) *= std::sqrt(0.5 * dsegs_x[j]); + u_data(i, j, k) *= sqrt_impl(T(0.5) * dsegs_x[j]); } } } @@ -265,7 +265,7 @@ class SamplingSVE : public AbstractSVE { { for (int i = 0; i < v_data.dimension(0); ++i) { - v_data(i, j, k) *= std::sqrt(0.5 * dsegs_y[j]); + v_data(i, j, k) *= sqrt_impl(T(0.5) * dsegs_y[j]); } } } @@ -349,80 +349,6 @@ class CentrosymmSVE : public AbstractSVE { auto mats_odd = odd.matrices(); return {mats_even[0], mats_odd[0]}; } - /* - SVEResult postprocess(const std::vector>& u_list, - const std::vector>& s_list, - const std::vector>& v_list) const override { - SVEResult result_even = even.postprocess({ u_list[0] }, { s_list[0] }, { v_list[0] }); - SVEResult result_odd = odd.postprocess({ u_list[1] }, { s_list[1] }, { v_list[1] }); - - - // Merge results - auto u = result_even.u; - u.insert(u.end(), result_odd.u.begin(), result_odd.u.end()); - auto s = result_even.s; - s.insert(s.end(), result_odd.s.begin(), result_odd.s.end()); - auto v = result_even.v; - v.insert(v.end(), result_odd.v.begin(), result_odd.v.end()); - - std::vector signs(result_even.s.size(), +1); - signs.insert(signs.end(), result_odd.s.size(), -1); - - // Sort singular values and associated vectors - std::vector indices(s.size()); - std::iota(indices.begin(), indices.end(), 0); - std::stable_sort(indices.begin(), indices.end(), [&](size_t i1, size_t i2) { - return s[i1] > s[i2]; - }); - - std::vector u_sorted(u.size()); - std::vector s_sorted(s.size()); - std::vector v_sorted(v.size()); - std::vector signs_sorted(signs.size()); - for (size_t i = 0; i < indices.size(); ++i) { - u_sorted[i] = u[indices[i]]; - s_sorted[i] = s[indices[i]]; - v_sorted[i] = v[indices[i]]; - signs_sorted[i] = signs[indices[i]]; - } - - // Extend to negative side - // Assuming definitions of necessary functions and data structures - auto full_hints = sve_hints(kernel, epsilon); - auto segs_x_full = segments_x(full_hints); - auto segs_y_full = segments_y(full_hints); - - std::vector u_complete(u_sorted.size()); - std::vector v_complete(v_sorted.size()); - - Eigen::Array poly_flip_x = Eigen::Array::LinSpaced(u_sorted[0].data.rows(), 0, u_sorted[0].data.rows() - 1); - poly_flip_x = poly_flip_x.unaryExpr([](T x) { return std::pow(-1, x); }); - - for (size_t i = 0; i < u_sorted.size(); ++i) { - Eigen::MatrixX u_pos_data = u_sorted[i].data / std::sqrt(T(2)); - Eigen::MatrixX v_pos_data = v_sorted[i].data / std::sqrt(T(2)); - - Eigen::MatrixX u_neg_data = u_pos_data.rowwise().reverse().array().colwise() * poly_flip_x.array() * T(signs_sorted[i]); - Eigen::MatrixX v_neg_data = v_pos_data.rowwise().reverse().array().colwise() * poly_flip_x.array() * T(signs_sorted[i]); - - Eigen::MatrixX u_data_full(u_pos_data.rows(), u_pos_data.cols() + u_neg_data.cols()); - u_data_full << u_neg_data, u_pos_data; - - Eigen::MatrixX v_data_full(v_pos_data.rows(), v_pos_data.cols() + v_neg_data.cols()); - v_data_full << v_neg_data, v_pos_data; - - u_complete[i] = PiecewiseLegendrePoly(u_data_full, segs_x_full, static_cast(i), signs_sorted[i]); - v_complete[i] = PiecewiseLegendrePoly(v_data_full, segs_y_full, static_cast(i), signs_sorted[i]); - } - - - // Convert std::vector to Eigen::VectorXd - Eigen::Map s_eigen(s_sorted.data(), s_sorted.size()); - - // Then use the converted Eigen vector in the constructor - return SVEResult(u_complete, s_eigen, v_complete, kernel, epsilon); - } - */ // Replace the vector merging code with Eigen operations SVEResult postprocess(const std::vector>& u_list, @@ -487,15 +413,16 @@ std::shared_ptr> determine_sve(const K& kernel, double safe_ep } // Function to truncate singular values +template inline void truncate_singular_values( - std::vector &u_list, - std::vector &s_list, - std::vector &v_list, - double rtol, + std::vector> &u_list, + std::vector> &s_list, + std::vector> &v_list, + T rtol, int lmax) { // Collect all singular values - std::vector all_singular_values; + std::vector all_singular_values; for (const auto &s : s_list) { for (int i = 0; i < s.size(); ++i) @@ -503,10 +430,10 @@ inline void truncate_singular_values( all_singular_values.push_back(s(i)); } } - std::sort(all_singular_values.begin(), all_singular_values.end(), std::greater()); + std::sort(all_singular_values.begin(), all_singular_values.end(), std::greater()); // Determine cutoff - double cutoff = rtol * all_singular_values.front(); + T cutoff = rtol * all_singular_values.front(); if (lmax < static_cast(all_singular_values.size())) { cutoff = std::max(cutoff, all_singular_values[lmax - 1]); @@ -543,7 +470,7 @@ auto pre_postprocess(K &kernel, double safe_epsilon, int n_gauss, double cutoff { auto sve = determine_sve(kernel, safe_epsilon, n_gauss); // Compute SVDs - std::vector> matrices = sve.matrices(); + std::vector> matrices = sve->matrices(); std::vector>> svds; for (const auto& mat : matrices) { Eigen::BDCSVD> svd(mat, Eigen::ComputeThinU | Eigen::ComputeThinV); @@ -560,18 +487,18 @@ auto pre_postprocess(K &kernel, double safe_epsilon, int n_gauss, double cutoff } // Apply cutoff and lmax - T cutoff_actual = std::isnan(cutoff) ? 2 * std::numeric_limits::epsilon() : cutoff; - //truncate_singular_values(u_list, s_list, v_list, cutoff_actual, lmax); + T cutoff_actual = std::isnan(cutoff) ? 2 * T(std::numeric_limits::epsilon()) : T(cutoff); + truncate_singular_values(u_list, s_list, v_list, cutoff_actual, lmax); // Postprocess to get the SVEResult - return sve.postprocess(u_list, s_list, v_list); + return sve->postprocess(u_list, s_list, v_list); } // Function to compute SVE result template auto compute_sve(K kernel, - std::string Twork = "Floatt64", - double cutoff = std::numeric_limits::quiet_NaN(), double epsilon = std::numeric_limits::quiet_NaN(), + double cutoff = std::numeric_limits::quiet_NaN(), + std::string Twork = "Floatt64", int lmax = std::numeric_limits::max(), int n_gauss = -1, const std::string& svd_strat = "auto") { diff --git a/test/sve.cxx b/test/sve.cxx index 407a0a7..03185ab 100644 --- a/test/sve.cxx +++ b/test/sve.cxx @@ -54,7 +54,7 @@ void check_smooth(const std::function& u, const std::vector(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); + 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("CentrosymmSVE", "[CentrosymmSVE]"){ +TEST_CASE("CentrosymmSVE", "[CentrosymmSVE]") { sparseir::LogisticKernel lk(10.0); auto hints = sparseir::sve_hints(lk, 1e-6); int nsvals_hint = hints.nsvals(); @@ -93,23 +93,26 @@ TEST_CASE("CentrosymmSVE", "[CentrosymmSVE]"){ sparseir::Rule rule = sparseir::legendre(n_gauss); + // TODO: resole std::bad_cast //auto sve = sparseir::CentrosymmSVE(lk, 1e-6); //auto sve_double = sparseir::CentrosymmSVE(lk, 1e-6); //auto sve_ddouble = sparseir::CentrosymmSVE(lk, 1e-6); } -TEST_CASE("sve.cpp", "[sve]") -{ +TEST_CASE("sve.cpp", "[compute_sve]") { + //auto sve = sparseir::compute_sve(sparseir::LogisticKernel(10.0)); // Define a map to store SVEResult objects //auto sve_logistic = std::map < int, sparseir::SVEResult>{ // {10, sparseir::compute_sve(sparseir::LogisticKernel(10.0))}, // {42, sparseir::compute_sve(sparseir::LogisticKernel(42.0))}, // {10000, sparseir::compute_sve(sparseir::LogisticKernel(10000.0))}, - // {100000000, sparseir::compute_sve(sparseir::LogisticKernel(10000.0), 1e-12)}}; + // //{100000000, sparseir::compute_sve(sparseir::LogisticKernel(10000.0), 1e-12)}, + // }; SECTION("smooth with Λ =") { for (int Lambda : {10, 42, 10000}) { + REQUIRE(true); // sparseir::FiniteTempBasis basis(1, Lambda, sve_logistic[Lambda]); // TODO: Check that the maximum implementation is defined // check_smooth(basis.u, basis.s, 2 * sparseir::maximum(basis.u(1)), 24); @@ -130,7 +133,7 @@ TEST_CASE("sve.cpp", "[sve]") */ /* - SECTION("num roots û with stat =, Λ =") { + SECTION("num roots û with stat =, Λ =") { for (const auto& stat : {Fermionic(), Bosonic()}) { for (int Lambda : {10, 42, 10000}) { FiniteTempBasis basis(stat, 1, Lambda, sparseir::sve_logistic[Lambda]); @@ -155,25 +158,25 @@ TEST_CASE("sve.cpp", "[sve]") } } */ +} - SECTION("choose_accuracy") { - REQUIRE(sparseir::choose_accuracy(nullptr, nullptr) == std::make_tuple(2.2204460492503131e-16, "Float64x2", "default")); - REQUIRE(sparseir::choose_accuracy(nullptr, "Float64") == std::make_tuple(1.4901161193847656e-8, "Float64", "default")); - REQUIRE(sparseir::choose_accuracy(nullptr, "Float64x2") == std::make_tuple(2.2204460492503131e-16, "Float64x2", "default")); +TEST_CASE("sve.cpp", "[choose_accuracy]") { + REQUIRE(sparseir::choose_accuracy(nullptr, nullptr) == std::make_tuple(2.2204460492503131e-16, "Float64x2", "default")); + REQUIRE(sparseir::choose_accuracy(nullptr, "Float64") == std::make_tuple(1.4901161193847656e-8, "Float64", "default")); + REQUIRE(sparseir::choose_accuracy(nullptr, "Float64x2") == std::make_tuple(2.2204460492503131e-16, "Float64x2", "default")); - REQUIRE(sparseir::choose_accuracy(1e-6, nullptr) == std::make_tuple(1.0e-6, "Float64", "default")); - // Note: Catch2 doesn't have a built-in way to capture logs. - // You might need to implement a custom logger or use a library that supports log capturing. - // Add debug output to see the actual return value - REQUIRE(sparseir::choose_accuracy(1e-8, nullptr) == std::make_tuple(1.0e-8, "Float64x2", "default")); - REQUIRE(sparseir::choose_accuracy(1e-20, nullptr) == std::make_tuple(1.0e-20, "Float64x2", "default")); + REQUIRE(sparseir::choose_accuracy(1e-6, nullptr) == std::make_tuple(1.0e-6, "Float64", "default")); + // Note: Catch2 doesn't have a built-in way to capture logs. + // You might need to implement a custom logger or use a library that supports log capturing. + // Add debug output to see the actual return value + REQUIRE(sparseir::choose_accuracy(1e-8, nullptr) == std::make_tuple(1.0e-8, "Float64x2", "default")); + REQUIRE(sparseir::choose_accuracy(1e-20, nullptr) == std::make_tuple(1.0e-20, "Float64x2", "default")); - REQUIRE(sparseir::choose_accuracy(1e-10, "Float64") == std::make_tuple(1.0e-10, "Float64", "accurate")); + REQUIRE(sparseir::choose_accuracy(1e-10, "Float64") == std::make_tuple(1.0e-10, "Float64", "accurate")); - REQUIRE(sparseir::choose_accuracy(1e-6, "Float64") == std::make_tuple(1.0e-6, "Float64", "default")); - REQUIRE(sparseir::choose_accuracy(1e-6, "Float64", "auto") == std::make_tuple(1.0e-6, "Float64", "default")); - REQUIRE(sparseir::choose_accuracy(1e-6, "Float64", "accurate") == std::make_tuple(1.0e-6, "Float64", "accurate")); - } + REQUIRE(sparseir::choose_accuracy(1e-6, "Float64") == std::make_tuple(1.0e-6, "Float64", "default")); + REQUIRE(sparseir::choose_accuracy(1e-6, "Float64", "auto") == std::make_tuple(1.0e-6, "Float64", "default")); + REQUIRE(sparseir::choose_accuracy(1e-6, "Float64", "accurate") == std::make_tuple(1.0e-6, "Float64", "accurate")); /* SECTION("truncate") {