sign-canonicalize dlr and if frequencies for determinsitic grid selection

Author: Wentzell

c++/cppdlr/dlr_build.cpp

@@ -17,6 +17,7 @@
 #include "dlr_build.hpp"
 #include "utils.hpp"
 #include "dlr_kernels.hpp"
+#include <algorithm>
 #include <numbers>
 
 using namespace std;
@@ -280,6 +281,14 @@ namespace cppdlr {
     auto omega = nda::vector<double>(r);
     for (int i = 0; i < r; ++i) { omega(i) = om(piv(i)); }
 
+    // Sign-canonicalize for deterministic grid selection: orient the grid so that the
+    // most positive frequency is the one with the largest absolute value.
+    // This helps with reproducibility across different BLAS backends.
+    if (std::abs(nda::min_element(omega)) > std::abs(nda::max_element(omega))) {
+      omega *= -1;
+      std::ranges::reverse(omega);
+    }
+
     return omega;
   }
 

c++/cppdlr/dlr_imfreq.cpp

@@ -39,6 +39,19 @@ namespace cppdlr {
     // Pivoted Gram-Schmidt to obtain DLR imaginary frequency nodes
     auto [q, norms, piv] = (symmetrize ? pivrgs_sym(kmat, niom) : pivrgs(kmat, 1e-100));
     std::sort(piv.begin(), piv.end()); // Sort pivots in ascending order
+
+    // Sign-canonicalize for deterministic grid selection (see build_dlr_rf): orient the node
+    // set so its most positive Matsubara frequency has the largest absolute value, flipping
+    // nu_n -> -nu_n (row p -> max_piv - p) otherwise. No effect on grid quality.
+    long min_if = nda::min_element(piv) - nmax; // index n of the most negative selected node
+    long max_if = nda::max_element(piv) - nmax; // index n of the most positive selected node
+    long nu_off = (statistic == Fermion);       // Matsubara frequency is 2n+1 (Fermion) / 2n (Boson)
+    if (std::abs(2 * max_if + nu_off) < std::abs(2 * min_if + nu_off)) {
+      int max_piv = static_cast<int>(kmat.extent(0)) - 1; // mirror row of p under nu_n -> -nu_n
+      for (auto &p : piv) p = max_piv - p;                // ascending -> descending,
+      std::ranges::reverse(piv);                          // so reverse restores ascending order
+    }
+
     for (int i = 0; i < niom; ++i) { dlr_if(i) = piv(i) - nmax; }
 
     // Obtain coefficients to imaginary frequency values transformation matrix