Include self-symmetric fixed points in symmetrized DLR grids

c++/cppdlr/dlr_build.cpp

@@ -42,48 +42,61 @@ namespace cppdlr {
     if (p <= 0) throw std::runtime_error("Choose p > 0.");
   }
 
-  nda::vector<double> build_rf_fine(fineparams &fine) {
+  nda::vector<double> build_rf_fine(fineparams const &fine, bool symmetrize) {
 
     int p    = fine.p;
     int npom = fine.npom;
 
     auto bc = barycheb(p);             // Get barycheb object for Chebyshev nodes
     auto xc = (bc.getnodes() + 1) / 2; // Cheb nodes on [0,1]
 
-    // Real frequency grid points
+    // Real frequency grid points, mirror-symmetric about zero: om(j) = -om(n-1-j).
+    // The symmetric case additionally includes omega=0, which is its own mirror.
 
-    auto om = nda::vector<double>(fine.nom);
+    int nhalf = npom * p;                       // # points on (0,lambda)
+    int upper = symmetrize ? nhalf + 1 : nhalf; // first index of the positive half
+    int ntot  = upper + nhalf;                  // Total # points
+    auto om   = nda::vector<double>(ntot);
 
     double a = 0, b = 0;
 
     // Points on (0,lambda)
 
-    a = 0.0;
     for (int i = 0; i < npom; ++i) {
       b                                             = fine.lambda / pow(2.0, npom - i - 1);
-      om(range((npom + i) * p, (npom + i + 1) * p)) = a + (b - a) * xc;
+      om(range(upper + i * p, upper + (i + 1) * p)) = a + (b - a) * xc;
       a                                             = b;
     }
 
-    // Points on (-lambda,0)
+    // Central point omega=0, its own mirror
 
-    om(range(0, npom * p)) = -om(range(2 * npom * p - 1, npom * p - 1, -1));
+    if (symmetrize) om(nhalf) = 0.0;
+
+    // Points on (-lambda,0): mirror of the positive half
+
+    om(range(0, nhalf)) = -om(range(ntot - 1, upper - 1, -1));
 
     return om;
   }
 
-  std::tuple<nda::vector<double>, nda::vector<double>> build_it_fine(fineparams &fine) {
+  std::tuple<nda::vector<double>, nda::vector<double>> build_it_fine(fineparams const &fine, bool symmetrize) {
 
     int p   = fine.p;
     int npt = fine.npt;
 
     auto [xgl, wgl] = gaussquad(p);  // Gauss-Legendre nodes and weights on [-1,1]
     xgl             = (xgl + 1) / 2; // Transform to [0,1]
 
-    // Imaginary time grid points
+    // Imaginary time grid points, mirror-symmetric about beta/2: t(i) = -t(n-1-i) in
+    // relative format. The symmetric case additionally includes tau=beta/2, its own
+    // mirror since k_it(0.5,om) = k_it(-0.5,om), with zero quadrature weight so the
+    // quadrature rule is unchanged.
 
-    auto t = nda::vector<double>(fine.nt);
-    auto w = nda::vector<double>(fine.nt);
+    int nhalf = npt * p;                        // # points on (0,1/2)
+    int upper = symmetrize ? nhalf + 1 : nhalf; // first index of the (1/2,1) half
+    int ntot  = upper + nhalf;                  // Total # points
+    auto t    = nda::vector<double>(ntot);
+    auto w    = nda::vector<double>(ntot);
 
     double a = 0, b = 0;
 
@@ -96,10 +109,17 @@ namespace cppdlr {
       a                            = b;
     }
 
-    // Points on (1/2,1) in relative format
+    // Central point tau=beta/2, carrying no quadrature weight
+
+    if (symmetrize) {
+      t(nhalf) = 0.5;
+      w(nhalf) = 0.0;
+    }
+
+    // Points on (1/2,1) in relative format: mirror of the first half
 
-    t(range(npt * p, 2 * npt * p)) = -t(range(npt * p - 1, -1, -1));
-    w(range(npt * p, 2 * npt * p)) = w(range(npt * p - 1, -1, -1));
+    t(range(upper, ntot)) = -t(range(nhalf - 1, -1, -1));
+    w(range(upper, ntot)) = w(range(nhalf - 1, -1, -1));
 
     return {t, w};
   }
@@ -162,7 +182,7 @@ namespace cppdlr {
     return kvec;
   }
 
-  std::tuple<double, double> geterr_k_it(fineparams &fine, nda::vector_const_view<double> t, nda::vector_const_view<double> om,
+  std::tuple<double, double> geterr_k_it(fineparams const &fine, nda::vector_const_view<double> t, nda::vector_const_view<double> om,
                                          nda::matrix_const_view<double> kmat) {
 
     int nt   = fine.nt;
@@ -263,8 +283,8 @@ namespace cppdlr {
     auto fine = fineparams(lambda);
 
     // Get fine grids in frequency and imaginary time
-    auto [t, w] = build_it_fine(fine);
-    auto om     = build_rf_fine(fine);
+    auto [t, w] = build_it_fine(fine, symmetrize);
+    auto om     = build_rf_fine(fine, symmetrize);
 
     // Get discretization of analytic continuation kernel on fine grids (the K
     // matrix), weighted in the time variable by the square root of the

c++/cppdlr/dlr_build.hpp

@@ -49,23 +49,27 @@ namespace cppdlr {
   * @brief Build fine composite Chebyshev grid in real frequency 
   *
   * @param[in] fine Fine grid parameters
+  * @param[in] symmetrize If true, include omega=0, making the grid odd-sized and
+  * mirror-symmetric about it.
   *
   * @return Fine real frequency grid
   */
-  nda::vector<double> build_rf_fine(fineparams &fine);
+  nda::vector<double> build_rf_fine(fineparams const &fine, bool symmetrize = false);
 
   /**
   * @brief Get fine composite Legendre grid in imaginary time and corresponding
   * square root quadrature weights
   *
   * @param[in] fine Fine grid parameters
+  * @param[in] symmetrize If true, include tau=beta/2 (relative t=0.5) with zero
+  * quadrature weight, making the grid odd-sized and mirror-symmetric about it.
   *
   * @return Tuple containing fine imaginary time grid and corresponding square
   * root quadrature weights
   *
   * \note Fine imaginary time grid is given in relative format 
   */
-  std::tuple<nda::vector<double>, nda::vector<double>> build_it_fine(fineparams &fine);
+  std::tuple<nda::vector<double>, nda::vector<double>> build_it_fine(fineparams const &fine, bool symmetrize = false);
 
   /**
   * @brief Get imaginary time discretization of analytic continuation kernel
@@ -134,7 +138,7 @@ namespace cppdlr {
   * \note \p kmat should be computed using the function get_kfine with composite
   * Chebyshev grids produced by get_omfine and get_tfine
   */
-  std::tuple<double, double> geterr_k_it(fineparams &fine, nda::vector_const_view<double> t, nda::vector_const_view<double> om,
+  std::tuple<double, double> geterr_k_it(fineparams const &fine, nda::vector_const_view<double> t, nda::vector_const_view<double> om,
                                          nda::matrix_const_view<double> kmat);
 
   /**
@@ -152,6 +156,10 @@ namespace cppdlr {
   *
   * @return DLR frequencies
   *
+  * @note With SYM the DLR frequencies are mirror-symmetric about omega=0, which is
+  * itself included as a self-paired node, so the rank is odd. Without it a
+  * symmetric grid needs a near-degenerate pole pair straddling zero, which
+  * ill-conditions the DLR fit.
   */
   nda::vector<double> build_dlr_rf(double lambda, double eps, bool symmetrize);
 

c++/cppdlr/dlr_imfreq.cpp

@@ -25,9 +25,11 @@ namespace cppdlr {
   imfreq_ops::imfreq_ops(double lambda, nda::vector_const_view<double> dlr_rf, statistic_t statistic, bool symmetrize)
      : lambda_(lambda), statistic(statistic), r(dlr_rf.size()), dlr_rf(dlr_rf) {
 
-    // Get # DLR imaginary frequency nodes; for symmetrized bosonic case, this
-    // is DLR rank + 1, otherwise it is DLR rank
-    niom   = (statistic == Boson && symmetrize) ? r + 1 : r;
+    // # DLR imaginary frequency nodes. The symmetrized grid is mirror-symmetric
+    // about i*nu=0. The bosonic grid contains the self-paired node n=0, so niom = r,
+    // which is odd. The fermionic grid has no self-paired node, so it consists of
+    // mirror pairs only and needs the even niom = r + 1.
+    niom   = (symmetrize && statistic == Fermion) ? r + 1 : r;
     dlr_if = nda::vector<int>(niom);
     cf2if  = nda::matrix<dcomplex>(niom, r);
 
@@ -55,9 +57,9 @@ namespace cppdlr {
       for (int j = 0; j < r; ++j) { cf2if(i, j) = kmat(piv(i), j); }
     }
 
-    if (!(symmetrize && statistic == Boson)) {
-      // Prepare imaginary time values to coefficients transformation by computing
-      // LU factors of coefficient to imaginary time matrix
+    if (niom == r) {
+      // Square case: LU factors invert cf2if. The over-determined case stores no
+      // factors and uses a least-squares solve instead (see vals2coefs).
       if2cf.lu  = nda::matrix<dcomplex>(cf2if);
       if2cf.piv = nda::vector<int>(r);
       lapack::getrf(if2cf.lu, if2cf.piv);

c++/cppdlr/dlr_imfreq.hpp

@@ -47,8 +47,9 @@ namespace cppdlr {
     * @param[in] symmetrize NONSYM or false for non-symmetrized DLR frequencies,
     * SYM or true for symmetrized
     *
-    * @note In case Boson and SYM options are selected, we enforce that i*nu_n=0 is
-    * chosen as a DLR imaginary frequency node.
+    * @note With SYM the DLR rank r is odd. Boson includes the self-symmetric node
+    * i*nu_n=0, giving a square niom = r. Fermion has no such node, so it uses
+    * niom = r + 1 mirror pairs and an over-determined map, inverted by least squares.
     */
     imfreq_ops(double lambda, nda::vector_const_view<double> dlr_rf, statistic_t statistic, bool symmetrize = false);
 
@@ -245,7 +246,7 @@ namespace cppdlr {
     double lambda_;                   ///< Energy cutoff divided by temperature
     statistic_t statistic;            ///< Particle statistic: Fermion or Boson
     int r;                            ///< DLR rank
-    int niom;                         ///< # DLR imaginary freq nodes (different from r in symmetrized bosonic case)
+    int niom;                         ///< # DLR imaginary freq nodes (= r + 1 in the symmetrized fermionic case, else = r)
     nda::vector<double> dlr_rf;       ///< DLR frequencies
     nda::vector<int> dlr_if;          ///< DLR imaginary frequency nodes
     nda::matrix<nda::dcomplex> cf2if; /// Transformation matrix from DLR coefficients to values at DLR imaginary frequency nodes

c++/cppdlr/dlr_imtime.cpp

@@ -27,6 +27,10 @@ namespace cppdlr {
 
   imtime_ops::imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, bool symmetrize) : lambda_(lambda), r(dlr_rf.size()), dlr_rf(dlr_rf) {
 
+    // The symmetrized selection always takes the self-paired tau=beta/2 node plus
+    // mirror pairs, so it can only produce an odd number of nodes.
+    if (symmetrize && r % 2 == 0) throw std::runtime_error("Symmetrized DLR frequency grid must have odd rank.");
+
     dlr_it    = nda::vector<double>(r);
     cf2it     = nda::matrix<double>(r, r);
     it2cf.lu  = nda::matrix<double>(r, r);
@@ -35,7 +39,7 @@ namespace cppdlr {
     // Get discretization of analytic continuation kernel on fine grid in
     // imaginary time, at DLR frequencies
     auto fine   = fineparams(lambda);
-    auto [t, w] = build_it_fine(fine);
+    auto [t, w] = build_it_fine(fine, symmetrize);
     auto kmat   = build_k_it(t, dlr_rf);
 
     // Pivoted Gram-Schmidt to obtain DLR imaginary time nodes

c++/cppdlr/dlr_imtime.hpp

@@ -51,6 +51,10 @@ namespace cppdlr {
     * @param[in] dlr_rf DLR frequencies
     * @param[in] symmetrize NONSYM or false for non-symmetrized DLR frequencies,
     * SYM or true for symmetrized
+    *
+    * @note With SYM the nodes are mirror-symmetric about tau=beta/2, which is
+    * itself included as a self-paired node. SYM therefore requires @p dlr_rf of odd
+    * rank, as built by build_dlr_rf with SYM, and throws otherwise.
     */
     imtime_ops(double lambda, nda::vector_const_view<double> dlr_rf, bool symmetrize);