[Bug report] Symmetrized DLR grids develop a near-null ±ω mode and diverge under repeated nonlinear transforms

### Description

With `symmetrize=true`, the DLR collocation matrices develop a nearly singular mirror-odd mode dominated by the central real-frequency pair `±δ`. In a simple self-consistent scalar GF2 iteration, roundoff injected into this mode first appears as rapidly growing, cancelling DLR coefficients and then causes the reconstructed Green's function/self-energy to diverge. The same calculation with `symmetrize=false` converges to machine precision.

This looks like the cppdlr-side root cause of [TRIQS/tprf#50: Eliashberg solver gives incorrect results for symmetrized DLR meshes](https://github.com/TRIQS/tprf/issues/50). It also agrees with the failing reproducer and analysis already present in [commit 53b6dd9](https://github.com/flatironinstitute/cppdlr/commit/53b6dd92afca5ab061f50d2c0148d2740e59fd6c) on the `tprf_issue` branch. Moreover, the same issue affects the CoQui self-consistent GW calculation with metallic compounds, which is not affected if using IR basis.

### Why GF2, GW, and Eliashberg expose the same cppdlr failure mode

These should not be viewed as three unrelated application-level failures. GF2, self-consistent GW, and the linearized Eliashberg solver all contain the same numerically critical sequence:

1. values on a symmetrized DLR grid are converted to coefficients through the ill-conditioned solve `c = K_sym^{-1} v`;
2. the coefficients are evaluated in another representation or on another grid;
3. pointwise products are formed;
4. the result is projected back to DLR and the operation is repeated.

The near-null mirror-odd mode associated with the central `±δ` pole pair is therefore common to all three calculations. Tiny floating-point components along this mode become large opposite DLR coefficients. They cancel at the original interpolation nodes, but the cancellation is lost after evaluation on another grid and multiplication.

The application-specific products are different,

- GF2: `Sigma(tau) ~ G(tau)^2 G(beta-tau)`;
- GW: `P(tau) ~ G(tau)G(-tau)` and `Sigma(tau) ~ G(tau)W(tau)`;
- Eliashberg: the kernel applies a `G Delta G` product,

but they amplify the same cppdlr projection error. In GF2 and self-consistent GW this appears as growth of DLR coefficients followed by divergence of the self-consistency loop. In the Eliashberg solver it appears as loss of floating-point linearity on non-spectral Krylov vectors, producing spurious and non-reproducible eigenvalues.

This interpretation is consistent with [TRIQS/tprf#50](https://github.com/TRIQS/tprf/issues/50). It is also supported by two independent controls: the metallic CoQui GW calculation remains stable with the IR basis, and the same GF2 equations remain stable with Lehmann.jl's numerically symmetrized `:sym` basis. Thus the common discriminator is the current cppdlr symmetrized-grid construction, not the GF2, GW, or Eliashberg physics.

### Steps to reproduce

Run the following with TRIQS built against cppdlr 1.3.0 (Julia comparison in attachment):

```python
import numpy as np
from triqs.gfs import (
Gf, MeshDLRImFreq, MeshDLRImTime,
make_gf_dlr, make_gf_dlr_imfreq, make_gf_imfreq,
)

BETA, LAMBDA, EPS = 100.0, 1.0e4, 1.0e-13
U, MU, MIXING, NITER = 1.0, 0.15, 0.70, 100

def bethe(z):
root = np.sqrt(z * z - 1.0)
root[np.imag(z) * np.imag(root) < 0.0] *= -1.0
return 2.0 / (z + root)

def run(symmetrize):
iw_mesh = MeshDLRImFreq(
BETA, "Fermion", LAMBDA / BETA, EPS, symmetrize
)
tau_mesh = MeshDLRImTime(
BETA, "Fermion", LAMBDA / BETA, EPS, symmetrize
)
iw = np.array([complex(point) for point in iw_mesh])
tau = np.array([float(point) for point in tau_mesh])
sigma = Gf(mesh=iw_mesh, target_shape=[1, 1])
sigma.data[:] = 0.0

for iteration in range(1, NITER + 1):
g_iw = Gf(mesh=iw_mesh, target_shape=[1, 1])
g_iw.data[:, 0, 0] = bethe(
iw + MU - sigma.data[:, 0, 0]
)
g = make_gf_dlr(g_iw)
g_tau = np.array([g(t)[0, 0] for t in tau])
g_reflected = np.array([g(BETA - t)[0, 0] for t in tau])

sigma_tau = Gf(mesh=tau_mesh, target_shape=[1, 1])
sigma_tau.data[:, 0, 0] = U**2 * g_tau**2 * g_reflected
sigma_new = make_gf_dlr_imfreq(make_gf_dlr(sigma_tau))

old = make_gf_imfreq(make_gf_dlr(sigma), n_iw=256).data[:, 0, 0]
new = make_gf_imfreq(make_gf_dlr(sigma_new), n_iw=256).data[:, 0, 0]
residual = np.max(np.abs(new - old))
sigma.data[:] = (
(1.0 - MIXING) * sigma.data + MIXING * sigma_new.data
)

dense = make_gf_imfreq(make_gf_dlr(sigma), n_iw=256).data[:, 0, 0]
coefficient = np.max(np.abs(make_gf_dlr(sigma).data))
hermiticity = np.max(
np.abs(dense - np.conjugate(dense[::-1]))
)
if iteration in (20, 40, 60, 100):
print(symmetrize, iteration, residual, coefficient, hermiticity)

for symmetric_dlr in (False, True):
run(symmetric_dlr)
```

Observed at iteration 100:

| grid | rank | residual | max DLR coefficient | Hermiticity error |
|---|---:|---:|---:|---:|
| non-symmetrized | 88 | `5.66e-16` | `2.38e-2` | `9.95e-15` |
| symmetrized | 88 | `4.38e3` | `2.34e14` | `2.53e3` |

For the symmetrized run, the coefficient already reaches `9.13e2` at iteration 20 and `3.02e8` at iteration 40, while the non-symmetrized run remains stable.

The failure is not caused by particle-hole asymmetry: repeating the test with `MU=0` also diverges.

### Cross-implementation check: Lehmann.jl

I ported the same GF2 equations and parameters to Julia using [Lehmann.jl](https://github.com/numericalEFT/Lehmann.jl) 0.2.7. The like-for-like comparison is Lehmann's `symmetry=:sym`, which uses a numerically symmetrized grid/basis without imposing physical particle-hole symmetry. Lehmann's `:ph` mode is a different, physically constrained kernel and was not used here (the test has `MU=0.15`).

After 100 iterations:

| Lehmann.jl grid | rank | residual | max DLR coefficient | Hermiticity error |
|---|---:|---:|---:|---:|
| `:none` | 89 | `1.10e-15` | `1.32e-1` | `1.30e-13` |
| `:sym` | 94 | `2.41e-15` | `1.30e-1` | `8.89e-16` |

Neither Lehmann.jl run shows coefficient growth or divergence. The early fixed-point residuals agree with the TRIQS/cppdlr calculation through iteration 15 to the displayed precision, which checks that the GF2 equation, the `G(beta-tau)` sign, and the Bethe square-root branch were ported consistently. This cross-check indicates that the instability is specific to cppdlr's current symmetric-grid construction/conditioning, rather than an unavoidable property of a symmetrized DLR representation.
_see attached file for Julia reproduce_

### Conditioning evidence

For the parameters above, the coefficient-to-imaginary-time matrix has:

| grid | `cond(K_tau)` | smallest singular value |
|---|---:|---:|
| non-symmetrized | `5.02e13` | `4.04e-13` |
| symmetrized | `1.25e16` | `1.62e-15` |

The smallest symmetric singular vector is mirror-odd and dominated by the near-zero pole pair `±0.00130680955`. Large opposite coefficients cancel at the DLR nodes but lose cancellation after evaluation in another representation and nonlinear multiplication.

The effect is not restricted to the tight GF2 parameters. At the smaller parameters used by the Eliashberg reproducer (`beta=2, lambda=10, eps=1e-8`), I obtain `cond(K_tau)=8.75e9` for the symmetric grid versus `2.16e8` for the non-symmetric grid.

A possibly contributing source-level issue is that the tolerance-based `pivrgs_sym` checks the residual norm of the first member of a mirror pair, but normalizes the second member without an equivalent small-norm check. This can retain a nearly dependent second direction. The resulting matrices are then inverted with an unregularized LU solve.

### Expected behavior

Symmetrizing the numerical node set should not introduce a nearly singular mode that destroys repeated DLR transforms or changes a converged physical result.

### Actual behavior

The symmetric basis can hide enormous mirror-odd coefficients behind cancellation at its nodes. Repeated transforms and nonlinear operations amplify this direction until residuals, Hermiticity, and the physical solution diverge.

### Proposed direction

The fixed-point construction described in the `tprf_issue` branch looks appropriate: include each symmetry fixed point as a single node (`omega=0`, `tau=beta/2`, and bosonic `n=0`). This makes the symmetric real-frequency rank odd; fermionic Matsubara nodes can then use `r+1` mirror-paired samples with a least-squares solve.

### Versions

- cppdlr 1.3.0, git hash `5f9adf91d7620ce145ead8e4e906151493038af3`
- TRIQS 4.0.0, git hash `86c3b4d6842a706b040dc5a364671097d8ce2473`
- Python 3.14.6
- macOS 15.7.3, arm64

[gf2_sym_vs_nonsym_lehmann.zip](https://github.com/user-attachments/files/30312417/gf2_sym_vs_nonsym_lehmann.zip)

### Additional Information
This bug detection and comparison is also suggested by [Xiansheng Cai](https://github.com/iintSjds)

Run the provided Python GF2 reproducer against the stated cppdlr version and compare the symmetric and non-symmetric residual, coefficient, and Hermiticity results. Then inspect the `pivrgs_sym` construction and the reproducer analysis on commit 53b6dd9 or the `tprf_issue` branch; done means repeated symmetric DLR transforms remain stable without the near-null mirror-odd mode.

Written by the indexing model from the issue text.