From f2c9cefe11fd6a36dd9affa8f1e8b6e6b54886dc Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Thu, 22 Oct 2026 09:43:28 -0400 Subject: [PATCH 1/2] RDCON - ENHANCEMENT - Galerkin LU: free each element matrix once it is assembled The element matrices (cell%mat, 256 mpert^2 bytes per cell) and the global band (gal%mat) were both resident during assembly, and the band was zeroed up front. With LU, each element matrix is now freed right after it is added to the band, and the band is zeroed column by column only as assembly reaches it, so peak memory is about max(element matrices, band) instead of their sum. Cholesky keeps the element matrices, since its LU fallback reassembles from them; the fallback itself then frees them as above. Delta' is bit-identical. The saving needs an allocator that returns freed blocks to the system. glibc keeps freed ~1 MB blocks once its dynamic mmap threshold has risen above them; MALLOC_MMAP_THRESHOLD_=131072 in rdcon's environment keeps them mmapped. Co-Authored-By: Claude Opus 5.5 --- rdcon/gal.f | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/rdcon/gal.f b/rdcon/gal.f index ab1d5f24..8ee1d7e6 100644 --- a/rdcon/gal.f +++ b/rdcon/gal.f @@ -184,7 +184,7 @@ SUBROUTINE gal_alloc(gal,nx,nq,dx0,dx1,dx2,pfac) $ gal%mat(gal%ldab,gal%ndim),gal%ipiv(gal%ndim)) gal%rhs=0 gal%sol=0 - gal%mat=0 +c mat is zeroed in gal_assemble_mat as assembly reaches each column. c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- @@ -211,7 +211,8 @@ SUBROUTINE gal_dealloc(gal) intvl => gal%intvl(ising) DO ix=1,gal%nx cell => intvl%cell(ix) - DEALLOCATE(cell%map,cell%mat) + DEALLOCATE(cell%map) + IF(ALLOCATED(cell%mat))DEALLOCATE(cell%mat) IF(ALLOCATED(cell%emat))DEALLOCATE(cell%emat) IF(ALLOCATED(cell%rhs))DEALLOCATE(cell%rhs) ENDDO @@ -699,7 +700,7 @@ SUBROUTINE gal_assemble_mat(gal) TYPE(gal_type), INTENT(INOUT) :: gal - INTEGER :: ising,ix,ip,ipert,i,jp,jpert,j,offset + INTEGER :: ising,ix,ip,ipert,i,jp,jpert,j,offset,jzero TYPE(interval_type), POINTER :: intvl TYPE(cell_type), POINTER :: cell c----------------------------------------------------------------------- @@ -710,10 +711,18 @@ SUBROUTINE gal_assemble_mat(gal) ELSEIF (solver == "cholesky") THEN offset=1 ENDIF + jzero=0 DO ising=0,msing intvl => gal%intvl(ising) DO ix=1,gal%nx cell => intvl%cell(ix) +c----------------------------------------------------------------------- +c zero mat (allocated untouched) up to this cell's last column, so +c its memory grows as the element matrices are freed. +c----------------------------------------------------------------------- + j=MAX(MAXVAL(cell%map),cell%emap) + gal%mat(:,jzero+1:j)=0 + jzero=MAX(jzero,j) DO ip=0,np DO ipert=1,mpert i=cell%map(ipert,ip) @@ -762,6 +771,11 @@ SUBROUTINE gal_assemble_mat(gal) c----------------------------------------------------------------------- ENDDO ENDDO +c----------------------------------------------------------------------- +c LU never reassembles, so free the element matrix now (cholesky +c keeps it for its LU fallback). +c----------------------------------------------------------------------- + IF (solver == "LU") DEALLOCATE(cell%mat) IF (cell%extra /= "none") THEN IF (cell%etype == "ext1".OR.cell%etype == "ext2") CYCLE j=cell%emap @@ -769,6 +783,7 @@ SUBROUTINE gal_assemble_mat(gal) ENDIF ENDDO ENDDO + gal%mat(:,jzero+1:)=0 c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- @@ -1435,7 +1450,6 @@ SUBROUTINE gal_solve gal%ldab=2*gal%kl+gal%ku+1 DEALLOCATE(gal%mat) ALLOCATE(gal%mat(gal%ldab,gal%ndim)) - gal%mat=0 CALL gal_assemble_mat(gal) ELSE WRITE(*,*)"Calculating Galerkin matrix solution" From ba52126e4252758c7c30feefe6ab471f131609ed Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 16:26:59 +0200 Subject: [PATCH 2/2] RDCON - ENHANCEMENT - Set the glibc mmap threshold so freed Galerkin element matrices return to the system --- install/RULES.inc | 4 ++++ rdcon/gal.f | 5 +++++ rdcon/gal_mmap.F | 43 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+) create mode 100644 rdcon/gal_mmap.F diff --git a/install/RULES.inc b/install/RULES.inc index be30ba9b..5cf18717 100644 --- a/install/RULES.inc +++ b/install/RULES.inc @@ -80,6 +80,10 @@ FCDCON = $(FC) $(FFLAGS) $(OMPFLAG) $(RECURSFLAG) $(LEGACYFLAG) -I$(MATHINC) -I$ # RDCON FCRDCON = $(FC) $(FFLAGS) $(OMPFLAG) $(LEGACYFLAG) -I$(MATHINC) -I$(NETCDFINC) -I../equil -I../vacuum +# glibc's mallopt, used by rdcon/gal_mmap.F; macOS has none +ifeq ($(PLATFORM),Linux) +FCRDCON += -DHAVE_MALLOPT +endif ../rdcon/%.o: ../rdcon/%.f $(FCRDCON) $(MODFLAG) ../rdcon -c $< -o $@ ../rdcon/%.o: ../rdcon/%.F diff --git a/rdcon/gal.f b/rdcon/gal.f index 8ee1d7e6..36cc8736 100644 --- a/rdcon/gal.f +++ b/rdcon/gal.f @@ -1427,6 +1427,11 @@ SUBROUTINE gal_solve WRITE(*,*)"Starting Galerkin method" msol_old=msol msol=2*mpert +c----------------------------------------------------------------------- +c let the element matrices freed in gal_assemble_mat go back to the +c system (see gal_mmap.F). +c----------------------------------------------------------------------- + CALL gal_mmap_threshold CALL gal_alloc(gal,nx,nq,dx0,dx1,dx2,pfac) CALL gal_make_arrays(gal) c----------------------------------------------------------------------- diff --git a/rdcon/gal_mmap.F b/rdcon/gal_mmap.F new file mode 100644 index 00000000..ae037528 --- /dev/null +++ b/rdcon/gal_mmap.F @@ -0,0 +1,43 @@ +c----------------------------------------------------------------------- +c file gal_mmap.F. +c lets glibc return the freed Galerkin element matrices to the +c system. +c----------------------------------------------------------------------- +c gal_assemble_mat frees each element matrix once it is in the +c band. glibc raises its mmap threshold as large blocks are freed, +c so later blocks of that size come from the heap, and their memory +c is kept after they are freed. fixing the threshold at 128 kB, +c below the size of an element matrix, keeps each element matrix in +c its own mapping, which is returned to the system when it is +c freed. this is the same as setting MALLOC_MMAP_THRESHOLD_=131072. +c the build defines HAVE_MALLOPT on Linux. macOS has no mallopt and +c already returns such blocks, so there this does nothing. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE gal_mmap_threshold +#ifdef HAVE_MALLOPT + USE, INTRINSIC :: iso_c_binding, ONLY: c_int + IMPLICIT NONE + + INTERFACE + FUNCTION mallopt(param,value) BIND(C,name="mallopt") + IMPORT :: c_int + INTEGER(c_int), VALUE :: param,value + INTEGER(c_int) :: mallopt + END FUNCTION mallopt + END INTERFACE + + INTEGER(c_int), PARAMETER :: m_mmap_threshold=-3 + INTEGER(c_int) :: info +c----------------------------------------------------------------------- +c set the threshold. +c----------------------------------------------------------------------- + info=mallopt(m_mmap_threshold,131072_c_int) +#endif +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE gal_mmap_threshold