[imfreq] select DLR Matsubara nodes on a dyadic fine grid

c++/cppdlr/dlr_build.cpp

@@ -18,6 +18,8 @@
 #include "utils.hpp"
 #include "dlr_kernels.hpp"
 #include <algorithm>
+#include <cmath>
+#include <limits>
 #include <numbers>
 
 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<int>::max();
+      double n              = std::ceil(lambda);
+      if (!(n < static_cast<double>(int_max))) return int_max;
+      return std::max(static_cast<int>(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<int>(max(ceil(log(lambda) / log(2.0)), 1.0))),
        npt(static_cast<int>(max(ceil(log(lambda) / log(2.0)) - 2, 1.0))),
        nom(2 * p * npom),
@@ -273,6 +290,47 @@ namespace cppdlr {
     }
   }
 
+  nda::matrix<dcomplex> build_k_if(nda::vector_const_view<int> n, nda::vector_const_view<double> om, statistic_t statistic) {
+
+    auto kmat = nda::matrix<dcomplex>(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<int> 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<int>{};
+    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<int>(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<int>(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<double> build_dlr_rf(double lambda, double eps, bool symmetrize) {
 
     if (eps < 1e-14) {

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
@@ -119,6 +119,41 @@ namespace cppdlr {
   */
   nda::matrix<dcomplex> build_k_if(int nmax, nda::vector_const_view<double> 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<dcomplex> build_k_if(nda::vector_const_view<int> n, nda::vector_const_view<double> 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<int> 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

c++/cppdlr/dlr_imfreq.cpp

@@ -34,9 +34,11 @@ namespace cppdlr {
     cf2if  = nda::matrix<dcomplex>(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) {

c++/cppdlr/dlr_kernels.cpp

@@ -42,13 +42,13 @@ namespace cppdlr {
     }
   }
 
-  std::complex<double> k_if_fermion(int n, double om) { return 1.0 / ((2 * n + 1) * pi * 1i - om); }
+  std::complex<double> k_if_fermion(int n, double om) { return 1.0 / ((2.0 * n + 1) * pi * 1i - om); }
 
   std::complex<double> 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);
     }
   }
 

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 - |<r_i, r_i'>|^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;
         }
       }
 

test/c++/dlr_build.cpp

@@ -18,6 +18,11 @@
 #include <cppdlr/dlr_build.hpp>
 #include <cppdlr/dlr_kernels.hpp>
 #include <fmt/format.h>
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <complex>
+#include <numbers>
 
 using namespace cppdlr;
 
@@ -101,3 +106,90 @@ TEST(dlr_build, symmetric_fixed_points) {
   EXPECT_EQ(t_non.size(), fine.nt);
   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<int>::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<double>(int_max));
+    EXPECT_EQ(fineparams(lambda).nmax, static_cast<int>(expected)) << "lambda = " << lambda;
+  }
+}
+
+/**
+* @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<double>{-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<double>(-om, (2.0 * n + 1) * pi);
+    auto kb   = std::tanh(0.5 * om) / std::complex<double>(-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));
+  }
+}