From a936c6fde022824963f6ac73f7a0469582666fcb Mon Sep 17 00:00:00 2001 From: Nils <32747859+schoeneberg@users.noreply.github.com> Date: Mon, 5 Feb 2024 11:31:32 +0100 Subject: [PATCH 1/3] Removing OpenMP from CLASS (#122) * Removing OpenMP from CLASS * Suggested changes according to Thomas Tram, as well as a few missed cases in lensing.c --------- Co-authored-by: schoeneberg --- .gitignore | 2 +- Makefile | 39 ++- include/common.h | 64 +++-- include/evolver_rkck.h | 4 +- include/parallel.h | 223 ++++++++++++++++ python/setup.py | 5 +- source/harmonic.c | 168 +++++------- source/lensing.c | 579 +++++++++++++++++++++++------------------ source/perturbations.c | 246 +++++------------ source/primordial.c | 76 +----- source/transfer.c | 193 +++++--------- tools/arrays.c | 69 ++--- tools/hyperspherical.c | 197 +++++++------- 13 files changed, 957 insertions(+), 908 deletions(-) create mode 100644 include/parallel.h diff --git a/.gitignore b/.gitignore index bf5fbf37..3a53eef1 100644 --- a/.gitignore +++ b/.gitignore @@ -6,7 +6,7 @@ class *.ini libclass.a *out -python/classy.c +python/classy.cpp *.pyc python/build/ *.dSYM diff --git a/Makefile b/Makefile index 9e4a8a77..10f3126e 100644 --- a/Makefile +++ b/Makefile @@ -11,6 +11,7 @@ WRKDIR = $(MDIR)/build vpath %.c source:tools:main:test vpath %.o build +vpath %.opp build vpath .base build ######################################################## @@ -21,6 +22,7 @@ vpath .base build CC = gcc #CC = icc #CC = pgcc +CPP = g++ --std=c++11 -fpermissive -Wno-write-strings # your tool for creating static libraries: AR = ar rv @@ -38,7 +40,7 @@ OPTFLAG = -O3 #OPTFLAG = -fast # your openmp flag (comment for compiling without openmp) -OMPFLAG = -fopenmp +OMPFLAG = -pthread #-fopenmp #OMPFLAG = -mp -mp=nonuma -mp=allcores -g #OMPFLAG = -openmp @@ -91,9 +93,12 @@ 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 +%.opp: %.c .base $(HEADERFILES) + cd $(WRKDIR);$(CPP) $(OPTFLAG) $(OMPFLAG) $(CCFLAG) $(INCLUDES) -c ../$< -o $*.opp -SOURCE = input.o background.o thermodynamics.o perturbations.o primordial.o fourier.o transfer.o harmonic.o lensing.o distortions.o +TOOLS = growTable.o dei_rkck.o sparse.o evolver_rkck.o evolver_ndf15.o arrays.opp parser.o quadrature.o hyperspherical.opp common.o trigonometric_integrals.o + +SOURCE = input.o background.o thermodynamics.o perturbations.opp primordial.opp fourier.o transfer.opp harmonic.opp lensing.opp distortions.o INPUT = input.o @@ -156,51 +161,43 @@ libclass.a: $(TOOLS) $(SOURCE) $(EXTERNAL) $(AR) $@ $(addprefix build/, $(TOOLS) $(SOURCE) $(EXTERNAL)) class: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(CLASS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o class $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o class $(addprefix build/,$(notdir $^)) -lm test_loops: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_loops_omp: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS_OMP) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_harmonic: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_HARMONIC) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_transfer: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_TRANSFER) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_fourier: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_FOURIER) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_perturbations: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_PERTURBATIONS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_thermodynamics: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_THERMODYNAMICS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_background: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_BACKGROUND) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + $(CPP) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm test_hyperspherical: $(TOOLS) $(TEST_HYPERSPHERICAL) $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o test_hyperspherical $(addprefix build/,$(notdir $^)) -lm - tar: $(C_ALL) $(C_TEST) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) classy: libclass.a python/classy.pyx python/cclassy.pxd -ifdef OMPFLAG - cp python/setup.py python/autosetup.py -else - grep -v "lgomp" python/setup.py > python/autosetup.py -endif - cd python; export CC=$(CC); $(PYTHON) autosetup.py install || $(PYTHON) autosetup.py install --user - rm python/autosetup.py + cd python; export CC=$(CC); $(PYTHON) setup.py install || $(PYTHON) setup.py install --user clean: .base rm -rf $(WRKDIR); rm -f libclass.a rm -f $(MDIR)/python/classy.c rm -rf $(MDIR)/python/build - rm -f python/autosetup.py diff --git a/include/common.h b/include/common.h index 6d0b6a49..efda051b 100644 --- a/include/common.h +++ b/include/common.h @@ -79,19 +79,24 @@ typedef char FileName[_FILENAMESIZE_]; /* @endcond */ -/* needed because of weird openmp bug on macosx lion... */ -void class_protect_sprintf(char* dest, char* tpl,...); -void class_protect_fprintf(FILE* dest, char* tpl,...); -void* class_protect_memcpy(void* dest, void* from, size_t sz); - -/* some general functions */ - -int get_number_of_titles(char * titlestring); -int file_exists(const char *fname); -int compare_doubles(const void * a, - const void * b); -int string_begins_with(char* thestring, char beginchar); +#ifdef __cplusplus +extern "C" { +#endif + /* needed because of weird openmp bug on macosx lion... */ + void class_protect_sprintf(char* dest, char* tpl, ...); + void class_protect_fprintf(FILE* dest, char* tpl, ...); + void* class_protect_memcpy(void* dest, void* from, size_t sz); + + /* some general functions */ + int get_number_of_titles(char * titlestring); + int file_exists(const char *fname); + int compare_doubles(const void * a, + const void * b); + int string_begins_with(char* thestring, char beginchar); +#ifdef __cplusplus +} +#endif /* general CLASS macros */ @@ -128,15 +133,15 @@ int string_begins_with(char* thestring, char beginchar); #define class_call(function, error_message_from_function, error_message_output) \ class_call_except(function, error_message_from_function,error_message_output,) -/* same in parallel region */ -#define class_call_parallel(function, error_message_from_function, error_message_output) { \ - if (abort == _FALSE_) { \ +/* same in parallel region -- UNUSED NOW */ +/*#define class_call_parallel(function, error_message_from_function, error_message_output) { \ + if (abort_now == _FALSE_) { \ if (function == _FAILURE_) { \ class_call_message(error_message_output,#function,error_message_from_function); \ - abort=_TRUE_; \ + abort_now=_TRUE_; \ } \ } \ -} +}*/ @@ -147,7 +152,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=(__typeof__(pointer))malloc(size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -156,23 +161,24 @@ int string_begins_with(char* thestring, char beginchar); } \ } -/* same inside parallel structure */ + +/* same inside parallel structure -- UNUSED NOW #define class_alloc_parallel(pointer, size, error_message_output) { \ pointer=NULL; \ - if (abort == _FALSE_) { \ - pointer=malloc(size); \ + if (abort_now == _FALSE_) { \ + pointer=(__typeof__(pointer))malloc(size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ class_alloc_message(error_message_output,#pointer, size_int); \ - abort=_TRUE_; \ + abort_now=_TRUE_; \ } \ } \ -} +}*/ /* 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=(__typeof__(pointer))calloc(init,size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -183,7 +189,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=(__typeof__(pointer))realloc(newname,size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -218,14 +224,15 @@ int string_begins_with(char* thestring, char beginchar); } \ } +/* UNUSED NOW #define class_test_parallel(condition, error_message_output, args...) { \ - if (abort == _FALSE_) { \ + if (abort_now == _FALSE_) { \ if (condition) { \ class_test_message(error_message_output,#condition, args); \ - abort=_TRUE_; \ + abort_now=_TRUE_; \ } \ } \ -} +}*/ /* macro for returning error message; args is a variable list of optional arguments, e.g.: args="x=%d",x @@ -404,4 +411,5 @@ struct precision }; + #endif diff --git a/include/evolver_rkck.h b/include/evolver_rkck.h index 9da76c0f..1106a384 100644 --- a/include/evolver_rkck.h +++ b/include/evolver_rkck.h @@ -1,5 +1,5 @@ -#ifndef __EVO__ -#define __EVO__ +#ifndef __EVO_RK__ +#define __EVO_RK__ #include "dei_rkck.h" diff --git a/include/parallel.h b/include/parallel.h new file mode 100644 index 00000000..da18ee78 --- /dev/null +++ b/include/parallel.h @@ -0,0 +1,223 @@ +// 'MAGIC' parallelization module written by Nils Schöneberg, adapted from the one created by Thomas Tram + +//To be used for declaring which arguments a parallel clause 'captures' from the outside +//This would be equivalent to the 'firstprivate' in OPENMP +#define with_arguments(...) __VA_ARGS__ +//To declare a list of arguments a parallel clause only has internally +//(necessary since comma-separated lists are interpreted as separate arguments +// except when surrounded by backets as forced by this macro) +#define declare_list_of_variables_inside_parallel_region(...) __VA_ARGS__ + +// To be called WITHIN a given parallel loop. First argument is the parameters to 'capture', +// corresponding to the 'firstprivate' in OPENMP. +// Usually either '=' to capture all necessary variables, or a list of arguments +// declared within the 'with_arguments' macro for a more precise list +// The second argument is the actual code to execute. That code should be formulated like +// the body of a normal function, e.g. return _SUCCESS_ at the end. +// 'private' arguments in OPENMP can be emulated by just making them local to the function. +// Special care: A list of private arguments needs to be surrounded by the special macro +// 'declare_list_of_variables-inside_parallel_region' due to the peculiarities +// of the C++ preprocessor macros. See examples e.g. +// in source/perturbations.c, source/lensing.c, or tools/hypershperical.c +#define class_run_parallel(arg1, arg2) future_output.push_back(task_system.AsyncTask([arg1] () {arg2})); + +// To be called ONLY ONCE without arguments before the intended parallel loop(s). +#define class_setup_parallel() \ +Tools::TaskSystem task_system; \ +std::vector> future_output; + +// To be called without arguments AFTER ANY parallel loop in order to actually execute the jobs. +// NEEDS TO BE CALLED BEFORE USING THE RESULTS! +#define class_finish_parallel() \ +for (std::future& future : future_output) { \ + if(future.get()!=_SUCCESS_) return _FAILURE_; \ +} \ +future_output.clear(); + +// +// thread_pool.h +// ppCLASS +// +// Created by Thomas Tram on 02/03/2020. +// Copyright 2020 Aarhus University. All rights reserved. +// +#ifndef THREAD_POOL_H +#define THREAD_POOL_H +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Tools { + +class NotificationQueue { +public: + bool TryPop(std::function& x) { + std::unique_lock lock(mutex_, std::try_to_lock); + if (!lock || queue_.empty()) { + return false; + } + x = std::move(queue_.front()); + queue_.pop_front(); + return true; + } + + bool Pop(std::function& x) { + std::unique_lock lock(mutex_); + while (queue_.empty() && !done_) { + ready_.wait(lock); + } + if (queue_.empty()) { + return false; + } + x = std::move(queue_.front()); + queue_.pop_front(); + return true; + } + + template + bool TryPush(F&& f) { + { + std::unique_lock lock(mutex_, std::try_to_lock); + if (!lock) { + return false; + } + queue_.emplace_back(std::forward(f)); + } + ready_.notify_one(); + return true; + } + + template + void Push(F&& f) { + { + std::unique_lock lock(mutex_); + queue_.emplace_back(std::forward(f)); + } + ready_.notify_one(); + } + + void Done() { + { + std::unique_lock lock(mutex_); + done_ = true; + } + ready_.notify_all(); + } +private: + std::deque> queue_; + bool done_ = false; + std::mutex mutex_; + std::condition_variable ready_; +}; + +class TaskSystem { +public: + TaskSystem(unsigned int count = GetNumThreads()) + : count_(count) + , index_(0) + , queues_{count_} { + for (unsigned int n = 0; n < count_; ++n) { + threads_.emplace_back([&, n]{ Run(n); }); + } + } + + ~TaskSystem() { + for (auto& e : queues_) e.Done(); + for (auto& e : threads_) e.join(); + } + + static unsigned int GetNumThreads() { + unsigned int number_of_threads = std::thread::hardware_concurrency(); + for (const std::string& env_var_name : {"OMP_NUM_THREADS", "SLURM_CPUS_PER_TASK"}) { + if (char* s = std::getenv(env_var_name.c_str())) { + int threads = std::atoi(s); + if ((threads > 0) && (threads <= 8192)) { + number_of_threads = threads; + break; + } + } + } + return number_of_threads; + } + + template + void Async(F&& f) { + auto i = index_++; + for (unsigned int n = 0; n < count_; ++n) { + if (queues_[(i + n) % count_].TryPush(std::forward(f))) { + return; + } + } + queues_[i % count_].Push(std::forward(f)); + } + + template + std::future::type> AsyncTask(F&& f) { + using return_type = typename std::result_of::type; + auto task = std::make_shared>(f); + std::future res = task->get_future(); + + auto work = [task](){ (*task)(); }; + unsigned int i = index_++; + for(unsigned int n = 0; n < count_; ++n) { + if(queues_[(i + n) % count_].TryPush(work)){ + return res; + } + } + queues_[i % count_].Push(work); + + return res; + } + + template + std::future::type> AsyncTask(F&& f, Args&&... args) { + using return_type = typename std::result_of::type; + auto task = std::make_shared>(std::bind(std::forward(f), std::forward(args)...)); + std::future res = task->get_future(); + + auto work = [task](){ (*task)(); }; + unsigned int i = index_++; + for(unsigned int n = 0; n < count_; ++n) { + if(queues_[(i + n) % count_].TryPush(work)){ + return res; + } + } + queues_[i % count_].Push(work); + + return res; + } + + unsigned int get_num_threads(){ + return count_; + } + +private: + void Run(unsigned int i) { + while (true) { + std::function f; + for (unsigned n = 0; n != count_; ++n) { + if (queues_[(i + n) % count_].TryPop(f)) { + break; + } + } + if (!f && !queues_[i].Pop(f)) { + break; + } + f(); + } + } + + const unsigned int count_; + std::vector threads_; + std::vector queues_; + std::atomic index_; +}; + +} +#endif diff --git a/python/setup.py b/python/setup.py index 0d0d2b9f..5689a20a 100644 --- a/python/setup.py +++ b/python/setup.py @@ -41,7 +41,9 @@ include_dirs=[nm.get_include(), include_folder, heat_folder, recfast_folder, hyrec_folder], libraries=liblist, library_dirs=[root_folder, GCCPATH], - extra_link_args=['-lgomp'] + #extra_link_args=['-lgomp'], + language="c++", + extra_compile_args=["-std=c++11"] ) import sys classy_ext.cython_directives = {'language_level': "3" if sys.version_info.major>=3 else "2"} @@ -53,5 +55,4 @@ url='http://www.class-code.net', cmdclass={'build_ext': build_ext}, ext_modules=[classy_ext], - #data_files=[('bbn', ['../bbn/sBBN.dat'])] ) diff --git a/source/harmonic.c b/source/harmonic.c index a29d41f6..670834e9 100644 --- a/source/harmonic.c +++ b/source/harmonic.c @@ -13,6 +13,7 @@ */ #include "harmonic.h" +#include "parallel.h" /** * Anisotropy power spectra \f$ C_l\f$'s for all types, modes and initial conditions. @@ -678,25 +679,6 @@ int harmonic_cls( int index_ct; int cl_integrand_num_columns; - double * cl_integrand; /* array with argument cl_integrand[index_k*cl_integrand_num_columns+1+phr->index_ct] */ - double * cl_integrand_limber; /* similar array with same columns but different number of lines (less k values) */ - - double * transfer_ic1; /* array with argument transfer_ic1[index_tt] */ - double * transfer_ic2; /* idem */ - double * primordial_pk; /* array with argument primordial_pk[index_ic_ic]*/ - - /* This code can be optionally compiled with the openmp option for parallel computation. - Inside parallel regions, the use of the command "return" is forbidden. - For error management, instead of "return _FAILURE_", we will set the variable below - to "abort = _TRUE_". This will lead to a "return _FAILURE_" jus after leaving the - parallel region. */ - int abort; - -#ifdef _OPENMP - /* instrumentation times */ - double tstart, tstop; -#endif - /** - allocate pointers to arrays where results will be stored */ class_alloc(phr->l_size,sizeof(int)*phr->md_size,phr->error_message); @@ -727,6 +709,8 @@ int harmonic_cls( /** - --> (c) loop over initial conditions */ + class_setup_parallel(); + for (index_ic1 = 0; index_ic1 < phr->ic_size[index_md]; index_ic1++) { for (index_ic2 = index_ic1; index_ic2 < phr->ic_size[index_md]; index_ic2++) { index_ic1_ic2 = index_symmetric_matrix(index_ic1,index_ic2,phr->ic_size[index_md]); @@ -734,97 +718,75 @@ int harmonic_cls( /* non-diagonal coefficients should be computed only if non-zero correlation */ if (phr->is_non_zero[index_md][index_ic1_ic2] == _TRUE_) { - /* initialize error management flag */ - abort = _FALSE_; - - /* beginning of parallel region */ - -#pragma omp parallel \ - shared(ppr,ptr,ppm,index_md,phr,ppt,cl_integrand_num_columns,index_ic1,index_ic2,abort) \ - private(tstart,cl_integrand,cl_integrand_limber,primordial_pk,transfer_ic1,transfer_ic2,index_l,tstop) - - { + /** - ---> loop over l values defined in the transfer module. + For each l, compute the \f$ C_l\f$'s for all types (TT, TE, ...) + by convolving primordial spectra with transfer functions. + This elementary task is assigned to harmonic_compute_cl() */ -#ifdef _OPENMP - tstart = omp_get_wtime(); -#endif - - class_alloc_parallel(cl_integrand, - ptr->q_size*cl_integrand_num_columns*sizeof(double), - phr->error_message); + for (index_l=0; index_l < ptr->l_size[index_md]; index_l++) { - cl_integrand_limber = NULL; - if (ptr->do_lcmb_full_limber == _TRUE_) { - class_alloc_parallel(cl_integrand_limber, - ptr->q_size_limber*cl_integrand_num_columns*sizeof(double), - phr->error_message); - } + class_run_parallel(=, - class_alloc_parallel(primordial_pk, - phr->ic_ic_size[index_md]*sizeof(double), - phr->error_message); - - class_alloc_parallel(transfer_ic1, - ptr->tt_size[index_md]*sizeof(double), - phr->error_message); - - class_alloc_parallel(transfer_ic2, - ptr->tt_size[index_md]*sizeof(double), - phr->error_message); - -#pragma omp for schedule (dynamic) - - /** - ---> loop over l values defined in the transfer module. - For each l, compute the \f$ C_l\f$'s for all types (TT, TE, ...) - by convolving primordial spectra with transfer functions. - This elementary task is assigned to harmonic_compute_cl() */ - - for (index_l=0; index_l < ptr->l_size[index_md]; index_l++) { - -#pragma omp flush(abort) - - class_call_parallel(harmonic_compute_cl(ppr, - pba, - ppt, - ptr, - ppm, - phr, - index_md, - index_ic1, - index_ic2, - index_l, - cl_integrand_num_columns, - cl_integrand, - cl_integrand_limber, - primordial_pk, - transfer_ic1, - transfer_ic2), - phr->error_message, - phr->error_message); - - } /* end of loop over l */ - -#ifdef _OPENMP - tstop = omp_get_wtime(); - if (phr->harmonic_verbose > 1) - 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); - - if (ptr->do_lcmb_full_limber == _TRUE_) { - free(cl_integrand_limber); - } + double * cl_integrand; /* array with argument cl_integrand[index_k*cl_integrand_num_columns+1+phr->index_ct] */ + double * cl_integrand_limber; /* similar array with same columns but different number of lines (less k values) */ + double * transfer_ic1; /* array with argument transfer_ic1[index_tt] */ + double * transfer_ic2; /* idem */ + double * primordial_pk; /* array with argument primordial_pk[index_ic_ic]*/ - free(primordial_pk); - free(transfer_ic1); + class_alloc(cl_integrand, + ptr->q_size*cl_integrand_num_columns*sizeof(double), + phr->error_message); - free(transfer_ic2); + cl_integrand_limber = NULL; + if (ptr->do_lcmb_full_limber == _TRUE_) { + class_alloc(cl_integrand_limber, + ptr->q_size_limber*cl_integrand_num_columns*sizeof(double), + phr->error_message); + } - } /* end of parallel region */ + class_alloc(primordial_pk, + phr->ic_ic_size[index_md]*sizeof(double), + phr->error_message); + + class_alloc(transfer_ic1, + ptr->tt_size[index_md]*sizeof(double), + phr->error_message); + + class_alloc(transfer_ic2, + ptr->tt_size[index_md]*sizeof(double), + phr->error_message); + + class_call(harmonic_compute_cl(ppr, + pba, + ppt, + ptr, + ppm, + phr, + index_md, + index_ic1, + index_ic2, + index_l, + cl_integrand_num_columns, + cl_integrand, + cl_integrand_limber, + primordial_pk, + transfer_ic1, + transfer_ic2), + phr->error_message, + phr->error_message); + + free(cl_integrand); + if (ptr->do_lcmb_full_limber == _TRUE_) { + free(cl_integrand_limber); + } + free(primordial_pk); + free(transfer_ic1); + free(transfer_ic2); - if (abort == _TRUE_) return _FAILURE_; + return _SUCCESS_; + ); + } /* end of loop over l */ } else { @@ -842,6 +804,8 @@ int harmonic_cls( } } + class_finish_parallel(); + /** - --> (d) now that for a given mode, all possible \f$ C_l\f$'s have been computed, compute second derivative of the array in which they are stored, in view of spline interpolation. */ diff --git a/source/lensing.c b/source/lensing.c index 29b01d5c..c87ca9bc 100644 --- a/source/lensing.c +++ b/source/lensing.c @@ -17,6 +17,7 @@ #include "lensing.h" #include +#include "parallel.h" /** * Anisotropy power spectra \f$ C_l\f$'s for all types, modes and initial conditions. @@ -495,30 +496,37 @@ int lensing_init( /** - Compute sigma2\f$(\mu)\f$ and Cgl2(\f$\mu\f$) **/ - //debut = omp_get_wtime(); -#pragma omp parallel for \ - private (index_mu,l) \ - schedule (static) + class_setup_parallel(); + for (index_mu=0; index_mul_unlensed_max; + class_run_parallel(with_arguments(index_mu,l_unlensed_max,Cgl,Cgl2,cl_pp,d11,d1m1), + int l; - for (l=2; l<=ple->l_unlensed_max; l++) { + Cgl[index_mu]=0; + Cgl2[index_mu]=0; - Cgl[index_mu] += (2.*l+1.)*l*(l+1.)* - cl_pp[l]*d11[index_mu][l]; + for (l=2; l<=l_unlensed_max; l++) { - Cgl2[index_mu] += (2.*l+1.)*l*(l+1.)* - cl_pp[l]*d1m1[index_mu][l]; + Cgl[index_mu] += (2.*l+1.)*l*(l+1.)* + cl_pp[l]*d11[index_mu][l]; - } + Cgl2[index_mu] += (2.*l+1.)*l*(l+1.)* + cl_pp[l]*d1m1[index_mu][l]; - Cgl[index_mu] /= 4.*_PI_; - Cgl2[index_mu] /= 4.*_PI_; + } + + Cgl[index_mu] /= 4.*_PI_; + Cgl2[index_mu] /= 4.*_PI_; + return _SUCCESS_; + ); } + class_finish_parallel(); + for (index_mu=0; index_mul_unlensed_max;l++) { ll = (double)l; @@ -687,7 +694,12 @@ int lensing_init( ksim[index_mu] += resm; } } + return _SUCCESS_; + + ); } + + class_finish_parallel(); //fin = omp_get_wtime(); //cpu_time = (fin-debut); //printf("time in ksi=%4.3f s\n",cpu_time); @@ -1051,23 +1063,26 @@ int lensing_lensed_cl_tt( struct lensing * ple ) { - double cle; - int imu; int index_l; /** Integration by Gauss-Legendre quadrature. **/ -#pragma omp parallel for \ - private (imu,index_l,cle) \ - schedule (static) + class_setup_parallel(); for (index_l=0; index_ll_size; index_l++){ - cle=0; - for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ - } - ple->cl_lens[index_l*ple->lt_size+ple->index_lt_tt]=cle*2.0*_PI_; + class_run_parallel(=, + double cle; + int imu; + cle=0; + for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ + } + ple->cl_lens[index_l*ple->lt_size+ple->index_lt_tt]=cle*2.0*_PI_; + return _SUCCESS_; + ); } + class_finish_parallel(); + return _SUCCESS_; } @@ -1114,23 +1129,25 @@ int lensing_lensed_cl_te( struct lensing * ple ) { - double clte; - int imu; int index_l; /** Integration by Gauss-Legendre quadrature. **/ -#pragma omp parallel for \ - private (imu,index_l,clte) \ - schedule (static) + class_setup_parallel(); for (index_l=0; index_l < ple->l_size; index_l++){ - clte=0; - for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ - } - ple->cl_lens[index_l*ple->lt_size+ple->index_lt_te]=clte*2.0*_PI_; + class_run_parallel(=, + double clte; + int imu; + clte=0; + for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ + } + ple->cl_lens[index_l*ple->lt_size+ple->index_lt_te]=clte*2.0*_PI_; + return _SUCCESS_; + ); } + class_finish_parallel(); return _SUCCESS_; } @@ -1181,24 +1198,26 @@ int lensing_lensed_cl_ee_bb( struct lensing * ple ) { - double clp, clm; - int imu; int index_l; + class_setup_parallel(); /** Integration by Gauss-Legendre quadrature. **/ -#pragma omp parallel for \ - private (imu,index_l,clp,clm) \ - schedule (static) - for (index_l=0; index_l < ple->l_size; index_l++){ - clp=0; clm=0; - for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ - clm += ksim[imu]*d2m2[imu][(int)ple->l[index_l]]*w8[imu]; /* loop could be optimized */ - } - ple->cl_lens[index_l*ple->lt_size+ple->index_lt_ee]=(clp+clm)*_PI_; - ple->cl_lens[index_l*ple->lt_size+ple->index_lt_bb]=(clp-clm)*_PI_; + class_run_parallel(=, + double clp; + double clm; + int imu; + clp=0; clm=0; + for (imu=0;imul[index_l]]*w8[imu]; /* loop could be optimized */ + clm += ksim[imu]*d2m2[imu][(int)ple->l[index_l]]*w8[imu]; /* loop could be optimized */ + } + ple->cl_lens[index_l*ple->lt_size+ple->index_lt_ee]=(clp+clm)*_PI_; + ple->cl_lens[index_l*ple->lt_size+ple->index_lt_bb]=(clp-clm)*_PI_; + return _SUCCESS_; + ); } + class_finish_parallel(); return _SUCCESS_; } @@ -1249,7 +1268,7 @@ int lensing_d00( int lmax, double ** d00 ) { - double ll, dlm1, dl, dlp1; + double ll; int index_mu, l; double *fac1, *fac2, *fac3; ErrorMsg erreur; @@ -1264,24 +1283,27 @@ int lensing_d00( fac3[l] = sqrt(2./(2*ll+3)); } -#pragma omp parallel for \ - private (index_mu,dlm1,dl,dlp1,l,ll) \ - schedule (static) - + class_setup_parallel(); for (index_mu=0;index_muerror_message, ppt->error_message); - /** - create an array of workspaces in multi-thread case */ - -#ifdef _OPENMP - -#pragma omp parallel - { - number_of_threads = omp_get_num_threads(); - } -#endif - - class_alloc(pppw,number_of_threads * sizeof(struct perturbations_workspace *),ppt->error_message); - + /* Setup task system */ + class_setup_parallel(); /** - loop over modes (scalar, tensors, etc). For each mode: */ - for (index_md = 0; index_md < ppt->md_size; index_md++) { if (ppt->perturbations_verbose > 1) printf("Evolving mode %d/%d\n",index_md+1,ppt->md_size); - abort = _FALSE_; - - sz = sizeof(struct perturbations_workspace); - -#pragma omp parallel \ - shared(pppw,ppr,pba,pth,ppt,index_md,abort,number_of_threads) \ - private(thread) \ - num_threads(number_of_threads) - - { - -#ifdef _OPENMP - thread=omp_get_thread_num(); -#endif - - /** - --> (a) create a workspace (one per thread in multi-thread case) */ - - class_alloc_parallel(pppw[thread],sz,ppt->error_message); - - /** - --> (b) initialize indices of vectors of perturbations with perturbations_indices_of_current_vectors() */ - - class_call_parallel(perturbations_workspace_init(ppr, - pba, - pth, - ppt, - index_md, - pppw[thread]), - ppt->error_message, - ppt->error_message); - - } /* end of parallel region */ - - if (abort == _TRUE_) return _FAILURE_; - - /** - --> (c) loop over initial conditions and wavenumbers; for each of them, evolve perturbations and compute source functions with perturbations_solve() */ + /** - --> loop over initial conditions and wavenumbers; for each of them, evolve perturbations and compute source functions with perturbations_solve() */ for (index_ic = 0; index_ic < ppt->ic_size[index_md]; index_ic++) { @@ -973,136 +910,88 @@ int perturbations_init( printf("evolving %d wavenumbers\n",ppt->k_size[index_md]); } - abort = _FALSE_; - -#pragma omp parallel \ - shared(pppw,ppr,pba,pth,ppt,index_md,index_ic,abort,number_of_threads) \ - private(index_k,thread,tstart,tstop,tspent) \ - num_threads(number_of_threads) - - { + /* integrating backwards is slightly more optimal for parallel runs */ + //for (index_k = 0; index_k < ppt->k_size; index_k++) { + for (index_k = ppt->k_size[index_md]-1; index_k >=0; index_k--) { -#ifdef _OPENMP - thread=omp_get_thread_num(); - tspent=0.; -#endif + class_run_parallel(with_arguments(ppr,pba,pth,ppt,index_md,index_ic,index_k), -#pragma omp for schedule (dynamic) - - /* integrating backwards is slightly more optimal for parallel runs */ - //for (index_k = 0; index_k < ppt->k_size; index_k++) { - for (index_k = ppt->k_size[index_md]-1; index_k >=0; index_k--) { - - if ((ppt->perturbations_verbose > 2) && (abort == _FALSE_)) { + if (ppt->perturbations_verbose > 2) { printf("evolving mode k=%e /Mpc (%d/%d)",ppt->k[index_md][index_k],index_k+1,ppt->k_size[index_md]); if (pba->sgnK != 0) printf(" (for scalar modes, corresponds to nu=%e)",sqrt(ppt->k[index_md][index_k]*ppt->k[index_md][index_k]+pba->K)/sqrt(pba->sgnK*pba->K)); printf("\n"); } -#ifdef _OPENMP - tstart = omp_get_wtime(); -#endif - - class_call_parallel(perturbations_solve(ppr, + struct perturbations_workspace pw; + class_call(perturbations_workspace_init(ppr, pba, pth, ppt, index_md, - index_ic, - index_k, - pppw[thread]), - ppt->error_message, - ppt->error_message); - -#ifdef _OPENMP - tstop = omp_get_wtime(); - - tspent += tstop-tstart; -#endif - -#pragma omp flush(abort) + &pw), + ppt->error_message, + ppt->error_message); - } /* end of loop over wavenumbers */ + class_call(perturbations_solve(ppr, + pba, + pth, + ppt, + index_md, + index_ic, + index_k, + &pw), + ppt->error_message, + ppt->error_message); -#ifdef _OPENMP - if (ppt->perturbations_verbose>2) - printf("In %s: time spent in parallel region (loop over k's) = %e s for thread %d\n", - __func__,tspent,omp_get_thread_num()); -#endif + class_call(perturbations_workspace_free(ppt,index_md,&pw), + ppt->error_message, + ppt->error_message); + return _SUCCESS_; - } /* end of parallel region */ + ); - if (abort == _TRUE_) return _FAILURE_; + } /* end of loop over wavenumbers */ } /* end of loop over initial conditions */ - abort = _FALSE_; - -#pragma omp parallel \ - shared(pppw,ppt,index_md,abort,number_of_threads) \ - private(thread) \ - num_threads(number_of_threads) - - { - -#ifdef _OPENMP - thread=omp_get_thread_num(); -#endif - - class_call_parallel(perturbations_workspace_free(ppt,index_md,pppw[thread]), - ppt->error_message, - ppt->error_message); - - } /* end of parallel region */ - - if (abort == _TRUE_) return _FAILURE_; - } /* end loop over modes */ - free(pppw); + class_finish_parallel(); + /** - spline the source array with respect to the time variable */ if (ppt->ln_tau_size > 1) { + class_setup_parallel(); + for (index_md = 0; index_md < ppt->md_size; index_md++) { for (index_ic = 0; index_ic < ppt->ic_size[index_md]; index_ic++) { - abort = _FALSE_; - -#pragma omp parallel \ - shared(ppt,index_md,index_ic,abort,number_of_threads) \ - private(index_tp) \ - num_threads(number_of_threads) - - { - -#pragma omp for schedule (dynamic) - - for (index_tp = 0; index_tp < ppt->tp_size[index_md]; index_tp++) { - - class_call_parallel(array_spline_table_lines(ppt->ln_tau, - ppt->ln_tau_size, - ppt->late_sources[index_md][index_ic * ppt->tp_size[index_md] + index_tp], - ppt->k_size[index_md], - ppt->ddlate_sources[index_md][index_ic*ppt->tp_size[index_md] + index_tp], - _SPLINE_EST_DERIV_, - ppt->error_message), - ppt->error_message, - ppt->error_message); - - } + for (index_tp = 0; index_tp < ppt->tp_size[index_md]; index_tp++) { - } /* end of parallel region */ + class_run_parallel(with_arguments(ppt, index_md, index_ic, index_tp), + class_call(array_spline_table_lines(ppt->ln_tau, + ppt->ln_tau_size, + ppt->late_sources[index_md][index_ic * ppt->tp_size[index_md] + index_tp], + ppt->k_size[index_md], + ppt->ddlate_sources[index_md][index_ic*ppt->tp_size[index_md] + index_tp], + _SPLINE_EST_DERIV_, + ppt->error_message), + ppt->error_message, + ppt->error_message); + return _SUCCESS_; + ); - if (abort == _TRUE_) return _FAILURE_; + } /* end of loop over type of source function*/ } /* end of loop over initial condition */ } /* end of loop over mode */ + class_finish_parallel(); } return _SUCCESS_; @@ -2978,8 +2867,6 @@ int perturbations_workspace_free ( } } - free(ppw); - return _SUCCESS_; } @@ -3068,13 +2955,12 @@ int perturbations_solve( /* function pointer to ODE evolver and names of possible evolvers */ - extern int evolver_rk(); - extern int evolver_ndf15(); - int (*generic_evolver)(); + + auto generic_evolver = &(evolver_ndf15); /* Related to the perturbation output */ - int (*perhaps_print_variables)(); + int (*perhaps_print_variables)(double, double*, double*, void*, char*); int index_ikout; /** - initialize indices relevant for back/thermo tables search */ @@ -6163,7 +6049,7 @@ int perturbations_approximations( /** - evaluate background quantities with background_at_tau() and Hubble time scale \f$ \tau_h = a/a' \f$ */ - class_call(background_at_tau(pba,tau, normal_info, ppw->inter_mode, &(ppw->last_index_back), ppw->pvecback), + class_call(background_at_tau(pba,tau, normal_info, (interpolation_method)ppw->inter_mode, &(ppw->last_index_back), ppw->pvecback), pba->error_message, ppt->error_message); @@ -6182,7 +6068,7 @@ int perturbations_approximations( class_call(thermodynamics_at_z(pba, pth, 1./ppw->pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */ - ppw->inter_mode, + (interpolation_method)ppw->inter_mode, &(ppw->last_index_thermo), ppw->pvecback, ppw->pvecthermo), @@ -6338,7 +6224,7 @@ int perturbations_approximations( class_call(thermodynamics_at_z(pba, pth, 1./ppw->pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */ - ppw->inter_mode, + (interpolation_method)ppw->inter_mode, &(ppw->last_index_thermo), ppw->pvecback, ppw->pvecthermo), @@ -6431,7 +6317,7 @@ int perturbations_timescale( double * pvecthermo; /** - extract the fields of the parameter_and_workspace input structure */ - pppaw = parameters_and_workspace; + pppaw = (struct perturbations_parameters_and_workspace *)parameters_and_workspace; pba = pppaw->pba; pth = pppaw->pth; ppt = pppaw->ppt; @@ -6450,7 +6336,7 @@ int perturbations_timescale( /** - evaluate background quantities with background_at_tau() and Hubble time scale \f$ \tau_h = a/a' \f$ */ - class_call(background_at_tau(pba,tau, normal_info, ppw->inter_mode, &(ppw->last_index_back), pvecback), + class_call(background_at_tau(pba,tau, normal_info, (interpolation_method)ppw->inter_mode, &(ppw->last_index_back), pvecback), pba->error_message, error_message); @@ -6474,7 +6360,7 @@ int perturbations_timescale( class_call(thermodynamics_at_z(pba, pth, 1./pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */ - ppw->inter_mode, + (interpolation_method)ppw->inter_mode, &(ppw->last_index_thermo), pvecback, pvecthermo), @@ -6505,7 +6391,7 @@ int perturbations_timescale( class_call(thermodynamics_at_z(pba, pth, 1./pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */ - ppw->inter_mode, + (interpolation_method)ppw->inter_mode, &(ppw->last_index_thermo), pvecback, pvecthermo), @@ -6535,7 +6421,7 @@ int perturbations_timescale( class_call(thermodynamics_at_z(pba, pth, 1./pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */ - ppw->inter_mode, + (interpolation_method)ppw->inter_mode, &(ppw->last_index_thermo), pvecback, pvecthermo), @@ -7512,7 +7398,7 @@ int perturbations_sources( double dmu_idm_g = 0., ddmu_idm_g = 0., exp_mu_idm_g = 0.; /** - rename structure fields (just to avoid heavy notations) */ - pppaw = parameters_and_workspace; + pppaw = (struct perturbations_parameters_and_workspace *)parameters_and_workspace; ppr = pppaw->ppr; pba = pppaw->pba; pth = pppaw->pth; @@ -8189,7 +8075,7 @@ int perturbations_print_variables(double tau, /** - rename structure fields (just to avoid heavy notations) */ - pppaw = parameters_and_workspace; + pppaw = (struct perturbations_parameters_and_workspace *)parameters_and_workspace; k = pppaw->k; index_md = pppaw->index_md; @@ -8525,7 +8411,7 @@ int perturbations_print_variables(double tau, } else{ ppt->scalar_perturbations_data[ppw->index_ikout] = - realloc(ppt->scalar_perturbations_data[ppw->index_ikout], + (double*)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; @@ -8635,7 +8521,7 @@ int perturbations_print_variables(double tau, } else{ ppt->tensor_perturbations_data[ppw->index_ikout] = - realloc(ppt->tensor_perturbations_data[ppw->index_ikout], + (double*)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; @@ -8816,7 +8702,7 @@ int perturbations_derivs(double tau, /** - rename the fields of the input structure (just to avoid heavy notations) */ - pppaw = parameters_and_workspace; + pppaw = (struct perturbations_parameters_and_workspace *)parameters_and_workspace; k = pppaw->k; k2=k*k; @@ -10032,7 +9918,7 @@ int perturbations_tca_slip_and_shear(double * y, /** - rename the fields of the input structure (just to avoid heavy notations) */ - pppaw = parameters_and_workspace; + pppaw = (struct perturbations_parameters_and_workspace *)parameters_and_workspace; k = pppaw->k; k2=k*k; diff --git a/source/primordial.c b/source/primordial.c index 181ad29c..a5a2ec17 100644 --- a/source/primordial.c +++ b/source/primordial.c @@ -15,6 +15,7 @@ */ #include "primordial.h" +#include "parallel.h" /** * Primordial spectra for arbitrary argument and for all initial conditions. @@ -1598,73 +1599,20 @@ int primordial_inflation_spectra( ) { int index_k; - /* number of threads (always one if no openmp) */ - int number_of_threads=1; - /* index of the thread (always 0 if no openmp) */ - int thread=0; - - - /* This code can be optionally compiled with the openmp option for parallel computation. - Inside parallel regions, the use of the command "return" is forbidden. - For error management, instead of "return _FAILURE_", we will set the variable below - to "abort = _TRUE_". This will lead to a "return _FAILURE_" just after leaving the - parallel region. */ - int abort; - -#ifdef _OPENMP - /* instrumentation times */ - double tstart, tstop, tspent; -#endif + class_setup_parallel(); + /* loop over Fourier wavenumbers */ + for (index_k=0; index_k < ppm->lnk_size; index_k++) { -#ifdef _OPENMP + class_run_parallel(with_arguments(ppt,ppm,ppr,y_ini,index_k), -#pragma omp parallel - { - number_of_threads = omp_get_num_threads(); + class_call(primordial_inflation_one_wavenumber(ppt,ppm,ppr,y_ini,index_k), + ppm->error_message, + ppm->error_message); + return _SUCCESS_; + ); } -#endif - - abort = _FALSE_; - -#pragma omp parallel shared(ppt,ppm,ppr,abort,y_ini) private(index_k,thread,tspent,tstart,tstop) num_threads(number_of_threads) - - { - -#ifdef _OPENMP - thread = omp_get_thread_num(); - tspent=0.; -#endif - -#pragma omp for schedule (dynamic) - - /* loop over Fourier wavenumbers */ - for (index_k=0; index_k < ppm->lnk_size; index_k++) { - -#ifdef _OPENMP - tstart = omp_get_wtime(); -#endif - - class_call_parallel(primordial_inflation_one_wavenumber(ppt,ppm,ppr,y_ini,index_k), - ppm->error_message, - ppm->error_message); - -#ifdef _OPENMP - tstop = omp_get_wtime(); - - tspent += tstop-tstart; -#endif - - } - -#ifdef _OPENMP - if (ppm->primordial_verbose>1) - printf("In %s: time spent in parallel region (loop over k's) = %e s for thread %d\n", - __func__,tspent,thread); -#endif - - } /* end of parallel zone */ - if (abort == _TRUE_) return _FAILURE_; + class_finish_parallel(); ppm->is_non_zero[ppt->index_md_scalars][ppt->index_ic_ad] = _TRUE_; ppm->is_non_zero[ppt->index_md_tensors][ppt->index_ic_ten] = _TRUE_; @@ -3073,7 +3021,7 @@ int primordial_inflation_derivs( struct primordial_inflation_parameters_and_workspace * ppipaw; struct primordial * ppm; - ppipaw = parameters_and_workspace; + ppipaw = (struct primordial_inflation_parameters_and_workspace *)parameters_and_workspace; ppm = ppipaw->ppm; // a2 diff --git a/source/transfer.c b/source/transfer.c index 5a2191f6..60d5330c 100644 --- a/source/transfer.c +++ b/source/transfer.c @@ -28,6 +28,7 @@ */ #include "transfer.h" +#include "parallel.h" /** * Transfer function \f$ \Delta_l^{X} (q) \f$ at a given wavenumber q. @@ -151,8 +152,6 @@ int transfer_init( */ double *** sources_spline; - /* pointer on workspace (one per thread if openmp) */ - struct transfer_workspace * ptw; /** - array with the correspondence between the index of sources in the perturbation module and in the transfer module, @@ -165,19 +164,6 @@ int transfer_init( HyperInterpStruct BIS; double xmax; - /* This code can be optionally compiled with the openmp option for parallel computation. - Inside parallel regions, the use of the command "return" is forbidden. - For error management, instead of "return _FAILURE_", we will set the variable below - to "abort = _TRUE_". This will lead to a "return _FAILURE_" just after leaving the - parallel region. */ - int abort; - -#ifdef _OPENMP - - /* instrumentation times */ - double tstart, tstop, tspent; - -#endif /** - check whether any spectrum in harmonic space (i.e., any \f$C_l\f$'s) is actually requested */ @@ -307,130 +293,91 @@ int transfer_init( ptr->error_message); } - /* (a.3.) workspace, allocated in a parallel zone since in openmp - version there is one workspace per thread */ - - /* initialize error management flag */ - abort = _FALSE_; - - /* beginning of parallel region */ -#pragma omp parallel \ - shared(tau_size_max,ptr,ppr,pba,ppt,tp_of_tt,tau_rec,sources_spline,abort,BIS,tau0) \ - private(ptw,index_q,tstart,tstop,tspent) - { - -#ifdef _OPENMP - tspent = 0.; -#endif - - /* allocate workspace */ - - ptw = NULL; - - class_call_parallel(transfer_workspace_init(ptr, - ppr, - &ptw, - ppt->tau_size, - tau_size_max, - pba->K, - pba->sgnK, - tau0-pth->tau_cut, - &BIS), - ptr->error_message, - ptr->error_message); - - /** - loop over all wavenumbers (parallelized).*/ - /* For each wavenumber: */ + class_setup_parallel(); -#pragma omp for schedule (dynamic) - - for (index_q = 0; index_q < MAX(ptr->q_size,ptr->q_size_limber); index_q++) { - -#ifdef _OPENMP - tstart = omp_get_wtime(); -#endif + /** - loop over all wavenumbers (parallelized).*/ + /* For each wavenumber: */ + for (index_q = 0; index_q < MAX(ptr->q_size,ptr->q_size_limber); index_q++) { + class_run_parallel(with_arguments(pba,pth,ppt,ptr,ppr,index_q,tau_rec,tp_of_tt,sources,sources_spline,tau_size_max,window,tau0,&BIS), /* compute the transfer functions in the normal case (not the full Limber one) */ + struct transfer_workspace * ptw = NULL; + class_call(transfer_workspace_init(ptr, + ppr, + &ptw, + ppt->tau_size, + tau_size_max, + pba->K, + pba->sgnK, + tau0-pth->tau_cut, + &BIS), + ptr->error_message, + ptr->error_message); + if (index_q < ptr->q_size) { if (ptr->transfer_verbose > 2) printf("Compute transfer for wavenumber [%d/%zu]\n",index_q,ptr->q_size-1); /* Update interpolation structure: */ - class_call_parallel(transfer_update_HIS(ppr, - ptr, - ptw, - index_q, - tau0), - ptr->error_message, - ptr->error_message); + class_call(transfer_update_HIS(ppr, + ptr, + ptw, + index_q, + tau0), + ptr->error_message, + ptr->error_message); - class_call_parallel(transfer_compute_for_each_q(ppr, - pba, - ppt, - ptr, - tp_of_tt, - index_q, - tau_size_max, - tau_rec, - sources, - sources_spline, - window, - ptw, - _FALSE_), - ptr->error_message, - ptr->error_message); + class_call(transfer_compute_for_each_q(ppr, + pba, + ppt, + ptr, + tp_of_tt, + index_q, + tau_size_max, + tau_rec, + sources, + sources_spline, + window, + ptw, + _FALSE_), + ptr->error_message, + ptr->error_message); } /* compute the transfer functions in the full Limber case (if - this case is not needed, ptr->q_size_limber=0 and the - condition is never met) */ + this case is not needed, ptr->q_size_limber=0 and the + condition is never met) */ if (index_q < ptr->q_size_limber) { - class_call_parallel(transfer_compute_for_each_q(ppr, - pba, - ppt, - ptr, - tp_of_tt, - index_q, - tau_size_max, - tau_rec, - sources, - sources_spline, - window, - ptw, - _TRUE_), + class_call(transfer_compute_for_each_q(ppr, + pba, + ppt, + ptr, + tp_of_tt, + index_q, + tau_size_max, + tau_rec, + sources, + sources_spline, + window, + ptw, + _TRUE_), ptr->error_message, ptr->error_message); } -#ifdef _OPENMP - tstop = omp_get_wtime(); - - tspent += tstop-tstart; -#endif - -#pragma omp flush(abort) - - } /* end of loop over wavenumber */ - - /* free workspace allocated inside parallel zone */ - class_call_parallel(transfer_workspace_free(ptr,ptw), - ptr->error_message, - ptr->error_message); - -#ifdef _OPENMP - if (ptr->transfer_verbose>1) - printf("In %s: time spent in parallel region (loop over k's) = %e s for thread %d\n", - __func__,tspent,omp_get_thread_num()); -#endif - - } /* end of parallel region */ + class_call(transfer_workspace_free(ptr,ptw), + ptr->error_message, + ptr->error_message); + return _SUCCESS_; + ); + } /* end of loop over wavenumber */ - if (abort == _TRUE_) return _FAILURE_; + class_finish_parallel(); /** - finally, free arrays allocated outside parallel zone */ free(window); @@ -4056,12 +4003,12 @@ int transfer_radial_function( double l = (double)ptr->l[index_l]; double rescale_argument; double rescale_amplitude; - double * rescale_function; - int (*interpolate_Phi)(); - int (*interpolate_dPhi)(); - int (*interpolate_Phid2Phi)(); - int (*interpolate_PhidPhi)(); - int (*interpolate_PhidPhid2Phi)(); + double* rescale_function; + int (*interpolate_Phi)(HyperInterpStruct*, int, int, double*, double*, char*); + int (*interpolate_dPhi)(HyperInterpStruct*, int, int, double*, double*, char*); + int (*interpolate_Phid2Phi)(HyperInterpStruct*, int, int, double*, double*, double*, char*); + int (*interpolate_PhidPhi)(HyperInterpStruct*, int, int, double*, double*, double*, char*); + int (*interpolate_PhidPhid2Phi)(HyperInterpStruct*, int, int, double*, double*, double*, double*, char*); enum Hermite_Interpolation_Order HIorder; K = ptw->K; @@ -4568,7 +4515,7 @@ int transfer_workspace_free( free(ptw->cscKgen); free(ptw->cotKgen); - free(ptw); + //free(ptw); return _SUCCESS_; } diff --git a/tools/arrays.c b/tools/arrays.c index 89837f31..13c82bfb 100644 --- a/tools/arrays.c +++ b/tools/arrays.c @@ -4,6 +4,7 @@ */ #include "arrays.h" +#include "parallel.h" /** * Called by thermodynamics_init(); perturbations_sources(). @@ -355,7 +356,7 @@ int array_spline( return _FAILURE_; } - u = malloc((n_lines-1) * sizeof(double)); + u = (double*)malloc((n_lines-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -469,7 +470,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 = (double*)malloc((n_lines-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -571,10 +572,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 = (double*)malloc((x_size-1) * y_size * sizeof(double)); + p = (double*)malloc(y_size * sizeof(double)); + qn = (double*)malloc(y_size * sizeof(double)); + un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); @@ -733,10 +734,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 = (double*)malloc((x_size-1) * y_size * sizeof(double)); + p = (double*)malloc(y_size * sizeof(double)); + qn = (double*)malloc(y_size * sizeof(double)); + un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -895,10 +896,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 = (double*)malloc((x_size-1) * y_size * sizeof(double)); + p = (double*)malloc(y_size * sizeof(double)); + qn = (double*)malloc(y_size * sizeof(double)); + un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1060,16 +1061,13 @@ int array_spline_table_columns2( double * qn; double * un; double * u; - double sig; - int index_x; + int index_y; - 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 = (double*)malloc((x_size-1) * y_size * sizeof(double)); + p = (double*)malloc(y_size * sizeof(double)); + qn = (double*)malloc(y_size * sizeof(double)); + un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1089,15 +1087,16 @@ int array_spline_table_columns2( if (x_size==2) spline_mode = _SPLINE_NATURAL_; // in the case of only 2 x-values, only the natural spline method is appropriate, for _SPLINE_EST_DERIV_ 3 x-values are needed. -#pragma omp parallel \ - shared(x,x_size,y_array,y_size,ddy_array,spline_mode,p,qn,un,u) \ - private(index_y,index_x,sig,dy_first,dy_last) - { + class_setup_parallel(); -#pragma omp for schedule (dynamic) + for (index_y=0; index_y < y_size; index_y++) { - for (index_y=0; index_y < y_size; index_y++) { + class_run_parallel(=, + double dy_first; + double dy_last; + int index_x; + double sig; if (spline_mode == _SPLINE_NATURAL_) { ddy_array[index_y*x_size+0] = 0.0; u[0*y_size+index_y] = 0.0; @@ -1174,8 +1173,12 @@ int array_spline_table_columns2( ddy_array[index_y*x_size+(index_x+1)] + u[index_x*y_size+index_y]; } - } + return _SUCCESS_; + ); } + + class_finish_parallel(); + free(qn); free(p); free(u); @@ -1205,7 +1208,7 @@ int array_spline_table_one_column( double dy_first; double dy_last; - u = malloc((x_size-1) * sizeof(double)); + u = (double*)malloc((x_size-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -1340,7 +1343,7 @@ int array_logspline_table_one_column( double dy_first; double dy_last; - u = malloc((x_stop-1) * sizeof(double)); + u = (double*)malloc((x_stop-1) * sizeof(double)); if (u == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; @@ -3131,7 +3134,7 @@ int array_smooth_trg(double * array, double weigth; double *coeff; - smooth=malloc(k_size*sizeof(double)); + smooth=(double*)malloc(k_size*sizeof(double)); if (smooth == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); return _FAILURE_; @@ -3261,7 +3264,7 @@ int array_smooth(double * array, int i,j,jmin,jmax; double weigth; - smooth=malloc(n_lines*sizeof(double)); + smooth=(double*)malloc(n_lines*sizeof(double)); if (smooth == NULL) { sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); return _FAILURE_; diff --git a/tools/hyperspherical.c b/tools/hyperspherical.c index 7a37bcd6..25e8fe01 100644 --- a/tools/hyperspherical.c +++ b/tools/hyperspherical.c @@ -6,6 +6,7 @@ */ #include "hyperspherical.h" +#include "parallel.h" int hyperspherical_HIS_create(int K, double beta, @@ -28,8 +29,6 @@ int hyperspherical_HIS_create(int K, double deltax, beta2, lambda, x, xfwd; double *sqrtK, *one_over_sqrtK,*PhiL; int j, k, l, nx, lmax, l_recurrence_max; - int abort; - int current_chunk, index_x; beta2 = beta*beta; lmax = lvec[nl-1]; @@ -132,107 +131,115 @@ int hyperspherical_HIS_create(int K, int xfwdidx = (xfwd-xmin)/deltax; //Calculate and assign Phi and dPhi values: - abort = _FALSE_; - -#pragma omp parallel \ - shared(nx,pHIS,xfwd,K,l_recurrence_max,beta,sqrtK,one_over_sqrtK,lvec,nl,xfwdidx,abort,error_message) \ - private(j,PhiL,k,l,current_chunk,index_x) \ - firstprivate(lmax) - { - class_alloc_parallel(PhiL,(lmax+2)*sizeof(double)*_HYPER_CHUNK_,error_message); - if ((K == 1) && ((int)(beta+0.2) == (lmax+1))) { - /** Take care of special case lmax = beta-1. - The routine below will try to compute - Phi_{lmax+1} which is not allowed. However, - the purpose is to calculate the derivative - Phi'_{lmax}, and the formula is correct if we set Phi_{lmax+1} = 0. - */ - for(j=0; j<(lmax+2)*_HYPER_CHUNK_; j++){ - PhiL[j] = 0.0; - } + int TNUM_i; + int TNUM_TOTAL; + class_setup_parallel(); + + /* Increase beyond threads of system to increase parallelization + * This will require more memory space to allocate the PhiL, + * but memory space is cheaper than computation time, usually. + * Default increase factor: 6 */ + TNUM_TOTAL = task_system.get_num_threads() * 6; + for(TNUM_i=0;TNUM_ix[j], - pHIS->sinK[j], - pHIS->cotK[j], - sqrtK, - one_over_sqrtK, - PhiL); - //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: - for (k=0; k<=index_recurrence_max; k++){ - l = lvec[k]; - pHIS->phi[k*nx+j] = PhiL[l]; - pHIS->dphi[k*nx+j] = l*pHIS->cotK[j]*PhiL[l]-sqrtK[l+1]*PhiL[l+1]; + for (j=0; jx[j], + pHIS->sinK[j], + pHIS->cotK[j], + sqrtK, + one_over_sqrtK, + PhiL); + //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: + for (k=0; k<=index_recurrence_max; k++){ + l = lvec[k]; + pHIS->phi[k*nx+j] = PhiL[l]; + pHIS->dphi[k*nx+j] = l*pHIS->cotK[j]*PhiL[l]-sqrtK[l+1]*PhiL[l+1]; + } } - } - /** - for (j=0; jx+j, - pHIS->sinK+j, - pHIS->cotK+j, - current_chunk, - sqrtK, - one_over_sqrtK, - PhiL); - //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: - for (k=0; k<=index_recurrence_max; k++){ - l = lvec[k]; - for (index_x=0; index_xphi[k*nx+j+index_x] = PhiL[l*current_chunk+index_x]; - pHIS->dphi[k*nx+j+index_x] = l*pHIS->cotK[j+index_x]* - PhiL[l*current_chunk+index_x]- - sqrtK[l+1]*PhiL[(l+1)*current_chunk+index_x]; + /** + for (j=0; jx+j, + pHIS->sinK+j, + pHIS->cotK+j, + current_chunk, + sqrtK, + one_over_sqrtK, + PhiL); + //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: + for (k=0; k<=index_recurrence_max; k++){ + l = lvec[k]; + for (index_x=0; index_xphi[k*nx+j+index_x] = PhiL[l*current_chunk+index_x]; + pHIS->dphi[k*nx+j+index_x] = l*pHIS->cotK[j+index_x]* + PhiL[l*current_chunk+index_x]- + sqrtK[l+1]*PhiL[(l+1)*current_chunk+index_x]; + } } } - } - */ - -#pragma omp for schedule (dynamic) \ - - for (j=xfwdidx; jx+j, - pHIS->sinK+j, - pHIS->cotK+j, - current_chunk, - sqrtK, - one_over_sqrtK, - PhiL); - - //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: - for (k=0; k<=index_recurrence_max; k++){ - l = lvec[k]; - for (index_x=0; index_xphi[k*nx+j+index_x] = PhiL[l*current_chunk+index_x]; - pHIS->dphi[k*nx+j+index_x] = l*pHIS->cotK[j+index_x]* - PhiL[l*current_chunk+index_x]- - sqrtK[l+1]*PhiL[(l+1)*current_chunk+index_x]; + */ + + for (j=xfwdidx; jx+j, + pHIS->sinK+j, + pHIS->cotK+j, + current_chunk, + sqrtK, + one_over_sqrtK, + PhiL); + + //We have now populated PhiL at x, assign Phi and dPhi for all l in lvec: + for (k=0; k<=index_recurrence_max; k++){ + l = lvec[k]; + for (index_x=0; index_xphi[k*nx+j+index_x] = PhiL[l*current_chunk+index_x]; + pHIS->dphi[k*nx+j+index_x] = l*pHIS->cotK[j+index_x]* + PhiL[l*current_chunk+index_x]- + sqrtK[l+1]*PhiL[(l+1)*current_chunk+index_x]; + } } } - } - free(PhiL); + free(PhiL); + return _SUCCESS_; + ); } - if (abort == _TRUE_) return _FAILURE_; + class_finish_parallel(); free(sqrtK); free(one_over_sqrtK); @@ -1195,7 +1202,7 @@ int hyperspherical_get_xmin_from_Airy(int K, double PhiWKB_minus_phiminabs(double x, void *param){ double phiwkb; - struct WKB_parameters *wkbparam = param; + struct WKB_parameters *wkbparam = (struct WKB_parameters *)param; hyperspherical_WKB(wkbparam->K,wkbparam->l,wkbparam->beta,x, &phiwkb); return(fabs(phiwkb)-wkbparam->phiminabs); } From adafaec44a8c137f935fe905c9eb8778d6e231d3 Mon Sep 17 00:00:00 2001 From: Nils <32747859+schoeneberg@users.noreply.github.com> Date: Thu, 8 Feb 2024 17:11:19 +0100 Subject: [PATCH 2/3] Fix Makefile pypi compatibility (#145) * Fixing 142 by using cubic interpolation of lnPk * Making makefile compatible with the pypi installation by promoting the CLASSDIR to an environment variable that can be passed during installation --------- Co-authored-by: schoeneberg --- Makefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 10f3126e..3fc4b8b5 100644 --- a/Makefile +++ b/Makefile @@ -59,7 +59,8 @@ HEATING = external/heating ######################################################## # pass current working directory to the code -CCFLAG += -D__CLASSDIR__='"$(MDIR)"' +CLASSDIR ?= $(MDIR) +CCFLAG += -D__CLASSDIR__='"$(CLASSDIR)"' # where to find include files *.h INCLUDES = -I../include From dcb8be02e8fc66e98331b540818fddb46e23aef4 Mon Sep 17 00:00:00 2001 From: Nils <32747859+schoeneberg@users.noreply.github.com> Date: Tue, 20 Feb 2024 13:41:22 +0100 Subject: [PATCH 3/3] fix compiler warnings (#150) * Fixed class_realloc issues, and simplified the macro * Removed unused parameters * Removed superflous brackets * Replaced sprintf (unsafe) with snprintf, making sure to calculate the size correctly of each string, and making it robust for uses in single-line if else satements * Removed unused variable * Removed unneccessary define that remained from previous ideas about snprintf * Fix evolver prototype * abs vs fabs fixes * Removed remaining compiler warnings. * changed version number to 3.2.3 --------- Co-authored-by: schoeneberg Co-authored-by: Thomas Tram Co-authored-by: lesgourg --- include/common.h | 27 ++++- include/parallel.h | 2 +- source/background.c | 19 ++-- source/distortions.c | 30 +++--- source/input.c | 48 ++++----- source/lensing.c | 71 +++++-------- source/output.c | 118 ++++++++++----------- source/perturbations.c | 33 +++--- source/primordial.c | 35 ++----- source/thermodynamics.c | 14 +-- source/transfer.c | 10 -- tools/arrays.c | 220 ++++++++++++++++++++-------------------- tools/hyperspherical.c | 5 +- tools/parser.c | 6 +- tools/quadrature.c | 6 +- 15 files changed, 305 insertions(+), 339 deletions(-) diff --git a/include/common.h b/include/common.h index efda051b..7dc96e57 100644 --- a/include/common.h +++ b/include/common.h @@ -15,7 +15,7 @@ #ifndef __COMMON__ #define __COMMON__ -#define _VERSION_ "v3.2.2" +#define _VERSION_ "v3.2.3" /* @cond INCLUDE_WITH_DOXYGEN */ @@ -100,6 +100,16 @@ extern "C" { /* general CLASS macros */ +//This macro receives additional 'do {' and '} while(0)' to safeguard +//in single-line if else clauses without '{' and '}' +//Also, careful: Since sprintf(NULL,0,x) returns the size of characters +//that are inside of the string x, then the buffer needs to be +//actually one character longer to hold also the null character '\0' +#define class_sprintf(string, format...) do { \ + int _buffer_size_sprintf = snprintf(NULL, 0, format); \ + snprintf(string, _buffer_size_sprintf+1, format); \ +} while (0) + #define class_build_error_string(dest,tmpl,...) { \ ErrorMsg FMsg; \ class_protect_sprintf(FMsg,tmpl,__VA_ARGS__); \ @@ -188,8 +198,8 @@ extern "C" { } /* macro for re-allocating memory, returning error if it failed */ -#define class_realloc(pointer, newname, size, error_message_output) { \ - pointer=(__typeof__(pointer))realloc(newname,size); \ +#define class_realloc(pointer, size, error_message_output) { \ + pointer=(__typeof__(pointer))realloc(pointer,size); \ if (pointer == NULL) { \ int size_int; \ size_int = size; \ @@ -335,6 +345,17 @@ extern "C" { #define class_print_species(name,type) \ printf("-> %-30s Omega = %-15g , omega = %-15g\n",name,pba->Omega0_##type,pba->Omega0_##type*pba->h*pba->h); +//Generic evolver prototype +#define EVOLVER_PROTOTYPE \ + int (*)(double, double *, double *, void *, ErrorMsg), \ + double, double, double *, int *, \ + int, void *, double, double, \ + int (*)(double, void *, double *, ErrorMsg), \ + double, double *, int, \ + int (*)(double, double *, double *, int, void *, ErrorMsg), \ + int (*)(double, double *, double *, void *, ErrorMsg), \ + ErrorMsg + /* Forward-Declare the structs of CLASS */ struct background; struct thermodynamics; diff --git a/include/parallel.h b/include/parallel.h index da18ee78..28fabf90 100644 --- a/include/parallel.h +++ b/include/parallel.h @@ -215,8 +215,8 @@ class TaskSystem { const unsigned int count_; std::vector threads_; - std::vector queues_; std::atomic index_; + std::vector queues_; }; } diff --git a/source/background.c b/source/background.c index 22246e23..65ea692b 100644 --- a/source/background.c +++ b/source/background.c @@ -1442,9 +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)); - + class_realloc(pba->q_ncdm[k],pba->q_size_ncdm[k]*sizeof(double), pba->error_message); + class_realloc(pba->w_ncdm[k],pba->q_size_ncdm[k]*sizeof(double), pba->error_message); if (pba->background_verbose > 0) { printf("ncdm species i=%d sampled with %d points for purpose of perturbation integration\n", @@ -1470,8 +1469,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)); + class_realloc(pba->q_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double), pba->error_message); + class_realloc(pba->w_ncdm_bg[k],pba->q_size_ncdm_bg[k]*sizeof(double), pba->error_message); /** - in verbose mode, inform user of number of sampled momenta for background quantities */ @@ -1886,9 +1885,9 @@ int background_solve( double conformal_distance; /* evolvers */ - extern int evolver_rk(); - extern int evolver_ndf15(); - int (*generic_evolver)() = evolver_ndf15; + extern int evolver_rk(EVOLVER_PROTOTYPE); + extern int evolver_ndf15(EVOLVER_PROTOTYPE); + int (*generic_evolver)(EVOLVER_PROTOTYPE) = evolver_ndf15; /* initial and final loga values */ double loga_ini, loga_final; @@ -2440,9 +2439,9 @@ int background_output_titles( class_store_columntitle(titles,"(.)rho_idm",pba->has_idm); if (pba->has_ncdm == _TRUE_) { for (n=0; nN_ncdm; n++) { - sprintf(tmp,"(.)rho_ncdm[%d]",n); + class_sprintf(tmp,"(.)rho_ncdm[%d]",n); class_store_columntitle(titles,tmp,_TRUE_); - sprintf(tmp,"(.)p_ncdm[%d]",n); + class_sprintf(tmp,"(.)p_ncdm[%d]",n); class_store_columntitle(titles,tmp,_TRUE_); } } diff --git a/source/distortions.c b/source/distortions.c index 507629d6..b198217e 100644 --- a/source/distortions.c +++ b/source/distortions.c @@ -161,8 +161,8 @@ int distortions_constants(struct precision * ppr, pow(pba->Omega0_b*pow(pba->h,2.)/0.02225,-2./5.)* pow(pba->T_cmb/2.726,1./5.); - sprintf(psd->sd_PCA_file_generator,"%s/%s",ppr->sd_external_path,"generate_PCA_files.py"); - sprintf(psd->sd_detector_list_file,"%s/%s",ppr->sd_external_path,"detectors_list.dat"); + class_sprintf(psd->sd_PCA_file_generator,"%s/%s",ppr->sd_external_path,"generate_PCA_files.py"); + class_sprintf(psd->sd_detector_list_file,"%s/%s",ppr->sd_external_path,"detectors_list.dat"); return _SUCCESS_; } @@ -223,12 +223,12 @@ int distortions_set_detector(struct precision * ppr, else { /* Generate a custom name for this custom detector, so we can check if it has already been defined */ if (psd->has_detector_file == _TRUE_) { - sprintf(psd->sd_detector_name, + class_sprintf(psd->sd_detector_name, "Custom__%.80s__Detector", psd->sd_detector_file_name); } else { - sprintf(psd->sd_detector_name, + class_sprintf(psd->sd_detector_name, "Custom__%7.2e_%7.2e_%7.2e_%i_%7.2e__Detector", psd->sd_detector_nu_min,psd->sd_detector_nu_max,psd->sd_detector_nu_delta,psd->sd_detector_bin_number,psd->sd_detector_delta_Ic); } @@ -282,7 +282,7 @@ int distortions_set_detector(struct precision * ppr, psd->error_message, "Detector property type disagrees between stored detector '%s' and input -> Userdefined (input) vs Noisefile (stored)", detector_name); - sprintf(psd->sd_detector_file_name, "%s", detector_noise_file_name); + class_sprintf(psd->sd_detector_file_name, "%s", detector_noise_file_name); psd->has_detector_file = has_detector_noise_file; } else { @@ -304,7 +304,7 @@ int distortions_set_detector(struct precision * ppr, psd->error_message, "Delta frequency (sd_detector_nu_delta) disagrees between stored detector '%s' and input -> %.10e (input) vs %.10e (stored)", detector_name,psd->sd_detector_nu_delta,nu_delta); - class_test(fabs(psd->sd_detector_bin_number-N_bins)>ppr->tol_sd_detector, + class_test(abs(psd->sd_detector_bin_number-N_bins)>ppr->tol_sd_detector, psd->error_message, "Number of bins (sd_detector_bin_number) disagrees between stored detector '%s' and input -> %i (input) vs %i (stored)", detector_name,psd->sd_detector_bin_number,N_bins); @@ -391,7 +391,7 @@ int distortions_generate_detector(struct precision * ppr, } if (psd->has_detector_file == _TRUE_) { - sprintf(temporary_string,"python %s %s %s %s %.10e %.10e %i %i %.10e %.10e %.10e", + class_sprintf(temporary_string,"python %s %s %s %s %.10e %.10e %i %i %.10e %.10e %.10e", psd->sd_PCA_file_generator, psd->sd_detector_name, ppr->sd_external_path, @@ -406,7 +406,7 @@ int distortions_generate_detector(struct precision * ppr, } else { - sprintf(temporary_string,"python %s %s %.10e %.10e %.10e %i %.10e %.10e %i %.10e %i %.10e %.10e %.10e", + class_sprintf(temporary_string,"python %s %s %.10e %.10e %.10e %i %.10e %.10e %i %.10e %i %.10e %.10e %.10e", psd->sd_PCA_file_generator, psd->sd_detector_name, psd->sd_detector_nu_min, @@ -452,7 +452,7 @@ int distortions_read_detector_noisefile(struct precision * ppr, int numcols; /** Open file */ - sprintf(psd->sd_detector_noise_file,"%s/%s",ppr->sd_external_path,psd->sd_detector_file_name); + class_sprintf(psd->sd_detector_noise_file,"%s/%s",ppr->sd_external_path,psd->sd_detector_file_name); class_open(infile, psd->sd_detector_noise_file, "r", psd->error_message); /** Read header */ @@ -1418,7 +1418,7 @@ int distortions_read_br_data(struct precision * ppr, int headlines = 0; /** Open file */ - sprintf(br_file,"%s/%s_branching_ratios.dat", ppr->sd_external_path, psd->sd_detector_name); + class_sprintf(br_file,"%s/%s_branching_ratios.dat", ppr->sd_external_path, psd->sd_detector_name); class_open(infile, br_file, "r", psd->error_message); /** Read header */ @@ -1651,7 +1651,7 @@ int distortions_read_sd_data(struct precision * ppr, int index_x,index_k; /** Open file */ - sprintf(sd_file,"%s/%s_distortions_shapes.dat",ppr->sd_external_path, psd->sd_detector_name); + class_sprintf(sd_file,"%s/%s_distortions_shapes.dat",ppr->sd_external_path, psd->sd_detector_name); class_open(infile, sd_file, "r", psd->error_message); /** Read header */ @@ -1931,16 +1931,16 @@ int distortions_output_sd_titles(struct distortions * psd, class_store_columntitle(titles,"SD_tot",_TRUE_); for (index_type=0;index_typetype_size;++index_type){ if (index_type==psd->index_type_g) { - sprintf(temp_title,"SD[g]"); + class_sprintf(temp_title,"SD[g]"); } if (index_type==psd->index_type_y) { - sprintf(temp_title,"SD[y]"); + class_sprintf(temp_title,"SD[y]"); } if (index_type==psd->index_type_mu) { - sprintf(temp_title,"SD[mu]"); + class_sprintf(temp_title,"SD[mu]"); } if (index_type>=psd->index_type_PCA) { - sprintf(temp_title,"SD[e_%i]",(index_type-psd->index_type_PCA)); + class_sprintf(temp_title,"SD[e_%i]",(index_type-psd->index_type_PCA)); } class_store_columntitle(titles,temp_title,_TRUE_); } diff --git a/source/input.c b/source/input.c index d4ca6fe2..fcdcb52d 100644 --- a/source/input.c +++ b/source/input.c @@ -277,7 +277,7 @@ int input_set_root(char* input_file, /* No file has been found yet */ found_filenum = _FALSE_; for (iextens = 0; iextens < n_extensions; ++iextens){ - sprintf(tmp_file,"%s%02d_%s", outfname, filenum, output_extensions[iextens]); + class_sprintf(tmp_file,"%s%02d_%s", outfname, filenum, output_extensions[iextens]); if (file_exists(tmp_file) == _TRUE_){ /* Found a file, the outer loop is forced to keep searching */ found_filenum = _TRUE_; @@ -295,8 +295,8 @@ int input_set_root(char* input_file, pfc->filename, errmsg), errmsg,errmsg); - sprintf(fc_root.name[0],"root"); - sprintf(fc_root.value[0],"%s%02d_",outfname,filenum); + class_sprintf(fc_root.name[0],"root"); + class_sprintf(fc_root.value[0],"%s%02d_",outfname,filenum); fc_root.read[0] = _FALSE_; class_call(parser_cat(pfc, &fc_root, @@ -314,7 +314,7 @@ int input_set_root(char* input_file, } /* If root was found, set the index in the fc_input struct */ else{ - sprintf(pfc->value[index_root_in_fc_input],"%s%02d_",outfname,filenum); + class_sprintf(pfc->value[index_root_in_fc_input],"%s%02d_",outfname,filenum); (*ppfc_input) = pfc; } } @@ -328,8 +328,8 @@ int input_set_root(char* input_file, pfc->filename, errmsg), errmsg,errmsg); - sprintf(fc_root.name[0],"root"); - sprintf(fc_root.value[0],"%s_",outfname); + class_sprintf(fc_root.name[0],"root"); + class_sprintf(fc_root.value[0],"%s_",outfname); fc_root.read[0] = _FALSE_; class_call(parser_cat(pfc, &fc_root, @@ -347,7 +347,7 @@ int input_set_root(char* input_file, } /* If root was found, set the index in the fc_input struct */ else{ - sprintf(pfc->value[index_root_in_fc_input],"%s_",outfname); + class_sprintf(pfc->value[index_root_in_fc_input],"%s_",outfname); (*ppfc_input) = pfc; } } @@ -658,7 +658,7 @@ int input_shooting(struct file_content * pfc, /* Store xzero */ // This needs to be done with enough accuracy. A standard double has a relative // precision of around 1e-16, so 1e-20 should be good enough for the shooting - sprintf(fzw.fc.value[fzw.unknown_parameters_index[0]],"%.20e",xzero); + class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[0]],"%.20e",xzero); if (input_verbose > 0) { fprintf(stdout," -> found '%s = %s'\n", fzw.fc.name[fzw.unknown_parameters_index[0]], @@ -705,7 +705,7 @@ int input_shooting(struct file_content * pfc, // This needs to be done with enough accuracy. A standard double has a relative // precision of around 1e-16, so 1e-20 should be good enough for the shooting for (counter = 0; counter < unknown_parameters_size; counter++){ - sprintf(fzw.fc.value[fzw.unknown_parameters_index[counter]], + class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[counter]], "%.20e",x_inout[counter]); if (input_verbose > 0) { fprintf(stdout," -> found '%s = %s'\n", @@ -838,7 +838,7 @@ int input_shooting(struct file_content * pfc, A_s = (fzw.target_value[0]/sigma8_or_S8) *(fzw.target_value[0]/sigma8_or_S8) * A_s; //(truesigma/sigma_for_guess)^2 *A_s_for_guess /* Store the derived value with high enough accuracy */ - sprintf(fzw.fc.value[pfc->size],"%.20e",A_s); + class_sprintf(fzw.fc.value[pfc->size],"%.20e",A_s); if (input_verbose > 0) { fprintf(stdout," -> found '%s = %s'\n", fzw.fc.name[pfc->size], @@ -1345,7 +1345,7 @@ int input_try_unknown_parameters(double * unknown_parameter, // This needs to be done with enough accuracy. A standard double has a relative // precision of around 1e-16, so 1e-20 should be good enough for the shooting for (i=0; i < unknown_parameters_size; i++) { - sprintf(pfzw->fc.value[pfzw->unknown_parameters_index[i]],"%.20e",unknown_parameter[i]); + class_sprintf(pfzw->fc.value[pfzw->unknown_parameters_index[i]],"%.20e",unknown_parameter[i]); } class_call(input_read_precisions(&(pfzw->fc),&pr,&ba,&th,&pt,&tr,&pm,&hr,&fo,&le,&sd,&op, @@ -3001,7 +3001,7 @@ int input_read_parameters_species(struct file_content * pfc, if (ppt->perturbations_verbose > 0) { printf("WARNING: only %i entries of alpha_idm_dr were provided for %i moments, filling up the rest with the last entry provided\n", entries_read, ppr->l_max_idr-1); } - class_realloc(ppt->alpha_idm_dr,ppt->alpha_idm_dr,(ppr->l_max_idr-1)*sizeof(double),errmsg); + class_realloc(ppt->alpha_idm_dr,(ppr->l_max_idr-1)*sizeof(double),errmsg); for (n=entries_read; n<(ppr->l_max_idr-1); n++) ppt->alpha_idm_dr[n] = ppt->alpha_idm_dr[entries_read-1]; } } @@ -3044,7 +3044,7 @@ int input_read_parameters_species(struct file_content * pfc, if (ppt->perturbations_verbose > 0) { printf("WARNING: only %i entries of beta_idr were provided for %i moments, filling up the rest with the last entry provided\n", entries_read, ppr->l_max_idr-1); } - class_realloc(ppt->beta_idr,ppt->beta_idr,(ppr->l_max_idr-1)*sizeof(double),errmsg); + class_realloc(ppt->beta_idr,(ppr->l_max_idr-1)*sizeof(double),errmsg); for (n=entries_read; n<(ppr->l_max_idr-1); n++) ppt->beta_idr[n] = ppt->beta_idr[entries_read-1]; } @@ -3133,7 +3133,7 @@ int input_read_parameters_species(struct file_content * pfc, /** 7.3) Final consistency checks for dark matter species */ - class_test(abs(f_cdm + f_idm - 1.) > 1e-10, + class_test(fabs(f_cdm + f_idm - 1.) > 1e-10, errmsg, "The dark matter species do not add up to the expected value"); @@ -3146,7 +3146,7 @@ int input_read_parameters_species(struct file_content * pfc, class_test((f_idm > 0.) && (pba->Omega0_cdm == 0.), errmsg, "If you want a fraction of interacting, to be consistent, you should not set the fraction of CDM to zero"); - class_test(abs(f_cdm + f_idm - 1.) > ppr->tol_fraction_accuracy, + class_test(fabs(f_cdm + f_idm - 1.) > ppr->tol_fraction_accuracy, errmsg, "The dark matter species do not add up to the expected value"); if ( f_idm > 0. ) @@ -4814,7 +4814,7 @@ int input_read_parameters_spectra(struct file_content * pfc, errmsg, errmsg); /* Complete set of parameters */ - if ((flag1 == _TRUE_)) { + if (flag1 == _TRUE_) { if ((strstr(string1,"analytic") != NULL)){ ptr->has_nz_analytic = _TRUE_; } @@ -4829,7 +4829,7 @@ int input_read_parameters_spectra(struct file_content * pfc, errmsg, errmsg); /* Complete set of parameters */ - if ((flag1 == _TRUE_)) { + if (flag1 == _TRUE_) { if ((strstr(string1,"analytic") != NULL)){ ptr->has_nz_evo_analytic = _TRUE_; } @@ -5535,7 +5535,7 @@ int input_write_info(struct file_content * pfc, /* Finally, since all variables are read, we can also print the parameters.ini and unused_parameters files */ if (flag1 == _TRUE_) { - sprintf(param_output_name,"%s%s",pop->root,"parameters.ini"); + class_sprintf(param_output_name,"%s%s",pop->root,"parameters.ini"); class_open(param_output,param_output_name,"w",errmsg); fprintf(param_output,"# List of input/precision parameters actually read\n"); fprintf(param_output,"# (all other parameters set to default values)\n"); @@ -5544,7 +5544,7 @@ int input_write_info(struct file_content * pfc, fprintf(param_output,"# This file can be used as the input file of another run\n"); fprintf(param_output,"#\n"); - sprintf(param_unused_name,"%s%s",pop->root,"unused_parameters"); + class_sprintf(param_unused_name,"%s%s",pop->root,"unused_parameters"); class_open(param_unused,param_unused_name,"w",errmsg); fprintf(param_unused,"# List of input/precision parameters passed\n"); fprintf(param_unused,"# but not used (just for info)\n"); @@ -5890,13 +5890,13 @@ int input_default_params(struct background *pba, /** 5) Injection efficiency */ pin->f_eff_type = f_eff_on_the_spot; pin->f_eff = 1.; - sprintf(pin->f_eff_file,"external/heating/example_f_eff_file.dat"); + class_sprintf(pin->f_eff_file,"external/heating/example_f_eff_file.dat"); /** 6) Deposition function */ pin->chi_type = chi_CK; /** 6.1) External file */ - sprintf(pin->chi_z_file,"external/heating/example_chiz_file.dat"); - sprintf(pin->chi_x_file,"external/heating/example_chix_file.dat"); + class_sprintf(pin->chi_z_file,"external/heating/example_chiz_file.dat"); + class_sprintf(pin->chi_x_file,"external/heating/example_chix_file.dat"); /** * Default to input_read_parameters_nonlinear @@ -6076,7 +6076,7 @@ int input_default_params(struct background *pba, /** 1.a.3) Detector name */ psd->has_user_defined_name = _FALSE_; psd->has_user_defined_detector = _FALSE_; - sprintf(psd->sd_detector_name,"PIXIE"); + class_sprintf(psd->sd_detector_name,"PIXIE"); /** 1.3.a.1) Detector nu min */ psd->sd_detector_nu_min = 30.; /** 1.3.a.2) Detector nu max */ @@ -6114,7 +6114,7 @@ int input_default_params(struct background *pba, /** 1) Output for external files */ /** 1.a) File name */ - sprintf(pop->root,"output/"); + class_sprintf(pop->root,"output/"); /** 1.b) Headers */ pop->write_header = _TRUE_; /** 1.c) Format */ diff --git a/source/lensing.c b/source/lensing.c index c87ca9bc..778af1f7 100644 --- a/source/lensing.c +++ b/source/lensing.c @@ -121,16 +121,6 @@ int lensing_init( double * ksip = NULL; /* ksip[index_mu] */ double * ksim = NULL; /* ksim[index_mu] */ - double fac,fac1; - double X_000; - double X_p000; - double X_220; - double X_022; - double X_p022; - double X_121; - double X_132; - double X_242; - int num_mu,index_mu,icount; int l; double ll; @@ -141,9 +131,6 @@ int lensing_init( double * cl_bb = NULL; /* unlensed cl, to be filled to avoid repeated calls to harmonic_cl_at_l */ double * cl_pp; /* potential cl, to be filled to avoid repeated calls to harmonic_cl_at_l */ - double res,resX,lens; - double resp, resm, lensp, lensm; - double * sqrt1; double * sqrt2; double * sqrt3; @@ -1286,14 +1273,13 @@ int lensing_d00( class_setup_parallel(); for (index_mu=0;index_muroot,"cl.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cl.dat"); class_call(output_open_cl_file(phr, pop, @@ -328,7 +328,7 @@ int output_cl( if (ple->has_lensed_cls == _TRUE_) { - sprintf(file_name,"%s%s",pop->root,"cl_lensed.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cl_lensed.dat"); class_call(output_open_cl_file(phr, pop, @@ -347,14 +347,14 @@ int output_cl( if (_scalars_) { - sprintf(file_name,"%s%s",pop->root,"cls.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar mode"); } if (_tensors_) { - sprintf(file_name,"%s%s",pop->root,"clt.dat"); + class_sprintf(file_name,"%s%s",pop->root,"clt.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for tensor mode"); } @@ -389,105 +389,105 @@ int output_cl( if ((ppt->has_ad == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_ad)) { - sprintf(file_name,"%s%s",pop->root,"cls_ad.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_ad.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar adiabatic (AD) mode"); } if ((ppt->has_bi == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_bi)) { - sprintf(file_name,"%s%s",pop->root,"cls_bi.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_bi.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar baryon isocurvature (BI) mode"); } if ((ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s",pop->root,"cls_cdi.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_cdi.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar CDM isocurvature (CDI) mode"); } if ((ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_nid) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s",pop->root,"cls_nid.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_nid.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar neutrino density isocurvature (NID) mode"); } if ((ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_niv) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s",pop->root,"cls_niv.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_niv.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar neutrino velocity isocurvature (NIV) mode"); } if ((ppt->has_ad == _TRUE_) && (ppt->has_bi == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_bi)) { - sprintf(file_name,"%s%s",pop->root,"cls_ad_bi.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_ad_bi.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross ADxBI mode"); } if ((ppt->has_ad == _TRUE_) && (ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s",pop->root,"cls_ad_cdi.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_ad_cdi.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross ADxCDI mode"); } if ((ppt->has_ad == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s",pop->root,"cls_ad_nid.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_ad_nid.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross ADxNID mode"); } if ((ppt->has_ad == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s",pop->root,"cls_ad_niv.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_ad_niv.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross ADxNIV mode"); } if ((ppt->has_bi == _TRUE_) && (ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s",pop->root,"cls_bi_cdi.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_bi_cdi.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross BIxCDI mode"); } if ((ppt->has_bi == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s",pop->root,"cls_bi_nid.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_bi_nid.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross BIxNID mode"); } if ((ppt->has_bi == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s",pop->root,"cls_bi_niv.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_bi_niv.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross BIxNIV mode"); } if ((ppt->has_cdi == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s",pop->root,"cls_cdi_nid.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_cdi_nid.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross CDIxNID mode"); } if ((ppt->has_cdi == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s",pop->root,"cls_cdi_niv.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_cdi_niv.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross CDIxNIV mode"); } if ((ppt->has_nid == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_nid) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s",pop->root,"cls_nid_niv.dat"); + class_sprintf(file_name,"%s%s",pop->root,"cls_nid_niv.dat"); strcpy(first_line,"[l(l+1)/2pi] C_l's for scalar cross NIDxNIV mode"); } @@ -684,15 +684,15 @@ int output_pk( if ((pfo->has_pk_m == _TRUE_) && (index_pk == pfo->index_pk_m)) { if (pk_output == pk_linear) - sprintf(type_suffix,"pk"); + class_sprintf(type_suffix,"pk"); else - sprintf(type_suffix,"pk_nl"); + class_sprintf(type_suffix,"pk_nl"); } if ((pfo->has_pk_cb == _TRUE_) && (index_pk == pfo->index_pk_cb)) { if (pk_output == pk_linear) - sprintf(type_suffix,"pk_cb"); + class_sprintf(type_suffix,"pk_cb"); else - sprintf(type_suffix,"pk_cb_nl"); + class_sprintf(type_suffix,"pk_cb_nl"); } /** - loop over z */ @@ -708,11 +708,11 @@ int output_pk( if (pop->z_pk_num == 1) redshift_suffix[0]='\0'; else - sprintf(redshift_suffix,"z%d_",index_z+1); + class_sprintf(redshift_suffix,"z%d_",index_z+1); /** - second, open only the relevant files and write a header in each of them */ - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,".dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,".dat"); class_call(output_open_pk_file(pba, pfo, @@ -732,77 +732,77 @@ int output_pk( for (index_ic2 = index_ic1; index_ic2 < pfo->ic_size; index_ic2++) { if ((ppt->has_ad == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_ad)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad.dat"); strcpy(first_line,"for adiabatic (AD) mode "); } if ((ppt->has_bi == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_bi)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi.dat"); strcpy(first_line,"for baryon isocurvature (BI) mode "); } if ((ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi.dat"); strcpy(first_line,"for CDM isocurvature (CDI) mode "); } if ((ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_nid) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_nid.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_nid.dat"); strcpy(first_line,"for neutrino density isocurvature (NID) mode "); } if ((ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_niv) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_niv.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_niv.dat"); strcpy(first_line,"for neutrino velocity isocurvature (NIV) mode "); } if ((ppt->has_ad == _TRUE_) && (ppt->has_bi == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_bi)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_bi.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_bi.dat"); strcpy(first_line,"for cross ADxBI mode "); } if ((ppt->has_ad == _TRUE_) && (ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_cdi.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_cdi.dat"); strcpy(first_line,"for cross ADxCDI mode "); } if ((ppt->has_ad == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_nid.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_nid.dat"); strcpy(first_line,"for scalar cross ADxNID mode "); } if ((ppt->has_ad == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_ad) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_niv.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_ad_niv.dat"); strcpy(first_line,"for cross ADxNIV mode "); } if ((ppt->has_bi == _TRUE_) && (ppt->has_cdi == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_cdi)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_cdi.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_cdi.dat"); strcpy(first_line,"for cross BIxCDI mode "); } if ((ppt->has_bi == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_nid.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_nid.dat"); strcpy(first_line,"for cross BIxNID mode "); } if ((ppt->has_bi == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_bi) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_niv.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_bi_niv.dat"); strcpy(first_line,"for cross BIxNIV mode "); } if ((ppt->has_cdi == _TRUE_) && (ppt->has_nid == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_nid)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi_nid.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi_nid.dat"); strcpy(first_line,"for cross CDIxNID mode "); } if ((ppt->has_cdi == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_cdi) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi_niv.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_cdi_niv.dat"); strcpy(first_line,"for cross CDIxNIV mode "); } if ((ppt->has_nid == _TRUE_) && (ppt->has_niv == _TRUE_) && (index_ic1 == ppt->index_ic_nid) && (index_ic2 == ppt->index_ic_niv)) { - sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_nid_niv.dat"); + class_sprintf(file_name,"%s%s%s%s",pop->root,redshift_suffix,type_suffix,"_nid_niv.dat"); strcpy(first_line,"for cross NIDxNIV mode "); } @@ -962,7 +962,7 @@ int output_tk( if (pop->z_pk_num == 1) redshift_suffix[0]='\0'; else - sprintf(redshift_suffix,"z%d_",index_z+1); + class_sprintf(redshift_suffix,"z%d_",index_z+1); /** - second, open only the relevant files, and write a heading in each of them */ @@ -982,9 +982,9 @@ int output_tk( ppt->error_message, pop->error_message); if ((ppt->has_ad == _TRUE_) && (ppt->ic_size[index_md] == 1) ) - sprintf(file_name,"%s%s%s",pop->root,redshift_suffix,"tk.dat"); + class_sprintf(file_name,"%s%s%s",pop->root,redshift_suffix,"tk.dat"); else - sprintf(file_name,"%s%s%s%s%s",pop->root,redshift_suffix,"tk_",ic_suffix,".dat"); + class_sprintf(file_name,"%s%s%s%s%s",pop->root,redshift_suffix,"tk_",ic_suffix,".dat"); class_open(tkfile, file_name, "w", pop->error_message); @@ -1061,7 +1061,7 @@ int output_background( pba->error_message, pop->error_message); - sprintf(file_name,"%s%s",pop->root,"background.dat"); + class_sprintf(file_name,"%s%s",pop->root,"background.dat"); class_open(backfile,file_name,"w",pop->error_message); if (pop->write_header == _TRUE_) { @@ -1113,7 +1113,7 @@ int output_thermodynamics( pth->error_message, pop->error_message); - sprintf(file_name,"%s%s",pop->root,"thermodynamics.dat"); + class_sprintf(file_name,"%s%s",pop->root,"thermodynamics.dat"); class_open(thermofile,file_name,"w",pop->error_message); if (pop->write_header == _TRUE_) { @@ -1170,7 +1170,7 @@ int output_perturbations( if (ppt->has_scalars == _TRUE_){ index_md = ppt->index_md_scalars; k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]]; - sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_s.dat"); + class_sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_s.dat"); class_open(out, file_name, "w", ppt->error_message); fprintf(out,"#scalar perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k); output_print_data(out, @@ -1183,7 +1183,7 @@ int output_perturbations( if (ppt->has_vectors == _TRUE_){ index_md = ppt->index_md_vectors; k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]]; - sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_v.dat"); + class_sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_v.dat"); class_open(out, file_name, "w", ppt->error_message); fprintf(out,"#vector perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k); output_print_data(out, @@ -1196,7 +1196,7 @@ int output_perturbations( if (ppt->has_tensors == _TRUE_){ index_md = ppt->index_md_tensors; k = ppt->k[index_md][ppt->index_k_output_values[index_md*ppt->k_output_values_num+index_ikout]]; - sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_t.dat"); + class_sprintf(file_name,"%s%s%d%s",pop->root,"perturbations_k",index_ikout,"_t.dat"); class_open(out, file_name, "w", ppt->error_message); fprintf(out,"#tensor perturbations for mode k = %.*e Mpc^(-1)\n",_OUTPUTPRECISION_,k); output_print_data(out, @@ -1224,7 +1224,7 @@ int output_primordial( double * data; int size_data, number_of_titles; - sprintf(file_name,"%s%s",pop->root,"primordial_Pk.dat"); + class_sprintf(file_name,"%s%s",pop->root,"primordial_Pk.dat"); class_call(primordial_output_titles(ppt,ppm,titles), ppm->error_message, @@ -1278,7 +1278,7 @@ int output_heating(struct injection* pin, struct noninjection* pni, struct outpu if (pop->write_exotic_injection == _TRUE_){ /* File name */ - sprintf(file_name_injection,"%s%s",pop->root,"exotic_injection.dat"); + class_sprintf(file_name_injection,"%s%s",pop->root,"exotic_injection.dat"); /* Titles */ class_call(injection_output_titles(pin,titles_injection), @@ -1320,7 +1320,7 @@ int output_heating(struct injection* pin, struct noninjection* pni, struct outpu if (pop->write_noninjection == _TRUE_){ /* File name */ - sprintf(file_name_noninjection,"%s%s",pop->root,"photon_noninjection.dat"); + class_sprintf(file_name_noninjection,"%s%s",pop->root,"photon_noninjection.dat"); /* Titles */ class_call(noninjection_output_titles(pni,titles_noninjection), @@ -1380,7 +1380,7 @@ int output_distortions( if (pop->write_distortions==_TRUE_ && psd->has_distortions == _TRUE_){ /* File name */ - sprintf(file_name_heat,"%s%s",pop->root,"sd_heating.dat"); + class_sprintf(file_name_heat,"%s%s",pop->root,"sd_heating.dat"); /* Titles */ class_call(distortions_output_heat_titles(psd,titles_heat), @@ -1419,7 +1419,7 @@ int output_distortions( fclose(out_heat); /* File name */ - sprintf(file_name_distortion,"%s%s",pop->root,"sd_distortions.dat"); + class_sprintf(file_name_distortion,"%s%s",pop->root,"sd_distortions.dat"); /* Titles */ class_call(distortions_output_sd_titles(psd,titles_distortion), @@ -1602,41 +1602,41 @@ int output_open_cl_file( if (phr->has_dd == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ for (index_d2=index_d1; index_d2<=MIN(index_d1+phr->non_diag,phr->d_size-1); index_d2++){ - sprintf(tmp,"dens[%d]-dens[%d]",index_d1+1,index_d2+1); + class_sprintf(tmp,"dens[%d]-dens[%d]",index_d1+1,index_d2+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } } if (phr->has_td == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ - sprintf(tmp,"T-dens[%d]",index_d1+1); + class_sprintf(tmp,"T-dens[%d]",index_d1+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } if (phr->has_pd == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ - sprintf(tmp,"phi-dens[%d]",index_d1+1); + class_sprintf(tmp,"phi-dens[%d]",index_d1+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } if (phr->has_ll == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ for (index_d2=index_d1; index_d2<=MIN(index_d1+phr->non_diag,phr->d_size-1); index_d2++){ - sprintf(tmp,"lens[%d]-lens[%d]",index_d1+1,index_d2+1); + class_sprintf(tmp,"lens[%d]-lens[%d]",index_d1+1,index_d2+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } } if (phr->has_tl == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ - sprintf(tmp,"T-lens[%d]",index_d1+1); + class_sprintf(tmp,"T-lens[%d]",index_d1+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } if (phr->has_dl == _TRUE_){ for (index_d1=0; index_d1d_size; index_d1++){ for (index_d2=MAX(index_d1-phr->non_diag,0); index_d2<=MIN(index_d1+phr->non_diag,phr->d_size-1); index_d2++) { - sprintf(tmp,"dens[%d]-lens[%d]",index_d1+1,index_d2+1); + class_sprintf(tmp,"dens[%d]-lens[%d]",index_d1+1,index_d2+1); class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum); } } diff --git a/source/perturbations.c b/source/perturbations.c index 0be0b7d4..77edbd9e 100644 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -562,7 +562,7 @@ int perturbations_output_titles( class_store_columntitle(titles,"d_idr",pba->has_idr); if (pba->has_ncdm == _TRUE_) { for (n_ncdm=0; n_ncdm < pba->N_ncdm; n_ncdm++) { - sprintf(tmp,"d_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"d_ncdm[%d]",n_ncdm); class_store_columntitle(titles,tmp,_TRUE_); } } @@ -591,7 +591,7 @@ int perturbations_output_titles( class_store_columntitle(titles,"t_idr",pba->has_idr); if (pba->has_ncdm == _TRUE_) { for (n_ncdm=0; n_ncdm < pba->N_ncdm; n_ncdm++) { - sprintf(tmp,"t_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"t_ncdm[%d]",n_ncdm); class_store_columntitle(titles,tmp,_TRUE_); } } @@ -708,8 +708,6 @@ int perturbations_init( int index_k; /* running index for type of perturbation */ int index_tp; - /* pointer to one struct perturbations_workspace per thread (one if no openmp) */ - struct perturbations_workspace ** pppw; /* background quantities */ double w_fld_ini, w_fld_0,dw_over_da_fld,integral_fld; @@ -2278,7 +2276,6 @@ int perturbations_get_k_list( ppt->k_size[ppt->index_md_scalars] = index_k; class_realloc(ppt->k[ppt->index_md_scalars], - ppt->k[ppt->index_md_scalars], ppt->k_size[ppt->index_md_scalars]*sizeof(double), ppt->error_message); } @@ -2412,7 +2409,6 @@ int perturbations_get_k_list( ppt->k_size[ppt->index_md_vectors] = index_k; class_realloc(ppt->k[ppt->index_md_vectors], - ppt->k[ppt->index_md_vectors], ppt->k_size[ppt->index_md_vectors]*sizeof(double), ppt->error_message); } @@ -2546,7 +2542,6 @@ int perturbations_get_k_list( ppt->k_size[ppt->index_md_tensors] = index_k; class_realloc(ppt->k[ppt->index_md_tensors], - ppt->k[ppt->index_md_tensors], ppt->k_size[ppt->index_md_tensors]*sizeof(double), ppt->error_message); } @@ -3338,13 +3333,13 @@ int perturbations_prepare_k_output(struct background * pba, /* Non-cold dark matter */ if ((pba->has_ncdm == _TRUE_) && ((ppt->has_density_transfers == _TRUE_) || (ppt->has_velocity_transfers == _TRUE_) || (ppt->has_source_delta_m == _TRUE_))) { for (n_ncdm=0; n_ncdm < pba->N_ncdm; n_ncdm++){ - sprintf(tmp,"delta_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"delta_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->scalar_titles,tmp,_TRUE_); - sprintf(tmp,"theta_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"theta_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->scalar_titles,tmp,_TRUE_); - sprintf(tmp,"shear_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"shear_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->scalar_titles,tmp,_TRUE_); - sprintf(tmp,"cs2_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"cs2_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->scalar_titles,tmp,_TRUE_); } } @@ -3386,11 +3381,11 @@ int perturbations_prepare_k_output(struct background * pba, if (ppt->evolve_tensor_ncdm == _TRUE_) { for (n_ncdm=0; n_ncdm < pba->N_ncdm; n_ncdm++){ - sprintf(tmp,"delta_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"delta_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->tensor_titles,tmp,_TRUE_); - sprintf(tmp,"theta_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"theta_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->tensor_titles,tmp,_TRUE_); - sprintf(tmp,"shear_ncdm[%d]",n_ncdm); + class_sprintf(tmp,"shear_ncdm[%d]",n_ncdm); class_store_columntitle(ppt->tensor_titles,tmp,_TRUE_); } } @@ -8202,7 +8197,7 @@ int perturbations_print_variables(double tau, theta_idr = y[ppw->pv->index_pt_theta_idr]; if (ppt->idr_nature == idr_free_streaming){ - if ((ppw->approx[ppw->index_ap_tca_idm_dr] == (int)tca_idm_dr_on)){ + if (ppw->approx[ppw->index_ap_tca_idm_dr] == (int)tca_idm_dr_on){ shear_idr = ppw->tca_shear_idm_dr; } else{ @@ -8410,9 +8405,9 @@ int perturbations_print_variables(double tau, ppt->size_scalar_perturbation_data[ppw->index_ikout] = 0; } else{ - ppt->scalar_perturbations_data[ppw->index_ikout] = - (double*)realloc(ppt->scalar_perturbations_data[ppw->index_ikout], - sizeof(double)*(ppt->size_scalar_perturbation_data[ppw->index_ikout]+ppt->number_of_scalar_titles)); + class_realloc(ppt->scalar_perturbations_data[ppw->index_ikout], + (ppt->size_scalar_perturbation_data[ppw->index_ikout]+ppt->number_of_scalar_titles)*sizeof(double), + ppt->error_message); } storeidx = 0; dataptr = ppt->scalar_perturbations_data[ppw->index_ikout]+ @@ -8776,7 +8771,7 @@ int perturbations_derivs(double tau, photon_scattering_rate += pvecthermo[pth->index_th_dmu_idm_g]; S_idm_g = 4./3. * pvecback[pba->index_bg_rho_g] / pvecback[pba->index_bg_rho_idm]; } - if ((pth->has_idm_dr == _TRUE_)){ + if (pth->has_idm_dr == _TRUE_){ S_idm_dr = 4./3. * pvecback[pba->index_bg_rho_idr]/ pvecback[pba->index_bg_rho_idm]; dmu_idm_dr = pvecthermo[pth->index_th_dmu_idm_dr]; dmu_idr = pth->b_idr/pth->a_idm_dr*pba->Omega0_idr/pba->Omega0_idm*dmu_idm_dr; diff --git a/source/primordial.c b/source/primordial.c index a5a2ec17..4831d17e 100644 --- a/source/primordial.c +++ b/source/primordial.c @@ -3227,7 +3227,7 @@ int primordial_external_spectrum_init( char command_with_arguments[2*_ARGUMENT_LENGTH_MAX_]; FILE *process; int n_data_guess, n_data = 0; - double *k = NULL, *pks = NULL, *pkt = NULL, *tmp = NULL; + double *k = NULL, *pks = NULL, *pkt = NULL; double this_k, this_pks, this_pkt; int status; int index_k; @@ -3242,16 +3242,16 @@ int primordial_external_spectrum_init( /* Prepare the command */ /* If the command is just a "cat", no arguments need to be passed */ if (strncmp("cat ", ppm->command, 4) == 0) { - sprintf(arguments, " "); + class_sprintf(arguments, " "); } /* otherwise pass the list of arguments */ else { - sprintf(arguments, " %g %g %g %g %g %g %g %g %g %g", + class_sprintf(arguments, " %g %g %g %g %g %g %g %g %g %g", ppm->custom1, ppm->custom2, ppm->custom3, ppm->custom4, ppm->custom5, ppm->custom6, ppm->custom7, ppm->custom8, ppm->custom9, ppm->custom10); } /* write the actual command in a string */ - sprintf(command_with_arguments, "%s %s", ppm->command, arguments); + class_sprintf(command_with_arguments, "%s %s", ppm->command, arguments); if (ppm->primordial_verbose > 0) printf(" -> running: %s\n",command_with_arguments); @@ -3273,24 +3273,12 @@ 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)); - 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)); - class_test(tmp == NULL, - ppm->error_message, - "Error allocating memory to read the external spectrum.\n"); - pks = tmp; + class_realloc(k, n_data_guess*sizeof(double), ppm->error_message); + class_realloc(pks, n_data_guess*sizeof(double), ppm->error_message); if (ppt->has_tensors == _TRUE_) { - tmp = (double *)realloc(pkt, n_data_guess*sizeof(double)); - class_test(tmp == NULL, - ppm->error_message, - "Error allocating memory to read the external spectrum.\n"); - pkt = tmp; - }; - }; + class_realloc(pkt, n_data_guess*sizeof(double), ppm->error_message); + } + } /* Store */ k [n_data] = this_k; pks[n_data] = this_pks; @@ -3328,24 +3316,19 @@ int primordial_external_spectrum_init( ppm->lnk_size = n_data; /** - Make room */ class_realloc(ppm->lnk, - ppm->lnk, ppm->lnk_size*sizeof(double), ppm->error_message); class_realloc(ppm->lnpk[ppt->index_md_scalars], - ppm->lnpk[ppt->index_md_scalars], ppm->lnk_size*sizeof(double), ppm->error_message); class_realloc(ppm->ddlnpk[ppt->index_md_scalars], - ppm->ddlnpk[ppt->index_md_scalars], ppm->lnk_size*sizeof(double), ppm->error_message); if (ppt->has_tensors == _TRUE_) { class_realloc(ppm->lnpk[ppt->index_md_tensors], - ppm->lnpk[ppt->index_md_tensors], ppm->lnk_size*sizeof(double), ppm->error_message); class_realloc(ppm->ddlnpk[ppt->index_md_tensors], - ppm->ddlnpk[ppt->index_md_tensors], ppm->lnk_size*sizeof(double), ppm->error_message); }; diff --git a/source/thermodynamics.c b/source/thermodynamics.c index c1683f56..b5cda4a5 100644 --- a/source/thermodynamics.c +++ b/source/thermodynamics.c @@ -1553,9 +1553,9 @@ int thermodynamics_solve( struct thermodynamics_parameters_and_workspace tpaw; /* function pointer to ODE evolver and names of possible evolvers. */ - extern int evolver_rk(); - extern int evolver_ndf15(); - int (*generic_evolver)() = evolver_ndf15; + extern int evolver_rk(EVOLVER_PROTOTYPE); + extern int evolver_ndf15(EVOLVER_PROTOTYPE); + int (*generic_evolver)(EVOLVER_PROTOTYPE) = evolver_ndf15; /** - choose evolver */ switch (ppr->thermo_evolver) { @@ -2190,9 +2190,9 @@ int thermodynamics_reionization_evolve_with_tau( struct thermo_workspace * ptw; /* function pointer to ODE evolver and names of possible evolvers */ - extern int evolver_rk(); - extern int evolver_ndf15(); - int (*generic_evolver)() = evolver_ndf15; + extern int evolver_rk(EVOLVER_PROTOTYPE); + extern int evolver_ndf15(EVOLVER_PROTOTYPE); + int (*generic_evolver)(EVOLVER_PROTOTYPE) = evolver_ndf15; /* pointers towards two thermo vector stuctures (see below) */ @@ -4647,7 +4647,7 @@ int thermodynamics_obtain_z_ini( if (pth->thermodynamics_verbose > 3) printf("The decoupling redshift for idm_dr is z_idm_dec = %.5e\n", z_idm_dec); /* we need to be careful if idm is coupled to photons and idr at the same time */ - class_test(z_idm_dec_min != _HUGE_ && abs(pba->T_idr - pba->T_cmb) > 1e-2, + class_test(z_idm_dec_min != _HUGE_ && fabs(pba->T_idr - pba->T_cmb) > 1e-2, pth->error_message, "It seems that at early times idm is thermally coupled to both idr and photons (possibly through baryons).\nPlease set the initial temperatures equal or disable this error."); diff --git a/source/transfer.c b/source/transfer.c index 60d5330c..026b4554 100644 --- a/source/transfer.c +++ b/source/transfer.c @@ -1218,7 +1218,6 @@ int transfer_get_q_list( /* now, readjust array size */ class_realloc(ptr->q, - ptr->q, ptr->q_size*sizeof(double), ptr->error_message); @@ -4658,7 +4657,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, int fevals=0, index_l_mid; int multiplier; int right_boundary_checked = _FALSE_; - int hil=0,hir=0,bini=0; class_call(get_xmin_generic(sgnK, lvec[0], nu, @@ -4692,7 +4690,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, } /* Hunt for left boundary: */ for (multiplier=1; ;multiplier *= 5){ - hil++; class_call(get_xmin_generic(sgnK, lvec[*index_l_left], nu, @@ -4702,7 +4699,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, &fevals), error_message, error_message); - //printf("Hunt left, iter = %d, x_nonzero=%g\n",hil,x_nonzero); if (x_nonzero <= xmax){ //Boundary found break; @@ -4722,8 +4718,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, /* If not found, hunt for right boundary: */ if (right_boundary_checked == _FALSE_){ for (multiplier=1; ;multiplier *= 5){ - hir++; - //printf("right iteration %d,index_l_right:%d\n",hir,*index_l_right); class_call(get_xmin_generic(sgnK, lvec[*index_l_right], nu, @@ -4755,7 +4749,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, // printf("Do binary search in get_lmax. \n"); //printf("Region: [%d, %d]\n",*index_l_left,*index_l_right); while (((*index_l_right) - (*index_l_left)) > 1) { - bini++; index_l_mid= (int)(0.5*((*index_l_right)+(*index_l_left))); //printf("left:%d, mid=%d, right=%d\n",*index_l_left,index_l_mid,*index_l_right); class_call(get_xmin_generic(sgnK, @@ -4773,9 +4766,6 @@ int transfer_get_lmax(int (*get_xmin_generic)(int sgnK, *index_l_right=index_l_mid; } //printf("Done\n"); - /* printf("Hunt left iter=%d, hunt right iter=%d (fevals: %d). For binary search: %d (fevals: %d)\n", - hil,hir,fevalshunt,bini,fevals); - */ return _SUCCESS_; } diff --git a/tools/arrays.c b/tools/arrays.c index 13c82bfb..dca382dc 100644 --- a/tools/arrays.c +++ b/tools/arrays.c @@ -80,7 +80,7 @@ int array_derive_spline( h = x_array[i+1] - x_array[i]; if (h == 0) { - sprintf(errmsg,"%s(L:%d) h=0, stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) h=0, stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -133,7 +133,7 @@ int array_derive_spline_table_line_to_line( h = x_array[i+1] - x_array[i]; if (h == 0) { - sprintf(errmsg,"%s(L:%d) h=0, stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) h=0, stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -165,7 +165,7 @@ int array_derive1_order2_table_line_to_line( double dxp,dxm,dyp,dym; if (n_lines < 2) { - sprintf(errmsg,"%s(L:%d) routine called with n_lines=%d, should be at least 2",__func__,__LINE__,n_lines); + class_sprintf(errmsg,"%s(L:%d) routine called with n_lines=%d, should be at least 2",__func__,__LINE__,n_lines); return _FAILURE_; } @@ -175,7 +175,7 @@ int array_derive1_order2_table_line_to_line( dym = *(array+0*n_columns+index_y) - *(array+1*n_columns+index_y); if ((dxp*dxm*(dxm-dxp)) == 0.) { - sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -192,7 +192,7 @@ int array_derive1_order2_table_line_to_line( dym = *(array+(i-1)*n_columns+index_y) - *(array+i*n_columns+index_y); if ((dxp*dxm*(dxm-dxp)) == 0.) { - sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -228,7 +228,7 @@ int array_derive2_order2_table_line_to_line( dym = *(array+(i-1)*n_columns+index_y) - *(array+i*n_columns+index_y); if ((dxp*dxm*(dxm-dxp)) == 0.) { - sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -297,7 +297,7 @@ int array_derive_two( double dx1,dx2,dy1,dy2,weight1,weight2; if ((index_dydx == index_x) || (index_dydx == index_y)) { - sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d and %d",__func__,__LINE__,index_dydx,index_x,index_y); + class_sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d and %d",__func__,__LINE__,index_dydx,index_x,index_y); return _FAILURE_; } @@ -314,7 +314,7 @@ int array_derive_two( weight2 = dx1*dx1; if ((dx1 == 0.) && (dx2 == 0.)) { - sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) stop to avoid division by zero",__func__,__LINE__); return _FAILURE_; } @@ -352,13 +352,13 @@ int array_spline( double dy_last; if (n_lines < 3) { - sprintf(errmsg,"%s(L:%d) n_lines=%d, while routine needs n_lines >= 3",__func__,__LINE__,n_lines); + class_sprintf(errmsg,"%s(L:%d) n_lines=%d, while routine needs n_lines >= 3",__func__,__LINE__,n_lines); return _FAILURE_; } u = (double*)malloc((n_lines-1) * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } @@ -387,7 +387,7 @@ int array_spline( -dy_first); } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -433,7 +433,7 @@ int array_spline( (*(array+(n_lines-1)*n_columns+index_x) - *(array+(n_lines-2)*n_columns+index_x))); } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -472,7 +472,7 @@ int array_spline_table_line_to_line( u = (double*)malloc((n_lines-1) * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } @@ -494,7 +494,7 @@ int array_spline_table_line_to_line( (x[1] - x[0])-dy_first); } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -534,7 +534,7 @@ int array_spline_table_line_to_line( (x[n_lines-1] - x[n_lines-2])); } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -578,19 +578,19 @@ int array_spline_table_lines( un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } if (p == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); return _FAILURE_; } if (qn == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); return _FAILURE_; } if (un == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); return _FAILURE_; } @@ -625,7 +625,7 @@ int array_spline_table_lines( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -683,7 +683,7 @@ int array_spline_table_lines( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -739,19 +739,19 @@ int array_logspline_table_lines( qn = (double*)malloc(y_size * sizeof(double)); un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } if (p == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); return _FAILURE_; } if (qn == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); return _FAILURE_; } if (un == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); return _FAILURE_; } @@ -786,7 +786,7 @@ int array_logspline_table_lines( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -844,7 +844,7 @@ int array_logspline_table_lines( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -901,19 +901,19 @@ int array_spline_table_columns( qn = (double*)malloc(y_size * sizeof(double)); un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } if (p == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); return _FAILURE_; } if (qn == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); return _FAILURE_; } if (un == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); return _FAILURE_; } @@ -959,7 +959,7 @@ int array_spline_table_columns( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1016,7 +1016,7 @@ int array_spline_table_columns( } } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1069,19 +1069,19 @@ int array_spline_table_columns2( qn = (double*)malloc(y_size * sizeof(double)); un = (double*)malloc(y_size * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } if (p == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate p",__func__,__LINE__); return _FAILURE_; } if (qn == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate qn",__func__,__LINE__); return _FAILURE_; } if (un == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate un",__func__,__LINE__); return _FAILURE_; } @@ -1210,7 +1210,7 @@ int array_spline_table_one_column( u = (double*)malloc((x_size-1) * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } @@ -1243,7 +1243,7 @@ int array_spline_table_one_column( } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1296,7 +1296,7 @@ int array_spline_table_one_column( } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1345,7 +1345,7 @@ int array_logspline_table_one_column( u = (double*)malloc((x_stop-1) * sizeof(double)); if (u == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate u",__func__,__LINE__); return _FAILURE_; } @@ -1378,7 +1378,7 @@ int array_logspline_table_one_column( } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1432,7 +1432,7 @@ int array_logspline_table_one_column( } else { - sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); + class_sprintf(errmsg,"%s(L:%d) Spline mode not identified: %d",__func__,__LINE__,spline_mode); return _FAILURE_; } } @@ -1528,7 +1528,7 @@ int array_integrate_all_trapzd_or_spline( double h; if ((index_start_spline<0) || (index_start_spline>=n_lines)) { - sprintf(errmsg,"%s(L:%d) index_start_spline outside of range",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) index_start_spline outside of range",__func__,__LINE__); return _FAILURE_; } @@ -1576,7 +1576,7 @@ int array_integrate( double sum; if ((index_int_y_dx == index_x) || (index_int_y_dx == index_y)) { - sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d and %d",__func__,__LINE__,index_int_y_dx,index_x,index_y); + class_sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d and %d",__func__,__LINE__,index_int_y_dx,index_x,index_y); return _FAILURE_; } @@ -1612,7 +1612,7 @@ int array_integrate_ratio( double sum; if ((index_int_y1_over_y2_dx == index_x) || (index_int_y1_over_y2_dx == index_y1) || (index_int_y1_over_y2_dx == index_y2)) { - sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d, %d and %d",__func__,__LINE__,index_int_y1_over_y2_dx,index_x,index_y1,index_y2); + class_sprintf(errmsg,"%s(L:%d) : Output column %d must differ from input columns %d, %d and %d",__func__,__LINE__,index_int_y1_over_y2_dx,index_x,index_y1,index_y2); return _FAILURE_; } @@ -1658,12 +1658,12 @@ int array_interpolate( if (*(array+inf*n_columns+index_x) < *(array+sup*n_columns+index_x)){ if (x < *(array+inf*n_columns+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array+inf*n_columns+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array+inf*n_columns+index_x)); return _FAILURE_; } if (x > *(array+sup*n_columns+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array+sup*n_columns+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array+sup*n_columns+index_x)); return _FAILURE_; } @@ -1680,12 +1680,12 @@ int array_interpolate( else { if (x < *(array+sup*n_columns+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array+sup*n_columns+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array+sup*n_columns+index_x)); return _FAILURE_; } if (x > *(array+inf*n_columns+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array+inf*n_columns+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array+inf*n_columns+index_x)); return _FAILURE_; } @@ -1738,12 +1738,12 @@ int array_interpolate_spline_transposed(double * array, if (array[inf*y_size+index_x] < array[sup*y_size+index_x]){ if (x < array[inf*y_size+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array[inf*y_size+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array[inf*y_size+index_x]); return _FAILURE_; } if (x > array[sup*y_size+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array[sup*y_size+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array[sup*y_size+index_x]); return _FAILURE_; } @@ -1760,12 +1760,12 @@ int array_interpolate_spline_transposed(double * array, else { if (x < array[sup*y_size+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array[sup*y_size+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array[sup*y_size+index_x]); return _FAILURE_; } if (x > array[inf*y_size+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array[inf*y_size+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array[inf*y_size+index_x]); return _FAILURE_; } @@ -1817,12 +1817,12 @@ int array_interpolate_spline( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -1839,12 +1839,12 @@ int array_interpolate_spline( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -1894,12 +1894,12 @@ int array_search_bisect( if (array[inf] < array[sup]){ if (c < array[inf]) { - sprintf(errmsg,"%s(L:%d) : c=%e < y_min=%e",__func__,__LINE__,c,array[inf]); + class_sprintf(errmsg,"%s(L:%d) : c=%e < y_min=%e",__func__,__LINE__,c,array[inf]); return _FAILURE_; } if (c > array[sup]) { - sprintf(errmsg,"%s(L:%d) : c=%e > y_max=%e",__func__,__LINE__,c,array[sup]); + class_sprintf(errmsg,"%s(L:%d) : c=%e > y_max=%e",__func__,__LINE__,c,array[sup]); return _FAILURE_; } @@ -1916,12 +1916,12 @@ int array_search_bisect( else { if (c < array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,c,array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,c,array[sup]); return _FAILURE_; } if (c > array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,c,array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,c,array[inf]); return _FAILURE_; } @@ -1965,12 +1965,12 @@ int array_interpolate_linear( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -1987,12 +1987,12 @@ int array_interpolate_linear( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -2047,12 +2047,12 @@ int array_interpolate_logspline( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -2069,12 +2069,12 @@ int array_interpolate_logspline( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -2132,12 +2132,12 @@ int array_interpolate_spline_one_column( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -2154,12 +2154,12 @@ int array_interpolate_spline_one_column( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -2232,12 +2232,12 @@ int array_interpolate_extrapolate_spline_one_column( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -2254,12 +2254,12 @@ int array_interpolate_extrapolate_spline_one_column( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -2340,12 +2340,12 @@ int array_interpolate_extrapolate_logspline_loglinear_one_column( if (x_array[inf] < x_array[sup]){ if (x < x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } if (x > x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } @@ -2362,12 +2362,12 @@ int array_interpolate_extrapolate_logspline_loglinear_one_column( else { if (x < x_array[sup]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,x_array[sup]); return _FAILURE_; } if (x > x_array[inf]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,x_array[inf]); return _FAILURE_; } @@ -2420,7 +2420,7 @@ int array_interpolate_growing_closeby( while (x < *(array+inf*n_columns+index_x)) { inf--; if (inf < 0) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, x,array[index_x]); return _FAILURE_; } @@ -2429,7 +2429,7 @@ int array_interpolate_growing_closeby( while (x > *(array+sup*n_columns+index_x)) { sup++; if (sup > (n_lines-1)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, x,array[(n_lines-1)*n_columns+index_x]); return _FAILURE_; } @@ -2474,7 +2474,7 @@ int array_interpolate_one_growing_closeby( while (x < *(array+inf*n_columns+index_x)) { inf--; if (inf < 0) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, x,array[index_x]); return _FAILURE_; } @@ -2483,7 +2483,7 @@ int array_interpolate_one_growing_closeby( while (x > *(array+sup*n_columns+index_x)) { sup++; if (sup > (n_lines-1)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, x,array[(n_lines-1)*n_columns+index_x]); return _FAILURE_; } @@ -2521,11 +2521,11 @@ int array_interpolate_spline_growing_closeby( /* if (*last_index < 0) { - sprintf(errmsg,"%s(L:%d) problem with last_index =%d < 0",__func__,__LINE__,*last_index); + class_sprintf(errmsg,"%s(L:%d) problem with last_index =%d < 0",__func__,__LINE__,*last_index); return _FAILURE_; } if (*last_index > (n_lines-1)) { - sprintf(errmsg,"%s(L:%d) problem with last_index =%d > %d",__func__,__LINE__,*last_index,n_lines-1); + class_sprintf(errmsg,"%s(L:%d) problem with last_index =%d > %d",__func__,__LINE__,*last_index,n_lines-1); return _FAILURE_; } */ @@ -2537,7 +2537,7 @@ int array_interpolate_spline_growing_closeby( while (x < x_array[inf]) { inf--; if (inf < 0) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, x,x_array[0]); return _FAILURE_; } @@ -2546,7 +2546,7 @@ int array_interpolate_spline_growing_closeby( while (x > x_array[sup]) { sup++; if (sup > (n_lines-1)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, x,x_array[n_lines-1]); return _FAILURE_; } @@ -2593,7 +2593,7 @@ int array_interpolate_spline_growing_hunt( if (x >= x_array[*last_index]) { if (x > x_array[n_lines-1]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, x,x_array[n_lines-1]); return _FAILURE_; } @@ -2620,7 +2620,7 @@ int array_interpolate_spline_growing_hunt( } else { if (x < x_array[0]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__, x,x_array[0]); return _FAILURE_; } @@ -2688,7 +2688,7 @@ int array_spline_hunt(double* x_array, inc=1; if (x >= x_array[last_index]) { if (x > x_array[x_size-1]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__, x,x_array[x_size-1]); return _FAILURE_; } @@ -2715,7 +2715,7 @@ int array_spline_hunt(double* x_array, } else { if (x < x_array[0]) { - sprintf(errmsg,"%s(L:%d) : x=%.20e < x_min=%.20e",__func__,__LINE__, + class_sprintf(errmsg,"%s(L:%d) : x=%.20e < x_min=%.20e",__func__,__LINE__, x,x_array[0]); return _FAILURE_; } @@ -2778,12 +2778,12 @@ int array_interpolate_two( if (x < array_x[inf*n_columns_x+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array_x[inf*n_columns_x+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array_x[inf*n_columns_x+index_x]); return _FAILURE_; } if (x > array_x[sup*n_columns_x+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array_x[sup*n_columns_x+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array_x[sup*n_columns_x+index_x]); return _FAILURE_; } @@ -2800,12 +2800,12 @@ int array_interpolate_two( else { if (x < *(array_x+sup*n_columns_x+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array_x+sup*n_columns_x+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array_x+sup*n_columns_x+index_x)); return _FAILURE_; } if (x > *(array_x+inf*n_columns_x+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array_x+inf*n_columns_x+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array_x+inf*n_columns_x+index_x)); return _FAILURE_; } @@ -2853,12 +2853,12 @@ int array_interpolate_two_bis( if (x < array_x[inf*n_columns_x+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array_x[inf*n_columns_x+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,array_x[inf*n_columns_x+index_x]); return _FAILURE_; } if (x > array_x[sup*n_columns_x+index_x]) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array_x[sup*n_columns_x+index_x]); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,array_x[sup*n_columns_x+index_x]); return _FAILURE_; } @@ -2875,12 +2875,12 @@ int array_interpolate_two_bis( else { if (x < *(array_x+sup*n_columns_x+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array_x+sup*n_columns_x+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e < x_min=%e",__func__,__LINE__,x,*(array_x+sup*n_columns_x+index_x)); return _FAILURE_; } if (x > *(array_x+inf*n_columns_x+index_x)) { - sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array_x+inf*n_columns_x+index_x)); + class_sprintf(errmsg,"%s(L:%d) : x=%e > x_max=%e",__func__,__LINE__,x,*(array_x+inf*n_columns_x+index_x)); return _FAILURE_; } @@ -2991,12 +2991,12 @@ int array_interpolate_equal( double x_step,x_minus,weight; if (x < x_min) { - sprintf(errmsg,"%s(L:%d) : x out of bounds: x=%e,x_min=%e",__func__,__LINE__,x,x_min); + class_sprintf(errmsg,"%s(L:%d) : x out of bounds: x=%e,x_min=%e",__func__,__LINE__,x,x_min); return _FAILURE_; } if (x > x_max) { - sprintf(errmsg,"%s(L:%d) : x out of bounds: x=%e,x_max=%e",__func__,__LINE__,x,x_max); + class_sprintf(errmsg,"%s(L:%d) : x out of bounds: x=%e,x_max=%e",__func__,__LINE__,x,x_max); return _FAILURE_; } @@ -3136,7 +3136,7 @@ int array_smooth_trg(double * array, smooth=(double*)malloc(k_size*sizeof(double)); if (smooth == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); return _FAILURE_; } @@ -3266,7 +3266,7 @@ int array_smooth(double * array, smooth=(double*)malloc(n_lines*sizeof(double)); if (smooth == NULL) { - sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) Cannot allocate smooth",__func__,__LINE__); return _FAILURE_; } @@ -3489,11 +3489,11 @@ int array_hunt_descending( /* checks */ if (array[i_inf] < array[i_sup]) { - sprintf(errmsg,"%s(L:%d) array is not in descending order (checked only the boundaries)",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) array is not in descending order (checked only the boundaries)",__func__,__LINE__); return _FAILURE_; } if ((value > array[i_inf]) || (value < array[i_sup])) { - sprintf(errmsg,"%s(L:%d) %e is outside the range [%e, %e]",__func__,__LINE__,value,array[size-1],array[0]); + class_sprintf(errmsg,"%s(L:%d) %e is outside the range [%e, %e]",__func__,__LINE__,value,array[size-1],array[0]); return _FAILURE_; } @@ -3542,11 +3542,11 @@ int array_hunt_ascending( /* checks */ if (array[i_inf] > array[i_sup]) { - sprintf(errmsg,"%s(L:%d) array is not in ascending order (checked only the boundaries)",__func__,__LINE__); + class_sprintf(errmsg,"%s(L:%d) array is not in ascending order (checked only the boundaries)",__func__,__LINE__); return _FAILURE_; } if ((value < array[i_inf]) || (value > array[i_sup])) { - sprintf(errmsg,"%s(L:%d) %e is outside the range [%e, %e]",__func__,__LINE__,value,array[0],array[size-1]); + class_sprintf(errmsg,"%s(L:%d) %e is outside the range [%e, %e]",__func__,__LINE__,value,array[0],array[size-1]); return _FAILURE_; } diff --git a/tools/hyperspherical.c b/tools/hyperspherical.c index 25e8fe01..1a679741 100644 --- a/tools/hyperspherical.c +++ b/tools/hyperspherical.c @@ -27,7 +27,7 @@ int hyperspherical_HIS_create(int K, then relative to ppHIS. */ double deltax, beta2, lambda, x, xfwd; - double *sqrtK, *one_over_sqrtK,*PhiL; + double *sqrtK, *one_over_sqrtK; int j, k, l, nx, lmax, l_recurrence_max; beta2 = beta*beta; @@ -881,6 +881,9 @@ int hyperspherical_WKB(int K,int l,double beta,double y, double *Phi){ airy_sign = 1; } } + else{ + return _FAILURE_; + } argu = 3.0*S/(2.0*e); Q = CscK*CscK-alpha2; C = 0.5*sqrt(alpha)/beta; diff --git a/tools/parser.c b/tools/parser.c index 602e2b72..66dabae7 100644 --- a/tools/parser.c +++ b/tools/parser.c @@ -632,15 +632,15 @@ int parser_cat(struct file_content * pfc1, if (pfc1->size == 0) { class_alloc(pfc3->filename,(strlen(pfc2->filename)+1)*sizeof(char),errmsg); - sprintf(pfc3->filename,"%s",pfc2->filename); + class_sprintf(pfc3->filename,"%s",pfc2->filename); } if (pfc2->size == 0) { class_alloc(pfc3->filename,(strlen(pfc1->filename)+1)*sizeof(char),errmsg); - sprintf(pfc3->filename,"%s",pfc1->filename); + class_sprintf(pfc3->filename,"%s",pfc1->filename); } if ((pfc1->size !=0) && (pfc2->size != 0)) { class_alloc(pfc3->filename,(strlen(pfc1->filename)+strlen(pfc2->filename)+5)*sizeof(char),errmsg); - sprintf(pfc3->filename,"%s or %s",pfc1->filename,pfc2->filename); + class_sprintf(pfc3->filename,"%s or %s",pfc1->filename,pfc2->filename); } pfc3->size = pfc1->size + pfc2->size; diff --git a/tools/quadrature.c b/tools/quadrature.c index 053cd040..41921e4a 100644 --- a/tools/quadrature.c +++ b/tools/quadrature.c @@ -304,17 +304,17 @@ int get_qsampling(double *x, } //printf("N_adapt=%d, N_combined=%d at level=%d, Nlag=%d\n",Nadapt,N_comb,level,NLag); if (adapt_converging==_TRUE_){ - sprintf(method_chosen,"Adaptive Gauss-Kronrod Quadrature"); + class_sprintf(method_chosen,"Adaptive Gauss-Kronrod Quadrature"); /* Gather weights and xvalues from tree: */ i = Nadapt-1; get_leaf_x_and_w(root,&i,x,w,_TRUE_); } else if (Laguerre_converging==_TRUE_){ - sprintf(method_chosen,"Gauss-Laguerre Quadrature"); + class_sprintf(method_chosen,"Gauss-Laguerre Quadrature"); /* x and w is already populated in this case. */ } else if (combined_converging == _TRUE_){ - sprintf(method_chosen,"Combined Quadrature"); + class_sprintf(method_chosen,"Combined Quadrature"); for(i=0; i