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 ab1d5f24..36cc8736 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----------------------------------------------------------------------- @@ -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" 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