Performance and Stability Optimisation

Performance

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.

Compile flags

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.

Stability

Rejecting Unphysical Scalar Field Parameter Combinations

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:

1. $V\left(\phi\right)=\infty$ or NaN, or any of its derivatives
2. Scalar field values causing the integration to produce $\infty$ or NaN
3. Scalar field momentum (perturbation) integration producing $\infty$ or NaN
4. General check in evolver for $\infty$ or NaN results 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 $\infty$ and 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 $\infty$, yet, for a 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);
    }
  }

Adding New Quantities to Background Table

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); \
    }                                                                                      \
  }