From 0a82cfb1a2be84d1552c28481644bd9ae15a32d6 Mon Sep 17 00:00:00 2001 From: "Hugo U. R. Strand" Date: Fri, 6 Nov 2026 11:14:04 +0100 Subject: [PATCH 1/3] [imfreq] select DLR Matsubara nodes on a dyadic fine grid imfreq_ops built the analytic continuation kernel at every Matsubara index up to the cutoff nmax = max(ceil(lambda), 20) before pivoting, so memory and time grew linearly in lambda: 6.8-9.4 GB and 54-104 s at lambda = 1e6. Pivot instead on a dyadic grid (build_if_fine): panels of 2r consecutive indices whose spacing doubles from panel to panel, up to the cutoff, which is always included, mirrored to negative frequencies as pivrgs_sym and flip_to_mirror expect. When 2r >= nmax this is the full grid and the non-symmetrized nodes are unchanged. At lambda = 1e6 construction now takes 18-32 MB and 0.15-0.36 s. On this grid the sum-of-norms pair score in the rank-r pivrgs_sym can pick a mirror pair whose second row is parallel to the first after projection, leaving an ill-conditioned node set (up to 9.5 eps at eps = 1e-10). Score pairs by the squared volume of their residuals instead, and keep the sum for a pair that exceeds the column dimension, where every volume vanishes. Also evaluate 2n+1 and 2n in double in k_if, which overflowed int for |n| >= 2^30, now reachable at lambda ~ 1e9. DLR Matsubara nodes change for symmetrized grids and, where 2r < nmax, for non-symmetrized grids. Benchmarked against dense node selection on lambda in [10, 1e7], eps 1e-6 to 1e-14, both statistics, with and without symmetrization: errors within a few eps of dense everywhere dense is affordable. Co-Authored-By: Claude Opus 5.5 --- c++/cppdlr/dlr_build.cpp | 41 ++++++++++++++++++++++ c++/cppdlr/dlr_build.hpp | 35 ++++++++++++++++++ c++/cppdlr/dlr_imfreq.cpp | 8 +++-- c++/cppdlr/dlr_kernels.cpp | 4 +-- c++/cppdlr/utils.hpp | 23 ++++++++++-- test/c++/dlr_build.cpp | 72 ++++++++++++++++++++++++++++++++++++++ test/c++/imfreq_ops.cpp | 47 +++++++++++++++++++++++++ test/c++/pivrgs.cpp | 32 +++++++++++++++++ 8 files changed, 254 insertions(+), 8 deletions(-) diff --git a/c++/cppdlr/dlr_build.cpp b/c++/cppdlr/dlr_build.cpp index 36d96fa..436cb55 100644 --- a/c++/cppdlr/dlr_build.cpp +++ b/c++/cppdlr/dlr_build.cpp @@ -273,6 +273,47 @@ namespace cppdlr { } } + nda::matrix build_k_if(nda::vector_const_view n, nda::vector_const_view om, statistic_t statistic) { + + auto kmat = nda::matrix(n.size(), om.size()); + for (int i = 0; i < n.size(); ++i) { + for (int j = 0; j < om.size(); ++j) { kmat(i, j) = k_if(n(i), om(j), statistic); } + } + return kmat; + } + + nda::vector build_if_fine(int nmax, statistic_t statistic, int npt) { + + if (nmax < 1) throw std::runtime_error("Choose nmax >= 1."); + if (npt < 1) throw std::runtime_error("Choose npt >= 1."); + + // Non-negative half of the grid: n = 0, ..., nmax - 1 (fermionic case) or + // n = 1, ..., nmax (bosonic case), in panels of npt indices with spacing 2^j + int nfirst = (statistic == Fermion ? 0 : 1); + int nlast = (statistic == Fermion ? nmax - 1 : nmax); + + auto half = std::vector{}; + long nn = nfirst; + for (long step = 1; nn < nlast; step *= 2) { + for (int i = 0; i < npt && nn < nlast; ++i, nn += step) half.push_back(static_cast(nn)); + } + half.push_back(nlast); + + // Mirror to negative frequencies, n -> -n-1 (fermionic case) or n -> -n + // (bosonic case), and include n = 0 in the bosonic case + int nh = half.size(); + int shift = (statistic == Fermion ? 1 : 0); + int mid = (statistic == Fermion ? 0 : 1); + auto n = nda::vector(2 * nh + mid); + for (int i = 0; i < nh; ++i) { + n(nh - 1 - i) = -half[i] - shift; + n(nh + mid + i) = half[i]; + } + if (mid == 1) n(nh) = 0; + + return n; + } + nda::vector build_dlr_rf(double lambda, double eps, bool symmetrize) { if (eps < 1e-14) { diff --git a/c++/cppdlr/dlr_build.hpp b/c++/cppdlr/dlr_build.hpp index 87ea6a7..7bb394d 100644 --- a/c++/cppdlr/dlr_build.hpp +++ b/c++/cppdlr/dlr_build.hpp @@ -119,6 +119,41 @@ namespace cppdlr { */ nda::matrix build_k_if(int nmax, nda::vector_const_view om, statistic_t statistic); + /** + * @brief Get imaginary frequency discretization of analytic continuation kernel + * on a given set of imaginary frequency indices + * + * @param[in] n Imaginary frequency indices + * @param[in] om Real frequency grid + * @param[in] statistic Particle Statistic: Boson or Fermion + * + * @return Discretization of analytic continuation kernel, with row i + * corresponding to imaginary frequency index n(i) + */ + nda::matrix build_k_if(nda::vector_const_view n, nda::vector_const_view om, statistic_t statistic); + + /** + * @brief Get dyadic fine grid of imaginary frequency indices + * + * The non-negative half of the grid consists of consecutive panels of @p npt + * indices each, where the spacing in panel j is 2^j. The grid is truncated + * at the imaginary frequency cutoff @p nmax, which is always included, and + * is mirrored to negative frequencies. + * + * @param[in] nmax Imaginary frequency cutoff + * @param[in] statistic Particle Statistic: Boson or Fermion + * @param[in] npt Number of indices per dyadic panel + * + * @return Sorted imaginary frequency indices, covering -nmax <= n < nmax + * (fermionic case) or -nmax <= n <= nmax (bosonic case) + * + * \note The grid is symmetric in the sense required by pivrgs_sym: if n is in + * entry i, then -n-1 (fermionic case) or -n (bosonic case) is in entry m-i-1, + * where m is the grid size. For @p npt >= @p nmax the grid contains all + * indices up to the cutoff. + */ + nda::vector build_if_fine(int nmax, statistic_t statistic, int npt); + /** * @brief Get error of fine composite Chebyshev discretization of analytic * continuation kernel in imaginary time diff --git a/c++/cppdlr/dlr_imfreq.cpp b/c++/cppdlr/dlr_imfreq.cpp index 50b2923..5eda594 100644 --- a/c++/cppdlr/dlr_imfreq.cpp +++ b/c++/cppdlr/dlr_imfreq.cpp @@ -34,9 +34,11 @@ namespace cppdlr { cf2if = nda::matrix(niom, r); // Get analytic continuation kernel at DLR frequencies, up to imaginary - // frequency cutoff + // frequency cutoff, on dyadic fine grid with 2r indices per panel auto nmax = fineparams(lambda).nmax; - auto kmat = build_k_if(nmax, dlr_rf, statistic); + auto nf = build_if_fine(nmax, statistic, 2 * r); + if (nf.size() < niom) throw std::runtime_error("Fine imaginary frequency grid smaller than # DLR imaginary frequency nodes."); + auto kmat = build_k_if(nf, dlr_rf, statistic); // Pivoted Gram-Schmidt to obtain DLR imaginary frequency nodes auto [q, norms, piv] = (symmetrize ? pivrgs_sym(kmat, niom) : pivrgs(kmat, 1e-100)); @@ -50,7 +52,7 @@ namespace cppdlr { std::ranges::reverse(piv); // so reverse restores ascending order } - for (int i = 0; i < niom; ++i) { dlr_if(i) = piv(i) - nmax; } + for (int i = 0; i < niom; ++i) { dlr_if(i) = nf(piv(i)); } // Obtain coefficients to imaginary frequency values transformation matrix for (int i = 0; i < niom; ++i) { diff --git a/c++/cppdlr/dlr_kernels.cpp b/c++/cppdlr/dlr_kernels.cpp index 8011e12..c8a5323 100644 --- a/c++/cppdlr/dlr_kernels.cpp +++ b/c++/cppdlr/dlr_kernels.cpp @@ -42,13 +42,13 @@ namespace cppdlr { } } - std::complex k_if_fermion(int n, double om) { return 1.0 / ((2 * n + 1) * pi * 1i - om); } + std::complex k_if_fermion(int n, double om) { return 1.0 / ((2.0 * n + 1) * pi * 1i - om); } std::complex k_if_boson(int n, double om) { if (n == 0 && om == 0.0) { return -0.5; } else { - return std::tanh(0.5 * om) / (2 * n * pi * 1i - om); + return std::tanh(0.5 * om) / (2.0 * n * pi * 1i - om); } } diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index 52b27c9..b1b6d2c 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -357,6 +357,11 @@ namespace cppdlr { * vector and its entry of the returned norms vector are meaningless; only the * pivots should be used. The symmetrized fermionic Matsubara grid is constructed * this way; see imfreq_ops. + * + * \note Pairs of pivots are chosen to maximize the squared volume + * |r_i|^2 |r_i'|^2 - ||^2 spanned by the residuals r_i, r_i' of + * the pair, rather than the sum of their squared norms, as long as the + * number of selected rows does not exceed the column dimension. */ // Type T must be scalar-valued rank 2 array/array_view or matrix/matrix_view @@ -410,13 +415,25 @@ namespace cppdlr { for (int j = firstpair; j < r; j += 2) { + // Score candidate row pairs by the squared volume they add to the + // selected subspace, which penalizes pairs whose rows are nearly + // parallel after orthogonalization. If the pair exceeds the column + // dimension, the volume vanishes for all pairs, so use the sum of + // squared norms instead. + bool usevol = (j + 1 < n); + auto pairnrm = [&](int k) { + if (!usevol) return norms(k) + norms(k + 1); + return norms(k) * norms(k + 1) - std::norm(nda::blas::dotc(aa(k, _), aa(k + 1, _))); + }; + // Find next pair of pivots jpiv = j; - nrm = norms(j) + norms(j + 1); + nrm = pairnrm(j); for (int k = j + 2; k < m; k += 2) { - if (norms(k) + norms(k + 1) > nrm) { + double nk = pairnrm(k); + if (nk > nrm) { jpiv = k; - nrm = norms(k) + norms(k + 1); + nrm = nk; } } diff --git a/test/c++/dlr_build.cpp b/test/c++/dlr_build.cpp index e8f6f6f..c2b25fe 100644 --- a/test/c++/dlr_build.cpp +++ b/test/c++/dlr_build.cpp @@ -18,6 +18,9 @@ #include #include #include +#include +#include +#include using namespace cppdlr; @@ -101,3 +104,72 @@ TEST(dlr_build, symmetric_fixed_points) { EXPECT_EQ(t_non.size(), fine.nt); EXPECT_EQ(w_non.size(), fine.nt); } + +/** +* @brief Test structure of dyadic fine imaginary frequency grid +*/ +TEST(dlr_build, build_if_fine) { + + for (auto statistic : {Fermion, Boson}) { + for (int nmax : {1, 20, 1000, 1000000}) { + for (int npt : {1, 7, 128}) { + + auto n = build_if_fine(nmax, statistic, npt); + int m = n.size(); + + // Covers cutoff range, sorted, symmetric in the sense of pivrgs_sym + int shift = (statistic == Fermion ? 1 : 0); + EXPECT_EQ(n(0), -nmax); + EXPECT_EQ(n(m - 1), nmax - shift); + EXPECT_EQ(m % 2, 1 - shift); + for (int i = 0; i + 1 < m; ++i) { EXPECT_LT(n(i), n(i + 1)); } + for (int i = 0; i < m; ++i) { EXPECT_EQ(n(i) + n(m - 1 - i), -shift); } + + // Spacing on the half with n >= 1 - shift doubles every npt points, except for + // the last step to the cutoff + int i0 = (m + 1) / 2; + for (int i = i0; i + 2 < m; ++i) { EXPECT_EQ(n(i + 1) - n(i), 1 << ((i - i0) / npt)); } + if (m - i0 >= 2) { EXPECT_LE(n(m - 1) - n(m - 2), 1 << ((m - 2 - i0) / npt)); } + + // Full grid if panels cover cutoff range + if (npt >= nmax) { EXPECT_EQ(m, 2 * nmax + 1 - shift); } + } + } + } + + EXPECT_THROW(build_if_fine(0, Fermion, 128), std::runtime_error); + EXPECT_THROW(build_if_fine(100, Fermion, 0), std::runtime_error); +} + +/** +* @brief Test that kernel on given imaginary frequency indices matches kernel +* on full grid +*/ +TEST(dlr_build, build_k_if_grid) { + + int nmax = 300; + auto om = nda::vector{-250.0, -1.0, 0.0, 0.5, 300.0}; + + for (auto statistic : {Fermion, Boson}) { + auto n = build_if_fine(nmax, statistic, 16); + auto kfull = build_k_if(nmax, om, statistic); + auto kgrid = build_k_if(n, om, statistic); + for (int i = 0; i < n.size(); ++i) { EXPECT_EQ(nda::max_element(nda::abs(kgrid(i, nda::range::all) - kfull(n(i) + nmax, nda::range::all))), 0); } + } +} + +/** +* @brief Test analytic continuation kernel in imaginary frequency for indices +* where 2n+1 overflows int +*/ +TEST(dlr_build, k_if_large_n) { + + using std::numbers::pi; + for (int n : {(1 << 30) + 5, -(1 << 30) - 7, 2147483000, -2147483000}) { + double om = 3.0; + auto kf = 1.0 / std::complex(-om, (2.0 * n + 1) * pi); + auto kb = std::tanh(0.5 * om) / std::complex(-om, 2.0 * n * pi); + EXPECT_LE(std::abs(k_if(n, om, Fermion) - kf), 1e-15 * std::abs(kf)); + EXPECT_LE(std::abs(k_if(n, om, Boson) - kb), 1e-15 * std::abs(kb)); + } +} diff --git a/test/c++/imfreq_ops.cpp b/test/c++/imfreq_ops.cpp index 9ca1c19..76c7fef 100644 --- a/test/c++/imfreq_ops.cpp +++ b/test/c++/imfreq_ops.cpp @@ -437,6 +437,53 @@ static void check_h5_roundtrip(imfreq_ops const &ifops, std::string const &filen EXPECT_EQ_ARRAY(ifops.get_if2cf_piv(), ifops_ref.get_if2cf_piv()); } +/** +* @brief Test DLR interpolation in imaginary frequency for large Lambda, for +* all combinations of statistics and symmetrization +*/ +TEST(imfreq_ops, interp_large_lambda) { + + double eps = 1e-10; // DLR tolerance + + for (double lambda : {31623.0, 3162278.0}) { + + // Test poles, spread over [-0.7 lambda, lambda] on a logarithmic scale + auto om = nda::vector(83); + om(0) = 0; + for (int l = 0; l < 41; ++l) { + om(1 + l) = lambda * pow(10.0, -6.0 * l / 40); + om(42 + l) = -0.7 * om(1 + l); + } + + // Test indices: all |n| <= 1000, and logarithmically spaced up to 10 * lambda + auto ntst = std::vector{}; + for (int n = -1001; n <= 1000; ++n) ntst.push_back(n); + for (int i = 0; i <= 2000; ++i) { + int n = static_cast(1000 * pow(10 * lambda / 1000, i / 2000.0)); + ntst.push_back(n); + ntst.push_back(-n - 1); + } + auto nt = nda::vector(ntst.size()); + for (int i = 0; i < nt.size(); ++i) nt(i) = ntst[i]; + + for (bool symmetrize : {false, true}) { + auto dlr_rf = build_dlr_rf(lambda, eps, symmetrize); + for (auto statistic : {Fermion, Boson}) { + + auto ifops = imfreq_ops(lambda, dlr_rf, statistic, symmetrize); + + // Fit poles from values at DLR nodes, and evaluate at test indices + auto gc = ifops.vals2coefs(build_k_if(nda::vector(ifops.get_ifnodes()), om, statistic)); + auto err = nda::max_element(nda::abs(build_k_if(nt, dlr_rf, statistic) * gc - build_k_if(nt, om, statistic))); + + std::cout << fmt::format("Lambda = {:e}, statistic = {}, symmetrize = {}: L^inf err = {:e}\n", lambda, + (statistic == Fermion ? "Fermion" : "Boson"), symmetrize, err); + EXPECT_LT(err, 2 * eps); + } + } + } +} + TEST(dlr_imfreq, h5_rw) { double lambda = 1000; // DLR cutoff diff --git a/test/c++/pivrgs.cpp b/test/c++/pivrgs.cpp index 2940a2e..5e32453 100644 --- a/test/c++/pivrgs.cpp +++ b/test/c++/pivrgs.cpp @@ -231,3 +231,35 @@ TEST(pivrgs, pivrgs_cplx) { EXPECT_EQ(pivthin, arange(r)); EXPECT_LE(frobenius_norm(q - qthin), 1e-14); } + +/** +* @brief Test that rank-r pivrgs_sym chooses mirror pairs by the volume they +* span, not by the sum of their squared norms +*/ +TEST(pivrgs, pivrgs_sym_pair_volume) { + + // Pair (0,3) has the larger norms but nearly parallel rows, pair (1,2) spans + // two dimensions + auto a = nda::matrix{{1.0, 0.0}, {0.6, 0.0}, {0.0, 0.6}, {1.0, 1e-8}}; + + auto [q, norms, piv] = pivrgs_sym(a, 2); + std::sort(piv.begin(), piv.end()); + + EXPECT_EQ(piv, (nda::vector{1, 2})); +} + +/** +* @brief Test that rank-r pivrgs_sym chooses a mirror pair exceeding the column +* dimension by the sum of squared norms +*/ +TEST(pivrgs, pivrgs_sym_overcomplete_pair) { + + // With one column every pair spans zero volume, so the larger pair (1,2) + // must be chosen by its norms + auto a = nda::matrix{{0.5}, {1.0}, {1.0}, {0.5}}; + + auto [q, norms, piv] = pivrgs_sym(a, 2); + std::sort(piv.begin(), piv.end()); + + EXPECT_EQ(piv, (nda::vector{1, 2})); +} From 1101b42b47bee734fdab3485a866c64d56cb31bd Mon Sep 17 00:00:00 2001 From: "Hugo U. R. Strand" Date: Mon, 9 Nov 2026 08:50:53 +0100 Subject: [PATCH 2/3] [test] imaginary frequency cutoff for Lambda > 2^31 fineparams computes nmax = max((int)ceil(lambda), 20). For ceil(lambda) > INT_MAX the conversion to int is undefined behavior: on x86-64 it gives INT_MIN, hence nmax = 20, and imfreq_ops throws 'Fine imaginary frequency grid smaller than # DLR imaginary frequency nodes'; on ARM64 it saturates to INT_MAX and the construction succeeds. Add a test of nmax across 2^31 and a DLR interpolation test at lambda = 2.2e9 with test indices over the int range. Both fail on x86-64 before the fix. Co-Authored-By: Claude Opus 5.5 --- test/c++/dlr_build.cpp | 20 ++++++++++++++++++ test/c++/imfreq_ops.cpp | 45 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 65 insertions(+) diff --git a/test/c++/dlr_build.cpp b/test/c++/dlr_build.cpp index c2b25fe..2b225f2 100644 --- a/test/c++/dlr_build.cpp +++ b/test/c++/dlr_build.cpp @@ -18,7 +18,9 @@ #include #include #include +#include #include +#include #include #include @@ -105,6 +107,24 @@ TEST(dlr_build, symmetric_fixed_points) { EXPECT_EQ(w_non.size(), fine.nt); } +/** +* @brief Test the imaginary frequency cutoff of the fine grid parameters for +* Lambda beyond the range of int +* +* nmax = max(ceil(Lambda), 20) must not overflow. Converting ceil(Lambda) > +* INT_MAX to int is undefined behavior: it gives INT_MIN, i.e. nmax = 20, on +* x86-64, and saturates to INT_MAX on ARM64, so imfreq_ops failed for Lambda > +* 2^31 on x86-64 only. +*/ +TEST(dlr_build, fineparams_nmax_large_lambda) { + + constexpr int int_max = std::numeric_limits::max(); + for (double lambda : {1.0, 19.5, 1000.0, 2.0e9, 2147483646.5, 2147483647.0, 2147483648.0, 2.2e9, 1.0e10, 1.0e12}) { + double expected = std::min(std::max(std::ceil(lambda), 20.0), static_cast(int_max)); + EXPECT_EQ(fineparams(lambda).nmax, static_cast(expected)) << "lambda = " << lambda; + } +} + /** * @brief Test structure of dyadic fine imaginary frequency grid */ diff --git a/test/c++/imfreq_ops.cpp b/test/c++/imfreq_ops.cpp index 76c7fef..9e2f248 100644 --- a/test/c++/imfreq_ops.cpp +++ b/test/c++/imfreq_ops.cpp @@ -437,6 +437,51 @@ static void check_h5_roundtrip(imfreq_ops const &ifops, std::string const &filen EXPECT_EQ_ARRAY(ifops.get_if2cf_piv(), ifops_ref.get_if2cf_piv()); } +/** +* @brief Test DLR interpolation in imaginary frequency for Lambda > 2^31, where +* the imaginary frequency cutoff of the fine grid is clamped to the range of int +* (see dlr_build test fineparams_nmax_large_lambda). Test indices span the int range. +*/ +TEST(imfreq_ops, interp_lambda_above_int_max) { + + double eps = 1e-10; // DLR tolerance + double lambda = 2.2e9; + + // Test poles, spread over [-0.7 lambda, lambda] on a logarithmic scale + auto om = nda::vector(83); + om(0) = 0; + for (int l = 0; l < 41; ++l) { + om(1 + l) = lambda * pow(10.0, -9.0 * l / 40); + om(42 + l) = -0.7 * om(1 + l); + } + + // Test indices: all |n| <= 1000, and logarithmically spaced up to the largest int + auto ntst = std::vector{}; + for (int n = -1001; n <= 1000; ++n) ntst.push_back(n); + for (int i = 0; i <= 2000; ++i) { + int n = static_cast(std::min(1000 * pow(2147483646.0 / 1000, i / 2000.0), 2147483646.0)); + ntst.push_back(n); + ntst.push_back(-n - 1); + } + auto nt = nda::vector(ntst.size()); + for (int i = 0; i < nt.size(); ++i) nt(i) = ntst[i]; + + for (bool symmetrize : {false, true}) { + auto dlr_rf = build_dlr_rf(lambda, eps, symmetrize); + for (auto statistic : {Fermion, Boson}) { + + auto ifops = imfreq_ops(lambda, dlr_rf, statistic, symmetrize); + + auto gc = ifops.vals2coefs(build_k_if(nda::vector(ifops.get_ifnodes()), om, statistic)); + auto err = nda::max_element(nda::abs(build_k_if(nt, dlr_rf, statistic) * gc - build_k_if(nt, om, statistic))); + + std::cout << fmt::format("Lambda = {:e}, statistic = {}, symmetrize = {}: L^inf err = {:e}\n", lambda, + (statistic == Fermion ? "Fermion" : "Boson"), symmetrize, err); + EXPECT_LT(err, 2 * eps); + } + } +} + /** * @brief Test DLR interpolation in imaginary frequency for large Lambda, for * all combinations of statistics and symmetrization From 7aa2a64073982220abfddf09e711897d45419da5 Mon Sep 17 00:00:00 2001 From: "Hugo U. R. Strand" Date: Mon, 9 Nov 2026 08:52:41 +0100 Subject: [PATCH 3/3] [build] clamp the imaginary frequency cutoff to the range of int fineparams set nmax = max((int)ceil(lambda), 20). For lambda > INT_MAX the conversion is undefined behavior: INT_MIN on x86-64, so nmax = 20 and the dyadic fine grid of imfreq_ops had fewer points than DLR nodes, while ARM64 saturates to INT_MAX. Clamp nmax to INT_MAX explicitly, the ARM64 result. The cutoff nu_max = (2 INT_MAX - 1) pi ~ 1.35e10 still exceeds lambda up to ~1.35e10. Fixes dlr_build.fineparams_nmax_large_lambda and imfreq_ops.interp_lambda_above_int_max (L^inf err 2-3e-11 at lambda = 2.2e9, eps = 1e-10). All 137 cppdlr tests pass. Co-Authored-By: Claude Opus 5.5 --- c++/cppdlr/dlr_build.cpp | 19 ++++++++++++++++++- c++/cppdlr/dlr_build.hpp | 2 +- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/c++/cppdlr/dlr_build.cpp b/c++/cppdlr/dlr_build.cpp index 436cb55..18bbc94 100644 --- a/c++/cppdlr/dlr_build.cpp +++ b/c++/cppdlr/dlr_build.cpp @@ -18,6 +18,8 @@ #include "utils.hpp" #include "dlr_kernels.hpp" #include +#include +#include #include using namespace std; @@ -27,12 +29,27 @@ using std::numbers::pi; namespace cppdlr { + namespace { + + // Imaginary frequency cutoff max(ceil(lambda), 20), clamped to the range of int. Converting + // ceil(lambda) > INT_MAX to int is undefined behavior (INT_MIN, hence nmax = 20, on x86-64; + // INT_MAX on ARM64). The clamped cutoff, nu_max = (2 INT_MAX - 1) pi ~ 1.35e10 for fermions, + // still exceeds lambda for lambda up to ~1.35e10. + int nmax_of_lambda(double lambda) { + constexpr int int_max = std::numeric_limits::max(); + double n = std::ceil(lambda); + if (!(n < static_cast(int_max))) return int_max; + return std::max(static_cast(n), 20); + } + + } // namespace + fineparams::fineparams(double lambda, int p) : // TODO: make alt constructor in which iommax is a parameter lambda(lambda), p(p), - nmax(max((int)ceil(lambda), 20)), // TODO: is this a good choice? + nmax(nmax_of_lambda(lambda)), // TODO: is max(ceil(lambda), 20) a good choice? npom(static_cast(max(ceil(log(lambda) / log(2.0)), 1.0))), npt(static_cast(max(ceil(log(lambda) / log(2.0)) - 2, 1.0))), nom(2 * p * npom), diff --git a/c++/cppdlr/dlr_build.hpp b/c++/cppdlr/dlr_build.hpp index 7bb394d..6273a67 100644 --- a/c++/cppdlr/dlr_build.hpp +++ b/c++/cppdlr/dlr_build.hpp @@ -38,7 +38,7 @@ namespace cppdlr { const double lambda; ///< DLR cutoff parameter const int p; ///< Order of composite Chebyshev grid - const int nmax; ///< Imaginary frequency cutoff + const int nmax; ///< Imaginary frequency cutoff, max(ceil(lambda), 20), clamped to INT_MAX const int npom; ///< # fine real frequency grid panels const int npt; ///< # fine imaginary time grid panels const int nom; ///< Total # fine real frequency grid points