From 83b80e40954d278dc0813aef6b2339b5a9977705 Mon Sep 17 00:00:00 2001 From: Kevin Croker Date: Sun, 3 Sep 2023 23:53:28 -0700 Subject: [PATCH 1/3] ADDED memory tracking and garbage collection to CLASS v3.2.0 and Python classy. Works in parallel with OpenMP (default 100 threads, can be adjusted at compile time) FIXES memory leaks when CLASS terminates on _FAILURE_ when called from Python. --- Makefile | 2 +- default.ini | 2 +- external/HyRec2020/history.c | 33 +-- external/HyRec2020/hydrogen.c | 38 +-- external/HyRec2020/hyrectools.c | 13 +- external/HyRec2020/wrap_hyrec.c | 2 +- external/heating/injection.c | 32 +-- external/heating/noninjection.c | 22 +- include/alloc_track.h | 23 ++ include/common.h | 20 +- main/class.c | 97 +++----- python/cclassy.pxd | 10 +- python/classy.pyx | 21 +- source/background.c | 88 +++---- source/distortions.c | 70 +++--- source/fourier.c | 194 +++++++-------- source/harmonic.c | 38 +-- source/input.c | 40 ++-- source/lensing.c | 104 ++++---- source/output.c | 53 +++-- source/perturbations.c | 145 ++++++------ source/primordial.c | 88 +++---- source/thermodynamics.c | 56 ++--- source/transfer.c | 113 +++++---- tools/alloc_track.c | 406 ++++++++++++++++++++++++++++++++ tools/arrays.c | 90 +++---- tools/dei_rkck.c | 26 +- tools/evolver_ndf15.c | 112 ++++----- tools/evolver_rkck.c | 2 +- tools/growTable.c | 4 +- tools/hyperspherical.c | 20 +- tools/parser.c | 8 +- tools/quadrature.c | 16 +- tools/sparse.c | 26 +- 34 files changed, 1210 insertions(+), 804 deletions(-) create mode 100644 include/alloc_track.h create mode 100644 tools/alloc_track.c diff --git a/Makefile b/Makefile index 9e4a8a779..c0823b287 100644 --- a/Makefile +++ b/Makefile @@ -91,7 +91,7 @@ endif %.o: %.c .base $(HEADERFILES) cd $(WRKDIR);$(CC) $(OPTFLAG) $(OMPFLAG) $(CCFLAG) $(INCLUDES) -c ../$< -o $*.o -TOOLS = growTable.o dei_rkck.o sparse.o evolver_rkck.o evolver_ndf15.o arrays.o parser.o quadrature.o hyperspherical.o common.o trigonometric_integrals.o +TOOLS = growTable.o dei_rkck.o sparse.o evolver_rkck.o evolver_ndf15.o arrays.o parser.o quadrature.o hyperspherical.o common.o trigonometric_integrals.o alloc_track.o SOURCE = input.o background.o thermodynamics.o perturbations.o primordial.o fourier.o transfer.o harmonic.o lensing.o distortions.o diff --git a/default.ini b/default.ini index 4ee76f819..5ac037a3c 100644 --- a/default.ini +++ b/default.ini @@ -181,7 +181,7 @@ sd_detector_name = PIXIE # Name of the detector overwrite_root = no # Overwrite the output files? write_background = no # Write background parameter table write_thermodynamics = no # Write thermodynamics parameter table -#k_output_values = 1e-3,1e-2 # Write perturbations parameter table (at given k) +k_output_values = 1e-3,1e-2 # Write perturbations parameter table (at given k) write_primordial = no # Write primordial parameter table write_parameters = yeap # Write used/unused parameter files write_warnings = yes # Warn about forgotten/wrong inputs diff --git a/external/HyRec2020/history.c b/external/HyRec2020/history.c index 9710fb166..26ad6c76f 100644 --- a/external/HyRec2020/history.c +++ b/external/HyRec2020/history.c @@ -40,6 +40,7 @@ #include "history.h" #include "helium.h" +#include "alloc_track.h" /************************************************************************************* Hubble expansion rate in sec^-1. @@ -74,7 +75,7 @@ void hyrec_init() { rec_data.path_to_hyrec = ""; hyrec_allocate(&rec_data, zmax, zmin); chdir(buffer); - free(buffer); + tracked_free(buffer); } void rec_build_history_camb_(const double* OmegaC, const double* OmegaB, const double* h0inp, const double* tcmb, @@ -160,7 +161,7 @@ double hyrec_xe_(double* a, double *xe){ double loga = log(*a); int error; /* error and error_message are meaningless here */ char *error_message; - error_message = malloc(SIZE_ErrorM); + error_message = tracked_malloc(SIZE_ErrorM); if (loga < logstart) return xe[0]; return rec_interp1d(logstart, dlna, xe, Nz, loga, &error, error_message); } @@ -934,23 +935,23 @@ void hyrec_allocate(HYREC_DATA *data, double zmax, double zmin) { else DLNA = DLNA_HYREC; data->error = 0; - data->error_message=malloc(SIZE_ErrorM); + data->error_message=tracked_malloc(SIZE_ErrorM); sprintf(data->error_message, "\n**** ERROR HAS OCCURRED in HYREC-2 ****\n"); data->zmax = (zmax > 3000.? zmax : 3000.); data->zmin = zmin; - data->atomic = (HYREC_ATOMIC *) malloc(sizeof(HYREC_ATOMIC)); + data->atomic = (HYREC_ATOMIC *) tracked_malloc(sizeof(HYREC_ATOMIC)); allocate_and_read_atomic(data->atomic, &data->error, data->path_to_hyrec, data->error_message); - data->fit = (FIT_FUNC *) malloc(sizeof(FIT_FUNC)); + data->fit = (FIT_FUNC *) tracked_malloc(sizeof(FIT_FUNC)); allocate_and_read_fit(data->fit, &data->error, data->path_to_hyrec, data->error_message); - data->cosmo = (REC_COSMOPARAMS *) malloc(sizeof(REC_COSMOPARAMS)); - data->cosmo->inj_params = (INJ_PARAMS *) malloc(sizeof(INJ_PARAMS)); + data->cosmo = (REC_COSMOPARAMS *) tracked_malloc(sizeof(REC_COSMOPARAMS)); + data->cosmo->inj_params = (INJ_PARAMS *) tracked_malloc(sizeof(INJ_PARAMS)); data->Nz = (long int) (log((1.+zmax)/(1.+zmin))/DLNA) + 2; - data->rad = (RADIATION *) malloc(sizeof(RADIATION)); + data->rad = (RADIATION *) tracked_malloc(sizeof(RADIATION)); // For now assume that radiation field never needed over more than 1 decade in redshift // (typically from z ~ 1700 to 800 for recombination history) @@ -964,16 +965,16 @@ void hyrec_allocate(HYREC_DATA *data, double zmax, double zmin) { void hyrec_free(HYREC_DATA *data) { free_atomic(data->atomic); - free(data->cosmo->inj_params); - free(data->cosmo); - free(data->atomic); - free(data->xe_output); - free(data->Tm_output); - free(data->error_message); + tracked_free(data->cosmo->inj_params); + tracked_free(data->cosmo); + tracked_free(data->atomic); + tracked_free(data->xe_output); + tracked_free(data->Tm_output); + tracked_free(data->error_message); if (MODEL == FULL) free_radiation(data->rad); - free(data->rad); + tracked_free(data->rad); free_fit(data->fit); - free(data->fit); + tracked_free(data->fit); } /****************************************************************** diff --git a/external/HyRec2020/hydrogen.c b/external/HyRec2020/hydrogen.c index 233f174bf..63e279f1e 100644 --- a/external/HyRec2020/hydrogen.c +++ b/external/HyRec2020/hydrogen.c @@ -34,7 +34,7 @@ #include "hyrectools.h" #include "hydrogen.h" - +#include "alloc_track.h" /*********************************************************************************************************** Some constants appropriately rescaled for different values of the fine-structure constant and electron mass @@ -120,7 +120,7 @@ void allocate_and_read_atomic(HYREC_ATOMIC *atomic, int *error, char *path_to_hy char *alpha_file, *rr_file, *twog_file; char sub_message[128]; - alpha_file = malloc(SIZE_InputFile); + alpha_file = tracked_malloc(SIZE_InputFile); alpha_file[0] = 0; strcat(alpha_file, path_to_hyrec); strcat(alpha_file, ALPHA_FILE); @@ -133,7 +133,7 @@ void allocate_and_read_atomic(HYREC_ATOMIC *atomic, int *error, char *path_to_hy return; } - rr_file = malloc(SIZE_InputFile); + rr_file = tracked_malloc(SIZE_InputFile); rr_file[0] = 0; strcat(rr_file, path_to_hyrec); strcat(rr_file, RR_FILE); @@ -187,7 +187,7 @@ void allocate_and_read_atomic(HYREC_ATOMIC *atomic, int *error, char *path_to_hy double L2s1s_current; int fscanf_counter; - twog_file = malloc(SIZE_InputFile); + twog_file = tracked_malloc(SIZE_InputFile); twog_file[0] = 0; strcat(twog_file, path_to_hyrec); strcat(twog_file, TWOG_FILE); @@ -244,9 +244,9 @@ void allocate_and_read_atomic(HYREC_ATOMIC *atomic, int *error, char *path_to_hy for (b = 0; b < NVIRT; b++) atomic->A1s_tab[b] = 0; #endif - free(alpha_file); - free(rr_file); - free(twog_file); + tracked_free(alpha_file); + tracked_free(rr_file); + tracked_free(twog_file); } @@ -273,7 +273,7 @@ void allocate_and_read_fit(FIT_FUNC *fit, int *error, char *path_to_hyrec, char /*********** Effective rates *************/ char *fit_file; - fit_file = malloc(SIZE_InputFile); + fit_file = tracked_malloc(SIZE_InputFile); fit_file[0] = 0; strcat(fit_file, path_to_hyrec); strcat(fit_file, FIT_FILE); @@ -305,7 +305,7 @@ void allocate_and_read_fit(FIT_FUNC *fit, int *error, char *path_to_hyrec, char } } fclose(fA); - free(fit_file); + tracked_free(fit_file); } @@ -315,7 +315,7 @@ Free the memory for rate tables. void free_fit(FIT_FUNC *fit){ unsigned j; - for (j = 0; j < 5; j++) free(fit->swift_func[j]); + for (j = 0; j < 5; j++) tracked_free(fit->swift_func[j]); } /************************************************************************************************ @@ -728,8 +728,8 @@ void populateTS_2photon(double Trr[2][2], double *Trv[2], double *Tvr[2], double } - free(Aup); - free(Adn); + tracked_free(Aup); + tracked_free(Adn); } /********************************************************************* @@ -759,8 +759,8 @@ void solveTXeqB(double *diag, double *updiag, double *dndiag, double *X, double X[N-1] = gamma[N-1]; for (i = N-2; i >= 0; i--) X[i] = gamma[i] - alpha[i] * X[i+1]; - free(alpha); - free(gamma); + tracked_free(alpha); + tracked_free(gamma); } /************************************************************************************************************** @@ -818,8 +818,8 @@ void solve_real_virt(double xr[2], double xv[NVIRT], double Trr[2][2], double *T for (b = 0; b < NVIRT; b++) xv[b] = Tvv_inv_sv[b] - Tvv_inv_Tvr[0][b]*xr[0] - Tvv_inv_Tvr[1][b]*xr[1]; /** Free memory **/ - for (i = 0; i < 2; i++) free(Tvv_inv_Tvr[i]); - free(Tvv_inv_sv); + for (i = 0; i < 2; i++) tracked_free(Tvv_inv_Tvr[i]); + tracked_free(Tvv_inv_sv); } @@ -1016,9 +1016,9 @@ double rec_HMLA_2photon_dxHIIdlna(HYREC_DATA *data, double xe, double xHII, doub /* Average radiation field in each bin */ for (b = 0; b < NVIRT; b++) Dfnu_hist[b][iz] = xv[b]/x1s; - for (i = 0; i < 2; i++) free(Trv[i]); - for (i = 0; i < 2; i++) free(Tvr[i]); - for (i = 0; i < 3; i++) free(Tvv[i]); + for (i = 0; i < 2; i++) tracked_free(Trv[i]); + for (i = 0; i < 2; i++) tracked_free(Tvr[i]); + for (i = 0; i < 3; i++) tracked_free(Tvv[i]); return dxedlna; diff --git a/external/HyRec2020/hyrectools.c b/external/HyRec2020/hyrectools.c index 5e8ac170f..4dfe6d996 100644 --- a/external/HyRec2020/hyrectools.c +++ b/external/HyRec2020/hyrectools.c @@ -10,6 +10,7 @@ January 2015 - added cubic interpolation for non-evenly spaced table) #include #include "hyrectools.h" +#include "alloc_track.h" /****************************************************************************************************** @@ -30,7 +31,7 @@ Creates a [n1] array. double *create_1D_array(unsigned n1, int *error, char error_message[SIZE_ErrorM]){ - double *matrix = (double *) calloc(n1, sizeof(double)); + double *matrix = (double *) tracked_calloc(n1, sizeof(double)); char sub_message[128]; if (*error == 1) return matrix; @@ -49,7 +50,7 @@ Creates a [n1][n2] array. double **create_2D_array(unsigned n1, unsigned n2, int *error, char error_message[SIZE_ErrorM]){ unsigned i; - double **matrix = (double **) calloc(n1, sizeof(double *)); + double **matrix = (double **) tracked_calloc(n1, sizeof(double *)); char sub_message[128]; if (*error == 1) return matrix; @@ -70,9 +71,9 @@ Frees the memory of a [n1][] array. void free_2D_array(double **matrix, unsigned n1){ unsigned i; - for (i = 0; i < n1; i++) free(matrix[i]); + for (i = 0; i < n1; i++) tracked_free(matrix[i]); - free(matrix); + tracked_free(matrix); } /********************************************************************************* @@ -82,7 +83,7 @@ Creates a [n1][n2][n3] matrix. double ***create_3D_array(unsigned n1, unsigned n2, unsigned n3, int *error, char error_message[SIZE_ErrorM]){ unsigned i; - double ***matrix = (double ***) calloc(n1, sizeof(double **)); + double ***matrix = (double ***) tracked_calloc(n1, sizeof(double **)); char sub_message[128]; if (*error == 1) return matrix; @@ -104,7 +105,7 @@ void free_3D_array(double ***matrix, unsigned n1, unsigned n2) { unsigned i; for (i = 0; i < n1; i++) free_2D_array(matrix[i], n2); - free(matrix); + tracked_free(matrix); } /******************************************************************************************** diff --git a/external/HyRec2020/wrap_hyrec.c b/external/HyRec2020/wrap_hyrec.c index fa44c774c..d9d3989a3 100644 --- a/external/HyRec2020/wrap_hyrec.c +++ b/external/HyRec2020/wrap_hyrec.c @@ -78,7 +78,7 @@ int thermodynamics_hyrec_free(struct thermohyrec* phy){ /* We just need to free hyrec (without error management) */ hyrec_free(phy->data); - free(phy->data); + class_free(phy->data); return _SUCCESS_; } diff --git a/external/heating/injection.c b/external/heating/injection.c index 788d9385e..bf6e12fda 100644 --- a/external/heating/injection.c +++ b/external/heating/injection.c @@ -212,42 +212,42 @@ int injection_free(struct thermodynamics* pth){ int index_inj, index_dep; /* Redshift */ - free(pin->z_table); + class_free(pin->z_table); /* Energy injection */ for(index_inj=0;index_injinj_size;++index_inj){ - free(pin->injection_table[index_inj]); + class_free(pin->injection_table[index_inj]); } - free(pin->injection_table); + class_free(pin->injection_table); if(pin->has_PBH_eva==_TRUE_){ - free(pin->PBH_table_z); - free(pin->PBH_table_mass); - free(pin->PBH_table_mass_dd); - free(pin->PBH_table_F); - free(pin->PBH_table_F_dd); + class_free(pin->PBH_table_z); + class_free(pin->PBH_table_mass); + class_free(pin->PBH_table_mass_dd); + class_free(pin->PBH_table_F); + class_free(pin->PBH_table_F_dd); } /* Injection efficiency */ if(pin->f_eff_type == f_eff_from_file){ - free(pin->feff_table); + class_free(pin->feff_table); } if(pin->chi_type == chi_from_z_file){ - free(pin->chiz_table); + class_free(pin->chiz_table); } if(pin->chi_type == chi_from_x_file || pin->chi_type == chi_Galli_file){ - free(pin->chix_table); + class_free(pin->chix_table); } /* Deposition function */ - free(pin->chi); + class_free(pin->chi); /* Energy deposition */ for(index_dep=0;index_depdep_size;++index_dep){ - free(pin->deposition_table[index_dep]); + class_free(pin->deposition_table[index_dep]); } - free(pin->deposition_table); - free(pin->pvecdeposition); + class_free(pin->deposition_table); + class_free(pin->pvecdeposition); return _SUCCESS_; } @@ -874,7 +874,7 @@ int injection_rate_PBH_evaporation_mass_evolution(struct background * pba, } /** - Free local variables */ - free(pvecback_loop); + class_free(pvecback_loop); /** - Spline mass and F(M) evolution in z */ class_call(array_spline_table_lines(pin->PBH_table_z, diff --git a/external/heating/noninjection.c b/external/heating/noninjection.c index 1d1627da4..fbaacd1f6 100644 --- a/external/heating/noninjection.c +++ b/external/heating/noninjection.c @@ -241,8 +241,8 @@ int noninjection_init(struct precision* ppr, } /** - Free temporary variables */ - free(pvecback); - free(pvecthermo); + class_free(pvecback); + class_free(pvecthermo); return _SUCCESS_; } @@ -255,18 +255,18 @@ int noninjection_init(struct precision* ppr, */ int noninjection_free(struct noninjection* pni){ - free(pni->z_table); - free(pni->photon_dep_table); + class_free(pni->z_table); + class_free(pni->photon_dep_table); - free(pni->k); - free(pni->k_weights); - free(pni->pk_primordial_k); + class_free(pni->k); + class_free(pni->k_weights); + class_free(pni->pk_primordial_k); - free(pni->z_table_coarse); - free(pni->noninjection_table); - free(pni->ddnoninjection_table); + class_free(pni->z_table_coarse); + class_free(pni->noninjection_table); + class_free(pni->ddnoninjection_table); - free(pni->integrand_approx); + class_free(pni->integrand_approx); return _SUCCESS_; } diff --git a/include/alloc_track.h b/include/alloc_track.h new file mode 100644 index 000000000..5df9fd337 --- /dev/null +++ b/include/alloc_track.h @@ -0,0 +1,23 @@ +// Header guards +#ifndef ALLOC_TRACK_H +#define ALLOC_TRACK_H + +#include + +#ifdef _OPENMP +#include "omp.h" +#endif + +// Enable if you want to standalone test +// #define DEBUG_ALLOC_TRACK + +void * tracked_malloc(size_t size); +void * tracked_calloc(size_t nmemb, size_t size); +void tracked_free(void *ptr); +void * tracked_realloc(void *ptr, size_t size); +void tracked_free_all(void); +void walk_chunks(void); +int init_memory_tracking(void); +void shutdown_memory_tracking(void); + +#endif diff --git a/include/common.h b/include/common.h index 870081bc8..44d87a6bc 100644 --- a/include/common.h +++ b/include/common.h @@ -8,6 +8,13 @@ #include "svnversion.h" #include +#include "alloc_track.h" + +// Turn this on if you want to make sure none of the code +// (outside of cython) is using untracked malloc/free (or +// related functions +//#pragma GCC poison free malloc calloc realloc + #ifdef _OPENMP #include "omp.h" #endif @@ -147,7 +154,7 @@ int string_begins_with(char* thestring, char beginchar); /* macro for allocating memory and returning error if it failed */ #define class_alloc(pointer, size, error_message_output) { \ - pointer=malloc(size); \ + pointer=tracked_malloc(size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -160,7 +167,7 @@ int string_begins_with(char* thestring, char beginchar); #define class_alloc_parallel(pointer, size, error_message_output) { \ pointer=NULL; \ if (abort == _FALSE_) { \ - pointer=malloc(size); \ + pointer=tracked_malloc(size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -172,7 +179,7 @@ int string_begins_with(char* thestring, char beginchar); /* macro for allocating memory, initializing it with zeros/ and returning error if it failed */ #define class_calloc(pointer, init,size, error_message_output) { \ - pointer=calloc(init,size); \ + pointer=tracked_calloc(init,size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -183,7 +190,7 @@ int string_begins_with(char* thestring, char beginchar); /* macro for re-allocating memory, returning error if it failed */ #define class_realloc(pointer, newname, size, error_message_output) { \ - pointer=realloc(newname,size); \ + pointer=tracked_realloc(newname,size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -192,6 +199,11 @@ int string_begins_with(char* thestring, char beginchar); } \ } +/* macro for freeing memory */ +#define class_free(pointer) { \ + tracked_free(pointer); \ +} + // Testing #define class_test_message(err_out,extra,args...) { \ diff --git a/main/class.c b/main/class.c index a94ae6edd..3c89eb0ac 100644 --- a/main/class.c +++ b/main/class.c @@ -19,108 +19,69 @@ int main(int argc, char **argv) { struct output op; /* for output files */ ErrorMsg errmsg; /* for error messages */ + char failed = 0; + if (input_init(argc, argv,&pr,&ba,&th,&pt,&tr,&pm,&hr,&fo,&le,&sd,&op,errmsg) == _FAILURE_) { printf("\n\nError running input_init \n=>%s\n",errmsg); - return _FAILURE_; + failed = 1; } - if (background_init(&pr,&ba) == _FAILURE_) { + if (!failed && (background_init(&pr,&ba) == _FAILURE_)) { printf("\n\nError running background_init \n=>%s\n",ba.error_message); - return _FAILURE_; + failed = 1; } - if (thermodynamics_init(&pr,&ba,&th) == _FAILURE_) { + if (!failed && (thermodynamics_init(&pr,&ba,&th) == _FAILURE_)) { printf("\n\nError in thermodynamics_init \n=>%s\n",th.error_message); - return _FAILURE_; + failed = 1; } - if (perturbations_init(&pr,&ba,&th,&pt) == _FAILURE_) { + if (!failed && (perturbations_init(&pr,&ba,&th,&pt) == _FAILURE_)) { printf("\n\nError in perturbations_init \n=>%s\n",pt.error_message); - return _FAILURE_; + failed = 1; } - if (primordial_init(&pr,&pt,&pm) == _FAILURE_) { + if (!failed && (primordial_init(&pr,&pt,&pm) == _FAILURE_)) { printf("\n\nError in primordial_init \n=>%s\n",pm.error_message); - return _FAILURE_; + failed = 1; } - if (fourier_init(&pr,&ba,&th,&pt,&pm,&fo) == _FAILURE_) { + if (!failed && (fourier_init(&pr,&ba,&th,&pt,&pm,&fo) == _FAILURE_)) { printf("\n\nError in fourier_init \n=>%s\n",fo.error_message); - return _FAILURE_; + failed = 1; } - if (transfer_init(&pr,&ba,&th,&pt,&fo,&tr) == _FAILURE_) { + if (!failed && (transfer_init(&pr,&ba,&th,&pt,&fo,&tr) == _FAILURE_)) { printf("\n\nError in transfer_init \n=>%s\n",tr.error_message); - return _FAILURE_; + failed = 1; } - if (harmonic_init(&pr,&ba,&pt,&pm,&fo,&tr,&hr) == _FAILURE_) { + if (!failed && (harmonic_init(&pr,&ba,&pt,&pm,&fo,&tr,&hr) == _FAILURE_)) { printf("\n\nError in harmonic_init \n=>%s\n",hr.error_message); - return _FAILURE_; + failed = 1; } - if (lensing_init(&pr,&pt,&hr,&fo,&le) == _FAILURE_) { + if (!failed && (lensing_init(&pr,&pt,&hr,&fo,&le) == _FAILURE_)) { printf("\n\nError in lensing_init \n=>%s\n",le.error_message); - return _FAILURE_; + failed = 1; } - if (distortions_init(&pr,&ba,&th,&pt,&pm,&sd) == _FAILURE_) { + if (!failed && (distortions_init(&pr,&ba,&th,&pt,&pm,&sd) == _FAILURE_)) { printf("\n\nError in distortions_init \n=>%s\n",sd.error_message); - return _FAILURE_; + failed = 1; } - if (output_init(&ba,&th,&pt,&pm,&tr,&hr,&fo,&le,&sd,&op) == _FAILURE_) { + if (!failed && (output_init(&ba,&th,&pt,&pm,&tr,&hr,&fo,&le,&sd,&op) == _FAILURE_)) { printf("\n\nError in output_init \n=>%s\n",op.error_message); - return _FAILURE_; - } - - /****** all calculations done, now free the structures ******/ - - if (distortions_free(&sd) == _FAILURE_) { - printf("\n\nError in distortions_free \n=>%s\n",sd.error_message); - return _FAILURE_; + failed = 1; } - if (lensing_free(&le) == _FAILURE_) { - printf("\n\nError in lensing_free \n=>%s\n",le.error_message); - return _FAILURE_; - } + /****** all calculations done, now free any remaining open memory ******/ - if (harmonic_free(&hr) == _FAILURE_) { - printf("\n\nError in harmonic_free \n=>%s\n",hr.error_message); + tracked_free_all(); + + if(failed) return _FAILURE_; - } - - if (transfer_free(&tr) == _FAILURE_) { - printf("\n\nError in transfer_free \n=>%s\n",tr.error_message); - return _FAILURE_; - } - - if (fourier_free(&fo) == _FAILURE_) { - printf("\n\nError in fourier_free \n=>%s\n",fo.error_message); - return _FAILURE_; - } - - if (primordial_free(&pm) == _FAILURE_) { - printf("\n\nError in primordial_free \n=>%s\n",pm.error_message); - return _FAILURE_; - } - - if (perturbations_free(&pt) == _FAILURE_) { - printf("\n\nError in perturbations_free \n=>%s\n",pt.error_message); - return _FAILURE_; - } - - if (thermodynamics_free(&th) == _FAILURE_) { - printf("\n\nError in thermodynamics_free \n=>%s\n",th.error_message); - return _FAILURE_; - } - - if (background_free(&ba) == _FAILURE_) { - printf("\n\nError in background_free \n=>%s\n",ba.error_message); - return _FAILURE_; - } - - return _SUCCESS_; - + else + return _SUCCESS_; } diff --git a/python/cclassy.pxd b/python/cclassy.pxd index 1d9963f42..af646f365 100644 --- a/python/cclassy.pxd +++ b/python/cclassy.pxd @@ -431,15 +431,7 @@ cdef extern from "class.h": FileArg * value short * read - void lensing_free(void*) - void harmonic_free(void*) - void transfer_free(void*) - void primordial_free(void*) - void perturbations_free(void*) - void thermodynamics_free(void*) - void background_free(void*) - void fourier_free(void*) - void distortions_free(void*) + void tracked_free_all() cdef int _FAILURE_ cdef int _FALSE_ diff --git a/python/classy.pyx b/python/classy.pyx index 4d178b62f..6c2155efe 100644 --- a/python/classy.pyx +++ b/python/classy.pyx @@ -205,24 +205,9 @@ cdef class Class: def struct_cleanup(self): if(self.allocated != True): return - if "distortions" in self.ncp: - distortions_free(&self.sd) - if "lensing" in self.ncp: - lensing_free(&self.le) - if "harmonic" in self.ncp: - harmonic_free(&self.hr) - if "transfer" in self.ncp: - transfer_free(&self.tr) - if "fourier" in self.ncp: - fourier_free(&self.fo) - if "primordial" in self.ncp: - primordial_free(&self.pm) - if "perturb" in self.ncp: - perturbations_free(&self.pt) - if "thermodynamics" in self.ncp: - thermodynamics_free(&self.th) - if "background" in self.ncp: - background_free(&self.ba) + + tracked_free_all() + self.allocated = False self.computed = False diff --git a/source/background.c b/source/background.c index 1ef092478..472789a49 100644 --- a/source/background.c +++ b/source/background.c @@ -885,13 +885,13 @@ int background_free_noinput( struct background *pba ) { - free(pba->tau_table); - free(pba->z_table); - free(pba->loga_table); - free(pba->d2tau_dz2_table); - free(pba->d2z_dtau2_table); - free(pba->background_table); - free(pba->d2background_dloga2_table); + class_free(pba->tau_table); + class_free(pba->z_table); + class_free(pba->loga_table); + class_free(pba->d2tau_dz2_table); + class_free(pba->d2z_dtau2_table); + class_free(pba->background_table); + class_free(pba->d2background_dloga2_table); return _SUCCESS_; } @@ -911,40 +911,40 @@ int background_free_input( if (pba->Omega0_ncdm_tot != 0.) { for (k=0; kN_ncdm; k++) { - free(pba->q_ncdm[k]); - free(pba->w_ncdm[k]); - free(pba->q_ncdm_bg[k]); - free(pba->w_ncdm_bg[k]); - free(pba->dlnf0_dlnq_ncdm[k]); + class_free(pba->q_ncdm[k]); + class_free(pba->w_ncdm[k]); + class_free(pba->q_ncdm_bg[k]); + class_free(pba->w_ncdm_bg[k]); + class_free(pba->dlnf0_dlnq_ncdm[k]); } - free(pba->ncdm_quadrature_strategy); - free(pba->ncdm_input_q_size); - free(pba->ncdm_qmax); - free(pba->q_ncdm); - free(pba->w_ncdm); - free(pba->q_ncdm_bg); - free(pba->w_ncdm_bg); - free(pba->dlnf0_dlnq_ncdm); - free(pba->q_size_ncdm); - free(pba->q_size_ncdm_bg); - free(pba->M_ncdm); - free(pba->T_ncdm); - free(pba->ksi_ncdm); - free(pba->deg_ncdm); - free(pba->Omega0_ncdm); - free(pba->m_ncdm_in_eV); - free(pba->factor_ncdm); + class_free(pba->ncdm_quadrature_strategy); + class_free(pba->ncdm_input_q_size); + class_free(pba->ncdm_qmax); + class_free(pba->q_ncdm); + class_free(pba->w_ncdm); + class_free(pba->q_ncdm_bg); + class_free(pba->w_ncdm_bg); + class_free(pba->dlnf0_dlnq_ncdm); + class_free(pba->q_size_ncdm); + class_free(pba->q_size_ncdm_bg); + class_free(pba->M_ncdm); + class_free(pba->T_ncdm); + class_free(pba->ksi_ncdm); + class_free(pba->deg_ncdm); + class_free(pba->Omega0_ncdm); + class_free(pba->m_ncdm_in_eV); + class_free(pba->factor_ncdm); if (pba->got_files!=NULL) - free(pba->got_files); + class_free(pba->got_files); if (pba->ncdm_psd_files!=NULL) - free(pba->ncdm_psd_files); + class_free(pba->ncdm_psd_files); if (pba->ncdm_psd_parameters!=NULL) - free(pba->ncdm_psd_parameters); + class_free(pba->ncdm_psd_parameters); } if (pba->Omega0_scf != 0.) { if (pba->scf_parameters != NULL) - free(pba->scf_parameters); + class_free(pba->scf_parameters); } return _SUCCESS_; } @@ -1442,8 +1442,8 @@ int background_ncdm_init( pba->error_message), pba->error_message, pba->error_message); - pba->q_ncdm[k]=realloc(pba->q_ncdm[k],pba->q_size_ncdm[k]*sizeof(double)); - pba->w_ncdm[k]=realloc(pba->w_ncdm[k],pba->q_size_ncdm[k]*sizeof(double)); + pba->q_ncdm[k]=tracked_realloc(pba->q_ncdm[k],pba->q_size_ncdm[k]*sizeof(double)); + pba->w_ncdm[k]=tracked_realloc(pba->w_ncdm[k],pba->q_size_ncdm[k]*sizeof(double)); if (pba->background_verbose > 0) { @@ -1470,8 +1470,8 @@ int background_ncdm_init( pba->error_message, pba->error_message); - pba->q_ncdm_bg[k]=realloc(pba->q_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double)); - pba->w_ncdm_bg[k]=realloc(pba->w_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double)); + pba->q_ncdm_bg[k]=tracked_realloc(pba->q_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double)); + pba->w_ncdm_bg[k]=tracked_realloc(pba->w_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double)); /** - in verbose mode, inform user of number of sampled momenta for background quantities */ @@ -1564,9 +1564,9 @@ int background_ncdm_init( /* If allocated, deallocate interpolation table: */ if ((pba->got_files!=NULL)&&(pba->got_files[k]==_TRUE_)) { - free(pbadist.q); - free(pbadist.f0); - free(pbadist.d2f0); + class_free(pbadist.q); + class_free(pbadist.f0); + class_free(pbadist.d2f0); } } @@ -2106,9 +2106,9 @@ int background_solve( } } - free(pvecback); - free(pvecback_integration); - free(used_in_output); + class_free(pvecback); + class_free(pvecback_integration); + class_free(used_in_output); return _SUCCESS_; @@ -2401,7 +2401,7 @@ int background_find_equality( printf(" corresponding to conformal time = %f Mpc\n",pba->tau_eq); } - free(pvecback); + class_free(pvecback); return _SUCCESS_; diff --git a/source/distortions.c b/source/distortions.c index 507629d60..3ddb17e50 100644 --- a/source/distortions.c +++ b/source/distortions.c @@ -98,38 +98,38 @@ int distortions_free(struct distortions * psd) { if (psd->has_distortions == _TRUE_) { /** Delete lists */ - free(psd->z); - free(psd->z_weights); - free(psd->x); - free(psd->x_weights); + class_free(psd->z); + class_free(psd->z_weights); + class_free(psd->x); + class_free(psd->x_weights); /** Delete noise file */ if (psd->has_detector_file == _TRUE_) { - free(psd->delta_Ic_array); + class_free(psd->delta_Ic_array); } /** Delete branching ratios */ for (index_type=0;index_typetype_size;++index_type){ - free(psd->br_table[index_type]); + class_free(psd->br_table[index_type]); } - free(psd->br_table); + class_free(psd->br_table); /** Delete heating functions */ - free(psd->dQrho_dz_tot); + class_free(psd->dQrho_dz_tot); /** Delete distortion shapes */ for (index_type=0;index_typetype_size;++index_type){ - free(psd->sd_shape_table[index_type]); - free(psd->sd_table[index_type]); + class_free(psd->sd_shape_table[index_type]); + class_free(psd->sd_table[index_type]); } - free(psd->sd_shape_table); - free(psd->sd_table); + class_free(psd->sd_shape_table); + class_free(psd->sd_table); /** Delete distortion amplitudes */ - free(psd->sd_parameter_table); + class_free(psd->sd_parameter_table); /** Delete total distortion */ - free(psd->DI); + class_free(psd->DI); } return _SUCCESS_; @@ -759,7 +759,7 @@ int distortions_compute_branching_ratios(struct precision * ppr, class_call(distortions_free_br_data(psd), psd->error_message, psd->error_message); - free(f_E); + class_free(f_E); } @@ -859,7 +859,7 @@ int distortions_compute_heating_rate(struct precision* ppr, psd->dQrho_dz_tot[index_z] = heat*a/(H*rho_g); // [-] } - free(pvecback); + class_free(pvecback); if (psd->include_only_exotic == _FALSE_) { /** Update heating table with second order contributions */ @@ -988,7 +988,7 @@ int distortions_compute_spectral_shapes(struct precision * ppr, class_call(distortions_free_sd_data(psd), psd->error_message, psd->error_message); - free(S); + class_free(S); } /** Compute distortion amplitude for residual parameter epsilon */ @@ -1618,15 +1618,15 @@ int distortions_interpolate_br_data(struct distortions* psd, int distortions_free_br_data(struct distortions * psd){ - free(psd->br_exact_z); - free(psd->f_g_exact); - free(psd->ddf_g_exact); - free(psd->f_y_exact); - free(psd->ddf_y_exact); - free(psd->f_mu_exact); - free(psd->ddf_mu_exact); - free(psd->E_vec); - free(psd->ddE_vec); + class_free(psd->br_exact_z); + class_free(psd->f_g_exact); + class_free(psd->ddf_g_exact); + class_free(psd->f_y_exact); + class_free(psd->ddf_y_exact); + class_free(psd->f_mu_exact); + class_free(psd->ddf_mu_exact); + class_free(psd->E_vec); + class_free(psd->ddE_vec); return _SUCCESS_; } @@ -1858,15 +1858,15 @@ int distortions_interpolate_sd_data(struct distortions* psd, int distortions_free_sd_data(struct distortions * psd){ - free(psd->PCA_nu); - free(psd->PCA_G_T); - free(psd->ddPCA_G_T); - free(psd->PCA_Y_SZ); - free(psd->ddPCA_Y_SZ); - free(psd->PCA_M_mu); - free(psd->ddPCA_M_mu); - free(psd->S_vec); - free(psd->ddS_vec); + class_free(psd->PCA_nu); + class_free(psd->PCA_G_T); + class_free(psd->ddPCA_G_T); + class_free(psd->PCA_Y_SZ); + class_free(psd->ddPCA_Y_SZ); + class_free(psd->PCA_M_mu); + class_free(psd->ddPCA_M_mu); + class_free(psd->S_vec); + class_free(psd->ddS_vec); return _SUCCESS_; } diff --git a/source/fourier.c b/source/fourier.c index cd534ac8b..2c8f235dd 100644 --- a/source/fourier.c +++ b/source/fourier.c @@ -488,7 +488,7 @@ int fourier_pk_at_k_and_z( pfo->error_message, pfo->error_message); - free(ddout_pk_at_z); + class_free(ddout_pk_at_z); // uncomment this part if you prefer a linear interpolation @@ -541,7 +541,7 @@ int fourier_pk_at_k_and_z( pfo->error_message, pfo->error_message); - free(ddout_pk_ic_at_z); + class_free(ddout_pk_ic_at_z); /** --> convert each ic component from logarithmic to linear format */ @@ -640,13 +640,13 @@ int fourier_pk_at_k_and_z( } } - free(pk_primordial_k); - free(pk_primordial_kmin); + class_free(pk_primordial_k); + class_free(pk_primordial_kmin); } - free(out_pk_at_z); + class_free(out_pk_at_z); if (do_ic == _TRUE_) { - free(out_pk_ic_at_z); + class_free(out_pk_ic_at_z); } } @@ -925,14 +925,14 @@ int fourier_pks_at_kvec_and_zvec( index_kvec++; } - free(ln_kvec); + class_free(ln_kvec); if (pfo->has_pk_m == _TRUE_) { - free(ln_pk_table); - free(ddln_pk_table); + class_free(ln_pk_table); + class_free(ddln_pk_table); } if (pfo->has_pk_cb == _TRUE_) { - free(ln_pk_cb_table); - free(ddln_pk_cb_table); + class_free(ln_pk_cb_table); + class_free(ddln_pk_cb_table); } return _SUCCESS_; @@ -1086,8 +1086,8 @@ int fourier_sigmas_at_z( /** - free allocated arrays */ - free(out_pk); - free(ddout_pk); + class_free(out_pk); + class_free(ddout_pk); return _SUCCESS_; } @@ -1564,7 +1564,7 @@ int fourier_init( fprintf(stdout, " -> [WARNING:] Non-linear corrections could not be computed at redshift z=%5.2f and higher.\n This is because k_max is too small for the algorithm (Halofit or HMcode) to be able to compute the scale k_NL at this redshift.\n If non-linear corrections at such high redshift really matter for you,\n just try to increase the precision parameter nonlinear_min_k_max (currently at %e) until k_NL can be computed at the desired z.\n",z,ppr->nonlinear_min_k_max); - free(pvecback); + class_free(pvecback); } } } @@ -1615,14 +1615,14 @@ int fourier_init( /* --> free temporary arrays */ for (index_pk=0; index_pkpk_size; index_pk++){ - free(pk_nl[index_pk]); - free(lnpk_l[index_pk]); - free(ddlnpk_l[index_pk]); + class_free(pk_nl[index_pk]); + class_free(lnpk_l[index_pk]); + class_free(ddlnpk_l[index_pk]); } - free(pk_nl); - free(lnpk_l); - free(ddlnpk_l); + class_free(pk_nl); + class_free(lnpk_l); + class_free(ddlnpk_l); /** --> free the nonlinear workspace */ @@ -1658,52 +1658,52 @@ int fourier_free( if ((pfo->has_pk_matter == _TRUE_) || (pfo->method > nl_none)) { - free(pfo->k); - free(pfo->ln_k); + class_free(pfo->k); + class_free(pfo->ln_k); for (index_pk=0; index_pkpk_size; index_pk++) { - free(pfo->ln_pk_ic_l[index_pk]); - free(pfo->ln_pk_l[index_pk]); + class_free(pfo->ln_pk_ic_l[index_pk]); + class_free(pfo->ln_pk_l[index_pk]); if (pfo->ln_tau_size>1) { - free(pfo->ddln_pk_ic_l[index_pk]); - free(pfo->ddln_pk_l[index_pk]); + class_free(pfo->ddln_pk_ic_l[index_pk]); + class_free(pfo->ddln_pk_l[index_pk]); } } - free(pfo->ln_pk_ic_l); - free(pfo->ln_pk_l); + class_free(pfo->ln_pk_ic_l); + class_free(pfo->ln_pk_l); - free (pfo->sigma8); + class_free(pfo->sigma8); if (pfo->ln_tau_size>1) { - free(pfo->ddln_pk_ic_l); - free(pfo->ddln_pk_l); - free(pfo->ln_tau); + class_free(pfo->ddln_pk_ic_l); + class_free(pfo->ddln_pk_l); + class_free(pfo->ln_tau); } - free(pfo->is_non_zero); + class_free(pfo->is_non_zero); } if (pfo->method > nl_none) { - free(pfo->tau); + class_free(pfo->tau); for (index_pk=0;index_pkpk_size;index_pk++){ - free(pfo->nl_corr_density[index_pk]); - free(pfo->k_nl[index_pk]); - free(pfo->ln_pk_nl[index_pk]); + class_free(pfo->nl_corr_density[index_pk]); + class_free(pfo->k_nl[index_pk]); + class_free(pfo->ln_pk_nl[index_pk]); if (pfo->ln_tau_size > 1) - free(pfo->ddln_pk_nl[index_pk]); + class_free(pfo->ddln_pk_nl[index_pk]); } - free(pfo->nl_corr_density); - free(pfo->k_nl); - free(pfo->ln_pk_nl); + class_free(pfo->nl_corr_density); + class_free(pfo->k_nl); + class_free(pfo->ln_pk_nl); if (pfo->ln_tau_size > 1) - free(pfo->ddln_pk_nl); + class_free(pfo->ddln_pk_nl); } if (pfo->has_pk_eq == _TRUE_) { - free(pfo->pk_eq_tau); - free(pfo->pk_eq_w_and_Omega); - free(pfo->pk_eq_ddw_and_ddOmega); + class_free(pfo->pk_eq_tau); + class_free(pfo->pk_eq_w_and_Omega); + class_free(pfo->pk_eq_ddw_and_ddOmega); } return _SUCCESS_; @@ -2261,8 +2261,8 @@ int fourier_pk_linear( lnpk[index_k] = log(pk); } - free(primordial_pk); - free(pk_ic); + class_free(primordial_pk); + class_free(pk_ic); return _SUCCESS_; @@ -2435,7 +2435,7 @@ int fourier_sigmas( /** - free allocated array */ - free(array_for_sigma); + class_free(array_for_sigma); return _SUCCESS_; } @@ -2522,8 +2522,8 @@ int fourier_sigma_at_z( /** - free allocated arrays */ - free(out_pk); - free(ddout_pk); + class_free(out_pk); + class_free(ddout_pk); return _SUCCESS_; } @@ -2668,7 +2668,7 @@ int fourier_halofit( Omega_m = w_and_Omega[pfo->index_pk_eq_Omega_m]; Omega_v = 1.-Omega_m; - free(w_and_Omega); + class_free(w_and_Omega); } anorm = 1./(2*pow(_PI_,2)); @@ -2783,7 +2783,7 @@ int fourier_halofit( /* class_test_except(sigma < 1., pfo->error_message, - free(pvecback);free(integrand_array), + class_free(pvecback);free(integrand_array), "Your k_max=%g 1/Mpc is too small for Halofit to find the non-linearity scale z_nl at z=%g. Increase input parameter P_k_max_h/Mpc or P_k_max_1/Mpc", pfo->k[pfo->k_size-1], 1./pvecback[pba->index_bg_a]-1.); @@ -2791,8 +2791,8 @@ int fourier_halofit( if (sigma < 1.) { * nl_corr_not_computable_at_this_k = _TRUE_; - free(pvecback); - free(integrand_array); + class_free(pvecback); + class_free(integrand_array); return _SUCCESS_; } else { @@ -2874,7 +2874,7 @@ int fourier_halofit( /* class_test_except(counter > _MAX_IT_, pfo->error_message, - free(pvecback);free(integrand_array), + class_free(pvecback);free(integrand_array), "could not converge within maximum allowed number of iterations"); */ /* ... but in this situation it sounds better to make it stop and return an error! */ @@ -2995,8 +2995,8 @@ int fourier_halofit( } } - free(pvecback); - free(integrand_array); + class_free(pvecback); + class_free(integrand_array); return _SUCCESS_; } @@ -3203,7 +3203,7 @@ int fourier_hmcode( /* The number below is the critical density today, rho_c = 3 * H0^2 / 8*pi*G, in units of M_sun over Mpc^3 */ rho_crit_today_in_msun_mpc3 = 3.*pow(1.e5*pba->h, 2)/8./_PI_/_G_*_Mpc_over_m_/_M_SUN_; - free(pvecback); + class_free(pvecback); /** Test whether pk_cb has to be taken into account (only if we have massive neutrinos)*/ if (pba->has_ncdm==_TRUE_){ @@ -3325,12 +3325,12 @@ int fourier_hmcode( if (nu_min > nu_nl) { if (pfo->fourier_verbose>0) fprintf(stdout, " -> [WARNING:] the minimum mass in the mass-table is too large to find the nonlinear scale at this redshift.\n Decrease mmin_for_p1h_integral\n"); * nl_corr_not_computable_at_this_k = _TRUE_; - free(mass); - free(r_real); - free(r_virial); - free(sigma_r); - free(sigmaf_r); - free(nu_arr); + class_free(mass); + class_free(r_real); + class_free(r_virial); + class_free(sigma_r); + class_free(sigmaf_r); + class_free(nu_arr); return _SUCCESS_; } @@ -3392,12 +3392,12 @@ int fourier_hmcode( if (*k_nl > pfo->k[pfo->k_size-1]) { * nl_corr_not_computable_at_this_k = _TRUE_; - free(mass); - free(r_real); - free(r_virial); - free(sigma_r); - free(sigmaf_r); - free(nu_arr); + class_free(mass); + class_free(r_real); + class_free(r_virial); + class_free(sigma_r); + class_free(sigmaf_r); + class_free(nu_arr); return _SUCCESS_; } else { @@ -3546,7 +3546,7 @@ int fourier_hmcode( if (pk_2h<0.) pk_2h=0.; pk_nl[index_k] = pow((pow(pk_1h, alpha) + pow(pk_2h, alpha)), (1./alpha))/pow(pfo->k[index_k],3)/anorm; //converted back to P_k - free(p1h_integrand); + class_free(p1h_integrand); } // print parameter values @@ -3577,13 +3577,13 @@ int fourier_hmcode( } - free(conc); - free(mass); - free(r_real); - free(r_virial); - free(sigma_r); - free(sigmaf_r); - free(nu_arr); + class_free(conc); + class_free(mass); + class_free(r_real); + class_free(r_virial); + class_free(sigma_r); + class_free(sigmaf_r); + class_free(nu_arr); return _SUCCESS_; } @@ -3656,25 +3656,25 @@ int fourier_hmcode_workspace_free( ) { int index_pk; - free(pnw->rtab); - free(pnw->stab); - free(pnw->ddstab); + class_free(pnw->rtab); + class_free(pnw->stab); + class_free(pnw->ddstab); - free(pnw->growtable); - free(pnw->ztable); - free(pnw->tautable); + class_free(pnw->growtable); + class_free(pnw->ztable); + class_free(pnw->tautable); for (index_pk=0; index_pkpk_size; index_pk++){ - free(pnw->sigma_8[index_pk]); - free(pnw->sigma_disp[index_pk]); - free(pnw->sigma_disp_100[index_pk]); - free(pnw->sigma_prime[index_pk]); + class_free(pnw->sigma_8[index_pk]); + class_free(pnw->sigma_disp[index_pk]); + class_free(pnw->sigma_disp_100[index_pk]); + class_free(pnw->sigma_prime[index_pk]); } - free(pnw->sigma_8); - free(pnw->sigma_disp); - free(pnw->sigma_disp_100); - free(pnw->sigma_prime); + class_free(pnw->sigma_8); + class_free(pnw->sigma_disp); + class_free(pnw->sigma_disp_100); + class_free(pnw->sigma_prime); return _SUCCESS_; } @@ -3731,7 +3731,7 @@ int fourier_hmcode_dark_energy_correction( pfo->error_message, pfo->error_message); - free(pvecback); + class_free(pvecback); pnw->dark_energy_correction = pow(g_wcdm/g_lcdm, 1.5); } @@ -3895,7 +3895,7 @@ int fourier_hmcode_fill_sigtab( } } - free(sigtab); + class_free(sigtab); return _SUCCESS_; } @@ -3955,7 +3955,7 @@ int fourier_hmcode_fill_growtab( } - free(pvecback); + class_free(pvecback); return _SUCCESS_; } @@ -4065,8 +4065,8 @@ int fourier_hmcode_growint( } //fprintf(stdout, "%e %e \n", a, *growth); - free(pvecback); - free(integrand); + class_free(pvecback); + class_free(integrand); return _SUCCESS_; } diff --git a/source/harmonic.c b/source/harmonic.c index a749723ac..a8b457120 100644 --- a/source/harmonic.c +++ b/source/harmonic.c @@ -337,24 +337,24 @@ int harmonic_free( if (phr->ct_size > 0) { for (index_md = 0; index_md < phr->md_size; index_md++) { - free(phr->l_max_ct[index_md]); - free(phr->cl[index_md]); - free(phr->ddcl[index_md]); + class_free(phr->l_max_ct[index_md]); + class_free(phr->cl[index_md]); + class_free(phr->ddcl[index_md]); } - free(phr->l); - free(phr->l_size); - free(phr->l_max_ct); - free(phr->l_max); - free(phr->cl); - free(phr->ddcl); + class_free(phr->l); + class_free(phr->l_size); + class_free(phr->l_max_ct); + class_free(phr->l_max); + class_free(phr->cl); + class_free(phr->ddcl); } for (index_md=0; index_md < phr->md_size; index_md++) - free(phr->is_non_zero[index_md]); + class_free(phr->is_non_zero[index_md]); - free(phr->is_non_zero); - free(phr->ic_size); - free(phr->ic_ic_size); + class_free(phr->is_non_zero); + class_free(phr->ic_size); + class_free(phr->ic_ic_size); } @@ -798,13 +798,13 @@ int harmonic_cls( printf("In %s: time spent in parallel region (loop over l's) = %e s for thread %d\n", __func__,tstop-tstart,omp_get_thread_num()); #endif - free(cl_integrand); + class_free(cl_integrand); - free(primordial_pk); + class_free(primordial_pk); - free(transfer_ic1); + class_free(transfer_ic1); - free(transfer_ic2); + class_free(transfer_ic2); } /* end of parallel region */ @@ -1268,8 +1268,8 @@ int harmonic_compute_cl( } if (ppt->has_cl_number_count == _TRUE_ && _scalars_) { - free(transfer_ic1_nc); - free(transfer_ic2_nc); + class_free(transfer_ic1_nc); + class_free(transfer_ic2_nc); } return _SUCCESS_; diff --git a/source/input.c b/source/input.c index 6f3e3c85b..e0d36ee17 100644 --- a/source/input.c +++ b/source/input.c @@ -712,8 +712,8 @@ int input_shooting(struct file_content * pfc, } /* Free local variables */ - free(x_inout); - free(dxdF); + class_free(x_inout); + class_free(dxdF); } if (input_verbose > 1) { @@ -745,10 +745,10 @@ int input_shooting(struct file_content * pfc, parser_free(&(fzw.fc)); /** Free arrays allocated */ - free(unknown_parameter); - free(fzw.unknown_parameters_index); - free(fzw.target_name); - free(fzw.target_value); + class_free(unknown_parameter); + class_free(fzw.unknown_parameters_index); + class_free(fzw.target_name); + class_free(fzw.target_value); } @@ -841,9 +841,9 @@ int input_shooting(struct file_content * pfc, parser_free(&(fzw.fc)); /** Free arrays allocated */ - free(fzw.unknown_parameters_index); - free(fzw.target_name); - free(fzw.target_value); + class_free(fzw.unknown_parameters_index); + class_free(fzw.target_name); + class_free(fzw.target_value); } return _SUCCESS_; @@ -2973,7 +2973,7 @@ int input_read_parameters_species(struct file_content * pfc, } /* If we don't have perturbations, we should free the arrays again if necessary */ else if (ppt->alpha_idm_dr != NULL) { - free(ppt->alpha_idm_dr); + class_free(ppt->alpha_idm_dr); } } @@ -3018,7 +3018,7 @@ int input_read_parameters_species(struct file_content * pfc, } /* If we don't have perturbations, we should free the arrays again if necessary */ else if (ppt->beta_idr != NULL) { - free(ppt->beta_idr); + class_free(ppt->beta_idr); } } @@ -3922,7 +3922,7 @@ int input_prepare_pk_eq(struct precision * ppr, pvecback), pba->error_message, errmsg); pfo->pk_eq_w_and_Omega[pfo->pk_eq_size*index_pk_eq_z+pfo->index_pk_eq_Omega_m] = pvecback[pba->index_bg_Omega_m]; - free(pvecback); + class_free(pvecback); class_call(background_free_noinput(pba), pba->error_message, @@ -3953,7 +3953,7 @@ int input_prepare_pk_eq(struct precision * ppr, pfo->pk_eq_w_and_Omega[pfo->pk_eq_size*index_pk_eq_z+pfo->index_pk_eq_Omega_m]); } } - free(z); + class_free(z); /** Spline the table for later interpolation */ class_call(array_spline_table_lines(pfo->pk_eq_tau, @@ -4534,7 +4534,7 @@ int input_read_parameters_primordial(struct file_content * pfc, errmsg, "You omitted to write a command for the external Pk"); /* Complete set of parameters */ - ppm->command = (char *) malloc (strlen(string1) + 1); + ppm->command = (char *) tracked_malloc (strlen(string1) + 1); strcpy(ppm->command, string1); /** 1.g.2) Command generating the table */ @@ -4670,7 +4670,7 @@ int input_read_parameters_spectra(struct file_content * pfc, /* Complete set of parameters */ ppt->selection_mean[i] = pointer1[i]; } - free(pointer1); + class_free(pointer1); for (i=1; iselection_mean[i]<=ppt->selection_mean[i-1], @@ -4702,7 +4702,7 @@ int input_read_parameters_spectra(struct file_content * pfc, class_stop(errmsg, "In input for selection function, you asked for %d bin centers and %d bin widths; number of bins unclear; you should pass either one bin width (common to all bins) or %d bin widths.",ppt->selection_num,int1,ppt->selection_num); } - free(pointer1); + class_free(pointer1); } /* Read */ @@ -4726,7 +4726,7 @@ int input_read_parameters_spectra(struct file_content * pfc, "In input for selection function, you asked for %d bin centers and %d bin biases; number of bins unclear; you should pass either one bin bias (common to all bins) or %d bin biases.", ppt->selection_num,int1,ppt->selection_num); } - free(pointer1); + class_free(pointer1); } /* Read */ @@ -4750,7 +4750,7 @@ int input_read_parameters_spectra(struct file_content * pfc, "In input for selection function, you asked for %d bin centers and %d bin biases; number of bins unclear; you should pass either one bin bias (common to all bins) or %d bin biases.", ppt->selection_num,int1,ppt->selection_num); } - free(pointer1); + class_free(pointer1); } } @@ -4857,7 +4857,7 @@ int input_read_parameters_spectra(struct file_content * pfc, for (i=0; iz_pk[i] = pointer1[i]; } - free(pointer1); + class_free(pointer1); } } @@ -5404,7 +5404,7 @@ int input_read_parameters_output(struct file_content * pfc, for (i=0; ik_output_values[i] = pointer1[i]; } - free(pointer1); + class_free(pointer1); qsort (ppt->k_output_values, ppt->k_output_values_num, sizeof(double), compare_doubles); // Sort the k_array using qsort ppt->store_perturbations = _TRUE_; pop->write_perturbations = _TRUE_; diff --git a/source/lensing.c b/source/lensing.c index 29b01d5c4..46fc20588 100644 --- a/source/lensing.c +++ b/source/lensing.c @@ -483,15 +483,15 @@ int lensing_init( for (index_md = 0; index_md < phr->md_size; index_md++) { if (phr->md_size > 1) - free(cl_md[index_md]); + class_free(cl_md[index_md]); if (phr->ic_size[index_md] > 1) - free(cl_md_ic[index_md]); + class_free(cl_md_ic[index_md]); } - free(cl_md_ic); - free(cl_md); + class_free(cl_md_ic); + class_free(cl_md); /** - Compute sigma2\f$(\mu)\f$ and Cgl2(\f$\mu\f$) **/ @@ -745,49 +745,49 @@ int lensing_init( ple->error_message); /** - Free lots of stuff **/ - free(buf_dxx); + class_free(buf_dxx); - free(d00); - free(d11); - free(d1m1); - free(d2m2); + class_free(d00); + class_free(d11); + class_free(d1m1); + class_free(d2m2); if (ple->has_te==_TRUE_) { - free(d20); - free(d3m1); - free(d4m2); + class_free(d20); + class_free(d3m1); + class_free(d4m2); } if (ple->has_ee==_TRUE_ || ple->has_bb==_TRUE_) { - free(d22); - free(d31); - free(d3m3); - free(d40); - free(d4m4); + class_free(d22); + class_free(d31); + class_free(d3m3); + class_free(d40); + class_free(d4m4); } if (ple->has_tt==_TRUE_) - free(ksi); + class_free(ksi); if (ple->has_te==_TRUE_) - free(ksiX); + class_free(ksiX); if (ple->has_ee==_TRUE_ || ple->has_bb==_TRUE_) { - free(ksip); - free(ksim); + class_free(ksip); + class_free(ksim); } - free(Cgl); - free(Cgl2); - free(sigma2); + class_free(Cgl); + class_free(Cgl2); + class_free(sigma2); - free(mu); - free(w8); + class_free(mu); + class_free(w8); - free(cl_unlensed); - free(cl_tt); + class_free(cl_unlensed); + class_free(cl_tt); if (ple->has_te==_TRUE_) - free(cl_te); + class_free(cl_te); if (ple->has_ee==_TRUE_ || ple->has_bb==_TRUE_) { - free(cl_ee); - free(cl_bb); + class_free(cl_ee); + class_free(cl_bb); } - free(cl_pp); + class_free(cl_pp); /** - Exit **/ return _SUCCESS_; @@ -810,10 +810,10 @@ int lensing_free( if (ple->has_lensed_cls == _TRUE_) { - free(ple->l); - free(ple->cl_lens); - free(ple->ddcl_lens); - free(ple->l_max_lt); + class_free(ple->l); + class_free(ple->cl_lens); + class_free(ple->ddcl_lens); + class_free(ple->l_max_lt); } @@ -997,15 +997,15 @@ int lensing_indices( for (index_md = 0; index_md < phr->md_size; index_md++) { if (phr->md_size > 1) - free(cl_md[index_md]); + class_free(cl_md[index_md]); if (phr->ic_size[index_md] > 1) - free(cl_md_ic[index_md]); + class_free(cl_md_ic[index_md]); } - free(cl_md_ic); - free(cl_md); + class_free(cl_md_ic); + class_free(cl_md); /* we want to output Cl_lensed up to the same l_max as Cl_unlensed (even if a number delta_l_max of extra values of l have been used @@ -1282,7 +1282,7 @@ int lensing_d00( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); + class_free(fac1); class_free(fac2); class_free(fac3); return _SUCCESS_; } @@ -1340,7 +1340,7 @@ int lensing_d11( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1397,7 +1397,7 @@ int lensing_d1m1( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1454,7 +1454,7 @@ int lensing_d2m2( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1511,7 +1511,7 @@ int lensing_d22( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1566,7 +1566,7 @@ int lensing_d20( dl = dlp1; } } - free(fac1); free(fac3); free(fac4); + class_free(fac1); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1624,7 +1624,7 @@ int lensing_d31( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1682,7 +1682,7 @@ int lensing_d3m1( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1740,7 +1740,7 @@ int lensing_d3m3( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1797,7 +1797,7 @@ int lensing_d40( dl = dlp1; } } - free(fac1); free(fac3); free(fac4); + class_free(fac1); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1856,7 +1856,7 @@ int lensing_d4m2( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } @@ -1915,6 +1915,6 @@ int lensing_d4m4( dl = dlp1; } } - free(fac1); free(fac2); free(fac3); free(fac4); + class_free(fac1); class_free(fac2); class_free(fac3); class_free(fac4); return _SUCCESS_; } diff --git a/source/output.c b/source/output.c index 92cd751fc..83959b371 100644 --- a/source/output.c +++ b/source/output.c @@ -74,15 +74,15 @@ int output_total_cl_at_l( for (index_md = 0; index_md < phr->md_size; index_md++) { if (phr->md_size > 1) - free(cl_md[index_md]); + class_free(cl_md[index_md]); if (phr->ic_size[index_md] > 1) - free(cl_md_ic[index_md]); + class_free(cl_md_ic[index_md]); } - free(cl_md_ic); - free(cl_md); + class_free(cl_md_ic); + class_free(cl_md); } @@ -586,27 +586,27 @@ int output_cl( fclose(out_md_ic[index_md][index_ic1_ic2]); } } - free(cl_md_ic[index_md]); + class_free(cl_md_ic[index_md]); } } if (ppt->md_size > 1) { for (index_md = 0; index_md < ppt->md_size; index_md++) { fclose(out_md[index_md]); - free(cl_md[index_md]); + class_free(cl_md[index_md]); } } fclose(out); if (ple->has_lensed_cls == _TRUE_) { fclose(out_lensed); } - free(cl_tot); + class_free(cl_tot); for (index_md = 0; index_md < ppt->md_size; index_md++) { - free(out_md_ic[index_md]); + class_free(out_md_ic[index_md]); } - free(out_md_ic); - free(cl_md_ic); - free(out_md); - free(cl_md); + class_free(out_md_ic); + class_free(cl_md_ic); + class_free(out_md); + class_free(cl_md); return _SUCCESS_; @@ -883,10 +883,15 @@ int output_pk( } /* end loop over index_pk */ /* free arrays and pointers */ - free(ln_pk); - if (pk_output == pk_linear) { - free(ln_pk_ic); - free(out_pk_ic); + class_free(ln_pk); + + // KC 8/22/23 + // XXX This was the wrong condition for freeing. We free if do_ic is true, + // because that's when these things get allocated + //if (pk_output == pk_linear) { + if(do_ic == _TRUE_) { + class_free(ln_pk_ic); + class_free(out_pk_ic); } return _SUCCESS_; @@ -1031,7 +1036,7 @@ int output_tk( } - free(data); + class_free(data); return _SUCCESS_; @@ -1081,7 +1086,7 @@ int output_background( data, size_data); - free(data); + class_free(data); fclose(backfile); return _SUCCESS_; @@ -1146,7 +1151,7 @@ int output_thermodynamics( data, size_data); - free(data); + class_free(data); fclose(thermofile); return _SUCCESS_; @@ -1249,7 +1254,7 @@ int output_primordial( data, size_data); - free(data); + class_free(data); fclose(out); return _SUCCESS_; @@ -1312,7 +1317,7 @@ int output_heating(struct injection* pin, struct noninjection* pni, struct outpu titles_injection, data_injection, size_data_injection); - free(data_injection); + class_free(data_injection); fclose(out_injection); } @@ -1353,7 +1358,7 @@ int output_heating(struct injection* pin, struct noninjection* pni, struct outpu titles_noninjection, data_noninjection, size_data_noninjection); - free(data_noninjection); + class_free(data_noninjection); fclose(out_noninjection); } @@ -1415,7 +1420,7 @@ int output_distortions( titles_heat, data_heat, size_data_heat); - free(data_heat); + class_free(data_heat); fclose(out_heat); /* File name */ @@ -1455,7 +1460,7 @@ int output_distortions( titles_distortion, data_distortion, size_data_distortion); - free(data_distortion); + class_free(data_distortion); fclose(out_distortion); } diff --git a/source/perturbations.c b/source/perturbations.c index db8d7fea8..fa9b004cb 100644 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -325,7 +325,7 @@ int perturbations_output_data_at_z( } } } - free(pvecsources); + class_free(pvecsources); } /** - store data */ @@ -337,7 +337,7 @@ int perturbations_output_data_at_z( /** - free tkfull */ // condition necessary because the size could be zero (if ppt->tp_size is zero) if (tkfull != NULL) - free(tkfull); + class_free(tkfull); return _SUCCESS_; } @@ -403,7 +403,7 @@ int perturbations_output_data_at_index_tau( /** - free tkfull */ // condition necessary because the size could be zero (if ppt->tp_size is zero) if (tkfull != NULL) - free(tkfull); + class_free(tkfull); return _SUCCESS_; } @@ -1059,7 +1059,7 @@ int perturbations_init( } /* end loop over modes */ - free(pppw); + class_free(pppw); /** - spline the source array with respect to the time variable */ @@ -1119,9 +1119,9 @@ int perturbations_init( int perturbations_free_input(struct perturbations* ppt) { if (ppt->alpha_idm_dr != NULL) - free(ppt->alpha_idm_dr); + class_free(ppt->alpha_idm_dr); if (ppt->beta_idr != NULL) - free(ppt->beta_idr); + class_free(ppt->beta_idr); return _SUCCESS_; } @@ -1153,54 +1153,54 @@ int perturbations_free( for (index_tp = 0; index_tp < ppt->tp_size[index_md]; index_tp++) { - free(ppt->sources[index_md][index_ic*ppt->tp_size[index_md]+index_tp]); + class_free(ppt->sources[index_md][index_ic*ppt->tp_size[index_md]+index_tp]); if (ppt->ln_tau_size > 1) - free(ppt->ddlate_sources[index_md][index_ic*ppt->tp_size[index_md]+index_tp]); + class_free(ppt->ddlate_sources[index_md][index_ic*ppt->tp_size[index_md]+index_tp]); } } - free(ppt->sources[index_md]); - free(ppt->late_sources[index_md]); - free(ppt->ddlate_sources[index_md]); + class_free(ppt->sources[index_md]); + class_free(ppt->late_sources[index_md]); + class_free(ppt->ddlate_sources[index_md]); - free(ppt->k[index_md]); + class_free(ppt->k[index_md]); } - free(ppt->tau_sampling); + class_free(ppt->tau_sampling); if (ppt->ln_tau_size > 1) - free(ppt->ln_tau); + class_free(ppt->ln_tau); - free(ppt->tp_size); + class_free(ppt->tp_size); - free(ppt->ic_size); + class_free(ppt->ic_size); - free(ppt->k); + class_free(ppt->k); - free(ppt->k_size_cmb); + class_free(ppt->k_size_cmb); - free(ppt->k_size_cl); + class_free(ppt->k_size_cl); - free(ppt->k_size); + class_free(ppt->k_size); - free(ppt->sources); - free(ppt->late_sources); - free(ppt->ddlate_sources); + class_free(ppt->sources); + class_free(ppt->late_sources); + class_free(ppt->ddlate_sources); /** Stuff related to perturbations output: */ /** - Free non-NULL pointers */ if (ppt->k_output_values_num > 0 ) - free(ppt->index_k_output_values); + class_free(ppt->index_k_output_values); for (filenum = 0; filenum<_MAX_NUMBER_OF_K_FILES_; filenum++){ if (ppt->scalar_perturbations_data[filenum] != NULL) - free(ppt->scalar_perturbations_data[filenum]); + class_free(ppt->scalar_perturbations_data[filenum]); if (ppt->vector_perturbations_data[filenum] != NULL) - free(ppt->vector_perturbations_data[filenum]); + class_free(ppt->vector_perturbations_data[filenum]); if (ppt->tensor_perturbations_data[filenum] != NULL) - free(ppt->tensor_perturbations_data[filenum]); + class_free(ppt->tensor_perturbations_data[filenum]); } } @@ -1622,15 +1622,20 @@ int perturbations_indices( } - /* Allocate the titles and data sections for the output file */ - ppt->number_of_scalar_titles=0; - ppt->number_of_vector_titles=0; - ppt->number_of_tensor_titles=0; - for (filenum = 0; filenum<_MAX_NUMBER_OF_K_FILES_; filenum++){ - ppt->scalar_perturbations_data[filenum] = NULL; - ppt->vector_perturbations_data[filenum] = NULL; - ppt->tensor_perturbations_data[filenum] = NULL; - } + // KC 8/29/23 + // This was already done at the top of this function... + // Nothing here gets mutated above... + // So why are we doing this again? + + /* /\* Allocate the titles and data sections for the output file *\/ */ + /* ppt->number_of_scalar_titles=0; */ + /* ppt->number_of_vector_titles=0; */ + /* ppt->number_of_tensor_titles=0; */ + /* for (filenum = 0; filenum<_MAX_NUMBER_OF_K_FILES_; filenum++){ */ + /* ppt->scalar_perturbations_data[filenum] = NULL; */ + /* ppt->vector_perturbations_data[filenum] = NULL; */ + /* ppt->tensor_perturbations_data[filenum] = NULL; */ + /* } */ return _SUCCESS_; @@ -2001,8 +2006,8 @@ int perturbations_timesampling_for_sources( /** - last sampling point = exactly today */ ppt->tau_sampling[counter] = pba->conformal_age; - free(pvecback); - free(pvecthermo); + class_free(pvecback); + class_free(pvecthermo); /** - check the maximum redshift z_max_pk at which the Fourier transfer functions \f$ T_i(k,z)\f$ should be computable by @@ -2680,7 +2685,7 @@ int perturbations_get_k_list( } } - free(ppt->k[index_mode]); + class_free(ppt->k[index_mode]); ppt->k[index_mode] = tmp_k_list; ppt->k_size[index_mode] = newk_size; @@ -2729,8 +2734,8 @@ int perturbations_get_k_list( ppt->k_max = MAX(ppt->k_max,ppt->k[ppt->index_md_tensors][ppt->k_size[ppt->index_md_tensors]-1]); /* last value, inferred from perturbations structure */ } - free(k_max_cmb); - free(k_max_cl); + class_free(k_max_cmb); + class_free(k_max_cl); return _SUCCESS_; @@ -2937,23 +2942,23 @@ int perturbations_workspace_free ( struct perturbations_workspace * ppw ) { - free(ppw->s_l); - free(ppw->pvecback); - free(ppw->pvecthermo); - free(ppw->pvecmetric); + class_free(ppw->s_l); + class_free(ppw->pvecback); + class_free(ppw->pvecthermo); + class_free(ppw->pvecmetric); if (ppw->ap_size > 0) - free(ppw->approx); + class_free(ppw->approx); if (_scalars_) { if ((ppt->has_density_transfers == _TRUE_) || (ppt->has_velocity_transfers == _TRUE_) || (ppt->has_source_delta_m == _TRUE_)) { - free(ppw->delta_ncdm); - free(ppw->theta_ncdm); - free(ppw->shear_ncdm); + class_free(ppw->delta_ncdm); + class_free(ppw->theta_ncdm); + class_free(ppw->shear_ncdm); } } - free(ppw); + class_free(ppw); return _SUCCESS_; } @@ -3194,7 +3199,7 @@ int perturbations_solve( /** - find the number of intervals over which approximation scheme is constant */ class_alloc(interval_number_of,ppw->ap_size*sizeof(int),ppt->error_message); - + ppw->inter_mode = inter_normal; class_call(perturbations_find_approximation_number(ppr, @@ -3235,7 +3240,7 @@ int perturbations_solve( ppt->error_message, ppt->error_message); - free(interval_number_of); + class_free(interval_number_of); /** - fill the structure containing all fixed parameters, indices and workspaces needed by perturbations_derivs */ @@ -3359,11 +3364,11 @@ int perturbations_solve( ppt->error_message); for (index_interval=0; index_intervall_max_ncdm != NULL) free(pv->l_max_ncdm); - if (pv->q_size_ncdm != NULL) free(pv->q_size_ncdm); - free(pv->y); - free(pv->dy); - free(pv->used_in_sources); - free(pv); + if (pv->l_max_ncdm != NULL) class_free(pv->l_max_ncdm); + if (pv->q_size_ncdm != NULL) class_free(pv->q_size_ncdm); + class_free(pv->y); + class_free(pv->dy); + class_free(pv->used_in_sources); + class_free(pv); return _SUCCESS_; } @@ -8499,8 +8504,8 @@ int perturbations_print_variables(double tau, } else{ ppt->scalar_perturbations_data[ppw->index_ikout] = - realloc(ppt->scalar_perturbations_data[ppw->index_ikout], - sizeof(double)*(ppt->size_scalar_perturbation_data[ppw->index_ikout]+ppt->number_of_scalar_titles)); + tracked_realloc(ppt->scalar_perturbations_data[ppw->index_ikout], + sizeof(double)*(ppt->size_scalar_perturbation_data[ppw->index_ikout]+ppt->number_of_scalar_titles)); } storeidx = 0; dataptr = ppt->scalar_perturbations_data[ppw->index_ikout]+ @@ -8609,8 +8614,8 @@ int perturbations_print_variables(double tau, } else{ ppt->tensor_perturbations_data[ppw->index_ikout] = - realloc(ppt->tensor_perturbations_data[ppw->index_ikout], - sizeof(double)*(ppt->size_tensor_perturbation_data[ppw->index_ikout]+ppt->number_of_tensor_titles)); + tracked_realloc(ppt->tensor_perturbations_data[ppw->index_ikout], + sizeof(double)*(ppt->size_tensor_perturbation_data[ppw->index_ikout]+ppt->number_of_tensor_titles)); } storeidx = 0; dataptr = ppt->tensor_perturbations_data[ppw->index_ikout]+ @@ -8684,10 +8689,10 @@ int perturbations_print_variables(double tau, } if (pba->has_ncdm == _TRUE_){ - free(delta_ncdm); - free(theta_ncdm); - free(shear_ncdm); - free(delta_p_over_delta_rho_ncdm); + class_free(delta_ncdm); + class_free(theta_ncdm); + class_free(shear_ncdm); + class_free(delta_p_over_delta_rho_ncdm); } return _SUCCESS_; diff --git a/source/primordial.c b/source/primordial.c index 181ad29ca..10006e371 100644 --- a/source/primordial.c +++ b/source/primordial.c @@ -563,31 +563,31 @@ int primordial_free( if (ppm->primordial_spec_type == analytic_Pk) { for (index_md = 0; index_md < ppm->md_size; index_md++) { - free(ppm->amplitude[index_md]); - free(ppm->tilt[index_md]); - free(ppm->running[index_md]); + class_free(ppm->amplitude[index_md]); + class_free(ppm->tilt[index_md]); + class_free(ppm->running[index_md]); } - free(ppm->amplitude); - free(ppm->tilt); - free(ppm->running); + class_free(ppm->amplitude); + class_free(ppm->tilt); + class_free(ppm->running); } else if (ppm->primordial_spec_type == external_Pk) { - free(ppm->command); + class_free(ppm->command); } for (index_md = 0; index_md < ppm->md_size; index_md++) { - free(ppm->lnpk[index_md]); - free(ppm->ddlnpk[index_md]); - free(ppm->is_non_zero[index_md]); + class_free(ppm->lnpk[index_md]); + class_free(ppm->ddlnpk[index_md]); + class_free(ppm->is_non_zero[index_md]); } - free(ppm->lnpk); - free(ppm->ddlnpk); - free(ppm->is_non_zero); - free(ppm->ic_size); - free(ppm->ic_ic_size); + class_free(ppm->lnpk); + class_free(ppm->ddlnpk); + class_free(ppm->is_non_zero); + class_free(ppm->ic_size); + class_free(ppm->ic_ic_size); - free(ppm->lnk); + class_free(ppm->lnk); } @@ -1211,7 +1211,7 @@ int primordial_inflation_solve_inflation( &dphidt_pivot), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); break; case inflation_H: @@ -1227,11 +1227,11 @@ int primordial_inflation_solve_inflation( &dddH), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); break; default: - free(y);free(y_ini);free(dy); + class_free(y);class_free(y_ini);class_free(dy); class_stop(ppm->error_message,"ppm->primordial_spec_type=%d different from possible relevant cases",ppm->primordial_spec_type); break; } @@ -1268,7 +1268,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); /* we need to do the opposite: to check that there is an initial time such that k_min << (aH)_ini. A guess is made by integrating @@ -1310,7 +1310,7 @@ int primordial_inflation_solve_inflation( class_test_except(counter >= ppr->primordial_inflation_phi_ini_maxit, ppm->error_message, - free(y);free(y_ini);free(dy), + class_free(y);class_free(y_ini);class_free(dy), "when searching for an initial value of phi just before observable inflation takes place, could not converge after %d iterations. The potential does not allow eough inflationary e-folds before reaching the pivot scale", counter); @@ -1333,7 +1333,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); phi_try = y[ppm->index_in_phi]; @@ -1352,7 +1352,7 @@ int primordial_inflation_solve_inflation( &dphidt_try), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); /* we need to normalize a properly so that a=a_pivot when phi=phi_pivot. To do so, we evolve starting arbitrarily from @@ -1373,7 +1373,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); /* now impose the correct a_ini */ a_try = a_pivot/y[ppm->index_in_a]; @@ -1406,7 +1406,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); y_ini[ppm->index_in_a] = y[ppm->index_in_a]; y_ini[ppm->index_in_phi] = y[ppm->index_in_phi]; @@ -1414,7 +1414,7 @@ int primordial_inflation_solve_inflation( break; default: - free(y);free(y_ini);free(dy); + class_free(y);class_free(y_ini);class_free(dy); class_stop(ppm->error_message,"ppm->primordial_spec_type=%d different from possible relevant cases",ppm->primordial_spec_type); break; } @@ -1433,7 +1433,7 @@ int primordial_inflation_solve_inflation( y_ini), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); } else if (ppm->behavior == analytical) { @@ -1443,7 +1443,7 @@ int primordial_inflation_solve_inflation( y_ini), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); } else { class_stop(ppm->error_message,"Uncomprehensible value of the flag ppm->behavior=%d",ppm->behavior); @@ -1468,7 +1468,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); ppm->phi_min=y[ppm->index_in_phi]; @@ -1483,7 +1483,7 @@ int primordial_inflation_solve_inflation( conformal), ppm->error_message, ppm->error_message, - free(y);free(y_ini);free(dy)); + class_free(y);class_free(y_ini);class_free(dy)); ppm->phi_max=y[ppm->index_in_phi]; @@ -1494,9 +1494,9 @@ int primordial_inflation_solve_inflation( /** - finally, we can de-allocate */ - free(y); - free(y_ini); - free(dy); + class_free(y); + class_free(y_ini); + class_free(dy); return _SUCCESS_; } @@ -1738,8 +1738,8 @@ int primordial_inflation_one_wavenumber( ppm->error_message, ppm->error_message); - free(y); - free(dy); + class_free(y); + class_free(dy); class_test(curvature<=0., ppm->error_message, @@ -3287,10 +3287,10 @@ int primordial_external_spectrum_init( /** - Initialization */ /* Prepare the data (with some initial size) */ n_data_guess = 100; - k = (double *)malloc(n_data_guess*sizeof(double)); - pks = (double *)malloc(n_data_guess*sizeof(double)); + k = (double *)tracked_malloc(n_data_guess*sizeof(double)); + pks = (double *)tracked_malloc(n_data_guess*sizeof(double)); if (ppt->has_tensors == _TRUE_) - pkt = (double *)malloc(n_data_guess*sizeof(double)); + pkt = (double *)tracked_malloc(n_data_guess*sizeof(double)); /* Prepare the command */ /* If the command is just a "cat", no arguments need to be passed */ if (strncmp("cat ", ppm->command, 4) == 0) { @@ -3325,18 +3325,18 @@ int primordial_external_spectrum_init( /* (it is faster and safer that reallocating every new line) */ if ((n_data+1) > n_data_guess) { n_data_guess *= 2; - tmp = (double *)realloc(k, n_data_guess*sizeof(double)); + tmp = (double *)tracked_realloc(k, n_data_guess*sizeof(double)); class_test(tmp == NULL, ppm->error_message, "Error allocating memory to read the external spectrum.\n"); k = tmp; - tmp = (double *)realloc(pks, n_data_guess*sizeof(double)); + tmp = (double *)tracked_realloc(pks, n_data_guess*sizeof(double)); class_test(tmp == NULL, ppm->error_message, "Error allocating memory to read the external spectrum.\n"); pks = tmp; if (ppt->has_tensors == _TRUE_) { - tmp = (double *)realloc(pkt, n_data_guess*sizeof(double)); + tmp = (double *)tracked_realloc(pkt, n_data_guess*sizeof(double)); class_test(tmp == NULL, ppm->error_message, "Error allocating memory to read the external spectrum.\n"); @@ -3416,10 +3416,10 @@ int primordial_external_spectrum_init( */ }; /** - Release the memory used locally */ - free(k); - free(pks); + class_free(k); + class_free(pks); if (ppt->has_tensors == _TRUE_) - free(pkt); + class_free(pkt); /** - Tell CLASS that there are scalar (and tensor) modes */ ppm->is_non_zero[ppt->index_md_scalars][ppt->index_ic_ad] = _TRUE_; if (ppt->has_tensors == _TRUE_) diff --git a/source/thermodynamics.c b/source/thermodynamics.c index 42f29c257..ec10f3a85 100644 --- a/source/thermodynamics.c +++ b/source/thermodynamics.c @@ -424,7 +424,7 @@ int thermodynamics_init( pth->error_message, pth->error_message); - free(pvecback); + class_free(pvecback); return _SUCCESS_; } @@ -448,10 +448,10 @@ int thermodynamics_free( } /* Free thermodynamics-related functions */ - free(pth->z_table); - free(pth->tau_table); - free(pth->thermodynamics_table); - free(pth->d2thermodynamics_dz2_table); + class_free(pth->z_table); + class_free(pth->tau_table); + class_free(pth->thermodynamics_table); + class_free(pth->d2thermodynamics_dz2_table); return _SUCCESS_; } @@ -520,7 +520,7 @@ int thermodynamics_helium_from_bbn( -pvecback[pba->index_bg_rho_g]) /(7./8.*pow(4./11.,4./3.)*pvecback[pba->index_bg_rho_g]); - free(pvecback); + class_free(pvecback); // printf("Neff early = %g, Neff at bbn: %g\n",pba->Neff,Neff_bbn); @@ -669,12 +669,12 @@ int thermodynamics_helium_from_bbn( } /** - deallocate arrays */ - free(omegab); - free(deltaN); - free(YHe); - free(ddYHe); - free(YHe_at_deltaN); - free(ddYHe_at_deltaN); + class_free(omegab); + class_free(deltaN); + class_free(YHe); + class_free(ddYHe); + class_free(YHe_at_deltaN); + class_free(ddYHe_at_deltaN); return _SUCCESS_; } @@ -1654,8 +1654,8 @@ int thermodynamics_solve( pth->error_message); } - free(interval_limit); - free(mz_output); + class_free(interval_limit); + class_free(mz_output); return _SUCCESS_; @@ -1825,8 +1825,8 @@ int thermodynamics_workspace_free( struct thermo_workspace * ptw ) { - free(ptw->ptdw->ap_z_limits); - free(ptw->ptdw->ap_z_limits_delta); + class_free(ptw->ptdw->ap_z_limits); + class_free(ptw->ptdw->ap_z_limits_delta); switch (pth->recombination) { @@ -1834,19 +1834,19 @@ int thermodynamics_workspace_free( class_call(thermodynamics_hyrec_free(ptw->ptdw->phyrec), ptw->ptdw->phyrec->error_message, pth->error_message); - free(ptw->ptdw->phyrec); + class_free(ptw->ptdw->phyrec); break; case recfast: - free(ptw->ptdw->precfast); + class_free(ptw->ptdw->precfast); break; } - free(ptw->ptrp->reionization_parameters); - free(ptw->ptdw); - free(ptw->ptrp); + class_free(ptw->ptrp->reionization_parameters); + class_free(ptw->ptdw); + class_free(ptw->ptrp); - free(ptw); + class_free(ptw); return _SUCCESS_; } @@ -3088,10 +3088,10 @@ int thermodynamics_vector_free( struct thermo_vector * tv ) { - free(tv->y); - free(tv->dy); - free(tv->used_in_output); - free(tv); + class_free(tv->y); + class_free(tv->dy); + class_free(tv->used_in_output); + class_free(tv); return _SUCCESS_; } @@ -3266,7 +3266,7 @@ int thermodynamics_calculate_damping_scale( pth->error_message, pth->error_message); - free(tau_table_growing); + class_free(tau_table_growing); /* we could now write the result as r_d = 2pi * sqrt(integral), * but we will first better acount for the contribution frokm the tau_ini boundary. @@ -4714,7 +4714,7 @@ int thermodynamics_idm_initial_temperature( /* This formula (assuming alpha,beta,epsilon=const) approximates the steady-state solution of the IDM temperature evolution equation */ ptdw->T_idm = (alpha + beta + epsilon * pba->T_idr/pba->T_cmb)/(1.+epsilon+alpha+beta) * pba->T_cmb * (1.+z_ini); - free(pvecback); + class_free(pvecback); return _SUCCESS_; } diff --git a/source/transfer.c b/source/transfer.c index 25e68daef..d8ce6830c 100644 --- a/source/transfer.c +++ b/source/transfer.c @@ -295,6 +295,9 @@ int transfer_init( /** - precompute window function for integrated nCl/sCl quantities*/ double* window = NULL; if ((ppt->has_cl_lensing_potential == _TRUE_) || (ppt->has_cl_number_count == _TRUE_)) { + + // KC 8/22/23 + // This is sending a STACK address... class_call(transfer_precompute_selection(ppr, pba, ppt, @@ -402,7 +405,12 @@ int transfer_init( if (abort == _TRUE_) return _FAILURE_; /** - finally, free arrays allocated outside parallel zone */ - free(window); + + // KC 8/22/23 + // XXX? window is never allocated unless these conditions are true. So this was + // freeing a void pointer before. Surprisingly, free(0x0) does not die on arrival... + if ((ppt->has_cl_lensing_potential == _TRUE_) || (ppt->has_cl_number_count == _TRUE_)) + class_free(window); class_call(transfer_perturbation_sources_spline_free(ppt,ptr,sources_spline), ptr->error_message, @@ -441,30 +449,30 @@ int transfer_free( if (ptr->has_cls == _TRUE_) { for (index_md = 0; index_md < ptr->md_size; index_md++) { - free(ptr->l_size_tt[index_md]); - free(ptr->transfer[index_md]); - free(ptr->k[index_md]); + class_free(ptr->l_size_tt[index_md]); + class_free(ptr->transfer[index_md]); + class_free(ptr->k[index_md]); } - free(ptr->tt_size); - free(ptr->l_size_tt); - free(ptr->l_size); - free(ptr->l); - free(ptr->q); - free(ptr->k); - free(ptr->transfer); + class_free(ptr->tt_size); + class_free(ptr->l_size_tt); + class_free(ptr->l_size); + class_free(ptr->l); + class_free(ptr->q); + class_free(ptr->k); + class_free(ptr->transfer); if (ptr->nz_size > 0) { - free(ptr->nz_z); - free(ptr->nz_nz); - free(ptr->nz_ddnz); + class_free(ptr->nz_z); + class_free(ptr->nz_nz); + class_free(ptr->nz_ddnz); } if (ptr->nz_evo_size > 0) { - free(ptr->nz_evo_z); - free(ptr->nz_evo_nz); - free(ptr->nz_evo_dlog_nz); - free(ptr->nz_evo_dd_dlog_nz); + class_free(ptr->nz_evo_z); + class_free(ptr->nz_evo_nz); + class_free(ptr->nz_evo_dlog_nz); + class_free(ptr->nz_evo_dd_dlog_nz); } } @@ -774,13 +782,13 @@ int transfer_perturbation_sources_free( ((ppt->has_source_phi_plus_psi == _TRUE_) && (index_tp == ppt->index_tp_phi_plus_psi)) || ((ppt->has_source_psi == _TRUE_) && (index_tp == ppt->index_tp_psi)))) { - free(sources[index_md][index_ic * ppt->tp_size[index_md] + index_tp]); + class_free(sources[index_md][index_ic * ppt->tp_size[index_md] + index_tp]); } } } - free(sources[index_md]); + class_free(sources[index_md]); } - free(sources); + class_free(sources); return _SUCCESS_; } @@ -797,12 +805,12 @@ int transfer_perturbation_sources_spline_free( for (index_md = 0; index_md < ptr->md_size; index_md++) { for (index_ic = 0; index_ic < ppt->ic_size[index_md]; index_ic++) { for (index_tp = 0; index_tp < ppt->tp_size[index_md]; index_tp++) { - free(sources_spline[index_md][index_ic * ppt->tp_size[index_md] + index_tp]); + class_free(sources_spline[index_md][index_ic * ppt->tp_size[index_md] + index_tp]); } } - free(sources_spline[index_md]); + class_free(sources_spline[index_md]); } - free(sources_spline); + class_free(sources_spline); return _SUCCESS_; } @@ -1439,9 +1447,9 @@ int transfer_free_source_correspondence( int index_md; for (index_md = 0; index_md < ptr->md_size; index_md++) { - free(tp_of_tt[index_md]); + class_free(tp_of_tt[index_md]); } - free(tp_of_tt); + class_free(tp_of_tt); return _SUCCESS_; @@ -2774,7 +2782,7 @@ int transfer_source_resample( } /* deallocate the temporary array */ - free(source_at_tau); + class_free(source_at_tau); return _SUCCESS_; @@ -3340,7 +3348,7 @@ int transfer_integrate( } - free(radial_function); + class_free(radial_function); return _SUCCESS_; } @@ -4043,11 +4051,11 @@ int transfer_radial_function( break; } - free(Phi); - free(dPhi); - free(d2Phi); - free(chireverse); - free(rescale_function); + class_free(Phi); + class_free(dPhi); + class_free(d2Phi); + class_free(chireverse); + class_free(rescale_function); return _SUCCESS_; } @@ -4306,15 +4314,15 @@ int transfer_workspace_free( ptr->error_message, ptr->error_message); } - free(ptw->interpolated_sources); - free(ptw->sources); - free(ptw->tau0_minus_tau); - free(ptw->w_trapz); - free(ptw->chi); - free(ptw->cscKgen); - free(ptw->cotKgen); - - free(ptw); + class_free(ptw->interpolated_sources); + class_free(ptw->sources); + class_free(ptw->tau0_minus_tau); + class_free(ptw->w_trapz); + class_free(ptw->chi); + class_free(ptw->cscKgen); + class_free(ptw->cotKgen); + + class_free(ptw); return _SUCCESS_; } @@ -4665,7 +4673,14 @@ int transfer_precompute_selection( class_alloc(tau0_minus_tau,tau_size_max*sizeof(double),ptr->error_message); class_alloc(selection,tau_size_max*sizeof(double),ptr->error_message); class_alloc(w_trapz,tau_size_max*sizeof(double),ptr->error_message); - class_alloc((*window),tau_size_max*ptr->tt_size[index_md]*sizeof(double),ptr->error_message); + //class_alloc((*window),tau_size_max*ptr->tt_size[index_md]*sizeof(double),ptr->error_message); + + // KC 8/22/23 + // See if the macro substitution was going wild... + (*window) = tracked_malloc(tau_size_max*ptr->tt_size[index_md]*sizeof(double)); + + // Did this fail? + printf("I asked for a window, what did we get? Address: %x\n", (*window)); class_alloc(pvecback,pba->bg_size*sizeof(double),ptr->error_message); /* conformal time today */ @@ -5046,17 +5061,17 @@ int transfer_precompute_selection( } /* deallocate temporary arrays */ - free(tau0_minus_tau_lensing_sources); - free(w_trapz_lensing_sources); + class_free(tau0_minus_tau_lensing_sources); + class_free(w_trapz_lensing_sources); } /* End integrated contribution */ } /* deallocate temporary arrays */ - free(selection); - free(tau0_minus_tau); - free(w_trapz); - free(pvecback); + class_free(selection); + class_free(tau0_minus_tau); + class_free(w_trapz); + class_free(pvecback); return _SUCCESS_; } diff --git a/tools/alloc_track.c b/tools/alloc_track.c new file mode 100644 index 000000000..854203465 --- /dev/null +++ b/tools/alloc_track.c @@ -0,0 +1,406 @@ +#include +#include + +#include "alloc_track.h" + +struct chunk { + + void *block; + + // Double direction is easier to implememt. + struct chunk *prev; + struct chunk *next; +}; + +// The head always contains the most-recently +// allocated chunk. +// +// Make this simple: just space for 100 threads tops. +// Increase if you need it bigger. +// +#define MAX_TRACKING 100 +struct chunk *tracking_head[MAX_TRACKING]; +int allocations[MAX_TRACKING]; + +// We do this function pointer song and dance +// so that we only have one allocation implementation +// to deal with. +void * malloc_wrapper(size_t size, size_t unused) { + + return malloc(size); +} + +void * calloc_wrapper(size_t nmemb, size_t size) { + + return calloc(nmemb, size); +} + +// This method uses the requested allocator and tracks the +// returned pointer at the head of the list. +// +// We use the head of the list because the most frequent usage is to +// free memory in the reversed order from which you allocated it. +// So most of the time, we are landing right on the block we want +// to free. +void * tracked_alloc_core( void * (*allocation_wrapper)(size_t, size_t), + size_t arg1, + size_t arg2) { + + struct chunk *achunk; + void *ptr; + int me = 0; + + // Get a chunk to keep track of + achunk = (struct chunk *)malloc(sizeof(struct chunk)); + if(!achunk) + return 0x0; + + // We use the desired call here + ptr = allocation_wrapper(arg1, arg2); + + if(!ptr) { + + // It failed, so deallocate achunk and return 0x0 + free(achunk); + return 0x0; + } + + // Track it + achunk->block = ptr; + +#ifdef _OPENMP + me = omp_get_thread_num(); +#endif + + // Does the table exist yet? + if(!tracking_head[me]) { + + // Nope. This becomes the head + tracking_head[me] = achunk; + achunk->prev = NULL; + achunk->next = NULL; + } + else { + + // Table exists already. + // Attach new chunk ****AT THE HEAD**** + achunk->next = tracking_head[me]; + achunk->prev = NULL; + tracking_head[me]->prev = achunk; + tracking_head[me] = achunk; + } + ++allocations[me]; + + // Return the newly allocated block + return ptr; +} + +void * tracked_malloc(size_t size) { + + return tracked_alloc_core(malloc_wrapper, size, 0); +} + +void * tracked_calloc(size_t nmemb, size_t size) { + + return tracked_alloc_core(calloc_wrapper, nmemb, size); +} + +void tracked_free(void *ptr) { + + struct chunk *list_cur; + int pos = 0; + int me = 0; + + // Mimic POSIX free() behaviour on null pointeres + if(!ptr) + return; + +#ifdef _OPENMP + me = omp_get_thread_num(); +#endif + + list_cur = tracking_head[me]; + + // Check for empty list! + if(!list_cur) { + + // Force a segfault + fprintf(stderr, "tracked_free(): attempting to free unknown %x from an empty list\n", ptr); + walk_chunks(); + printf("%d", *(int *)(0x0)); + } + + // Search the list until we find it + while(list_cur->block != ptr) { + + // Move to the right if we can + if(list_cur->next != NULL) { + + ++pos; + + // Did we just run off the edge? + if(pos == allocations[me]) { + fprintf(stderr, "tracked_free(): table corruption! we just ran off the edge of the list (position %d)\n", pos); + walk_chunks(); + printf("%d", *(int *)(0x0)); + } + + list_cur = list_cur->next; + } + else { + // Force a segfault + fprintf(stderr, "tracked_free(): requested to free 0x%x, but this address is not tracked!\n", ptr); + walk_chunks(); + printf("%d", *(int *)(0x0)); + } + } + + // We have found the right chunk. + // Unlink it + // Three cases: + + if(list_cur == tracking_head[me]) { + + // We are at the head of the list. + + // SANITY CHECK + if(list_cur->prev != NULL) { + + fprintf(stderr, "tracked_free(): the head of the list has a non-null record to the left!\n"); + walk_chunks(); + printf("%d", *(int *)(0x0)); + } + + // Are we the only one? + if(list_cur->next == NULL) { + + // Yes. We only had one thing. + + // Free the memory + free(list_cur->block); + + // Free the accounting chunk associated to it + free(list_cur); + + // Set the tracking head to empty + tracking_head[me] = NULL; + + // Record that we freed memory + --allocations[me]; + } + else { + + // There is a subsequent one + + // Free the memory + free(list_cur->block); + + // Unlink it + list_cur->next->prev = NULL; + + // Set the head to its next + tracking_head[me] = list_cur->next; + + // Free the memory associated with the accounting structure + free(list_cur); + + // Record that we freed memory + --allocations[me]; + } + } + else { + + // We are either in the middle of the list, + // or at the end (and we are never in first position). + + // Are we at the end? + if(list_cur->next == NULL) { + + // Free its memory + free(list_cur->block); + + // We can just unlink it + list_cur->prev->next = NULL; + + // Now destroy the accounting structure + free(list_cur); + + // Record that we freed memory + --allocations[me]; + } + else { + + // We are in the middle. + + // Unlink me + list_cur->prev->next = list_cur->next; + list_cur->next->prev = list_cur->prev; + + // Free the allocated memory + free(list_cur->block); + + // Now free the accounting structure + free(list_cur); + + // Record that we freed memory + --allocations[me]; + } + } + return; +} + +void * tracked_realloc(void *ptr, size_t size) { + + void *newblock; + int me = 0; + struct chunk *list_cur; + +#ifdef _OPENMP + me = omp_get_thread_num(); +#endif + + list_cur = tracking_head[me]; + + // Accounting structures don't need to change, just + // what is being tracked is changing. + + // Search the list until we find it + while(list_cur->block != ptr) { + + // Move to the right if we can + if(list_cur->next != NULL) + list_cur = list_cur->next; + else { + fprintf(stderr, "tracked_realloc(): requested to reallocate 0x%x, but this address is not tracked!\n", ptr); + + // Walk the chunks + walk_chunks(); + + // Force a segfault + printf("%d", *(int *)(0x0)); + } + } + + // Call the underlying realloc() + newblock = realloc(ptr, size); + + // WTF was I doing here before? LOL + if(newblock != ptr) + list_cur->block = newblock; + + return newblock; +} + +void tracked_free_all(void) { + + struct chunk *list_next; + int i; + int counter; + + // Clear out all allocations in parallel + // Should probably just walk the arrays here + for(i = 0; i < MAX_TRACKING; ++i) { + + if(!tracking_head[i]) + continue; + + counter = 0; + + while(tracking_head[i]) { + + ++counter; + list_next = tracking_head[i]->next; + free(tracking_head[i]->block); + free(tracking_head[i]); + tracking_head[i] = list_next; + --allocations[i]; + } + //printf("[Thread %d] alloc_track: Cleaned up %d allocations\n", i, counter); + } +} + +void walk_chunks(void) { + + int i = 0; + struct chunk *achunk; + int me = 0; + + printf("Total allocations open: %d\n", allocations); + +#ifdef _OPENMP + me = omp_get_thread_num(); +#endif + + achunk = tracking_head[me]; + + // Walk forward + while(achunk) { + printf("(forward) Position %d: memory at %x, contains %d\n", i++, achunk->block, *(int *)(achunk->block)); + + if(achunk->next != NULL) + achunk = achunk->next; + else + break; + } + + // Walk backward + while(achunk) { + + printf("(backward) Position %d: memory at %x, contains %d\n", --i, achunk->block, *(int *)(achunk->block)); + + if(achunk->prev != NULL) + achunk = achunk->prev; + else + break; + } +} + +#ifdef DEBUG_ALLOC_TRACK +#define TEST_LEN 5 +int main(int argc, char **argv) { + + int i; + int *ptr; + void *addrs[TEST_LEN]; + + for(i=0; i < TEST_LEN; ++i) { + ptr = tracked_malloc(sizeof(int)); + addrs[i] = ptr; + *ptr = i; + } + + walk_chunks(); + + // Looks good. + + printf("Removing item in the (zero-indexed) 2nd position\n"); + tracked_free(addrs[2]); + walk_chunks(); + + // Looks good. + + printf("Removing item at the head of the list (most recent addition)\n"); + tracked_free(addrs[TEST_LEN-1]); + walk_chunks(); + + // Looks good. + + printf("Removing item at the tail of the list (oldest addition)\n"); + tracked_free(addrs[0]); + walk_chunks(); + + // Looks good. + + printf("Removing 2nd added item\n"); + tracked_free(addrs[1]); + walk_chunks(); + + printf("Removing 4th added item\n"); + tracked_free(addrs[3]); + walk_chunks(); + + tracked_free_all(); + + return 0; +} +#endif diff --git a/tools/arrays.c b/tools/arrays.c index 89837f316..8098f3db6 100644 --- a/tools/arrays.c +++ b/tools/arrays.c @@ -355,7 +355,7 @@ int array_spline( return _FAILURE_; } - u = malloc((n_lines-1) * sizeof(double)); + u = tracked_malloc((n_lines-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -444,7 +444,7 @@ int array_spline( *(array+k*n_columns+index_ddydx2) = *(array+k*n_columns+index_ddydx2) * *(array+(k+1)*n_columns+index_ddydx2) + u[k]; - free(u); + class_free(u); return _SUCCESS_; } @@ -469,7 +469,7 @@ int array_spline_table_line_to_line( errmsg, "no possible spline with less than three lines"); - u = malloc((n_lines-1) * sizeof(double)); + u = tracked_malloc((n_lines-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -545,7 +545,7 @@ int array_spline_table_line_to_line( *(array+k*n_columns+index_ddydx2) = *(array+k*n_columns+index_ddydx2) * *(array+(k+1)*n_columns+index_ddydx2) + u[k]; - free(u); + class_free(u); return _SUCCESS_; } @@ -571,10 +571,10 @@ int array_spline_table_lines( double dy_first; double dy_last; - u = malloc((x_size-1) * y_size * sizeof(double)); - p = malloc(y_size * sizeof(double)); - qn = malloc(y_size * sizeof(double)); - un = malloc(y_size * sizeof(double)); + u = tracked_malloc((x_size-1) * y_size * sizeof(double)); + p = tracked_malloc(y_size * sizeof(double)); + qn = tracked_malloc(y_size * sizeof(double)); + un = tracked_malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); @@ -704,10 +704,10 @@ int array_spline_table_lines( } } - free(qn); - free(un); - free(p); - free(u); + class_free(qn); + class_free(un); + class_free(p); + class_free(u); return _SUCCESS_; } @@ -733,10 +733,10 @@ int array_logspline_table_lines( double dy_first; double dy_last; - u = malloc((x_size-1) * y_size * sizeof(double)); - p = malloc(y_size * sizeof(double)); - qn = malloc(y_size * sizeof(double)); - un = malloc(y_size * sizeof(double)); + u = tracked_malloc((x_size-1) * y_size * sizeof(double)); + p = tracked_malloc(y_size * sizeof(double)); + qn = tracked_malloc(y_size * sizeof(double)); + un = tracked_malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -866,10 +866,10 @@ int array_logspline_table_lines( } } - free(qn); - free(un); - free(p); - free(u); + class_free(qn); + class_free(un); + class_free(p); + class_free(u); return _SUCCESS_; } @@ -895,10 +895,10 @@ int array_spline_table_columns( double dy_first; double dy_last; - u = malloc((x_size-1) * y_size * sizeof(double)); - p = malloc(y_size * sizeof(double)); - qn = malloc(y_size * sizeof(double)); - un = malloc(y_size * sizeof(double)); + u = tracked_malloc((x_size-1) * y_size * sizeof(double)); + p = tracked_malloc(y_size * sizeof(double)); + qn = tracked_malloc(y_size * sizeof(double)); + un = tracked_malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1037,10 +1037,10 @@ int array_spline_table_columns( } } - free(qn); - free(p); - free(u); - free(un); + class_free(qn); + class_free(p); + class_free(u); + class_free(un); return _SUCCESS_; } @@ -1066,10 +1066,10 @@ int array_spline_table_columns2( double dy_first; double dy_last; - u = malloc((x_size-1) * y_size * sizeof(double)); - p = malloc(y_size * sizeof(double)); - qn = malloc(y_size * sizeof(double)); - un = malloc(y_size * sizeof(double)); + u = tracked_malloc((x_size-1) * y_size * sizeof(double)); + p = tracked_malloc(y_size * sizeof(double)); + qn = tracked_malloc(y_size * sizeof(double)); + un = tracked_malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1176,10 +1176,10 @@ int array_spline_table_columns2( } } } - free(qn); - free(p); - free(u); - free(un); + class_free(qn); + class_free(p); + class_free(u); + class_free(un); return _SUCCESS_; } @@ -1205,7 +1205,7 @@ int array_spline_table_one_column( double dy_first; double dy_last; - u = malloc((x_size-1) * sizeof(double)); + u = tracked_malloc((x_size-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1313,7 +1313,7 @@ int array_spline_table_one_column( } - free(u); + class_free(u); return _SUCCESS_; } @@ -1340,7 +1340,7 @@ int array_logspline_table_one_column( double dy_first; double dy_last; - u = malloc((x_stop-1) * sizeof(double)); + u = tracked_malloc((x_stop-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1449,7 +1449,7 @@ int array_logspline_table_one_column( } - free(u); + class_free(u); return _SUCCESS_; } @@ -3131,7 +3131,7 @@ int array_smooth_trg(double * array, double weigth; double *coeff; - smooth=malloc(k_size*sizeof(double)); + smooth=tracked_malloc(k_size*sizeof(double)); if (smooth == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); return _FAILURE_; @@ -3243,8 +3243,8 @@ int array_smooth_trg(double * array, for (i=starting_k; iyscal); - free(pgi->y); - free(pgi->dydx); - - free(pgi->yerr); - free(pgi->ytempo); - - free(pgi->ak2); - free(pgi->ak3); - free(pgi->ak4); - free(pgi->ak5); - free(pgi->ak6); - free(pgi->ytemp); + class_free(pgi->yscal); + class_free(pgi->y); + class_free(pgi->dydx); + + class_free(pgi->yerr); + class_free(pgi->ytempo); + + class_free(pgi->ak2); + class_free(pgi->ak3); + class_free(pgi->ak4); + class_free(pgi->ak5); + class_free(pgi->ak6); + class_free(pgi->ytemp); return _SUCCESS_; } diff --git a/tools/evolver_ndf15.c b/tools/evolver_ndf15.c index 0d2937235..5c730934b 100644 --- a/tools/evolver_ndf15.c +++ b/tools/evolver_ndf15.c @@ -680,27 +680,27 @@ int evolver_ndf15( /** Deallocate memory */ - free(buffer); - - /* free(f0); */ - /* free(wt); */ - /* free(ddfddt); */ - /* free(pred); */ - /* free(y); */ - /* free(invwt); */ - /* free(rhs); */ - /* free(psi); */ - /* free(difkp1); */ - /* free(del); */ - /* free(yinterp); */ - /* free(ypinterp); */ - /* free(yppinterp); */ - /* free(tempvec1); */ - /* free(tempvec2); */ - - /* free(interpidx); */ - /* free(dif[1]); */ - /* free(dif); */ + class_free(buffer); + + /* class_free(f0); */ + /* class_free(wt); */ + /* class_free(ddfddt); */ + /* class_free(pred); */ + /* class_free(y); */ + /* class_free(invwt); */ + /* class_free(rhs); */ + /* class_free(psi); */ + /* class_free(difkp1); */ + /* class_free(del); */ + /* class_free(yinterp); */ + /* class_free(ypinterp); */ + /* class_free(yppinterp); */ + /* class_free(tempvec1); */ + /* class_free(tempvec2); */ + + /* class_free(interpidx); */ + /* class_free(dif[1]); */ + /* class_free(dif); */ uninitialize_jacobian(&jac); uninitialize_numjac_workspace(&nj_ws); @@ -1171,14 +1171,14 @@ int fzero_Newton(int (*func)(double *x, } } - free(p); - free(lu_work); - free(indx); - free(Fjac[1]); - free(Fjac); - free(F0); - free(delx); - free(Fdel); + class_free(p); + class_free(lu_work); + class_free(indx); + class_free(Fjac[1]); + class_free(Fjac); + class_free(F0); + class_free(delx); + class_free(Fdel); if (has_converged == _TRUE_){ return _SUCCESS_; @@ -1587,21 +1587,21 @@ int initialize_jacobian(struct jacobian *jac, int neq, ErrorMsg error_message){ } int uninitialize_jacobian(struct jacobian *jac){ - free(jac->dfdy[1]); - free(jac->dfdy); - free(jac->LU[1]); - free(jac->LU); + class_free(jac->dfdy[1]); + class_free(jac->dfdy); + class_free(jac->LU[1]); + class_free(jac->LU); - free(jac->luidx); - free(jac->LUw); - free(jac->jacvec); + class_free(jac->luidx); + class_free(jac->LUw); + class_free(jac->jacvec); if(jac->sparse_stuff_initialized){ - free(jac->xjac); - free(jac->col_wi); - free(jac->col_group); - free(jac->Cp); - free(jac->Ci); + class_free(jac->xjac); + class_free(jac->col_wi); + class_free(jac->col_group); + class_free(jac->Cp); + class_free(jac->Ci); sp_mat_free(jac->spJ); sp_num_free(jac->Numerical); } @@ -1637,20 +1637,20 @@ int initialize_numjac_workspace(struct numjac_workspace * nj_ws,int neq, ErrorMs int uninitialize_numjac_workspace(struct numjac_workspace * nj_ws){ /* Deallocate vectors and matrices: */ - free(nj_ws->yscale); - free(nj_ws->del); - free(nj_ws->Difmax); - free(nj_ws->absFdelRm); - free(nj_ws->absFvalue); - free(nj_ws->absFvalueRm); - free(nj_ws->Fscale); - free(nj_ws->ffdel); - free(nj_ws->yydel); - free(nj_ws->tmp); - - free(nj_ws->ydel_Fdel[1]); - free(nj_ws->ydel_Fdel); - free(nj_ws->logj); - free(nj_ws->Rowmax); + class_free(nj_ws->yscale); + class_free(nj_ws->del); + class_free(nj_ws->Difmax); + class_free(nj_ws->absFdelRm); + class_free(nj_ws->absFvalue); + class_free(nj_ws->absFvalueRm); + class_free(nj_ws->Fscale); + class_free(nj_ws->ffdel); + class_free(nj_ws->yydel); + class_free(nj_ws->tmp); + + class_free(nj_ws->ydel_Fdel[1]); + class_free(nj_ws->ydel_Fdel); + class_free(nj_ws->logj); + class_free(nj_ws->Rowmax); return _SUCCESS_; } diff --git a/tools/evolver_rkck.c b/tools/evolver_rkck.c index 7ecf412de..4f5dac72c 100644 --- a/tools/evolver_rkck.c +++ b/tools/evolver_rkck.c @@ -171,7 +171,7 @@ int evolver_rk(int (*derivs)(double x, gi.error_message, error_message); - free(dy); + class_free(dy); return _SUCCESS_; diff --git a/tools/growTable.c b/tools/growTable.c index 02ae9beea..fcf092407 100644 --- a/tools/growTable.c +++ b/tools/growTable.c @@ -54,7 +54,7 @@ int gt_add( if (ridx+sz>self->sz) { /** - test -> pass -> ok we need to grow */ - nbuffer=realloc(self->buffer,self->sz*_GT_FACTOR_); + nbuffer=tracked_realloc(self->buffer,self->sz*_GT_FACTOR_); class_test(nbuffer==NULL, self->error_message, "Cannot grow growTable"); @@ -151,7 +151,7 @@ int gt_getPtr( * Called by background_solve(). */ int gt_free(growTable* self) { - free(self->buffer); + class_free(self->buffer); self->csz=-1; self->sz=-1; self->freeze=_FALSE_; /**< This line added by JL */ diff --git a/tools/hyperspherical.c b/tools/hyperspherical.c index 5fd701080..20ab25bd9 100644 --- a/tools/hyperspherical.c +++ b/tools/hyperspherical.c @@ -228,12 +228,12 @@ int hyperspherical_HIS_create(int K, } } } - free(PhiL); + class_free(PhiL); } if (abort == _TRUE_) return _FAILURE_; - free(sqrtK); - free(one_over_sqrtK); + class_free(sqrtK); + class_free(one_over_sqrtK); for (k=0; kchi_at_phimin+k,NULL); @@ -269,13 +269,13 @@ int hyperspherical_update_pointers(HyperInterpStruct *pHIS_local, int hyperspherical_HIS_free(HyperInterpStruct *pHIS, ErrorMsg error_message){ /** Free the Hyperspherical Interpolation Structure. */ - free(pHIS->l); - free(pHIS->chi_at_phimin); - free(pHIS->x); - free(pHIS->sinK); - free(pHIS->cotK); - free(pHIS->phi); - free(pHIS->dphi); + class_free(pHIS->l); + class_free(pHIS->chi_at_phimin); + class_free(pHIS->x); + class_free(pHIS->sinK); + class_free(pHIS->cotK); + class_free(pHIS->phi); + class_free(pHIS->dphi); return _SUCCESS_; } diff --git a/tools/parser.c b/tools/parser.c index 602e2b729..ce43a5dcd 100644 --- a/tools/parser.c +++ b/tools/parser.c @@ -66,10 +66,10 @@ int parser_init(struct file_content * pfc, int parser_free(struct file_content * pfc) { if (pfc->size > 0) { - free(pfc->name); - free(pfc->value); - free(pfc->read); - free(pfc->filename); + class_free(pfc->name); + class_free(pfc->value); + class_free(pfc->read); + class_free(pfc->filename); } return _SUCCESS_; diff --git a/tools/quadrature.c b/tools/quadrature.c index 053cd0401..da2d214d8 100644 --- a/tools/quadrature.c +++ b/tools/quadrature.c @@ -31,8 +31,8 @@ int get_qsampling_manual(double *x, (*function)(params_for_function,x[i],&y); w[i] *= y; } - free(b); - free(c); + class_free(b); + class_free(c); return _SUCCESS_; case (qm_trapz) : for (i=0; ileft!=NULL) burn_tree(node->left); if (node->right!=NULL) burn_tree(node->right); - if (node->x!=NULL) free(node->x); - if (node->w!=NULL) free(node->w); - free(node); + if (node->x!=NULL) class_free(node->x); + if (node->w!=NULL) class_free(node->w); + class_free(node); } return _SUCCESS_; } diff --git a/tools/sparse.c b/tools/sparse.c index 2c8246fb8..ff0430ad9 100644 --- a/tools/sparse.c +++ b/tools/sparse.c @@ -27,10 +27,10 @@ int sp_mat_alloc(sp_mat** A, int ncols, int nrows, int maxnz, ErrorMsg error_mes } int sp_mat_free(sp_mat *A){ - free(A->Ax); - free(A->Ai); - free(A->Ap); - free(A); + class_free(A->Ax); + class_free(A->Ai); + class_free(A->Ap); + class_free(A); return _SUCCESS_; } @@ -62,15 +62,15 @@ int sp_num_alloc(sp_num** N, int n, ErrorMsg error_message){ int sp_num_free(sp_num *N){ sp_mat_free(N->L); sp_mat_free(N->U); - free(N->xi[0]); - free(N->xi); - free(N->topvec); - free(N->pinv); - free(N->p); - free(N->q); - free(N->w); - free(N->wamd); - free(N); + class_free(N->xi[0]); + class_free(N->xi); + class_free(N->topvec); + class_free(N->pinv); + class_free(N->p); + class_free(N->q); + class_free(N->w); + class_free(N->wamd); + class_free(N); return _SUCCESS_; } From 69a7934185995163bec29d42360a6e31fcbbc947 Mon Sep 17 00:00:00 2001 From: Operator Date: Thu, 18 Jul 2024 15:10:00 -0700 Subject: [PATCH 2/3] FIXED missing POSIX semantics in tracked_realloc() --- tools/alloc_track.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tools/alloc_track.c b/tools/alloc_track.c index 854203465..a93d266a2 100644 --- a/tools/alloc_track.c +++ b/tools/alloc_track.c @@ -254,6 +254,11 @@ void * tracked_realloc(void *ptr, size_t size) { void *newblock; int me = 0; struct chunk *list_cur; + + // Mimic POSIX semantics. + // realloc() is equivalent to malloc() if ptr is null + if(!ptr) + return tracked_alloc(size); #ifdef _OPENMP me = omp_get_thread_num(); From 9d99a7777aa5ab6ce2397607dd5d98c4cf104560 Mon Sep 17 00:00:00 2001 From: Operator Date: Thu, 18 Jul 2024 15:13:00 -0700 Subject: [PATCH 3/3] FIXED minor typo --- tools/alloc_track.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/alloc_track.c b/tools/alloc_track.c index a93d266a2..375e561b9 100644 --- a/tools/alloc_track.c +++ b/tools/alloc_track.c @@ -258,7 +258,7 @@ void * tracked_realloc(void *ptr, size_t size) { // Mimic POSIX semantics. // realloc() is equivalent to malloc() if ptr is null if(!ptr) - return tracked_alloc(size); + return tracked_malloc(size); #ifdef _OPENMP me = omp_get_thread_num();