Skip to content

Fix race condition in atolRtolFactor_() with per-thread tolerance arrays #983

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Copilot wants to merge 8 commits into
base: main
Choose a base branch
from
Draft

Fix race condition in atolRtolFactor_() with per-thread tolerance arrays #983

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5f0791b
Initial plan
Copilot Jan 28, 2026
24df908
Add OpenMP critical section to atolRtolFactor_() for thread safety
Copilot Jan 28, 2026
e15358b
Move pointer access inside critical section for complete thread safety
Copilot Jan 28, 2026
ba37bfc
Thread-safety fix complete
Copilot Jan 28, 2026
086c548
Remove accidentally added symlink
Copilot Jan 28, 2026
3a6aa7b
Refactor to use per-thread tolerance arrays instead of critical sections
Copilot Jan 29, 2026
239be04
Remove modification of shared op->ATOL and op->RTOL to prevent race c…
Copilot Jan 29, 2026
f20d540
Co-authored-by: mattfidler <514778+mattfidler@users.noreply.github.com>
Copilot Jan 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions _codeql_detected_source_root
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
.
7 changes: 4 additions & 3 deletions src/expm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,12 +310,13 @@ int meOnly(int cSub, double *yc_, double *yp_, double tp, double tf, double tcov
//' inductive linearization
//' @name rxIndLin_
//' @noRd
extern "C" int indLin(int cSub, rx_solving_options *op, double tp, double *yp_, double tf,
extern "C" int indLin(int cSub, rx_solving_options *op, rx_solving_options_ind *ind, double tp, double *yp_, double tf,
double *InfusionRate_, int *on_,
t_ME ME, t_IndF IndF){
int neq = op->neq;
double *rtol=op->rtol2;
double *atol=op->atol2;
// Use per-individual tolerance arrays if available
double *rtol = (ind != NULL && ind->rtol2 != NULL) ? ind->rtol2 : op->rtol2;
double *atol = (ind != NULL && ind->atol2 != NULL) ? ind->atol2 : op->atol2;
int maxsteps=op->mxstep;
int doIndLin=op->doIndLin;
// int indLinPerterb=10;
Expand Down
13 changes: 10 additions & 3 deletions src/par_solve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -732,7 +732,7 @@ static const char *err_msg_ls[] =
//dummy solout fn
extern "C" void solout(long int nr, double t_old, double t, double *y, int *nptr, int *irtrn){}

extern "C" int indLin(int cSub, rx_solving_options *op, double tp, double *yp_, double tf,
extern "C" int indLin(int cSub, rx_solving_options *op, rx_solving_options_ind *ind, double tp, double *yp_, double tf,
double *InfusionRate_, int *on_,
t_ME ME, t_IndF IndF);

Expand Down Expand Up @@ -760,7 +760,7 @@ static inline void solveWith1Pt(int *neq,
case 3:
if (!isSameTime(xout, xp)) {
preSolve(op, ind, xp, xout, yp);
idid = indLin(ind->id, op, xp, yp, xout, ind->InfusionRate, ind->on,
idid = indLin(ind->id, op, ind, xp, yp, xout, ind->InfusionRate, ind->on,
ME, IndF);
}
if (idid <= 0) {
Expand Down Expand Up @@ -2251,7 +2251,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
badSolveExit(i);
} else {
preSolve(op, ind, xoutp, xout, yp);
idid = indLin(solveid, op, xoutp, yp, xout, ind->InfusionRate, ind->on,
idid = indLin(solveid, op, ind, xoutp, yp, xout, ind->InfusionRate, ind->on,
ME, IndF);
xoutp=xout;
postSolve(neq, &idid, rc, &i, yp, NULL, 0, true, ind, op, rx);
Expand Down Expand Up @@ -2368,6 +2368,13 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
ctx = NULL;
return;
}

// Use per-individual tolerance arrays if available
if (ind->rtol2 != NULL && ind->atol2 != NULL) {
opt.rtol = ind->rtol2;
opt.atol = ind->atol2;
}

nx = ind->n_all_times;
BadDose = ind->BadDose;
InfusionRate = ind->InfusionRate;
Expand Down
57 changes: 49 additions & 8 deletions src/rxData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1540,6 +1540,12 @@ struct rx_globals {
// time per thread
double *timeThread = NULL;

// Per-thread tolerance arrays for thread-safe atolRtolFactor_
double *gatol2Thread = NULL;
double *grtol2Thread = NULL;
double *gssAtolThread = NULL;
double *gssRtolThread = NULL;

bool alloc=false;
};

Expand Down Expand Up @@ -1633,6 +1639,12 @@ extern "C" void _setIndPointersByThread(rx_solving_options_ind *ind) {
ind->timeThread = _globals.timeThread + rx->maxAllTimes*omp_get_thread_num();
ind->llikSave = _globals.gLlikSave + op->nLlik*rxLlikSaveSize*omp_get_thread_num();
ind->lhs = _globals.glhs+op->nlhs*omp_get_thread_num();

// Assign per-thread tolerance arrays
ind->atol2 = _globals.gatol2Thread + op->neq*omp_get_thread_num();
ind->rtol2 = _globals.grtol2Thread + op->neq*omp_get_thread_num();
ind->ssAtol = _globals.gssAtolThread + op->neq*omp_get_thread_num();
ind->ssRtol = _globals.gssRtolThread + op->neq*omp_get_thread_num();
}

extern "C" void setZeroMatrix(int which) {
Expand All @@ -1655,12 +1667,23 @@ double maxAtolRtolFactor = 0.1;
void atolRtolFactor_(double factor){
rx_solve* rx = getRxSolve_();
rx_solving_options* op = rx->op;
for (int i = op->neq;i--;){
_globals.grtol2[i] = min2(_globals.grtol2[i]*factor, maxAtolRtolFactor);
_globals.gatol2[i] = min2(_globals.gatol2[i]*factor, maxAtolRtolFactor);
}
op->ATOL = min2(op->ATOL*factor, maxAtolRtolFactor);
op->RTOL = min2(op->RTOL*factor, maxAtolRtolFactor);

// Get thread-specific tolerance arrays
int threadId = omp_get_thread_num();
double *atol2 = _globals.gatol2Thread + op->neq*threadId;
double *rtol2 = _globals.grtol2Thread + op->neq*threadId;
double *ssAtol = _globals.gssAtolThread + op->neq*threadId;
double *ssRtol = _globals.gssRtolThread + op->neq*threadId;

// Modify thread-specific tolerances
for (int i = op->neq; i--;){
atol2[i] = min2(atol2[i]*factor, maxAtolRtolFactor);
rtol2[i] = min2(rtol2[i]*factor, maxAtolRtolFactor);
ssAtol[i] = min2(ssAtol[i]*factor, maxAtolRtolFactor);
ssRtol[i] = min2(ssRtol[i]*factor, maxAtolRtolFactor);
}
// Note: op->ATOL and op->RTOL are not modified to avoid race conditions
// The per-thread tolerance arrays are used by the solvers
}

extern "C" double * getAol(int n, double atol){
Expand Down Expand Up @@ -5256,10 +5279,11 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
int n8 = rx->maxAllTimes*op->cores;
int n9 = (op->numLinSens+op->numLin)*op->cores;
int n10 = (op->neq)*op->cores;
int n11 = 4*op->neq*op->cores; // Per-thread tolerance arrays (atol2, rtol2, ssAtol, ssRtol)
int nlin = (rx->linB)* 7* rx->nsub * rx->nsim;
if (_globals.gsolve != NULL) free(_globals.gsolve);
_globals.gsolve = (double*)calloc(nlin+n0+3*nsave+n2+ n4+n5_c+n6+ n7 + n8 +
n9 + n10 +
n9 + n10 + n11 +
5*op->neq + 4*n3a_c + nllik_c,
sizeof(double));// [n0]
#ifdef rxSolveT
Expand Down Expand Up @@ -5298,7 +5322,11 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
_globals.gIndSim = _globals.gCurDoseS + n3a_c;// [n7]
_globals.gLinSave = _globals.gIndSim + n7; // [n9]
_globals.gLinDummy = _globals.gLinSave + n9; // [n10]
_globals.timeThread = _globals.gLinDummy + n10;
_globals.timeThread = _globals.gLinDummy + n10; // [n8]
_globals.gatol2Thread = _globals.timeThread + n8; // [op->neq * op->cores]
_globals.grtol2Thread = _globals.gatol2Thread + op->neq*op->cores; // [op->neq * op->cores]
_globals.gssAtolThread = _globals.grtol2Thread + op->neq*op->cores; // [op->neq * op->cores]
_globals.gssRtolThread = _globals.gssAtolThread + op->neq*op->cores; // [op->neq * op->cores]
std::fill_n(rx->ypNA, op->neq, NA_REAL);

std::fill_n(&_globals.gatol2[0],op->neq, atolNV[0]);
Expand All @@ -5310,6 +5338,19 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
std::fill_n(&_globals.gssRtol[0],op->neq, rtolNVss[0]);
std::copy(atolNVss.begin(), atolNVss.begin() + min2(op->neq, atolNVss.size()), &_globals.gssAtol[0]);
std::copy(rtolNVss.begin(), rtolNVss.begin() + min2(op->neq, rtolNVss.size()), &_globals.gssRtol[0]);

// Initialize per-thread tolerance arrays from global values
for (int core = 0; core < op->cores; core++) {
std::copy(&_globals.gatol2[0], &_globals.gatol2[0] + op->neq,
_globals.gatol2Thread + core*op->neq);
std::copy(&_globals.grtol2[0], &_globals.grtol2[0] + op->neq,
_globals.grtol2Thread + core*op->neq);
std::copy(&_globals.gssAtol[0], &_globals.gssAtol[0] + op->neq,
_globals.gssAtolThread + core*op->neq);
std::copy(&_globals.gssRtol[0], &_globals.gssRtol[0] + op->neq,
_globals.gssRtolThread + core*op->neq);
}

op->atol2 = &_globals.gatol2[0];
op->rtol2 = &_globals.grtol2[0];
op->ssAtol = _globals.gssAtol;
Expand Down
Loading