-
Notifications
You must be signed in to change notification settings - Fork 0
Performance and Stability Optimisation
CLASS tends to perform many computations like
pow(pba->H0, 2) / pow(a, 4)that require fetching data from a struct and calling a function like pow() that are less optimisable by the compiler than a simple multiplication (as those function also provide additional checks and safeguards). To slightly speed up those computations, we defined new variables like
double a2 = a * a;
double a3 = a2 * a;
double a4 = a3 * a;
double H0_sq = pba->H0 * pba->H0;and substituted those in the code. Besides this general substitutions, we have many function- or loop-specific substitutions to accelerate the computation, e.g. by writing for-loops in a way that should make it easy for the compiler to unroll the loop and vectorise the computation, using modern CPU instructions sets, or by first gathering parameter values and caching powers of those values that are repeatedly used. We mention such implementations in the code and the corresponding Wiki sections. An extensive example of the former is the scalar field potential. An example of the latter is the scalar field coupling to dark matter.
Besides the CLASS default compile flag -O3, we also activated the following flags to speed up the binary:
# your optimization flag
OPTFLAG = -O3
OPTFLAG += -funroll-loops -ftree-vectorize -ftree-slp-vectorize -flto=auto -fPIC
OPTFLAG += -march=native -mtune=native
OPTFLAG += -fno-math-errno -fno-trapping-math # CLASS does not check for math errors, and these flags allow the compiler to vectorize more code, which can lead to significant speed-ups. See https://gcc.gnu.org/wiki/Vectorization#Floating_Point_Math_Errors for more details.This should encourage the compiler to unroll loops and vectorise whatever is possible, to make use of the modern CPU instruction set. The latter is also encouraged by the -march=native -mtune=native flags, which optimise for the specific CPU the compute node is using (AMD EPYC CPUs).
Additionally, we use profile-guided optimisation. First, the code is compiled with the OPTFLAG and LDFLAG compile flags -fprofile-generate. Then, when the code runs under a typical workload, it creates files in a pre-defined folder that inform the compiler about the hot paths of an actual production run. This information is then used in a second compilation with the flags -fprofile-use to optimise the binary with the typical use case in mind. This further reduces the runtime by 12.3%.
Some scalar field parameters can lead to segfaults when some derived quantities are passed to Cobaya. The MCMC run then needs to be resumed. Not only can resuming be avoided, but the entire computation within CLASS can sometimes be avoided if the input parameters are checked early on. We implemented the following checks:
-
$V\left(\phi\right)=\infty$ orNaN, or any of its derivatives - Scalar field values causing the integration to produce
$\infty$ orNaN - Scalar field momentum (perturbation) integration producing
$\infty$ orNaN - General check in evolver for
$\infty$ orNaNresults during integration These checks are implemented as follows:
(1.) First (directly with the initial scalar field values, before any expensive ODE integration is done) in the background_checks() function in background.c
class_test(isnan(V_check) || isinf(V_check) ||
isnan(dV_check) || isinf(dV_check) ||
isnan(ddV_check) || isinf(ddV_check),
pba->error_message,
"Scalar field potential (type %d) evaluates to NaN or Inf at phi_ini=%e "
"(V=%e, dV=%e, ddV=%e). This typically happens when exp(), cosh() or "
"similar functions overflow for extreme parameter combinations — "
"skipping this parameter point.",
pba->scf_potential, phi_ini_check, V_check, dV_check, ddV_check);and then during the cosmological evolution in the background_functions() (also in background.c)
class_test(isnan(V_phi) || isinf(V_phi) || isnan(dV_p) || isinf(dV_p) || isnan(ddV_val) || isinf(ddV_val),
pba->error_message,
"NaN or Inf in scalar field potential at phi=%e (V=%e, dV=%e, ddV=%e). "
"The scalar field parameter combination is likely unphysical.",
phi, V_phi, dV_p, ddV_val);
class_test(isnan(pvecback[pba->index_bg_dV_scf]) || isinf(pvecback[pba->index_bg_dV_scf]),
pba->error_message,
"NaN or Inf in effective scalar field derivative dV_eff at phi=%e. "
"The coupling or scalar field parameters may be unphysical.",
phi);This assures the scalar field potential is always well-defined.
(2.) Check in background_derivs() if the scalar field integration was successful and ready for the next time-step. If not, the computation is aborted and the MCMC can move on from this combination.
if (pba->has_scf == _TRUE_)
{
/** - Scalar field equation: \f$ \phi'' + 2 a H \phi' + a^2 (dV + rho') = 0 \f$ (note H is wrt cosmological time)
written as \f$ d\phi/dlna = phi' / (aH) \f$ and \f$ d\phi'/dlna = -2*phi' - (a/H) (dV + rho_prime) \f$ */
dy[pba->index_bi_phi_scf] = y[pba->index_bi_phi_prime_scf] / a / H;
/** - Use the pre-computed effective derivative pvecback[index_bg_dV_scf]
(= V'_pure + coupling) from background_functions() above,
avoiding a redundant dV_p_scf() switch + coupling_scf() evaluation. */
if (pba->model_cdm == 2)
{
dy[pba->index_bi_phi_prime_scf] = -2 * y[pba->index_bi_phi_prime_scf] - a * pvecback[pba->index_bg_dV_scf] / H - a * rho_cdm_prime(pba, y[pba->index_bi_phi_scf], pvecback) / H; // rho_cdm_prime is part of the Klein–Gordon equation
}
else
{
dy[pba->index_bi_phi_prime_scf] = -2 * y[pba->index_bi_phi_prime_scf] - a * pvecback[pba->index_bg_dV_scf] / H;
}
/** - Check for NaN/Inf in scalar field derivatives (can happen for extreme parameter combinations) */
class_test(isnan(dy[pba->index_bi_phi_scf]) || isinf(dy[pba->index_bi_phi_scf]) ||
isnan(dy[pba->index_bi_phi_prime_scf]) || isinf(dy[pba->index_bi_phi_prime_scf]),
error_message,
"NaN or Inf in scalar field background derivatives at a=%e, phi=%e, phi_prime=%e. "
"This typically indicates extreme/unphysical scalar field parameters.",
a, y[pba->index_bi_phi_scf], y[pba->index_bi_phi_prime_scf]);
}(3.) The scalar field momentum in perturbations_derivs() from perturbations.c is also checked for being well-defined:
/** - ----> catch NaN/Inf in scalar field perturbation derivatives
before they propagate to numjac and cause a segfault */
class_test(isnan(dy[pv->index_pt_phi_scf]) || isinf(dy[pv->index_pt_phi_scf]) ||
isnan(dy[pv->index_pt_mom_scf]) || isinf(dy[pv->index_pt_mom_scf]),
error_message,
"NaN or Inf in scalar field perturbation derivatives (tau=%e, k=%e). "
"Background quantities: phi=%e, phi_prime=%e, V=%e, dV=%e, ddV=%e. "
"The scalar field parameter combination is likely unphysical.",
tau, k,
pvecback[pba->index_bg_phi_scf],
pvecback[pba->index_bg_phi_prime_scf],
pvecback[pba->index_bg_V_scf],
pvecback[pba->index_bg_dV_scf],
pvecback[pba->index_bg_ddV_scf]);(4.) To generally ascertain that the evolver returns well-defined output, we implemented a brief additional check directly in the evolver_ndf15.c. On the one hand, before the main loop, the 'storage' of the integrator is checked for well-defined values. This assures that we start our integration from a proper basis.
/* Catch NaN/Inf in initial derivative evaluation before entering the main loop */
for (ii = 1; ii <= neq; ii++)
{
class_test(isnan(f0[ii]) || isinf(f0[ii]),
error_message,
"NaN or Inf in initial derivative evaluation (component %d of %d at t=%e). "
"The physical model produced unphysical values for these input parameters.",
ii, neq, t0);
}On the other hand, there is a per-column check for the Jacobian. It is sufficient to check the entire column, instead of every element individually, as in C99 NaN propagate through summation. Checking every element individually would cause a too large overhead for little gain. We also considered if it was possible that a sum of well-defined values could reach double we consider this highly unlikely. The check is implemented as follows:
/* The next section should work regardless of sparse...*/
/* Evaluate the function at y+delta vectors:*/
for (j = 1; j <= colmax; j++)
{
for (i = 1; i <= neq; i++)
{
nj_ws->yydel[i] = nj_ws->ydel_Fdel[i][j];
}
class_call((*derivs)(t, nj_ws->yydel + 1, nj_ws->ffdel + 1,
parameters_and_workspace_for_derivs, error_message),
error_message, error_message);
*nfe += 1;
{
double finite_check = 0.0;
for (i = 1; i <= neq; i++)
{
nj_ws->ydel_Fdel[i][j] = nj_ws->ffdel[i];
finite_check += nj_ws->ffdel[i];
}
/* A single isfinite() check on the accumulated sum catches any NaN or Inf
component, since NaN/Inf propagate through addition (IEEE 754). */
class_test(!isfinite(finite_check),
error_message,
"numjac: NaN or Inf detected in derivative vector at t=%e. "
"This usually means the physical model produced unphysical values.",
t);
}
}When adding new quantities to the background table, space must be reserved for the columns with the additional data. This can easily be missed. To implement a little safeguard against this developer-error, we changed strcpy to strncpy in the output_print_data() function in output.c:
strncpy(thetitle, titles, _MAXTITLESTRINGLENGTH_ - 1);
thetitle[_MAXTITLESTRINGLENGTH_ - 1] = '\0';
pch = strtok(thetitle, _DELIMITER_);We did the same for the class_store_columntitle() function from common.h:
#define class_store_columntitle(titlestring, \
title, \
condition) \
{ \
if (condition == _TRUE_) \
{ \
strncat(titlestring, title, _MAXTITLESTRINGLENGTH_ - strlen(titlestring) - 1); \
strncat(titlestring, _DELIMITER_, _MAXTITLESTRINGLENGTH_ - strlen(titlestring) - 1); \
} \
}evolver_ndf15 allocates a single contiguous buffer and carves it into sub-arrays for doubles, an int array (interpidx), a double-pointer array (dif), and a second double region for backward differences. The original layout placed the int array (neqp × sizeof(int)) between the double arrays and the double-pointer array:
[15×neqp doubles] [neqp ints] [neqp double*] [(7×neq+1) doubles]
When neqp is odd, neqp × 4 bytes leaves the subsequent double* and double regions at a 4-byte boundary rather than the 8-byte alignment required for these types. This constitutes undefined behaviour in C11 and is flagged at runtime by the -fsanitize=undefined CCFLAG / LDFLAG as dozens of misaligned load/store errors during thermodynamics integration.
The solution is to move the int array to the end of the buffer, after all double and double* regions:
[15×neqp doubles] [neqp double*] [(7×neq+1) doubles] [neqp ints]
Since malloc returns memory aligned to at least 8 bytes, and all regions before the int array are multiples of 8 bytes in size, every double and double* access is now aligned after the fix. The int array at the tail has no alignment requirement beyond 4 bytes, so placing it last is always safe.
This is implemented in evolver_ndf15.c as
void *buffer;
struct jacobian jac;
struct numjac_workspace nj_ws;
int neqp = neq + 1;
int buffer_size =
15 * neqp * sizeof(double) + neqp * sizeof(double *) + (7 * neq + 1) * sizeof(double) + neqp * sizeof(int);
int status;There was a tiny chance that line 50 in input.h
if (flag_temp == _TRUE_)
{
strcpy(destination, string_temp);
}might cause a buffer overflow, as the string could be 1024 bit long, yet the buffer reserved for it was set to store 1000 bit. It's very unlikely to encounter this, but no downside in preventing it:
#define _BASEPATHSIZE_ 1024 /**< allowed size of the base path */in common.h, instead of _BASEPATHSIZE_ 1000.
The pNG potential sometimes causes CLASS to hang. Only killing the process can remedy the situation, which sometimes even kills the entire terminal. It seems that the ndf15 solver encounters a stiff system of equations and just hangs forever. Typically, 5k evolver steps are taken at the background level, and 1~3M steps at the perturbation level. I implemented class_tests that trigger if more than 100k steps are taken at the background (20x margin) or 60M at the evolver level (20x margin from observed perturbations maximum). This will kill truly infinite loops and (hopefully) not trigger prematurely.
In background.c a test is run in background_derivs() (#define _BACKGROUND_MAX_FEVALS_ 100000):
pbpaw->derivs_call_count++;
class_test(pbpaw->derivs_call_count > _BACKGROUND_MAX_FEVALS_,
error_message,
"Background integration exceeded %d derivative evaluations "
"(loga ~ %g). This indicates extreme stiffness, likely from "
"pathological scalar field parameters.",
_BACKGROUND_MAX_FEVALS_, loga);In evolver_ndf15_inner() in evolver_ndf15.c, another test is run in the main loop of the evolver (#define _NDF15_MAX_STEPS_ 60000000):
class_test(stepstat[0] > _NDF15_MAX_STEPS_, error_message,
"ndf15: exceeded %d steps (fevals=%d) in interval [%g:%g], "
"current h=%g. The ODE system is too stiff to solve.",
_NDF15_MAX_STEPS_, stepstat[2], t0, tfinal, absh);During the MCMC, I observed segfaults on certain parameter combinations, sometimes followed by excessive use of server resources, which causes the run to be aborted. A possible source of that are memory leaks. There are three basic strategies used:
-
Replace class_call with class_call_except, with the exception causing that allocated memory is freed.
-
The evolvers are wrapped in such a way that other functions calling the evolver don't see a difference: the outer layer kept its name and returns the same numerical values and error messages, but the actual evolution is done in the inner layer. If the inner layer fails, the outer layer frees the allocated memory and passes the error messages.
-
Tables are NULL-initialised: if anything happens that triggered the freeing of a table before the table was actually initialised, this assures that it's defined behaviour.
The implementation of this required many changes that best viewed in the diff of commit dcbc1ca.