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..375e561b9 --- /dev/null +++ b/tools/alloc_track.c @@ -0,0 +1,411 @@ +#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; + + // Mimic POSIX semantics. + // realloc() is equivalent to malloc() if ptr is null + if(!ptr) + return tracked_malloc(size); + +#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_; }