[StuartBenjamin]

RDCON's Galerkin solve keeps every element matrix (cell%mat, about 256 mpert² bytes per cell) in memory while it builds the banded global matrix, and it zeroes the whole band up front. Both are resident at once, so the peak is roughly element matrices + band + right-hand sides. With many rational surfaces this reaches tens of GB: the DIII-D resistive example at n = 8 peaks at 34 GiB.

With solver = "LU", this PR frees each element matrix as soon as it has been added to the band. It also zeroes the band column by column only as assembly reaches it, so the element matrices drain while the band fills and the peak drops to about max(element matrices, band) + right-hand sides. Δ′ is bit-identical. With solver = "cholesky" the element matrices are kept, because the LU fallback reassembles from them; if the fallback is taken, it frees them in the same way. The first commit changes rdcon/gal.f only.

The memory only goes back to the system if the allocator returns freed blocks. glibc keeps freed blocks of this size once its dynamic mmap threshold has grown past them, so the second commit has rdcon fix that threshold itself on Linux: rdcon/gal_mmap.F calls glibc's mallopt(M_MMAP_THRESHOLD, 131072) at the start of gal_solve. This is the same as setting MALLOC_MMAP_THRESHOLD_=131072, so no environment variable is needed. The build defines HAVE_MALLOPT for rdcon when uname -s is Linux (install/RULES.inc). macOS has no mallopt and already returns these blocks, so there the call is compiled out. The "without the MALLOC setting" column below is the first commit alone on glibc.

Checked with the second commit and no environment variable, DIII-D resistive example, develop → this PR, all outputs (delta_gw, pest3_data) byte-identical:

platform	n	peak memory
Linux, glibc 2.39 (Ubuntu 24.04 container)	1	847 MB → 636 MB (RSS)
Linux, glibc 2.39	4	10.45 GiB → 8.09 GiB (RSS)
macOS (Apple silicon)	1	887 MB → 656 MB (peak footprint; RSS stays high because macOS counts freed reusable pages in it)

GPEC examples, same machine, develop → this PR, both with MALLOC_MMAP_THRESHOLD_=131072:

example	n	Galerkin DOF	peak memory [GiB]	without the MALLOC setting	wall [s]
solovev_resistive	4	82 k	0.96 → 0.64	0.95	31 → 31
DIIID_resistive	2	164 k	1.9 → 1.3	1.9	28 → 27
DIIID_resistive	4	451 k	7.3 → 4.9	7.3	90 → 84
DIIID_resistive	6	832 k	17.2 → 11.7	17.2	248 → 242
DIIID_resistive	8	1.36 M	34.4 → 23.4	34.4	585 → 561

The delta_gw and pest3_data outputs were byte-identical in all 12 example runs. On 32 reconstructed DIII-D equilibria at n = 1–4 (our GPEC branch with and without this commit), Δ′ was bit-identical in all 32, and the largest peak went from 37.9 to 24.8 GiB. The zip holds the script that runs the examples (run_demo.sh DEVELOP_DIR PR_DIR OUT), the plotting script, both tables and the plot.

rdcon_free_mats_demo.zip

🤖 Generated with Claude Code

[StuartBenjamin]

A note on solver choice: with this PR, solver = "LU" and solver = "cholesky" need about the same memory. Cholesky's band is about a third the size of LU's, but Cholesky keeps the element matrices for its LU fallback, and element matrices plus the Cholesky band come to about the size of the LU band. On the 29 DIII-D runs where Cholesky held, its peak was within 1 % of this PR's LU peak (median; at most 9 %), whereas before this PR LU needed about 1.5× Cholesky's memory. The 3 runs that fell back to LU peaked like LU, since the fallback frees the element matrices as it reassembles. Freeing them in the Cholesky path as well would mean rebuilding them whenever Cholesky fails, and a job would still have to budget for that LU fallback, so that is left out of this PR.

[matt-pharr]

This is handy and great for doing scans. I made a tweak with an extra fortran file that fixes the behavior of this addition on mac. Darwin handles memory differently in compiled fortran applications so we have to make it explicit by hacking in this C symbol, but it works and is definitely an improvement. Thanks Stuart