Odd interaction between `ode_adams` and threading

Description

I have noticed some odd behavior using ode_adams and STAN_THREADS. This does not recur when using ode_bdf or ode_rk45.

Example

Using a simple SIR model, courtesy of @charlesm93:

functions {
  vector sir(real t, vector y, real beta, real gamma, int N) {
    real S = y[1];
    real I = y[2];
    // real R = y[3];

    real dS_dt = - beta * I * S / N;
    real dI_dt = beta * I * S / N - gamma * I;
    real dR_dt = gamma * I;

    return [dS_dt, dI_dt, dR_dt]';
  }

}

data {
  int<lower=1> n_days;
  vector[3] y0;
  real t0;
  array[n_days] real ts;
  int N;
  array[n_days] int cases;
}

parameters {
  real<lower=0> gamma;
  real<lower=0> beta;
  real<lower=0> phi_inv;
}

transformed parameters {
  real phi = 1. / phi_inv;
  array[n_days] vector[3] y = ode_adams(sir, y0, t0, ts, beta, gamma, N);
}

model {
  //priors
  beta ~ normal(2, 1);
  gamma ~ normal(0.4, 0.5);
  phi_inv ~ exponential(5);

  cases ~ neg_binomial_2(y[,2], phi);
  // cases ~ poisson(y[,2]);
}

generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  array[n_days] real pred_cases = neg_binomial_2_rng(y[,2], phi);
  // array[n_days] real pred_cases = poisson_rng(y[,2]);
}

data:

{
    "N": 763,
    "cases": [3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29, 14, 4],
    "n_days": 14,
    "ts": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14],
    "t0": 0,
    "y0": [762, 1, 0]
}

inits

{
    "beta": 1.00247,
    "gamma": 1.15014,
    "phi_inv": 0.008622444
}

When compiled via cmdstan with STAN_THREADS=1, I then run with

./sir sample num_chains=4 data file=sir.data.json init=sir.init.json num_threads=4

It seems to vary a bit between runs, but when run repeatedly the most common outcome is that two chains run/print their progress, until completion, and then the process hangs. I have recreated this on two separate linux systems.

strace reports a lot of using of sched_yield from tbb, but otherwise there is nothing happening. gdb seems to implicate set_zero_all_adjoints in a backtrace, but not in any way I can track down.

Expected Output

Expect all four chains to run simultaneously and the program to not hang.

Current Version:

v4.7.0

Reproduce the hang with the supplied SIR model, cmdstan command, and STAN_THREADS configuration, comparing ode_adams with ode_bdf and ode_rk45. Start by tracing the ode_adams threading path and the reported set_zero_all_adjoints and TBB activity; done means four chains complete repeatedly without hanging.