PE - BUG - compute_clebsch_displacements does an ideal A-inverse solve (and a Cholesky) on the kinetic path

[matt-pharr]

Summary

compute_clebsch_displacements recomputes the regularized Clebsch component
ξ^α as -A⁻¹(B·xmp1 + C·xsp) for every run, and factorizes A with a
Cholesky. On a kinetic run the A/B/C it reads are the kinetic-modified
matrices, and the reference implementation does something entirely different
there: it performs no linear solve at all.

Surfaced during review of #383. The defect predates that PR — before the
ideal/kinetic split, the same line read the same kinetic-overwritten amats —
so #383 neither introduces nor fixes it. #383 only makes the branch explicit.

Where

src/PerturbedEquilibrium/FieldReconstruction.jl, inside the
Threads.@threads :static loop of compute_clebsch_displacements
(develop: loop at L377, matrix block at L400-411).

ffit.amats(view(amat, :), psi_norm; hint=hint)   # kinetic A on a kinetic run
ffit.bmats(view(bmat, :), psi_norm; hint=hint)
ffit.cmats(view(cmat_buf, :), psi_norm; hint=hint)
...
# amat is positive-definite by construction (Newcomb kinetic-energy form), so cholesky is safe.
amat_fact = cholesky!(Hermitian(amat, :L))

What the reference does

GPEC_fortran/gpec/gpeq.f, gpeq_sol (L45-140). Two things differ.

1. On the kinetic path there is no A solve — ξ_s is read from its stored spline and simply scaled.

IF (kin_flag) THEN
   CALL cspline_eval(u3,psi,0)
   CALL cspline_eval(u4,psi,0)
   xsp1_mn=u3%f
   xss_mn=u4%f                ! stored kinetic ξ_s — A is never touched
ELSE
   ...
   CALL zhetrf('L',mpert,amat,mpert,ipiva,work,mpert*mpert,info)
   CALL zhetrs('L',mpert,mpert,amat,mpert,ipiva,bmat,mpert,info)
   CALL zhetrs('L',mpert,mpert,amat,mpert,ipiva,cmat,mpert,info)
   xss_mn=-MATMUL(bmat,xsp1_mn)-MATMUL(cmat,xsp_mn)
ENDIF

and in the regularization block:

IF (reg_flag) THEN
   xmp1_mn=xsp1_mn*(singfac**2/(singfac**2+reg_spot**2))
   IF (kin_flag) THEN
      xms_mn=xss_mn*(singfac**2/(singfac**2+reg_spot**2))   ! same per-m factor, no solve
   ELSE
      xms_mn=-MATMUL(bmat,xmp1_mn)-MATMUL(cmat,xsp_mn)
   ENDIF

So the correct kinetic result is the stored ξ_s multiplied by the same
per-mode regularization factor already applied to xsp1 — not a fresh solve
against the kinetic A.

2. Even on the ideal path, the reference does not assume positive-definiteness.
It uses zhetrf/zhetrs (Bunch-Kaufman LDLᴴ, valid for indefinite Hermitian
matrices), not zpotrf/zpotrs. Our in-code justification — "amat is
positive-definite by construction (Newcomb kinetic-energy form), so cholesky is
safe" — is an assumption the reference implementation deliberately declines to
make. Cholesky may well be a valid optimization for the ideal A; that should
be established rather than asserted.

Failure modes

With kinetic_factor > 0 and reg_spot > 0, the kinetic A is non-Hermitian
(Hermitian kwmat + anti-Hermitian ktmat). Hermitian(amat, :L) reads only
the lower triangle, so the wrong-Hermiticity case does not error — it
silently symmetrizes and returns a wrong clebsch_alpha, which propagates into
bmt/bmz and the PENTRC Clebsch output. If the symmetrized matrix is also not
positive-definite, it instead throws PosDefException from inside a @threads
region (surfacing as a TaskFailedException).

Both are wrong; the silent one is worse.

Why nothing has caught it

The branch is currently unreachable from any shipped deck. It needs
kinetic_factor > 0 and a [PerturbedEquilibrium] section with
compute_response = true and reg_spot > 0 simultaneously:

deck	kinetic_factor	[PerturbedEquilibrium]
a10_kinetic_example	1.0	absent
Solovev_kinetic_calculated_example	1.0	absent
Solovev_kinetic_NTV_example	0.0	present (reg_spot = 0.05)
all other decks	0.0	—

PE only runs when the section is present, so the two kinetic decks skip PE
entirely, and the one kinetic-NTV deck that runs PE is on the ideal path. This
is latent today and goes live the first time someone runs a genuinely kinetic
case through the perturbed-equilibrium stage — i.e. the intended purpose of the
kinetic NTV pipeline.

Note reg_spot defaults to 0.05, so the guard is not "nobody enables
regularization" — it is only the missing PE section.

Suggested fix

Branch on the model, mirroring gpeq_sol:

• kinetic (mats.kinetic !== nothing): drop the matrix evaluation, the
  factorization and the solve entirely; set
  clebsch_alpha[ipsi, ipert] = xi_s_modes[ipsi, ipert] * reg_factor / chi1
  using the reg_factor already computed for xmp1 in the same loop. xi_s_modes
  is already a parameter and is already the reg_spot == 0 return value, so this
  is close to free and removes per-thread amat/bmat/cmat buffers on that path.
• ideal: keep the solve. Either justify Cholesky properly, or match the
  reference and use a Bunch-Kaufman/bunchkaufman! factorization for safety.

Acceptance criteria

• This is a results-moving change for kinetic + PE runs, so it needs its own
  branch and its own regression report — it must not ride along inside a PR
  claiming inertness.
• Needs a new regression case that actually exercises the path: a deck with
  kinetic_factor > 0, a [PerturbedEquilibrium] section, compute_response = true
  and reg_spot > 0. No existing case covers it.
• Related harness debt noted in Repo - API! - Replace FourFitVars with immutable MatrixSplines split into ideal and kinetic #383: no case exercises kinetic_source = "fixed"
  with kinetic_factor > 0 either.