RDCON: free each Galerkin element matrix once it is assembled (LU)

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

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,13 +771,19 @@ 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
                gal%mat(offset,j)=gal%mat(offset,j)+cell%ediag
             ENDIF
          ENDDO
       ENDDO
+      gal%mat(:,jzero+1:)=0
 c-----------------------------------------------------------------------
 c     terminate.
 c-----------------------------------------------------------------------
@@ -1412,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-----------------------------------------------------------------------
@@ -1435,7 +1455,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"

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