GPEC, DCON, RDCON, STRIDE - BUGFIX - Fix LAPACK and ZVODE type and size errors

[matt-pharr]

Fixes a number of semantic errors according to the LAPACK documentation. Different implementations of LAPACK may or may not be unsafe to these errors, but some certainly caused undefined behavior on my laptop:

• Places where LWORK was set larger than the size of the ZHEEV work matrix, allowing it to write into undefined memory
• Places where real or complex arguments were complex or real
• Places where matrices used in ZHEEV were reused later (ZHEEV destroys the matrix input while computing eigenvalues)
• A wmat that was missed in a past variable renaming circa 2020 which messed up W_xe_eigenvalues in gpec

These instances are easily discernible from each other by looking at the small code diff. This does fix a number of behaviors, including the eigenvalues of the reluctance, the vacuum matrix, and eigenvalues of the inverse plasma energy matrix.

Here is Claude's list of what changed:

What changed (11 files, 56 insertions, 21 deletions):

#	Fix	Files
1	wmat typo, so W_xe_eigenvalue/W_xe_amp are real numbers again	gpec/gpout.f
2	permeabsv declared REAL(r8), so the s_P column has all singular values at full precision	gpec/gpglobal.F, gpec/gpout.f
3	rdcon ev is now real; the netcdf W_v_eigenvalue is 1-D like dcon's after #256	rdcon/free.f, rdcon/rdcon_netcdf.f, rdcon/dcon.f
4	evalsi computed from an independently built inverse response matrix	dcon/ode_output.f
5	riccati_eig copies X before ZHEEV	stride/riccati.f
6	gpdiag ev/evo declared real	gpec/gpdiag.f
7	lwork=2*mpert-1 reset before ZHEEV	dcon/free.f, rdcon/free.f
8	Complex zero passed as BETA to zgbmv/zhbmv	stride/sing.F
9	Return early when ZGETRF reports a singular matrix	dcon/ode_output.f, rdcon/ode_output.f

[matt-pharr]

Thanks to @StuartBenjamin for pointing out that this issue was more widespread than previously thought

[logan-nc]

I am surprised to see this needed here since I see no new zheev in the diff and even a Xtmp=X guard added. Is this fixing a bug introduced by another PR? If so, it needs to be clearly marked as urgent!

It seems like that style of guard is the better way to go rather than introducing repeat wp calculations that are hard to maintain.

Even the rdcon diff is different from the dcon one already. Inevitably someone will change the original wp calculation location and not realize they need to go fix these other locations as well.

[matt-pharr]

We could do instead an Xtmp=X type thing, yes, and I see your point why that one would be better. Though I am not sure I understand the rest of your comment about another PR. The ZHEEV statement that had the bug looks like it was added 10 years ago when you added RDCON.

[logan-nc]

@claude review this PR

[claude]

Claude finished @logan-nc's task in 46s —— View job

---

PR review: LAPACK and ZVODE type and size fixes

I read the full diff. I haven't compiled it or run any cases, so the review is from reading the code only. The fixes look correct, and I found nothing blocking. Two points below are worth acting on.

Reply to the inline question on dcon/ode_output.f:295

This is not a regression from another PR. In the base branch, ode_output_get_evals declares evalsi and indexi and uses them in the out_evals ASCII write (dcon/ode_output.f:319-320). Nothing in that routine ever assigned evalsi or indexi, so the eval and err columns printed uninitialized memory. The Xtmp=X guard in riccati.f is a different fix: it protects an INTENT(IN) argument. It does not replace this block.

This is not urgent. It only affects the diagnostic out_evals output, not the stability results.

Changes I checked

• gpec/gpout.f: wmat = MATMUL(MATMUL(ptof,wvecs),ptof) becomes wmat. This is clearly the intended fix. wvecs was stale from the earlier renaming.
• permeabsv is now REAL(r8) in gpglobal.F, and the REAL() wrapper is dropped in gpout.f. The two files are consistent. Please grep for any other use of permeabsv, such as netcdf writers or gpdiag, that still assumes complex.
• ev is now real in rdcon/free.f, rdcon_netcdf.f, and gpdiag.f. This matches the zheev signature. The netcdf W_v_eigenvalue is now 1-D and written as ev. The W_v_eigenvalue shape change is a user-visible output change. It matches dcon after Add W_v output eigenvalues fix #256, but it should be noted in docs/releases.rst.
• lwork=2*mpert-1 before zheev in the dcon and rdcon free.f files: correct. lwork had been reused from the earlier zgeev calls, which use a larger value.
• stride/riccati.f: the Xtemp copy before zheev is correct.
• stride/sing.F: zero is complex, and it was already used at lines 1175 and 1193 as BETA. The diff only adds the declaration, and no other lines needed to change.

Suggestions

1. Maintainability (the reviewer's concern). The inverse response matrix is now built by a second block that re-derives it from u. The dcon and rdcon copies already differ in how they handle early returns. Two options:

   • Keep a copy before zheev, for example wptmp=wp, and invert that copy with zgetrf/zgetri or zgetrs. This ties the inverse to the original wp and avoids a duplicated formula.
   • Factor the construction into a small shared helper that both get_evals routines call.

   ode_output_get_crit already repeats the same pattern at dcon/ode_output.f:441-446. A single helper would cover all three places.

2. Early RETURN on a singular zgetrf. IF(info /= 0)RETURN silently skips the evals ASCII and binary output for that step. That is safer than using a garbage factorization. Consider a WRITE/warning when out_evals is on, so missing rows are explained. The second RETURN is unreachable in practice if the first one passed. It is only reachable if u(:,:,2) is singular while u(:,:,1) is not, which is the case the guard is for. It is fine as written.

3. Release notes and validation. Add a note to docs/releases.rst about the changed W_v_eigenvalue shape and the corrected W_xe_*, s_P, reluctance, and inverse-energy eigenvalue outputs. Please also confirm that regression benchmarks (DIII-D example) were re-run, since several outputs change.

CLAUDE.md compliance: I did not check every changed .f line against the 72-column limit, so please run awk 'length > 72' on the changed .f files before merging. The new lines in dcon/ode_output.f and rdcon/ode_output.f look to be within the limit.

[logan-nc]

@matt-pharr ignore the claude item 3. Item 1 is overkill I think since we plan to stop actively developing fortran anyways and the wptmp approach still requires inversion (doesn't save many lines / opportunities for mistakes).

Item 2 would be nice if you have time. If not, then it's not blocking.

[logan-nc]

@matt-pharr I leave this to you to hit "merge" when happy