From 41dbede9edb001e8169c462a45282bb883b90f18 Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:31:55 -0400 Subject: [PATCH 1/9] Drop redundant re-initialization of panel edge in build_rf_fine a is already zero from its declaration two lines above. Assisted-by: Claude --- c++/cppdlr/dlr_build.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/c++/cppdlr/dlr_build.cpp b/c++/cppdlr/dlr_build.cpp index 69100fc..2668038 100644 --- a/c++/cppdlr/dlr_build.cpp +++ b/c++/cppdlr/dlr_build.cpp @@ -58,7 +58,6 @@ namespace cppdlr { // 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; From aeb0c3cb1f6203854ca4f9034d5b7ec39515b363 Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:31:09 -0400 Subject: [PATCH 2/9] Take fineparams by const reference in the fine grid builders build_rf_fine, build_it_fine and geterr_k_it only read their fineparams argument. Source-compatible for callers, but it does change the mangled names, so it is an ABI break for anything linking a prebuilt cppdlr. Assisted-by: Claude --- c++/cppdlr/dlr_build.cpp | 6 +++--- c++/cppdlr/dlr_build.hpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/c++/cppdlr/dlr_build.cpp b/c++/cppdlr/dlr_build.cpp index 2668038..a59235e 100644 --- a/c++/cppdlr/dlr_build.cpp +++ b/c++/cppdlr/dlr_build.cpp @@ -42,7 +42,7 @@ namespace cppdlr { if (p <= 0) throw std::runtime_error("Choose p > 0."); } - nda::vector build_rf_fine(fineparams &fine) { + nda::vector build_rf_fine(fineparams const &fine) { int p = fine.p; int npom = fine.npom; @@ -71,7 +71,7 @@ namespace cppdlr { return om; } - std::tuple, nda::vector> build_it_fine(fineparams &fine) { + std::tuple, nda::vector> build_it_fine(fineparams const &fine) { int p = fine.p; int npt = fine.npt; @@ -161,7 +161,7 @@ namespace cppdlr { return kvec; } - std::tuple geterr_k_it(fineparams &fine, nda::vector_const_view t, nda::vector_const_view om, + std::tuple geterr_k_it(fineparams const &fine, nda::vector_const_view t, nda::vector_const_view om, nda::matrix_const_view kmat) { int nt = fine.nt; diff --git a/c++/cppdlr/dlr_build.hpp b/c++/cppdlr/dlr_build.hpp index 16748b9..5d2ac9d 100644 --- a/c++/cppdlr/dlr_build.hpp +++ b/c++/cppdlr/dlr_build.hpp @@ -52,7 +52,7 @@ namespace cppdlr { * * @return Fine real frequency grid */ - nda::vector build_rf_fine(fineparams &fine); + nda::vector build_rf_fine(fineparams const &fine); /** * @brief Get fine composite Legendre grid in imaginary time and corresponding @@ -65,7 +65,7 @@ namespace cppdlr { * * \note Fine imaginary time grid is given in relative format */ - std::tuple, nda::vector> build_it_fine(fineparams &fine); + std::tuple, nda::vector> build_it_fine(fineparams const &fine); /** * @brief Get imaginary time discretization of analytic continuation kernel @@ -134,7 +134,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 geterr_k_it(fineparams &fine, nda::vector_const_view t, nda::vector_const_view om, + std::tuple geterr_k_it(fineparams const &fine, nda::vector_const_view t, nda::vector_const_view om, nda::matrix_const_view kmat); /** From 8c8a3854dc08ab9c0be42f31eeac85c12a97e08c Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:13:34 -0400 Subject: [PATCH 3/9] Rename jstrt to firstpair in fixed-rank pivrgs_sym Same quantity, same name as in the eps-cutoff overload: the row where the first mirror pair begins. Assisted-by: Claude --- c++/cppdlr/utils.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index bbdc2ab..94f37d7 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -392,9 +392,9 @@ namespace cppdlr { } } - int jstrt = (m % 2 == 1) ? 1 : 0; + int firstpair = (m % 2 == 1) ? 1 : 0; // Row where the first mirror pair begins // Then proceed with pivoted GS algorithm as normal - for (int j = jstrt; j < r; j += 2) { + for (int j = firstpair; j < r; j += 2) { // Find next pair of pivots jpiv = j; From 0982d76b8f9cf490eab20b2a80db55ce8e765020 Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:22:33 -0400 Subject: [PATCH 4/9] Drop middle-row-at-the-end leftovers in fixed-rank pivrgs_sym Two commented-out lines and two comments describe an abandoned layout in which the self-paired middle row was placed last; the code places it first. Assisted-by: Claude --- c++/cppdlr/utils.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index 94f37d7..3971211 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -349,7 +349,7 @@ namespace cppdlr { if (r == n + 1 && (n % 2 == 1 || n > m)) { throw std::runtime_error("If r = n+1, n must be even and less than or equal to m."); } // Copy input data, re-ordering rows to make symmetric rows adjacent. If m - // odd, put middle row at the end. + // odd, put middle row first. auto aa = typename T::regular_type(m, n); if (m % 2 == 0) { aa(nda::range(0, m, 2), _) = a(nda::range(0, m / 2), _); @@ -358,7 +358,6 @@ namespace cppdlr { aa(0, _) = a((m - 1) / 2, _); aa(nda::range(1, m, 2), _) = a(nda::range(0, (m - 1) / 2), _); aa(nda::range(2, m, 2), _) = a(nda::range(m - 1, (m - 1) / 2, -1), _); - //aa(m - 1, _) = a((m - 1) / 2, _); } // Compute norms of rows of input matrix @@ -376,10 +375,9 @@ namespace cppdlr { piv(0) = (m - 1) / 2; piv(nda::range(1, m, 2)) = nda::arange(0, (m - 1) / 2); piv(nda::range(2, m, 2)) = nda::arange(m - 1, (m - 1) / 2, -1); - //piv(m - 1) = (m - 1) / 2; } - // If m odd, first choose middle row (now last row) as first pivot + // If m odd, first choose middle row (now the first row) as first pivot if (m % 2 == 1) { // Normalize From 2a07d41b8c90c4be9cd34e188b48726a6be9a77c Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Tue, 4 Aug 2026 15:43:52 -0400 Subject: [PATCH 5/9] symcompare: report the symmetrized rank on the symmetrized line Both symcompare examples printed the unsymmetrized rank r on the "Symmetrized DLR" line, understating or overstating the symmetrized rank whenever the two differ. symcompare_it already computed rsym and used it to size its arrays, so only the format argument was wrong; symcompare_if did not compute it at all. Assisted-by: Claude --- test/c++/symcompare_if.cpp | 5 +++-- test/c++/symcompare_it.cpp | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/test/c++/symcompare_if.cpp b/test/c++/symcompare_if.cpp index 54919e1..ae18be5 100644 --- a/test/c++/symcompare_if.cpp +++ b/test/c++/symcompare_if.cpp @@ -67,7 +67,8 @@ int main() { auto dlr_rf = build_dlr_rf(lambda, eps); auto dlr_rf_sym = build_dlr_rf(lambda, eps, SYM); - int r = dlr_rf.size(); + int r = dlr_rf.size(); + int rsym = dlr_rf_sym.size(); // Get DLR imaginary frequency object auto ifops = imfreq_ops(lambda, dlr_rf, statistic); @@ -106,5 +107,5 @@ int main() { // Print results std::cout << fmt::format("Unsymmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", r, errl2, errlinf); - std::cout << fmt::format("Symmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", r, errl2_sym, errlinf_sym); + std::cout << fmt::format("Symmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", rsym, errl2_sym, errlinf_sym); } \ No newline at end of file diff --git a/test/c++/symcompare_it.cpp b/test/c++/symcompare_it.cpp index 2a3f9cc..39e7c9e 100644 --- a/test/c++/symcompare_it.cpp +++ b/test/c++/symcompare_it.cpp @@ -106,5 +106,5 @@ int main() { // Print results std::cout << fmt::format("Unsymmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", r, errl2, errlinf); - std::cout << fmt::format("Symmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", r, errl2_sym, errlinf_sym); + std::cout << fmt::format("Symmetrized DLR: rank = {}, L^2(tau) err = {:e}, L^inf(tau) err = {:e}\n", rsym, errl2_sym, errlinf_sym); } \ No newline at end of file From b81f8d101d07df9a78af24e00d9ca3e0f8ae4bd8 Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Tue, 4 Aug 2026 15:29:21 -0400 Subject: [PATCH 6/9] pivrgs_sym: support symmetrized selection with a self-symmetric row The symmetrized DLR grid construction needs to select the self-symmetric fixed point of a mirror-symmetric fine grid (omega=0, tau=beta/2) as a single self-paired pivot. Extend pivrgs_sym accordingly: * The eps-cutoff overload accepted only an even number of rows. Allow an odd row count, selecting the middle row (m-1)/2 first as a self-paired pivot and then proceeding in mirror pairs, for an odd epsilon-rank. This mirrors what the fixed-rank overload already does for odd m. The pivot order is built parity-agnostically and the input rows are copied in that order, so the two cannot disagree; results are unchanged for even m. * The fixed-rank overload rejected r = n+1 with n odd. The symmetrized fermionic Matsubara grid needs exactly this: it has no self-symmetric frequency, so it reaches an odd rank n only via n+1 mirror-paired nodes. The remaining checks (r <= min(m, n+1), matching parity of r and m) already cover the meaningful constraints. Assisted-by: Claude --- c++/cppdlr/utils.hpp | 70 +++++++++++++++++++++++++++++++------------- test/c++/pivrgs.cpp | 56 +++++++++++++++++++++++++++++++++++ 2 files changed, 106 insertions(+), 20 deletions(-) diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index 3971211..234000d 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -218,7 +218,11 @@ namespace cppdlr { * * \note The symmetrization condition is that if A(i,:), the ith row of A, is * selected as a pivot, then A(m-i-1,:) is also selected as a pivot. Here, m - * is the row dimension of A, and A is zero-indexed. m must be even. + * is the row dimension of A, and A is zero-indexed. If m is odd, the middle row + * (index (m-1)/2) is its own symmetric partner, and is automatically selected as + * the first pivot, so that the epsilon-rank is odd. Since that pivot is forced, + * the middle row must itself be above the @p eps tolerance; otherwise there is no + * symmetrized basis of odd rank, and an exception is thrown. */ // Type T must be scalar-valued rank 2 array/array_view or matrix/matrix_view @@ -226,15 +230,22 @@ namespace cppdlr { std::tuple, nda::vector> pivrgs_sym(T const &a, double eps) { // Get matrix dimensions - auto [m, n] = a.shape(); - int maxrnk = std::min(m, n); - - if (m % 2 != 0) { throw std::runtime_error("Input matrix must have even number of rows."); } - - // Copy input data, re-ordering rows to make symmetric rows adjacent. - auto aa = typename T::regular_type(m, n); - aa(nda::range(0, m, 2), _) = a(nda::range(0, m / 2), _); - aa(nda::range(1, m, 2), _) = a(nda::range(m - 1, m / 2 - 1, -1), _); + auto [m, n] = a.shape(); + int maxrnk = std::min(m, n); + int nprs = m / 2; // Number of mirror pairs of rows + bool hasmid = (m % 2 == 1); // If m is odd, the middle row is its own mirror + int firstpair = hasmid ? 1 : 0; // Row where the first mirror pair begins + + // Order rows to make symmetric rows adjacent, putting the middle row first + // if m is odd: piv = [(m-1)/2,] 0, m-1, 1, m-2, ... + auto piv = nda::arange(0, m); + if (hasmid) { piv(0) = (m - 1) / 2; } + piv(nda::range(firstpair, m, 2)) = nda::arange(0, nprs); + piv(nda::range(firstpair + 1, m, 2)) = nda::arange(m - 1, m - 1 - nprs, -1); + + // Copy input data in this order + auto aa = typename T::regular_type(m, n); + for (int j = 0; j < m; ++j) { aa(j, _) = a(piv(j), _); } // Compute norms of rows of input matrix, and rescale eps tolerance auto norms = nda::vector(m); @@ -242,17 +253,32 @@ namespace cppdlr { for (int j = 0; j < m; ++j) { norms(j) = normsq(aa(j, _)); } // Begin pivoted double Gram-Schmidt procedure - int jpiv = 0; - double nrm = 0; - auto piv = nda::arange(0, m); - piv(nda::range(0, m, 2)) = nda::arange(0, m / 2); // Re-order pivots to match re-ordered input matrix - piv(nda::range(1, m, 2)) = nda::arange(m - 1, m / 2 - 1, -1); - - if (maxrnk % 2 != 0) { // If n < m and n is odd, decrease maxrnk to maintain symmetry - maxrnk -= 1; + int jpiv = 0; + double nrm = 0; + + // maxrnk must have the same parity as m, since rows are selected in pairs, + // preceded by the middle row if m is odd + if (maxrnk % 2 != m % 2) { maxrnk -= 1; } + + // Choose the middle row (now the first row) as the first pivot. It is forced + // rather than chosen by norm, so a negligible middle row says nothing about the + // remaining rows: reporting an epsilon-rank of zero would be wrong, and no + // symmetrized basis of odd rank exists. + if (hasmid) { + if (norms(0) <= epssq) { throw std::runtime_error("Middle row of input matrix must be above the eps tolerance, since it is a forced pivot."); } + + // Normalize middle row + aa(0, _) /= sqrt(norms(0)); + + // Orthogonalize remaining rows against middle row + for (int k = 1; k < m; ++k) { + if (norms(k) <= epssq) { continue; } // Can skip rows with norm less than tolerance + aa(k, _) = aa(k, _) - aa(0, _) * nda::blas::dotc(aa(0, _), aa(k, _)); + norms(k) = normsq(aa(k, _)); + } } - for (int j = 0; j < maxrnk; j += 2) { + for (int j = firstpair; j < maxrnk; j += 2) { // Find next pair of pivots jpiv = j; @@ -334,6 +360,11 @@ namespace cppdlr { * is the row dimension of A, and A is zero-indexed. A can have an odd number of * rows if and only if @p r is odd, and in this case the middle row (index * (m-1)/2) of A is automatically selected as a pivot. + * + * \note @p r may be as large as n+1, in which case the last selected row is + * linearly dependent and only the pivots are meaningful. This is needed for a + * matrix with no self-symmetric row, which can only reach an odd rank n by + * selecting n+1 mirror-paired rows. */ // Type T must be scalar-valued rank 2 array/array_view or matrix/matrix_view @@ -346,7 +377,6 @@ namespace cppdlr { if (m % 2 == 1 && r % 2 == 0) { throw std::runtime_error("If input matrix has odd number of rows, r must be odd."); } if (r % 2 == 1 && m % 2 == 0) { throw std::runtime_error("If r is odd, input matrix must have odd number of rows."); } if (r > m || r > n + 1) { throw std::runtime_error("r must be less than or equal to min(m,n+1)."); } - if (r == n + 1 && (n % 2 == 1 || n > m)) { throw std::runtime_error("If r = n+1, n must be even and less than or equal to m."); } // Copy input data, re-ordering rows to make symmetric rows adjacent. If m // odd, put middle row first. diff --git a/test/c++/pivrgs.cpp b/test/c++/pivrgs.cpp index 4c972d1..2940a2e 100644 --- a/test/c++/pivrgs.cpp +++ b/test/c++/pivrgs.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include using namespace cppdlr; @@ -99,6 +100,61 @@ TEST(pivrgs, pivrgs_real) { EXPECT_LE(frobenius_norm(q - qthin), 1e-14); } +/** + * Test that the eps-cutoff symmetrized Gram-Schmidt pivrgs_sym handles an odd + * number of rows: no throw, the central row (m-1)/2 is selected as a + * self-paired pivot giving an odd rank, the remaining pivots come in mirror + * pairs, and the returned rows are orthonormal. + */ +TEST(pivrgs, pivrgs_sym_odd) { + + int m = 51; // Odd number of rows + int n = 40; + double eps = 1e-6; + + // Random matrix with rapidly decaying column norms, so numerically low-rank + auto a = nda::matrix::rand(m, n); + for (int j = 0; j < n; ++j) { a(_, j) *= pow(2.0, -j); } + + // Symmetrized pivoted Gram-Schmidt with eps cutoff (must not throw for odd m) + auto [q, norms, piv] = pivrgs_sym(a, eps); + int r = norms.size(); + + // Rank is odd, and below the number of columns + EXPECT_EQ(r % 2, 1); + EXPECT_LT(r, n); + + // The central row (m-1)/2, its own mirror, is selected as the first pivot + EXPECT_EQ(piv(0), (m - 1) / 2); + + // The remaining pivots come in mirror pairs: i selected => m-1-i selected + auto selected = std::vector(m, false); + for (auto p : piv) { selected[p] = true; } + for (auto p : piv) { EXPECT_TRUE(selected[m - 1 - p]); } + + // Returned rows are orthonormal to near double precision + EXPECT_LE(frobenius_norm(eye(r) - q * transpose(q)), 1e-13); +} + +/** + * Test that the eps-cutoff pivrgs_sym rejects an odd-row matrix whose middle row + * is negligible. That row is a forced pivot rather than one chosen by norm, so + * there is no symmetrized basis of odd rank, however large the other rows are. + */ +TEST(pivrgs, pivrgs_sym_negligible_middle_row) { + + int m = 51; // Odd number of rows + int n = 40; + double eps = 1e-6; + + auto a = nda::matrix::rand(m, n); + for (int j = 0; j < n; ++j) { a(_, j) *= pow(2.0, -j); } + EXPECT_NO_THROW(pivrgs_sym(a, eps)); + + a((m - 1) / 2, _) = 0; + EXPECT_THROW(pivrgs_sym(a, eps), std::runtime_error); +} + /** * Test pivoted reorthogonalized Gram-Schmidt function pivrgs for complex-valued * matrix From 450b92a1a7241bd7c87f0623f86f3a97596fa33a Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:24:30 -0400 Subject: [PATCH 7/9] Build the pivot order parity-agnostically in fixed-rank pivrgs_sym The row ordering was spelled out twice, once for the row copy and once for the pivot vector, each split into an even-m and an odd-m branch. Build the pivot order once and copy the rows through it, as the eps-cutoff overload already does, so the two cannot disagree. Bit-for-bit identical output. Assisted-by: Claude --- c++/cppdlr/utils.hpp | 35 ++++++++++++++--------------------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index 234000d..e285a03 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -378,17 +378,20 @@ namespace cppdlr { if (r % 2 == 1 && m % 2 == 0) { throw std::runtime_error("If r is odd, input matrix must have odd number of rows."); } if (r > m || r > n + 1) { throw std::runtime_error("r must be less than or equal to min(m,n+1)."); } - // Copy input data, re-ordering rows to make symmetric rows adjacent. If m - // odd, put middle row first. + int nprs = m / 2; // Number of mirror pairs of rows + bool hasmid = (m % 2 == 1); // If m is odd, the middle row is its own mirror + int firstpair = hasmid ? 1 : 0; // Row where the first mirror pair begins + + // Order rows to make symmetric rows adjacent, putting the middle row first + // if m is odd: piv = [(m-1)/2,] 0, m-1, 1, m-2, ... + auto piv = nda::arange(0, m); + if (hasmid) { piv(0) = (m - 1) / 2; } + piv(nda::range(firstpair, m, 2)) = nda::arange(0, nprs); + piv(nda::range(firstpair + 1, m, 2)) = nda::arange(m - 1, m - 1 - nprs, -1); + + // Copy input data in this order auto aa = typename T::regular_type(m, n); - if (m % 2 == 0) { - aa(nda::range(0, m, 2), _) = a(nda::range(0, m / 2), _); - aa(nda::range(1, m, 2), _) = a(nda::range(m - 1, m / 2 - 1, -1), _); - } else { - aa(0, _) = a((m - 1) / 2, _); - aa(nda::range(1, m, 2), _) = a(nda::range(0, (m - 1) / 2), _); - aa(nda::range(2, m, 2), _) = a(nda::range(m - 1, (m - 1) / 2, -1), _); - } + for (int j = 0; j < m; ++j) { aa(j, _) = a(piv(j), _); } // Compute norms of rows of input matrix auto norms = nda::vector(m); @@ -397,19 +400,10 @@ namespace cppdlr { // Begin pivoted double Gram-Schmidt procedure int jpiv = 0; double nrm = 0; - auto piv = nda::arange(0, m); - if (m % 2 == 0) { - piv(nda::range(0, m, 2)) = nda::arange(0, m / 2); // Re-order pivots to match re-ordered input matrix - piv(nda::range(1, m, 2)) = nda::arange(m - 1, m / 2 - 1, -1); - } else { - piv(0) = (m - 1) / 2; - piv(nda::range(1, m, 2)) = nda::arange(0, (m - 1) / 2); - piv(nda::range(2, m, 2)) = nda::arange(m - 1, (m - 1) / 2, -1); - } // If m odd, first choose middle row (now the first row) as first pivot - if (m % 2 == 1) { + if (hasmid) { // Normalize aa(0, _) /= sqrt(normsq(aa(0, _))); @@ -420,7 +414,6 @@ namespace cppdlr { } } - int firstpair = (m % 2 == 1) ? 1 : 0; // Row where the first mirror pair begins // Then proceed with pivoted GS algorithm as normal for (int j = firstpair; j < r; j += 2) { From 9435f7236c2a488e068707f76d2270223e2f64ab Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Wed, 5 Aug 2026 11:27:05 -0400 Subject: [PATCH 8/9] Extract the normalize-and-project step in both pivrgs_sym overloads Each overload spelled out "normalize row j, then orthogonalize every later row against it and update its norm" three times: once for the self-paired middle row and twice per selected mirror pair. Extract it as a lambda. Bit-for-bit identical output. Assisted-by: Claude --- c++/cppdlr/utils.hpp | 105 +++++++++++++++++-------------------------- 1 file changed, 40 insertions(+), 65 deletions(-) diff --git a/c++/cppdlr/utils.hpp b/c++/cppdlr/utils.hpp index e285a03..e1e90d8 100644 --- a/c++/cppdlr/utils.hpp +++ b/c++/cppdlr/utils.hpp @@ -229,8 +229,10 @@ namespace cppdlr { template T, nda::Scalar S = nda::get_value_t> std::tuple, nda::vector> pivrgs_sym(T const &a, double eps) { - // Get matrix dimensions - auto [m, n] = a.shape(); + // Get matrix dimensions. Not a structured binding, since clang before 19 cannot + // capture one in the normalize_and_project lambda below when OpenMP is enabled. + long m = a.extent(0); + long n = a.extent(1); int maxrnk = std::min(m, n); int nprs = m / 2; // Number of mirror pairs of rows bool hasmid = (m % 2 == 1); // If m is odd, the middle row is its own mirror @@ -252,9 +254,16 @@ namespace cppdlr { double epssq = eps * eps; for (int j = 0; j < m; ++j) { norms(j) = normsq(aa(j, _)); } - // Begin pivoted double Gram-Schmidt procedure - int jpiv = 0; - double nrm = 0; + // Normalize row j, given its squared norm, and orthogonalize the rows after it + // against it, updating their norms + auto normalize_and_project = [&](int j, double nrmsq) { + aa(j, _) /= sqrt(nrmsq); + for (int k = j + 1; k < m; ++k) { + if (norms(k) <= epssq) { continue; } // Can skip rows with norm less than tolerance + aa(k, _) = aa(k, _) - aa(j, _) * nda::blas::dotc(aa(j, _), aa(k, _)); + norms(k) = normsq(aa(k, _)); + } + }; // maxrnk must have the same parity as m, since rows are selected in pairs, // preceded by the middle row if m is odd @@ -266,18 +275,13 @@ namespace cppdlr { // symmetrized basis of odd rank exists. if (hasmid) { if (norms(0) <= epssq) { throw std::runtime_error("Middle row of input matrix must be above the eps tolerance, since it is a forced pivot."); } - - // Normalize middle row - aa(0, _) /= sqrt(norms(0)); - - // Orthogonalize remaining rows against middle row - for (int k = 1; k < m; ++k) { - if (norms(k) <= epssq) { continue; } // Can skip rows with norm less than tolerance - aa(k, _) = aa(k, _) - aa(0, _) * nda::blas::dotc(aa(0, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } + normalize_and_project(0, norms(0)); } + // Begin pivoted double Gram-Schmidt procedure + int jpiv = 0; + double nrm = 0; + for (int j = firstpair; j < maxrnk; j += 2) { // Find next pair of pivots @@ -309,29 +313,13 @@ namespace cppdlr { // (not including current row) if (nrm <= epssq) { return {aa(nda::range(0, j), _), norms(nda::range(0, j)), piv(nda::range(0, j))}; }; - // Normalize current row - aa(j, _) /= sqrt(nrm); - - // Orthogonalize remaining rows against current row - for (int k = j + 1; k < m; ++k) { - if (norms(k) <= epssq) { continue; } // Can skip rows with norm less than tolerance - aa(k, _) = aa(k, _) - aa(j, _) * nda::blas::dotc(aa(j, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } + normalize_and_project(j, nrm); // Orthogonalize current row (now the second chosen pivot row) against all // previously chosen rows for (int k = 0; k < j + 1; ++k) { aa(j + 1, _) = aa(j + 1, _) - aa(k, _) * nda::blas::dotc(aa(k, _), aa(j + 1, _)); } - // Normalize current row - aa(j + 1, _) /= sqrt(normsq(aa(j + 1, _))); - - // Orthogonalize remaining rows against current row - for (int k = j + 2; k < m; ++k) { - if (norms(k) <= epssq) { continue; } // Can skip rows with norm less than tolerance - aa(k, _) = aa(k, _) - aa(j + 1, _) * nda::blas::dotc(aa(j + 1, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } + normalize_and_project(j + 1, normsq(aa(j + 1, _))); } return {aa(nda::range(maxrnk), _), norms(nda::range(maxrnk)), piv(nda::range(maxrnk))}; @@ -371,8 +359,10 @@ namespace cppdlr { template T, nda::Scalar S = nda::get_value_t> std::tuple, nda::vector> pivrgs_sym(T const &a, int r) { - // Get matrix dimensions - auto [m, n] = a.shape(); + // Get matrix dimensions. Not a structured binding, since clang before 19 cannot + // capture one in the normalize_and_project lambda below when OpenMP is enabled. + long m = a.extent(0); + long n = a.extent(1); if (m % 2 == 1 && r % 2 == 0) { throw std::runtime_error("If input matrix has odd number of rows, r must be odd."); } if (r % 2 == 1 && m % 2 == 0) { throw std::runtime_error("If r is odd, input matrix must have odd number of rows."); } @@ -397,24 +387,23 @@ namespace cppdlr { auto norms = nda::vector(m); for (int j = 0; j < m; ++j) { norms(j) = normsq(aa(j, _)); } - // Begin pivoted double Gram-Schmidt procedure - int jpiv = 0; - double nrm = 0; + // Normalize row j and orthogonalize the rows after it against it, updating + // their norms + auto normalize_and_project = [&](int j) { + aa(j, _) /= sqrt(normsq(aa(j, _))); + for (int k = j + 1; k < m; ++k) { + aa(k, _) = aa(k, _) - aa(j, _) * nda::blas::dotc(aa(j, _), aa(k, _)); + norms(k) = normsq(aa(k, _)); + } + }; // If m odd, first choose middle row (now the first row) as first pivot + if (hasmid) { normalize_and_project(0); } - if (hasmid) { - // Normalize - aa(0, _) /= sqrt(normsq(aa(0, _))); - - // Orthogonalize remaining rows against current row - for (int k = 1; k < m; ++k) { - aa(k, _) = aa(k, _) - aa(0, _) * nda::blas::dotc(aa(0, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } - } + // Begin pivoted double Gram-Schmidt procedure + int jpiv = 0; + double nrm = 0; - // Then proceed with pivoted GS algorithm as normal for (int j = firstpair; j < r; j += 2) { // Find next pair of pivots @@ -439,27 +428,13 @@ namespace cppdlr { // previously chosen rows for (int k = 0; k < j; ++k) { aa(j, _) = aa(j, _) - aa(k, _) * nda::blas::dotc(aa(k, _), aa(j, _)); } - // Normalize current row - aa(j, _) /= sqrt(normsq(aa(j, _))); - - // Orthogonalize remaining rows against current row - for (int k = j + 1; k < m; ++k) { - aa(k, _) = aa(k, _) - aa(j, _) * nda::blas::dotc(aa(j, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } + normalize_and_project(j); // Orthogonalize current row (now the second chosen pivot row) against all // previously chosen rows for (int k = 0; k < j + 1; ++k) { aa(j + 1, _) = aa(j + 1, _) - aa(k, _) * nda::blas::dotc(aa(k, _), aa(j + 1, _)); } - // Normalize current row - aa(j + 1, _) /= sqrt(normsq(aa(j + 1, _))); - - // Orthogonalize remaining rows against current row - for (int k = j + 2; k < m; ++k) { - aa(k, _) = aa(k, _) - aa(j + 1, _) * nda::blas::dotc(aa(j + 1, _), aa(k, _)); - norms(k) = normsq(aa(k, _)); - } + normalize_and_project(j + 1); } return {aa(nda::range(r), _), norms(nda::range(r)), piv(nda::range(r))}; From 0c4dc711f40dc28e71916432b253a91fa1311e9f Mon Sep 17 00:00:00 2001 From: Nils Wentzell Date: Tue, 4 Aug 2026 12:26:52 -0400 Subject: [PATCH 9/9] Include self-symmetric fixed points in symmetrized DLR grids The symmetrized construction selected nodes only in exact mirror pairs, so the real-frequency grid had to take a near-degenerate pair +/-delta straddling omega=0, whose basis columns 1/(i*nu_n -/+ delta) are nearly parallel. This ill-conditioned the DLR fit matrices: at eps=1e-10, cond(cf2it) was 2.5e13 and 1.1e14 for lambda=10 and 1000, against ~1e10 unsymmetrized. Harmless on DLR-representable inputs, but workflows applying a DLR-built operator to non-spectral vectors (e.g. the Krylov iterates of an eigensolver) suffer O(1) floating-point error. Instead always include the self-symmetric fixed point as a single self-paired node: omega=0 in real frequency, tau=beta/2 in imaginary time, n=0 in the bosonic Matsubara grid. Conditioning becomes 4.4e10 and 6.8e10 for the same two cases, ranks stay within one of the unsymmetrized rank, and accuracy on spectral inputs is unchanged. Non-symmetrized grids are unchanged bit for bit. Consequences: - The symmetrized DLR rank r is now odd. - The symmetrized fermionic Matsubara grid has no self-symmetric frequency, so it decouples to an even niom = r + 1 of mirror pairs, inverted by least squares (as the symmetrized bosonic case already was). - imtime_ops now rejects an even-rank dlr_rf under SYM, which every symmetrized grid from cppdlr <= 1.3.0 has, so a persisted grid fed back with SYM hits it. Updates the symmetrized DLR grids section of doc/background.rst, which documented the previous pair-only behavior. Still pending: triqs_tprf's DLR Eliashberg tests remain pinned to non-symmetrized meshes (symmetrize=False in test/python/eliashberg/dlr_eliashberg_solver.py), since they were failing with the symmetrized ones. They should be revisited against this change. Assisted-by: Claude --- c++/cppdlr/dlr_build.cpp | 51 ++++++++++++++++++-------- c++/cppdlr/dlr_build.hpp | 12 +++++- c++/cppdlr/dlr_imfreq.cpp | 14 ++++--- c++/cppdlr/dlr_imfreq.hpp | 7 ++-- c++/cppdlr/dlr_imtime.cpp | 6 ++- c++/cppdlr/dlr_imtime.hpp | 4 ++ doc/background.rst | 29 +++++++++++---- test/c++/dlr_build.cpp | 41 +++++++++++++++++++++ test/c++/imfreq_ops.cpp | 77 +++++++++++++++++++++++++++------------ test/c++/imtime_ops.cpp | 53 ++++++++++++++++++++------- 10 files changed, 222 insertions(+), 72 deletions(-) diff --git a/c++/cppdlr/dlr_build.cpp b/c++/cppdlr/dlr_build.cpp index a59235e..36d96fa 100644 --- a/c++/cppdlr/dlr_build.cpp +++ b/c++/cppdlr/dlr_build.cpp @@ -42,7 +42,7 @@ namespace cppdlr { if (p <= 0) throw std::runtime_error("Choose p > 0."); } - nda::vector build_rf_fine(fineparams const &fine) { + nda::vector build_rf_fine(fineparams const &fine, bool symmetrize) { int p = fine.p; int npom = fine.npom; @@ -50,9 +50,13 @@ namespace cppdlr { 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(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(ntot); double a = 0, b = 0; @@ -60,18 +64,22 @@ namespace cppdlr { 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> build_it_fine(fineparams const &fine) { + std::tuple, nda::vector> build_it_fine(fineparams const &fine, bool symmetrize) { int p = fine.p; int npt = fine.npt; @@ -79,10 +87,16 @@ namespace cppdlr { 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(fine.nt); - auto w = nda::vector(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(ntot); + auto w = nda::vector(ntot); double a = 0, b = 0; @@ -95,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}; } @@ -262,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 diff --git a/c++/cppdlr/dlr_build.hpp b/c++/cppdlr/dlr_build.hpp index 5d2ac9d..20c8649 100644 --- a/c++/cppdlr/dlr_build.hpp +++ b/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 build_rf_fine(fineparams const &fine); + nda::vector 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> build_it_fine(fineparams const &fine); + std::tuple, nda::vector> build_it_fine(fineparams const &fine, bool symmetrize = false); /** * @brief Get imaginary time discretization of analytic continuation kernel @@ -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 build_dlr_rf(double lambda, double eps, bool symmetrize); diff --git a/c++/cppdlr/dlr_imfreq.cpp b/c++/cppdlr/dlr_imfreq.cpp index a76307a..50b2923 100644 --- a/c++/cppdlr/dlr_imfreq.cpp +++ b/c++/cppdlr/dlr_imfreq.cpp @@ -25,9 +25,11 @@ namespace cppdlr { imfreq_ops::imfreq_ops(double lambda, nda::vector_const_view 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(niom); cf2if = nda::matrix(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(cf2if); if2cf.piv = nda::vector(r); lapack::getrf(if2cf.lu, if2cf.piv); diff --git a/c++/cppdlr/dlr_imfreq.hpp b/c++/cppdlr/dlr_imfreq.hpp index 4b2a864..9da09d2 100644 --- a/c++/cppdlr/dlr_imfreq.hpp +++ b/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 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 dlr_rf; ///< DLR frequencies nda::vector dlr_if; ///< DLR imaginary frequency nodes nda::matrix cf2if; /// Transformation matrix from DLR coefficients to values at DLR imaginary frequency nodes diff --git a/c++/cppdlr/dlr_imtime.cpp b/c++/cppdlr/dlr_imtime.cpp index 4845129..4c89af1 100644 --- a/c++/cppdlr/dlr_imtime.cpp +++ b/c++/cppdlr/dlr_imtime.cpp @@ -27,6 +27,10 @@ namespace cppdlr { imtime_ops::imtime_ops(double lambda, nda::vector_const_view 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(r); cf2it = nda::matrix(r, r); it2cf.lu = nda::matrix(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 diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp index e44deec..29884c8 100644 --- a/c++/cppdlr/dlr_imtime.hpp +++ b/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 dlr_rf, bool symmetrize); diff --git a/doc/background.rst b/doc/background.rst index 950a240..cca3397 100644 --- a/doc/background.rst +++ b/doc/background.rst @@ -308,12 +308,25 @@ functionality via symmetrization flags. Please see the :ref:`list of other cppdlr capabilities` section on the :ref:`examples page` for a list of `cppdlr` tests which showcase this functionality. -We make a note about symmetrization for bosonic Green's functions. In this case, -we always include :math:`i \nu_n = 0` as a DLR imaginary frequency node. Since -we do not include the point :math:`\omega = 0` in the :math:`r` symmetrized DLR -frequencies, :math:`r` is even in this case, and the extra imaginary frequency -grid points makes the number of DLR imaginary frequency nodes :math:`r+1`. To -obtain the DLR coefficients from samples of the Green's function at the DLR -imaginary frequency nodes, we therefore solve a slightly overdetermined +A symmetrized grid includes its self-symmetric fixed point---the point which is +its own mirror image---as a single self-paired node: :math:`\omega = 0` among the +DLR frequencies, and :math:`\tau = \beta/2` among the DLR imaginary time nodes. +The symmetrized DLR rank :math:`r` is therefore odd. Selecting nodes in mirror +pairs alone would instead force the frequency grid onto a nearly degenerate pair +:math:`\pm \delta` straddling :math:`\omega = 0`, whose kernel columns are nearly +parallel, badly ill-conditioning the DLR fit matrices. + +We make a note about the symmetrized imaginary frequency grids. In the bosonic +case, :math:`i \nu_n = 0` is the self-symmetric Matsubara frequency, and is +included as the central node; there are then :math:`r` DLR imaginary frequency +nodes, and the values-to-coefficients map is a square :math:`r \times r` linear +system. In the fermionic case there is no self-symmetric Matsubara frequency, +since :math:`\nu_n = (2n+1) \pi / \beta` never vanishes, so a symmetrized +fermionic grid consists of mirror pairs only and we use :math:`r+1` imaginary +frequency nodes. To obtain the DLR coefficients from samples of the Green's +function at these nodes, we therefore solve a slightly overdetermined (:math:`(r+1) \times r`, rather than :math:`r \times r`) linear system using the -least squares method. \ No newline at end of file +least squares method. Consequently, evaluating the resulting expansion back at +the DLR imaginary frequency nodes recovers the samples to a relative accuracy of +order :math:`\epsilon` in general, rather than to machine precision as in the +square case. \ No newline at end of file diff --git a/test/c++/dlr_build.cpp b/test/c++/dlr_build.cpp index 68176c2..e8f6f6f 100644 --- a/test/c++/dlr_build.cpp +++ b/test/c++/dlr_build.cpp @@ -60,3 +60,44 @@ TEST(dlr_build, get_kfine) { std::cout << fmt::format("Max imag time err = {:e}, Max freq err = {:e}\n", errt, errom); } + +/** +* @brief Test that the symmetric DLR construction includes the self-symmetric fixed +* points, omega=0 and tau=beta/2, and leaves the non-symmetric grids unchanged. +*/ +TEST(dlr_build, symmetric_fixed_points) { + + double eps = 1e-10; + + for (double lambda : {10.0, 100.0, 1000.0}) { + auto dlr_rf = build_dlr_rf(lambda, eps, SYM); + long r = dlr_rf.size(); + + // Symmetric rank is odd, with omega=0 as the central self-paired node + EXPECT_EQ(r % 2, 1); + EXPECT_DOUBLE_EQ(dlr_rf((r - 1) / 2), 0.0); + + // Grid is mirror-symmetric about omega=0: omega_j = -omega_{r-1-j}, exactly, since + // the negative half of the fine grid is built by negating the positive half + for (int j = 0; j < r / 2; ++j) { EXPECT_EQ(dlr_rf(j), -dlr_rf(r - 1 - j)); } + } + + // Symmetric fine grids: omega=0 and tau=beta/2 (relative t=0.5) are present as + // central nodes, with one extra point relative to the even default grid. + fineparams fine(100.0); + + auto om_sym = build_rf_fine(fine, SYM); + EXPECT_EQ(om_sym.size(), fine.nom + 1); + EXPECT_DOUBLE_EQ(om_sym(fine.nom / 2), 0.0); + + auto [t_sym, w_sym] = build_it_fine(fine, SYM); + EXPECT_EQ(t_sym.size(), fine.nt + 1); + EXPECT_DOUBLE_EQ(t_sym(fine.nt / 2), 0.5); + EXPECT_DOUBLE_EQ(w_sym(fine.nt / 2), 0.0); + + // Non-symmetric grids are unchanged (even count, no central fixed-point node) + auto [t_non, w_non] = build_it_fine(fine); + EXPECT_EQ(build_rf_fine(fine).size(), fine.nom); + EXPECT_EQ(t_non.size(), fine.nt); + EXPECT_EQ(w_non.size(), fine.nt); +} diff --git a/test/c++/imfreq_ops.cpp b/test/c++/imfreq_ops.cpp index d0cc420..9ca1c19 100644 --- a/test/c++/imfreq_ops.cpp +++ b/test/c++/imfreq_ops.cpp @@ -262,26 +262,38 @@ TEST(imfreq_ops, interp_matrix_sym_fer) { auto dlr_rf = build_dlr_rf(lambda, eps, SYM); int r = dlr_rf.size(); - // Verify symmetry - EXPECT_EQ(max_element(abs(dlr_rf(range(r / 2)) + dlr_rf(range(r - 1, r / 2 - 1, -1)))), 0); + // Verify DLR rank is odd, with omega=0 as the central self-paired node and the + // remaining frequencies in mirror pairs about it + EXPECT_EQ(r % 2, 1); + EXPECT_EQ(dlr_rf((r - 1) / 2), 0.0); + for (int j = 0; j < r / 2; ++j) { EXPECT_EQ(dlr_rf(j), -dlr_rf(r - 1 - j)); } // Get DLR imaginary frequency object auto ifops = imfreq_ops(lambda, dlr_rf, statistic, SYM); // Sample Green's function G at DLR imaginary frequency nodes auto const &dlr_if = ifops.get_ifnodes(); + int niom = dlr_if.size(); - // Verify symmetry - EXPECT_EQ(max_element(abs(2 * dlr_if(range(r / 2)) + 1 + 2 * dlr_if(range(r - 1, r / 2 - 1, -1)) + 1)), 0); + // The fermionic Matsubara grid has no self-symmetric frequency, so it consists of + // mirror pairs only and uses the even niom = r + 1 rather than the odd r. + EXPECT_EQ(niom, r + 1); - auto g = nda::array(r, norb, norb); - for (int i = 0; i < r; ++i) { g(i, _, _) = gfun(norb, beta, dlr_if(i), statistic); } + // Verify symmetry: each mirror pair satisfies nu_n + nu_n' = 0, i.e. n + n' = -1 + for (int i = 0; i < niom / 2; ++i) { EXPECT_EQ(dlr_if(i) + dlr_if(niom - 1 - i), -1); } + + auto g = nda::array(niom, norb, norb); + for (int i = 0; i < niom; ++i) { g(i, _, _) = gfun(norb, beta, dlr_if(i), statistic); } // DLR coefficients of G auto gc = ifops.vals2coefs(beta, g); - // Check that G can be recovered at imaginary frequency nodes - EXPECT_LT(max_element(abs(ifops.coefs2vals(beta, gc) - g)), 3e-13); + // Check that G can be recovered at imaginary frequency nodes. The fermionic + // symmetric grid has no self-symmetric frequency, so it is over-determined + // (niom = r + 1) and vals2coefs is a least-squares fit rather than an exact + // interpolation; recovery is therefore accurate to the DLR tolerance eps, not + // machine precision. + EXPECT_LT(max_element(abs(ifops.coefs2vals(beta, gc) - g)), 10 * eps); // Compute error in imaginary frequency auto gtru = nda::matrix(norb, norb); @@ -400,21 +412,11 @@ TEST(imfreq_ops, interp_matrix_sym_bos) { std::cout << fmt::format("Imag time: L^2 err = {:e}, L^inf err = {:e}\n", errl2, errlinf); } -TEST(dlr_imfreq, h5_rw) { - - double lambda = 1000; // DLR cutoff - double eps = 1e-10; // DLR tolerance - auto statistic = Fermion; // Fermionic Green's function - - // Get DLR frequencies - auto dlr_rf = build_dlr_rf(lambda, eps); - - // Get DLR imaginary frequency object - auto ifops = imfreq_ops(lambda, dlr_rf, statistic); - - auto filename = "data_imfreq_ops_h5_rw.h5"; - auto name = "ifops"; +// Write ifops to an HDF5 file, read it back, and check that every stored member +// survives the round trip +static void check_h5_roundtrip(imfreq_ops const &ifops, std::string const &filename) { + auto name = "ifops"; { h5::file file(filename, 'w'); h5::write(file, name, ifops); @@ -426,7 +428,6 @@ TEST(dlr_imfreq, h5_rw) { h5::read(file, name, ifops_ref); } - // Check equal EXPECT_EQ(ifops.lambda(), ifops_ref.lambda()); EXPECT_EQ(ifops.rank(), ifops_ref.rank()); EXPECT_EQ_ARRAY(ifops.get_rfnodes(), ifops_ref.get_rfnodes()); @@ -436,6 +437,36 @@ TEST(dlr_imfreq, h5_rw) { EXPECT_EQ_ARRAY(ifops.get_if2cf_piv(), ifops_ref.get_if2cf_piv()); } +TEST(dlr_imfreq, h5_rw) { + + double lambda = 1000; // DLR cutoff + double eps = 1e-10; // DLR tolerance + + auto dlr_rf = build_dlr_rf(lambda, eps); + + check_h5_roundtrip(imfreq_ops(lambda, dlr_rf, Fermion), "data_imfreq_ops_h5_rw.h5"); +} + +/** +* @brief Test HDF5 round-trip for symmetrized DLR imaginary frequency objects. +* +* The symmetrized fermionic grid is over-determined (niom = r + 1), so no +* if2cf LU factorization is stored (get_if2cf_lu/piv are empty); this checks +* that the empty factors serialize and deserialize cleanly. +*/ +TEST(dlr_imfreq, h5_rw_sym) { + + double lambda = 1000; // DLR cutoff + double eps = 1e-10; // DLR tolerance + + auto dlr_rf = build_dlr_rf(lambda, eps, SYM); + + for (auto statistic : {Fermion, Boson}) { + SCOPED_TRACE(statistic == Fermion ? "Fermion" : "Boson"); + check_h5_roundtrip(imfreq_ops(lambda, dlr_rf, statistic, SYM), "data_imfreq_ops_h5_rw_sym.h5"); + } +} + /** * @brief Test direct application of if2it matrix against standard vals2coefs -> coefs2vals * pipeline, for all statistic/symmetrization combinations diff --git a/test/c++/imtime_ops.cpp b/test/c++/imtime_ops.cpp index c43818b..bd60333 100644 --- a/test/c++/imtime_ops.cpp +++ b/test/c++/imtime_ops.cpp @@ -870,6 +870,18 @@ TEST(imtime_ops, refl_matrix) { std::cout << fmt::format("Imag time: L^2 err = {:e}, L^inf err = {:e}\n", errl2, errlinf); } +// Check that a symmetrized DLR node set is mirror-symmetric about its center, the +// self-symmetric fixed point of the grid (@p center is 0 for frequencies, 0.5 for +// imaginary times). In both conventions the mirror of a node is its negation. +static void check_sym_nodes(nda::vector_const_view x, double center) { + + int r = x.size(); + + EXPECT_EQ(r % 2, 1); + EXPECT_EQ(x((r - 1) / 2), center); + for (int j = 0; j < r / 2; ++j) { EXPECT_EQ(x(j), -x(r - 1 - j)); } +} + /** * @brief Test symmetrized DLR interpolation and evaluation for fermionic * matrix-valued Green's function @@ -891,11 +903,8 @@ TEST(imtime_ops, interp_matrix_sym_fer) { auto dlr_rf = build_dlr_rf(lambda, eps, SYM); int r = dlr_rf.size(); - // Verify DLR rank is even - EXPECT_EQ(r % 2, 0); - - // Verify symmetry - EXPECT_EQ(max_element(abs(dlr_rf(range(r / 2)) + dlr_rf(range(r - 1, r / 2 - 1, -1)))), 0); + // Verify the DLR frequencies are mirror-symmetric about the omega=0 node + check_sym_nodes(dlr_rf, 0.0); // Get DLR imaginary time object auto itops = imtime_ops(lambda, dlr_rf, SYM); @@ -903,8 +912,8 @@ TEST(imtime_ops, interp_matrix_sym_fer) { // Sample Green's function G at DLR imaginary time nodes auto const &dlr_it = itops.get_itnodes(); - // Verify symmetry - EXPECT_EQ(max_element(abs(dlr_it(range(r / 2)) + dlr_it(range(r - 1, r / 2 - 1, -1)))), 0); + // Verify the DLR imaginary time nodes are mirror-symmetric about tau=beta/2 + check_sym_nodes(dlr_it, 0.5); auto g = nda::array(r, norb, norb); for (int i = 0; i < r; ++i) { g(i, _, _) = gfun(norb, beta, dlr_it(i)); } @@ -974,11 +983,8 @@ TEST(imtime_ops, interp_matrix_sym_bos) { auto dlr_rf = build_dlr_rf(lambda, eps, SYM); int r = dlr_rf.size(); - // Verify DLR rank is even - EXPECT_EQ(r % 2, 0); - - // Verify symmetry - EXPECT_EQ(max_element(abs(dlr_rf(range(r / 2)) + dlr_rf(range(r - 1, r / 2 - 1, -1)))), 0); + // Verify the DLR frequencies are mirror-symmetric about the omega=0 node + check_sym_nodes(dlr_rf, 0.0); // Get DLR imaginary time object auto itops = imtime_ops(lambda, dlr_rf, SYM); @@ -986,8 +992,8 @@ TEST(imtime_ops, interp_matrix_sym_bos) { // Obtain DLR imaginary time nodes auto const &dlr_it = itops.get_itnodes(); - // Verify symmetry - EXPECT_EQ(max_element(abs(dlr_it(range(r / 2)) + dlr_it(range(r - 1, r / 2 - 1, -1)))), 0); // r even + // Verify the DLR imaginary time nodes are mirror-symmetric about tau=beta/2 + check_sym_nodes(dlr_it, 0.5); // Sample Green's function at DLR nodes auto g = nda::array(r, norb, norb); @@ -1037,6 +1043,25 @@ TEST(imtime_ops, interp_matrix_sym_bos) { std::cout << fmt::format("Imag freq: l^2 err = {:e}, L^inf err = {:e}\n", errl2, errlinf); } +/** +* @brief Test that the SYM constructor rejects a DLR frequency grid of even rank +* +* The symmetrized selection takes the self-paired tau=beta/2 node plus mirror +* pairs, so an even-rank input has no valid symmetrized selection. +*/ +TEST(imtime_ops, sym_requires_odd_rank) { + + double lambda = 1000; // DLR cutoff + double eps = 1e-10; // DLR tolerance + + // Trim a non-symmetrized grid to even rank + auto dlr_rf = build_dlr_rf(lambda, eps); + int r = dlr_rf.size() - dlr_rf.size() % 2; + + EXPECT_THROW(imtime_ops(lambda, dlr_rf(range(r)), SYM), std::runtime_error); + EXPECT_NO_THROW(imtime_ops(lambda, dlr_rf(range(r)), NONSYM)); +} + /** * @brief Test inner product of two DLR expansions */