From e645f29b4be554e3a038fe53ac75b081b7e0a943 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 17:47:18 -0500
Subject: [PATCH 01/13] Replace shared op->badSolve reads with per-individual
 localBadSolve for thread safety

In parallel ODE solving, multiple threads reading/writing the shared
op->badSolve flag creates a race condition where one thread's failed
solve can incorrectly prevent other threads from continuing. This
introduces a localBadSolve variable in ind_indLin0, ind_liblsoda0,
and ind_lsoda0 that tracks bad-solve state per-individual, and removes
the shared op->badSolve check from iniSubject in par_solve.h.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/par_solve.cpp | 100 +++++++++++++++++++++++++++-------------------
 src/par_solve.h   |   5 ++-
 2 files changed, 62 insertions(+), 43 deletions(-)

diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index b23e6b052..c1469e868 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -2229,6 +2229,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
   double *yp;
   inits = op->inits;
   int idid = 0;
+  int localBadSolve = 0;
   ind = &(rx->subjects[neq[1]]);
   if (!iniSubject(neq[1], 0, ind, op, rx, u_inis)) return;
   nx = ind->n_all_times;
@@ -2249,6 +2250,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
         *rc = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         preSolve(op, ind, xoutp, xout, yp);
         idid = indLin(solveid, op, xoutp, yp, xout, ind->InfusionRate, ind->on,
@@ -2258,7 +2260,7 @@ extern "C" void ind_indLin0(rx_solve *rx, rx_solving_options *op, int solveid,
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &idid, u_inis);
@@ -2356,6 +2358,7 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
   double *yp;
   int neqOde = *neq - op->numLin - op->numLinSens;
   inits = op->inits;
+  int localBadSolve = 0;
   struct lsoda_context_t * ctx = lsoda_create_ctx();
   ctx->function = (_lsoda_f)dydt_liblsoda;
   ctx->data = neq;
@@ -2386,45 +2389,51 @@ extern "C" void ind_liblsoda0(rx_solve *rx, rx_solving_options *op, struct lsoda
         *rc = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         // REprintf("xp: %f xout: %f\n", xp, xout);
         if (handleExtraDose(neq, BadDose, InfusionRate, ind->dose, yp, xout,
                             xp, ind->id, &i, nx, &(ctx->state), op, ind, u_inis, ctx)) {
-          if (!isSameTime(ind->extraDoseNewXout, xp)) {
+          if (!localBadSolve && !isSameTime(ind->extraDoseNewXout, xp)) {
             preSolve(op, ind, xp, ind->extraDoseNewXout, yp);
             lsoda(ctx, yp, &xp, ind->extraDoseNewXout);
             copyLinCmt(neq, ind, op, yp);
             postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+            if (*rc < 0) localBadSolve = 1;
           }
-          int idx = ind->idx;
-          int ixds= ind->ixds;
-          int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
-          ind->idx = -1-trueIdx;
-          handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
-                      BadDose, InfusionRate, ind->dose, yp, xout, neq[1], ind);
-          ctx->state = 1;
-          ind->idx = idx;
-          ind->ixds = ixds;
-          ind->idxExtra++;
-          if (!isSameTime(xout, ind->extraDoseNewXout)) {
-            preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
-            lsoda(ctx,yp, &ind->extraDoseNewXout, xout);
-            copyLinCmt(neq, ind, op, yp);
-            postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+          if (!localBadSolve) {
+            int idx = ind->idx;
+            int ixds= ind->ixds;
+            int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
+            ind->idx = -1-trueIdx;
+            handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
+                        BadDose, InfusionRate, ind->dose, yp, xout, neq[1], ind);
+            ctx->state = 1;
+            ind->idx = idx;
+            ind->ixds = ixds;
+            ind->idxExtra++;
+            if (!isSameTime(xout, ind->extraDoseNewXout)) {
+              preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
+              lsoda(ctx,yp, &ind->extraDoseNewXout, xout);
+              copyLinCmt(neq, ind, op, yp);
+              postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+              if (*rc < 0) localBadSolve = 1;
+            }
+            xp =  ind->extraDoseNewXout;
           }
-          xp =  ind->extraDoseNewXout;
         }
-        if (!isSameTime(xout, xp)) {
+        if (!localBadSolve && !isSameTime(xout, xp)) {
           preSolve(op, ind, xp, xout, yp);
           lsoda(ctx, yp, &xp, xout);
           copyLinCmt(neq, ind, op, yp);
           postSolve(neq, &(ctx->state), rc, &i, yp, NULL, 0, false, ind, op, rx);
+          if (*rc < 0) localBadSolve = 1;
         }
         xp = xout;
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &(ctx->state), u_inis);
@@ -2816,6 +2825,7 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
 
   double xp = getAllTimes(ind, 0);
   double xout;
+  int localBadSolve = 0;
 
   if (!iniSubject(neq[1], 0, ind, op, rx, u_inis)) return;
   ind->solvedIdx = 0;
@@ -2829,10 +2839,11 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
         ind->rc[0] = -1000;
         // Bad Solve => NA
         badSolveExit(i);
+        localBadSolve = 1;
       } else {
         if (handleExtraDose(neq, ind->BadDose, ind->InfusionRate, ind->dose, yp, xout,
                             xp, ind->id, &i, ind->n_all_times, &istate, op, ind, u_inis, ctx)) {
-          if (!isSameTime(ind->extraDoseNewXout, xp)) {
+          if (!localBadSolve && !isSameTime(ind->extraDoseNewXout, xp)) {
             preSolve(op, ind, xp, ind->extraDoseNewXout, yp);
             neq[0] = op->neq - op->numLin - op->numLinSens;
             F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &xp, &ind->extraDoseNewXout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
@@ -2840,29 +2851,33 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
             neq[0] = op->neq;
             copyLinCmt(neq, ind, op, yp);
             postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+            if (*(ind->rc) < 0) localBadSolve = 1;
           }
-          int idx = ind->idx;
-          int ixds = ind->ixds;
-          int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
-          ind->idx = -1-trueIdx;
-          handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
-                      ind->BadDose, ind->InfusionRate, ind->dose, yp, xout, neq[1], ind);
-          istate = 1;
-          ind->ixds = ixds; // This is a fake dose, real dose stays in place
-          ind->idx = idx;
-          ind->idxExtra++;
-          if (!isSameTime(xout, ind->extraDoseNewXout)) {
-            preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
-            neq[0] = op->neq - op->numLin - op->numLinSens;
-            F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &ind->extraDoseNewXout, &xout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
-                             &istate, &giopt, rwork, &lrw, iwork, &liw, jdum, &jt);
-            neq[0] = op->neq;
-            copyLinCmt(neq, ind, op, yp);
-            postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+          if (!localBadSolve) {
+            int idx = ind->idx;
+            int ixds = ind->ixds;
+            int trueIdx = ind->extraDoseTimeIdx[ind->idxExtra];
+            ind->idx = -1-trueIdx;
+            handle_evid(ind->extraDoseEvid[trueIdx], neq[0],
+                        ind->BadDose, ind->InfusionRate, ind->dose, yp, xout, neq[1], ind);
+            istate = 1;
+            ind->ixds = ixds; // This is a fake dose, real dose stays in place
+            ind->idx = idx;
+            ind->idxExtra++;
+            if (!isSameTime(xout, ind->extraDoseNewXout)) {
+              preSolve(op, ind, ind->extraDoseNewXout, xout, yp);
+              neq[0] = op->neq - op->numLin - op->numLinSens;
+              F77_CALL(dlsoda)(dydt_lsoda, neq, yp, &ind->extraDoseNewXout, &xout, &gitol, &(op->RTOL), &(op->ATOL), &gitask,
+                               &istate, &giopt, rwork, &lrw, iwork, &liw, jdum, &jt);
+              neq[0] = op->neq;
+              copyLinCmt(neq, ind, op, yp);
+              postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+              if (*(ind->rc) < 0) localBadSolve = 1;
+            }
+            xp =  ind->extraDoseNewXout;
           }
-          xp =  ind->extraDoseNewXout;
         }
-        if (!isSameTime(xout, xp)) {
+        if (!localBadSolve && !isSameTime(xout, xp)) {
           preSolve(op, ind, xp, xout, yp);
           neq[0] = op->neq - op->numLin - op->numLinSens;
           F77_CALL(dlsoda)(dydt_lsoda, neq, yp,
@@ -2874,13 +2889,14 @@ extern "C" void ind_lsoda0(rx_solve *rx, rx_solving_options *op, int solveid, in
           neq[0] = op->neq;
           copyLinCmt(neq, ind, op, yp);
           postSolve(neq, &istate, ind->rc, &i, yp, err_msg_ls, 7, true, ind, op, rx);
+          if (*(ind->rc) < 0) localBadSolve = 1;
         }
         xp = xout;
         //dadt_counter = 0;
       }
     }
     ind->_newind = 2;
-    if (!op->badSolve){
+    if (!localBadSolve){
       ind->idx = i;
       if (getEvid(ind, ind->ix[i]) == 3) {
         handleEvid3(ind, op, rx, neq, &xp, &xout,  yp, &(istate), u_inis);
diff --git a/src/par_solve.h b/src/par_solve.h
index dad810589..a9363a2fb 100644
--- a/src/par_solve.h
+++ b/src/par_solve.h
@@ -78,7 +78,10 @@ extern "C" {
 			ind->solved = -1;
 		}
     sortInd(ind);
-    if (op->badSolve) return 0;
+    // Note: previously checked op->badSolve here, but that is a shared
+    // global flag — in parallel mode another thread's failed solve would
+    // prevent THIS individual from initializing.  The bad-solve state is
+    // now handled per-individual via localBadSolve / *rc in the caller.
 		ind->ixds=ind->idx=0;
     if (ncmt) ind->pendingDosesN[0] = 0;
 		return 1;

From dd0035b7b973ccced2a34e9ef1e5d38fad4f6b86 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:04:12 -0500
Subject: [PATCH 02/13] Use atomic writes for op->badSolve and op->naTime in
 parallel regions

These shared fields are written from multiple threads during parallel
ODE solving. Without synchronization this is undefined behavior in
C/C++. Use #pragma omp atomic write to ensure thread-safe stores.

Also refactor push*Dose allocation failure handling to use local flags
instead of writing op->badSolve inside critical sections then reading
it outside.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 inst/include/rxode2parseGetTime.h    | 40 +++++++++-------
 inst/include/rxode2parseHandleEvid.h | 72 ++++++++++++++--------------
 src/par_solve.cpp                    |  1 +
 3 files changed, 62 insertions(+), 51 deletions(-)

diff --git a/inst/include/rxode2parseGetTime.h b/inst/include/rxode2parseGetTime.h
index 7fd3b6dfc..68c9fe623 100644
--- a/inst/include/rxode2parseGetTime.h
+++ b/inst/include/rxode2parseGetTime.h
@@ -38,10 +38,12 @@ static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double
   }
   double ret = LAG(id, cmt, time);
   if (ISNA(ret)) {
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 1 + 10*cmt;
-    }
+    int newBadSolve = 1;
+    int newNaTime = 1 + 10*cmt;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+#pragma omp atomic write
+    op->naTime = newNaTime;
   }
   return ret;
 }
@@ -51,10 +53,12 @@ static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, doubl
   returnBadTime(t);
   double ret = RATE(id, cmt, dose, t);
   if (ISNA(ret)){
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 2 + 10*cmt;
-    }
+    int newBadSolve = 1;
+    int newNaTime = 2 + 10*cmt;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+#pragma omp atomic write
+    op->naTime = newNaTime;
   }
   return ret;
 }
@@ -65,10 +69,12 @@ static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double
   if (ISNA(t)) return t;
   double ret = DUR(id, cmt, dose, t);
   if (ISNA(ret)){
-    op->badSolve=1;
-    if (op->naTime == 0) {
-      op->naTime = 3 + 10*cmt;
-    }
+    int newBadSolve = 1;
+    int newNaTime = 3 + 10*cmt;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+#pragma omp atomic write
+    op->naTime = newNaTime;
   }
   return ret;
 }
@@ -340,10 +346,12 @@ static inline double handleInfusionItem(int idx, rx_solve *rx, rx_solving_option
     ind->idx = oIdx;
     if (ISNA(f)){
       rx_solving_options *op = &op_global;
-      op->badSolve=1;
-      if (op->naTime == 0) {
-        op->naTime = 4 + 10*ind->cmt;
-      }
+      int newBadSolve = 1;
+      int newNaTime = 4 + 10*ind->cmt;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+#pragma omp atomic write
+      op->naTime = newNaTime;
     }
     double durOld = (getAllTimes(ind, ind->idose[infEidx]) -
                      getAllTimes(ind, ind->idose[infBidx]));
diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h
index f90ed9143..e4679a707 100644
--- a/inst/include/rxode2parseHandleEvid.h
+++ b/inst/include/rxode2parseHandleEvid.h
@@ -285,8 +285,12 @@ static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt,
   double ret = AMT(id, cmt, dose, t, y);
   if (ISNA(ret)){
     rx_solving_options *op = &op_global;
-    op->badSolve=1;
-    op->naTime = 5 + 10*cmt;
+    int newBadSolve = 1;
+    int newNaTime = 5 + 10*cmt;
+#pragma omp atomic write
+    op->badSolve = newBadSolve;
+#pragma omp atomic write
+    op->naTime = newNaTime;
   }
   return ret;
 }
@@ -309,20 +313,22 @@ static inline int pushIgnoredDose(int doseIdx, rx_solving_options_ind *ind) {
   }
   if (ind->ignoredDosesN[0]+1 >= ind->ignoredDosesAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->ignoredDoses, (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = 1;
-        // return 0;
+        allocFailed = 1;
       } else {
         ind->ignoredDoses = tmpI;
         ind->ignoredDosesAllocN[0] = (ind->ignoredDosesN[0]+1+EVID_EXTRA_SIZE);
         re = 1;
       }
     }
-    if (op->badSolve) {
-      return 0; // don't continue if we have a bad solve.
+    if (allocFailed) {
+#pragma omp atomic write
+      op->badSolve = 1;
+      return 0;
     }
   }
   ind->ignoredDoses[ind->ignoredDosesN[0]] = doseIdx;
@@ -334,19 +340,21 @@ static inline int pushPendingDose(int doseIdx, rx_solving_options_ind *ind) {
   int re = 0;
   if (ind->pendingDosesN[0]+1 >= ind->pendingDosesAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->pendingDoses, (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = 1;
-        //return 0;
+        allocFailed = 1;
       } else {
         ind->pendingDoses = tmpI;
         ind->pendingDosesAllocN[0] = (ind->pendingDosesN[0]+1+EVID_EXTRA_SIZE);
         re = 1;
       }
     }
-    if (op->badSolve == 1)  {
+    if (allocFailed) {
+#pragma omp atomic write
+      op->badSolve = 1;
       return 0;
     }
   }
@@ -361,32 +369,29 @@ static inline int pushDosingEvent(double time, double amt, int evid,
   int re = 0;
   if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;  // 0=ok, 1=partial alloc, -1=first alloc failed
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = -1;
-        // return 0;
+        allocFailed = -1;
       } else {
         ind->extraDoseTimeIdx = tmpI;
 
         tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
         if (tmpI == NULL) {
-          op->badSolve = 1;
-          // return 1;
+          allocFailed = 1;
         } else {
           ind->extraDoseEvid = tmpI;
           double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
           if (tmpD == NULL) {
-            op->badSolve = 1;
-            //return 1;
+            allocFailed = 1;
           } else {
             ind->extraDoseTime = tmpD;
 
             tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
             if (tmpD == NULL) {
-              op->badSolve = 1;
-              //return 1;
+              allocFailed = 1;
             } else {
               ind->extraDoseDose = tmpD;
 
@@ -397,11 +402,11 @@ static inline int pushDosingEvent(double time, double amt, int evid,
         }
       }
     }
-    if (op->badSolve == 1) {
-      return 1;
-    } else if (op->badSolve == -1) {
-      op->badSolve = 1; // set to bad solve.
-      return 0;
+    if (allocFailed != 0) {
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return (allocFailed == 1) ? 1 : 0;
     }
   }
   ind->extraDoseTimeIdx[ind->extraDoseN[0]] = ind->extraDoseN[0];
@@ -419,33 +424,30 @@ static inline int pushUniqueDosingEvent(double time, double amt, int evid,
   int re = 0;
   if (ind->extraDoseN[0]+1 >= ind->extraDoseAllocN[0]) {
     rx_solving_options *op = &op_global;
+    int allocFailed = 0;  // 0=ok, 1=partial alloc, -1=first alloc failed
 #pragma omp critical
     {
       int *tmpI = (int*)realloc(ind->extraDoseTimeIdx, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
       if (tmpI == NULL) {
-        op->badSolve = -1;
-        // return 0;
+        allocFailed = -1;
       }  else {
         ind->extraDoseTimeIdx = tmpI;
 
         tmpI = (int*)realloc(ind->extraDoseEvid, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(int));
         if (tmpI == NULL) {
-          op->badSolve = 1;
-          // return 1;
+          allocFailed = 1;
         } else {
           ind->extraDoseEvid = tmpI;
 
           double * tmpD = (double*)realloc(ind->extraDoseTime, (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
           if (tmpD == NULL) {
-            op->badSolve = 1;
-            //return 1;
+            allocFailed = 1;
           } else {
             ind->extraDoseTime = tmpD;
 
             tmpD = (double*)realloc(ind->extraDoseDose,  (ind->extraDoseN[0]+1+EVID_EXTRA_SIZE)*sizeof(double));
             if (tmpD == NULL) {
-              op->badSolve = 1;
-              // return 1;
+              allocFailed = 1;
             } else {
               ind->extraDoseDose = tmpD;
 
@@ -455,11 +457,11 @@ static inline int pushUniqueDosingEvent(double time, double amt, int evid,
         }
       }
     }
-    if (op->badSolve == 1) {
-      return 1; // don't continue if we have a bad solve.
-    } else if (op->badSolve == -1) {
-      op->badSolve = 1; // set to bad solve.
-      return 0;
+    if (allocFailed != 0) {
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return (allocFailed == 1) ? 1 : 0;
     }
     re = 1;
   }
diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index c1469e868..9f35c2600 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -37,6 +37,7 @@ extern "C" {
 #define badSolveExit(i) for (int j = op->neq*(ind->n_all_times); j--;){ \
     ind->solve[j] = NA_REAL;                                           \
   }                                                                     \
+  _Pragma("omp atomic write")                                           \
   op->badSolve = 1;                                                     \
   i = ind->n_all_times-1; // Get out of here!
 // Yay easy parallel support

From c46d6174cb907a09e7b593cb4af96f0ad771e2fb Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:05:18 -0500
Subject: [PATCH 03/13] Replace Rf_errorcall in _getDur with NA_REAL return and
 badSolve flag

Rf_errorcall is not thread-safe and will crash or corrupt R state when
called from OpenMP worker threads. Replace with the same NA_REAL
return pattern already used for backward==2 case, plus setting the
badSolve flag atomically to signal the error after the parallel region.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/handle_evid.cpp | 27 +++++++++++++++++++++++----
 1 file changed, 23 insertions(+), 4 deletions(-)

diff --git a/src/handle_evid.cpp b/src/handle_evid.cpp
index 2850f5d26..585f48b90 100644
--- a/src/handle_evid.cpp
+++ b/src/handle_evid.cpp
@@ -31,24 +31,39 @@ extern "C" void handleTlast(double *time, rx_solving_options_ind *ind) {
 }
 
 // Linear compartment models/functions
+// Note: Rf_errorcall is not thread-safe and cannot be called from
+// within OpenMP parallel regions.  Return NA_REAL and set badSolve
+// instead so the error is handled after the parallel region completes.
 extern "C" double _getDur(int l, rx_solving_options_ind *ind, int backward, unsigned int *p) {
   double dose = getDoseNumber(ind, l);
   if (backward==1 && l != 0){
     if (l <= 0) {
-      Rf_errorcall(R_NilValue, _("could not find a start to the infusion #1"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     p[0] = l-1;
     while (p[0] > 0 && getDoseNumber(ind, p[0]) != -dose){
       p[0]--;
     }
     if (getDoseNumber(ind, p[0]) != -dose){
-      Rf_errorcall(R_NilValue, _("could not find a start to the infusion #2"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     return getAllTimes(ind, ind->idose[l]) - getAllTimes(ind, ind->idose[p[0]]);
   } else {
     if (l >= ind->ndoses) {
       if (backward==2) return(NA_REAL);
-      Rf_errorcall(R_NilValue, _("could not find an end to the infusion"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     p[0] = l+1;
     while (p[0] < ind->ndoses && getDoseNumber(ind, p[0]) != -dose){
@@ -56,7 +71,11 @@ extern "C" double _getDur(int l, rx_solving_options_ind *ind, int backward, unsi
     }
     if (getDoseNumber(ind, p[0]) != -dose){
       if (backward==2) return(NA_REAL);
-      Rf_errorcall(R_NilValue, _("could not find an end to the infusion"));
+      rx_solving_options *op = &op_global;
+      int newBadSolve = 1;
+#pragma omp atomic write
+      op->badSolve = newBadSolve;
+      return NA_REAL;
     }
     return getAllTimes(ind, ind->idose[p[0]]) - getAllTimes(ind, ind->idose[l]);
   }

From 39d9f445fa7daf1e9391687285b85a8b3cd2d298 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:07:00 -0500
Subject: [PATCH 04/13] Mark abort flag as volatile in parallel solve loops

The abort flag is read by all threads and written by the master thread
without synchronization. Using volatile prevents the compiler from
caching the value in a register, ensuring threads see updates promptly.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/par_solve.cpp | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/par_solve.cpp b/src/par_solve.cpp
index 9f35c2600..717c0a2f6 100644
--- a/src/par_solve.cpp
+++ b/src/par_solve.cpp
@@ -2304,7 +2304,8 @@ extern "C" void par_indLin(rx_solve *rx){
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   // FIXME parallel
   uint32_t seed0 = getRxSeed1(1);
   for (int solveid = 0; solveid < nsim*nsub; solveid++){
@@ -2504,7 +2505,8 @@ extern "C" void par_linCmt(rx_solve *rx) {
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(cores)
@@ -2581,7 +2583,8 @@ extern "C" void par_liblsodaR(rx_solve *rx) {
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(cores)
@@ -2656,7 +2659,8 @@ extern "C" void par_liblsoda(rx_solve *rx){
   // Breaking of of loop ideas came from http://www.thinkingparallel.com/2007/06/29/breaking-out-of-loops-in-openmp/
   // http://permalink.gmane.org/gmane.comp.lang.r.devel/27627
   // It was buggy due to Rprint.  Use REprint instead since Rprint calls the interrupt every so often....
-  int abort = 0;
+  // volatile ensures reads/writes are not cached in registers across threads
+  volatile int abort = 0;
   uint32_t seed0 = getRxSeed1(cores);
 #ifdef _OPENMP
 #pragma omp parallel for num_threads(op->cores)

From caf74981836dd301c961f1f56de9e92ac04c7f3e Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:07:34 -0500
Subject: [PATCH 05/13] Fix rx_get_thread off-by-one and fallback to prevent
 thread collisions

rx_get_thread(mx) used <= which allowed index==mx (one past the end
of per-thread arrays sized to mx). Changed to < for correct bounds.

Also changed the overflow fallback from returning 0 (which collides
with thread 0's slot) to clamping to mx-1.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/rxomp.h | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/rxomp.h b/src/rxomp.h
index 308f3630a..6d67bb410 100644
--- a/src/rxomp.h
+++ b/src/rxomp.h
@@ -8,8 +8,10 @@
 static inline int rx_get_thread(int mx) {
   int tn = omp_get_thread_num();
   if (tn < 0) return 0;
-  if (tn <= mx) return tn;
-  return 0;
+  if (tn < mx) return tn;
+  // Clamp to last valid index rather than collapsing to 0,
+  // which would cause two threads to share the same slot
+  return (mx > 0) ? mx - 1 : 0;
 }
 
 #else

From dd70fea659f7bc53b94323872ce86c0d44959221 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:14:36 -0500
Subject: [PATCH 06/13] Refactor linCmtB to per-thread vector for thread safety

linCmtB used a single global stan::math::linCmtStan object and
associated Eigen matrices, making it unsafe for parallel execution.
Refactored to use the same per-thread vector pattern as linCmtA:
- Created linB_t struct holding all per-thread state
- Changed __linCmtB from a single object to std::vector<linB_t>
- Added ensureLinCmtB(nCores) called during solve setup
- linCmtB() now indexes __linCmtB by omp_get_thread_num()
- R-facing utility functions use element 0 (main thread)

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/linCmt.cpp | 79 +++++++++++++++++++++++++++-----------------------
 src/rxData.cpp |  2 ++
 2 files changed, 44 insertions(+), 37 deletions(-)

diff --git a/src/linCmt.cpp b/src/linCmt.cpp
index cd0336d49..743c7d873 100644
--- a/src/linCmt.cpp
+++ b/src/linCmt.cpp
@@ -40,36 +40,40 @@ extern "C" void ensureLinCmtA(int nCores) {
 }
 
 // Global linear compartment B model object
-// Since this cannot be threaded, this is not a vector
-// object.  This is created once to reduce memory allocation
-// and deallocation time.
-stan::math::linCmtStan __linCmtB(0, 0, 0, true, 0, 0);
-// Maximum size can be 2*ncmt + 1;
-double __linCmtBdata[14];
-int __linCmtBnumSens=0;
+// Refactored to per-thread vector for thread safety, matching linCmtA pattern.
+typedef struct {
+  stan::math::linCmtStan lc;
+  double data[14];
+  int numSens;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> fx;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> yp;
+  Eigen::Matrix<double, Eigen::Dynamic, 2> g;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> J;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Js;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> Jg;
+} linB_t;
+
+std::vector<linB_t> __linCmtB;
+
+extern "C" void ensureLinCmtB(int nCores) {
+  if ((int)__linCmtB.size() < nCores) {
+    __linCmtB.resize(nCores);
+  }
+}
 
 #define linCmtBaddrTheta 0
 #define linCmtBaddrThetaSens 1
-static inline double * getLinCmtDoubleAddr(int type) {
+static inline double * getLinCmtDoubleAddr(linB_t &lcb, int type) {
   switch (type) {
   case linCmtBaddrTheta: // max 7
-    return __linCmtBdata;
+    return lcb.data;
   case linCmtBaddrThetaSens:  // max 7
-    return &__linCmtBdata[7];
+    return &lcb.data[7];
   // note fx needs cannot be a Map for use in stan :(
   }
   return NULL;
 }
 
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtBfx;
-
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtByp;
-Eigen::Matrix<double, Eigen::Dynamic, 2> __linCmtBg;
-
-Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> __linCmtBJ;
-Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> __linCmtBJs;
-Eigen::Matrix<double, Eigen::Dynamic, 1> __linCmtBJg;
-
 // [[Rcpp::export]]
 RObject linCmtModelDouble(double dt,
                           double p1, double v1, double p2,
@@ -383,16 +387,16 @@ extern "C" double linCmtA(rx_solve *rx, int id,
 }
 
 extern "C" double linCmtScaleInitPar(int which) {
-  return __linCmtB.initPar(which);
+  return __linCmtB[0].lc.initPar(which);
 }
 
 extern "C" double linCmtScaleInitN() {
-  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = __linCmtB.initPar();
+  Eigen::Matrix<double, Eigen::Dynamic, 1> theta = __linCmtB[0].lc.initPar();
   return theta.size();
 }
 
 extern "C" int linCmtZeroJac(int i) {
-  return __linCmtB.parDepV1(i);
+  return __linCmtB[0].lc.parDepV1(i);
 }
 
 
@@ -481,14 +485,15 @@ extern "C" double linCmtB(rx_solve *rx, int id,
                           double p4, double p5,
                           // Oral parameters
                           double ka) {
-#define fx        __linCmtBfx
-#define Jg        __linCmtBJg
-#define lc        __linCmtB
-#define AlastA    __linCmtBAlastA
-#define J         __linCmtBJ
-#define Js        __linCmtBJs
-#define yp        __linCmtByp
-#define g         __linCmtBg
+  // Per-thread linCmtB state (matching linCmtA pattern for thread safety)
+  linB_t &lcb = __linCmtB[omp_get_thread_num()];
+#define fx        lcb.fx
+#define Jg        lcb.Jg
+#define lc        lcb.lc
+#define J         lcb.J
+#define Js        lcb.Js
+#define yp        lcb.yp
+#define g         lcb.g
   rx_solving_options_ind *ind = &(rx->subjects[id]);
   rx_solving_options *op = rx->op;
   int idx = ind->idx;
@@ -523,8 +528,8 @@ extern "C" double linCmtB(rx_solve *rx, int id,
     // NA fill and resize
     J = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Constant(ncmt + oral0, npars, NA_REAL);
 
-    __linCmtBnumSens = lc.numSens();
-    Js = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(ncmt+oral0, __linCmtBnumSens);//(ncmt + oral0, 2*ncmt + oral0);
+    lcb.numSens = lc.numSens();
+    Js = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(ncmt+oral0, lcb.numSens);//(ncmt + oral0, 2*ncmt + oral0);
     // thetaSens.resize(numSens);
 
     // AlastA.resize(ncmt + oral0);
@@ -542,7 +547,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
   lc.setId(id);
 
   Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> >
-    theta(getLinCmtDoubleAddr(linCmtBaddrTheta), lc.getNpars());
+    theta(getLinCmtDoubleAddr(lcb, linCmtBaddrTheta), lc.getNpars());
 
   int sw = ncmt + 10*oral0;
   switch (sw) {
@@ -555,7 +560,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
   }
 
   Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 1> >
-    thetaSens(getLinCmtDoubleAddr(linCmtBaddrThetaSens), __linCmtBnumSens);
+    thetaSens(getLinCmtDoubleAddr(lcb, linCmtBaddrThetaSens), lcb.numSens);
 
   lc.sensTheta(theta, thetaSens, rx->sensType == 3, rx->linCmtScale);
   if (ind->linSS == linCmtSsInf) {
@@ -675,13 +680,13 @@ extern "C" double linCmtB(rx_solve *rx, int id,
       }
     }
   }
-  lc.getJacCp(__linCmtBJ, fx, theta, Jg);
+  lc.getJacCp(lcb.J, fx, theta, Jg);
   return lc.adjustF(fx, theta, ind->linCmtHV);
 #undef fx
 #undef J
 #undef Jg
 #undef lc
-#undef theta
-#undef AlastA
+#undef Js
 #undef yp
+#undef g
 }
diff --git a/src/rxData.cpp b/src/rxData.cpp
index 2b9aabc88..1c3f7d84d 100644
--- a/src/rxData.cpp
+++ b/src/rxData.cpp
@@ -43,6 +43,7 @@ using namespace arma;
 
 extern "C" void seedEng(int ncores);
 extern "C" void ensureLinCmtA(int nCores);
+extern "C" void ensureLinCmtB(int nCores);
 
 #include "cbindThetaOmega.h"
 #include "../inst/include/rxode2parseHandleEvid.h"
@@ -4994,6 +4995,7 @@ SEXP rxSolve_(const RObject &obj, const List &rxControl,
     }
     seedEng((int)(op->cores));
     ensureLinCmtA((int)op->cores);
+    ensureLinCmtB((int)op->cores);
 
     // Now set up events and parameters
     RObject par0 = params;

From 4422d42521a08610651e9daecb3b5c4e541fc699 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 20:26:04 -0500
Subject: [PATCH 07/13] Refactor DOP853 solver to use context struct for thread
 safety

Move all 16 file-scope static variables in dop853.c into a dop853_ctx_t
struct defined in dop853.h.  Each call to dop853() now allocates the
context on the stack, making the solver reentrant and safe for concurrent
use from multiple threads.

Changes:
- Add dop853_ctx_t struct to dop853.h containing nfcn/nstep/naccpt/nrejct,
  hout/xold/xout, nrds/indir, yy1/k1-k10, and rcont1-rcont8
- dopcor() takes dop853_ctx_t* as first parameter, uses local aliases
  for array pointers to minimize diff in the computation body
- contd8() becomes static and takes dop853_ctx_t* (never called externally
  since all callers use iout=0)
- Remove accessor functions (nfcnRead, nstepRead, etc.) which were never
  called from outside dop853.c
- Public dop853() API is unchanged - callers need no modifications

This enables future parallelization of par_dop() with OpenMP, which was
previously blocked by the comment "This part CAN be parallelized, if dop
is thread safe..."

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/dop853.c | 370 ++++++++++++++++++++++++---------------------------
 src/dop853.h |  23 ++--
 2 files changed, 183 insertions(+), 210 deletions(-)

diff --git a/src/dop853.c b/src/dop853.c
index 4716f5d89..f735d933e 100644
--- a/src/dop853.c
+++ b/src/dop853.c
@@ -15,54 +15,9 @@
 
 #define _(String) (String)
 
-static long int     nfcn, nstep, naccpt, nrejct;
-static double       hout, xold, xout;
-static int nrds, *indir;
-static double       *yy1, *k1, *k2, *k3, *k4, *k5, *k6, *k7, *k8, *k9, *k10;
-static double       *rcont1, *rcont2, *rcont3, *rcont4;
-static double       *rcont5, *rcont6, *rcont7, *rcont8;
-
-
-long int nfcnRead (void)
-{
-  return nfcn;
-
-} /* nfcnRead */
-
-
-long int nstepRead (void)
-{
-  return nstep;
-
-} /* stepRead */
-
-
-long int naccptRead (void)
-{
-  return naccpt;
-
-} /* naccptRead */
-
-
-long int nrejctRead (void)
-{
-  return nrejct;
-
-} /* nrejct */
-
-
-double hRead (void)
-{
-  return hout;
-
-} /* hRead */
-
-
-double xRead (void)
-{
-  return xout;
-
-} /* xRead */
+/* All former static variables are now in dop853_ctx_t (see dop853.h).
+   This makes the solver reentrant and thread-safe: each call to dop853()
+   uses its own local context on the stack. */
 
 
 static double sign (double a, double b)
@@ -162,8 +117,36 @@ static double hinit (int *nptr, FcnEqDiff fcn, double x, double* y,
 } /* hinit */
 
 
+/* dense output function */
+static double contd8 (dop853_ctx_t *ctx, int ii, double x)
+{
+  int i;
+  double       s, s1;
+
+  i = INT_MAX;
+
+  if (!ctx->indir)
+    i = ii;
+  else
+    i = ctx->indir[ii];
+
+  if (i == INT_MAX)
+    {
+      Rprintf (_("no dense output available for %uth component"), ii);
+      return 0.0;
+    }
+
+  s = (x - ctx->xold) / ctx->hout;
+  s1 = 1.0 - s;
+
+  return ctx->rcont1[i]+s*(ctx->rcont2[i]+s1*(ctx->rcont3[i]+s*(ctx->rcont4[i]+s1*(ctx->rcont5[i]+
+							       s*(ctx->rcont6[i]+s1*(ctx->rcont7[i]+s*ctx->rcont8[i]))))));
+
+} /* contd8 */
+
+
 /* core integrator */
-static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
+static int dopcor (dop853_ctx_t *ctx, int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                    double hmax, double h, double* rtoler, double* atoler,
                    int itoler, FILE* fileout, SolTrait solout, int iout,
                    long int nmax, double uround, int meth, long int nstiff, double safe,
@@ -191,6 +174,16 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   double   d61, d66, d67, d68, d69, d610, d611, d612, d613, d614, d615, d616;
   double   d71, d76, d77, d78, d79, d710, d711, d712, d713, d714, d715, d716;
 
+  /* local aliases for context arrays */
+  double *yy1 = ctx->yy1, *k1 = ctx->k1, *k2 = ctx->k2, *k3 = ctx->k3;
+  double *k4 = ctx->k4, *k5 = ctx->k5, *k6 = ctx->k6, *k7 = ctx->k7;
+  double *k8 = ctx->k8, *k9 = ctx->k9, *k10 = ctx->k10;
+  double *rcont1 = ctx->rcont1, *rcont2 = ctx->rcont2;
+  double *rcont3 = ctx->rcont3, *rcont4 = ctx->rcont4;
+  double *rcont5 = ctx->rcont5, *rcont6 = ctx->rcont6;
+  double *rcont7 = ctx->rcont7, *rcont8 = ctx->rcont8;
+  int nrds = ctx->nrds;
+
   /* initialisations */
   switch (meth)
     {
@@ -384,16 +377,16 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   iord = 8;
   if (h == 0.0)
     h = hinit (nptr, fcn, x, y, posneg, k1, k2, k3, iord, hmax, atoler, rtoler, itoler);
-  nfcn += 2;
+  ctx->nfcn += 2;
   reject = 0;
-  xold = x;
+  ctx->xold = x;
 
   if (iout)
     {
       irtrn = 1;
-      hout = 1.0;
-      xout = x;
-      solout (naccpt+1, xold, x, y, nptr, &irtrn);
+      ctx->hout = 1.0;
+      ctx->xout = x;
+      solout (ctx->naccpt+1, ctx->xold, x, y, nptr, &irtrn);
       if (irtrn < 0)
         {
           /* if (fileout) */
@@ -405,12 +398,12 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
   /* basic integration step */
   while (1)
     {
-      if (nstep > nmax)
+      if (ctx->nstep > nmax)
         {
           /* if (fileout) */
 	  Rprintf (_("exit of dop853 at x = %.16e, more than nmax = %li are needed\n"), x, nmax);
-          xout = x;
-          hout = h;
+          ctx->xout = x;
+          ctx->hout = h;
           return -2;
         }
 
@@ -418,8 +411,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
         {
           /* if (fileout) */
 	  Rprintf (_("exit of dop853 at x = %.16e, step size too small h = %.16e\n"), x, h);
-          xout = x;
-          hout = h;
+          ctx->xout = x;
+          ctx->hout = h;
           return -3;
         }
 
@@ -429,7 +422,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           last = 1;
         }
 
-      nstep++;
+      ctx->nstep++;
 
       /* the twelve stages */
       for (i = 0; i < n; i++)
@@ -472,7 +465,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                              a127*k7[i] + a128*k8[i] + a129*k9[i] +
                              a1210*k10[i] + a1211*k2[i]);
       fcn (nptr, xph, yy1, k3);
-      nfcn += 11;
+      ctx->nfcn += 11;
       for (i = 0; i < n; i++)
         {
           k4[i] = b1*k1[i] + b6*k6[i] + b7*k7[i] + b8*k8[i] + b9*k9[i] +
@@ -525,12 +518,12 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* step accepted */
 
           facold = max_d (err, 1.0E-4);
-          naccpt++;
+          ctx->naccpt++;
           fcn (nptr, xph, k5, k4);
-          nfcn++;
+          ctx->nfcn++;
 
           /* stiffness detection */
-          if (!(naccpt % nstiff) || (iasti > 0))
+          if (!(ctx->naccpt % nstiff) || (iasti > 0))
             {
               stnum = 0.0;
               stden = 0.0;
@@ -550,8 +543,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                   if (iasti == 15)
                     {
 		      Rprintf (_("the problem seems to become stiff at x = %.16e\n"), x);
-		      xout = x;
-		      hout = h;
+		      ctx->xout = x;
+		      ctx->hout = h;
 		      return -4;
                     }
                 }
@@ -621,7 +614,7 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
                                      a169*k9[i] + a1613*k4[i] + a1614*k10[i] +
                                      a1615*k2[i]);
               fcn (nptr, x+c16*h, yy1, k3);
-              nfcn += 3;
+              ctx->nfcn += 3;
 
               /* final preparation */
               if (nrds == n)
@@ -653,14 +646,14 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
 
           memcpy (k1, k4, n * sizeof(double));
           memcpy (y, k5, n * sizeof(double));
-          xold = x;
+          ctx->xold = x;
           x = xph;
 
           if (iout)
             {
-              hout = h;
-              xout = x;
-              solout (naccpt+1, xold, x, y, nptr, &irtrn);
+              ctx->hout = h;
+              ctx->xout = x;
+              solout (ctx->naccpt+1, ctx->xold, x, y, nptr, &irtrn);
               if (irtrn < 0)
                 {
                   /* if (fileout) */
@@ -672,8 +665,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* normal exit */
           if (last)
             {
-              hout=hnew;
-              xout = x;
+              ctx->hout=hnew;
+              ctx->xout = x;
               return 1;
             }
 
@@ -689,8 +682,8 @@ static int dopcor (int *nptr, FcnEqDiff fcn, double x, double* y, double xend,
           /* step rejected */
           hnew = h / min_d (facc1, fac11/safe);
           reject = 1;
-          if (naccpt >= 1)
-            nrejct=nrejct + 1;
+          if (ctx->naccpt >= 1)
+            ctx->nrejct=ctx->nrejct + 1;
           last = 0;
         }
 
@@ -707,14 +700,22 @@ int dop853
  double safe, double fac1, double fac2, double beta, double hmax, double h,
  long int nmax, int meth, long int nstiff, int nrdens, int* icont, int licont)
 {
+  dop853_ctx_t ctx;
   int          arret, idid;
   int i;
   int n = nptr[0];
 
-  /* initialisations */
-  nfcn = nstep = naccpt = nrejct = arret = 0;
-  rcont1 = rcont2 = rcont3 = rcont4 = rcont5 = rcont6 = rcont7 = rcont8 = NULL;
-  indir = NULL;
+  /* initialise context */
+  ctx.nfcn = ctx.nstep = ctx.naccpt = ctx.nrejct = 0;
+  ctx.hout = ctx.xold = ctx.xout = 0.0;
+  ctx.nrds = 0;
+  ctx.indir = NULL;
+  ctx.rcont1 = ctx.rcont2 = ctx.rcont3 = ctx.rcont4 = NULL;
+  ctx.rcont5 = ctx.rcont6 = ctx.rcont7 = ctx.rcont8 = NULL;
+  ctx.yy1 = ctx.k1 = ctx.k2 = ctx.k3 = ctx.k4 = ctx.k5 = NULL;
+  ctx.k6 = ctx.k7 = ctx.k8 = ctx.k9 = ctx.k10 = NULL;
+
+  arret = 0;
 
   /* n, the dimension of the system */
   if (n == INT_MAX)
@@ -768,19 +769,19 @@ int dop853
   else if (nrdens)
     {
       /* is there enough memory to allocate rcont12345678&indir ? */
-      rcont1 = R_Calloc(nrdens,double);
-      rcont2 = R_Calloc(nrdens,double);
-      rcont3 = R_Calloc(nrdens,double);
-      rcont4 = R_Calloc(nrdens,double);
-      rcont5 = R_Calloc(nrdens,double);
-      rcont6 = R_Calloc(nrdens,double);
-      rcont7 = R_Calloc(nrdens,double);
-      rcont8 = R_Calloc(nrdens,double);
+      ctx.rcont1 = R_Calloc(nrdens,double);
+      ctx.rcont2 = R_Calloc(nrdens,double);
+      ctx.rcont3 = R_Calloc(nrdens,double);
+      ctx.rcont4 = R_Calloc(nrdens,double);
+      ctx.rcont5 = R_Calloc(nrdens,double);
+      ctx.rcont6 = R_Calloc(nrdens,double);
+      ctx.rcont7 = R_Calloc(nrdens,double);
+      ctx.rcont8 = R_Calloc(nrdens,double);
       if (nrdens < n)
-        indir = R_Calloc(n,int);
+        ctx.indir = R_Calloc(n,int);
 
-      if (!rcont1 || !rcont2 || !rcont3 || !rcont4 || !rcont5 ||
-          !rcont6 || !rcont7 || !rcont8 || (!indir && (nrdens < n)))
+      if (!ctx.rcont1 || !ctx.rcont2 || !ctx.rcont3 || !ctx.rcont4 || !ctx.rcont5 ||
+          !ctx.rcont6 || !ctx.rcont7 || !ctx.rcont8 || (!ctx.indir && (nrdens < n)))
         {
           /* if (fileout) */
 	  Rprintf ( _("not enough free memory for rcont12345678&indir\n"));
@@ -792,7 +793,7 @@ int dop853
         {
           if (icont)
             Rprintf ( _("warning : when nrdens = n there is no need allocating memory for icont\n"));
-          nrds = n;
+          ctx.nrds = n;
         }
       else if (licont < nrdens)
         {
@@ -804,11 +805,11 @@ int dop853
         {
           if ((iout < 2) && fileout)
             fprintf (fileout, "Warning : put iout = 2 for dense output\n");
-          nrds = nrdens;
+          ctx.nrds = nrdens;
           for (i = 0; i < n; i++)
-            indir[i] = INT_MAX;
+            ctx.indir[i] = INT_MAX;
           for (i = 0; i < nrdens; i++)
-            indir[icont[i]] = i;
+            ctx.indir[icont[i]] = i;
         }
     }
 
@@ -855,19 +856,20 @@ int dop853
     hmax = xend - x;
 
   /* is there enough free memory for the method ? */
-  yy1 =R_Calloc(n,double);
-  k1 = R_Calloc(n,double);
-  k2 = R_Calloc(n,double);
-  k3 = R_Calloc(n,double);
-  k4 = R_Calloc(n,double);
-  k5 = R_Calloc(n,double);
-  k6 = R_Calloc(n,double);
-  k7 = R_Calloc(n,double);
-  k8 = R_Calloc(n,double);
-  k9 = R_Calloc(n,double);
-  k10 = R_Calloc(n,double);
-
-  if (!yy1 || !k1 || !k2 || !k3 || !k4 || !k5 || !k6 || !k7 || !k8 || !k9 || !k10)
+  ctx.yy1 =R_Calloc(n,double);
+  ctx.k1 = R_Calloc(n,double);
+  ctx.k2 = R_Calloc(n,double);
+  ctx.k3 = R_Calloc(n,double);
+  ctx.k4 = R_Calloc(n,double);
+  ctx.k5 = R_Calloc(n,double);
+  ctx.k6 = R_Calloc(n,double);
+  ctx.k7 = R_Calloc(n,double);
+  ctx.k8 = R_Calloc(n,double);
+  ctx.k9 = R_Calloc(n,double);
+  ctx.k10 = R_Calloc(n,double);
+
+  if (!ctx.yy1 || !ctx.k1 || !ctx.k2 || !ctx.k3 || !ctx.k4 || !ctx.k5 ||
+      !ctx.k6 || !ctx.k7 || !ctx.k8 || !ctx.k9 || !ctx.k10)
     {
       if (fileout)
         fprintf (fileout, "Not enough free memory for the method\n");
@@ -877,107 +879,79 @@ int dop853
   /* when a failure has occured, we return -1 */
   if (arret)
     {
-      if (k10)
-        R_Free (k10);
-      if (k9)
-        R_Free (k9);
-      if (k8)
-        R_Free (k8);
-      if (k7)
-        R_Free (k7);
-      if (k6)
-        R_Free (k6);
-      if (k5)
-        R_Free (k5);
-      if (k4)
-        R_Free (k4);
-      if (k3)
-        R_Free (k3);
-      if (k2)
-        R_Free (k2);
-      if (k1)
-        R_Free (k1);
-      if (yy1)
-        R_Free (yy1);
-      if (indir)
-        R_Free (indir);
-      if (rcont8)
-        R_Free (rcont8);
-      if (rcont7)
-        R_Free (rcont7);
-      if (rcont6)
-        R_Free (rcont6);
-      if (rcont5)
-        R_Free (rcont5);
-      if (rcont4)
-        R_Free (rcont4);
-      if (rcont3)
-        R_Free (rcont3);
-      if (rcont2)
-        R_Free (rcont2);
-      if (rcont1)
-        R_Free (rcont1);
+      if (ctx.k10)
+        R_Free (ctx.k10);
+      if (ctx.k9)
+        R_Free (ctx.k9);
+      if (ctx.k8)
+        R_Free (ctx.k8);
+      if (ctx.k7)
+        R_Free (ctx.k7);
+      if (ctx.k6)
+        R_Free (ctx.k6);
+      if (ctx.k5)
+        R_Free (ctx.k5);
+      if (ctx.k4)
+        R_Free (ctx.k4);
+      if (ctx.k3)
+        R_Free (ctx.k3);
+      if (ctx.k2)
+        R_Free (ctx.k2);
+      if (ctx.k1)
+        R_Free (ctx.k1);
+      if (ctx.yy1)
+        R_Free (ctx.yy1);
+      if (ctx.indir)
+        R_Free (ctx.indir);
+      if (ctx.rcont8)
+        R_Free (ctx.rcont8);
+      if (ctx.rcont7)
+        R_Free (ctx.rcont7);
+      if (ctx.rcont6)
+        R_Free (ctx.rcont6);
+      if (ctx.rcont5)
+        R_Free (ctx.rcont5);
+      if (ctx.rcont4)
+        R_Free (ctx.rcont4);
+      if (ctx.rcont3)
+        R_Free (ctx.rcont3);
+      if (ctx.rcont2)
+        R_Free (ctx.rcont2);
+      if (ctx.rcont1)
+        R_Free (ctx.rcont1);
 
       return -1;
     }
   else
     {
-      idid = dopcor (nptr, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout,
+      idid = dopcor (&ctx, nptr, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout,
                      solout, iout, nmax, uround, meth, nstiff, safe, beta, fac1, fac2, icont);
-      R_Free (k10);
-      R_Free (k9);
-      R_Free (k8);
-      R_Free (k7);
-      R_Free (k6);
-      R_Free (k5);    /* reverse order R_Freeing too increase chances */
-      R_Free (k4);    /* of efficient dynamic memory managing       */
-      R_Free (k3);
-      R_Free (k2);
-      R_Free (k1);
-      R_Free (yy1);
-      if (indir)
-        R_Free (indir);
-      if (rcont8)
+      R_Free (ctx.k10);
+      R_Free (ctx.k9);
+      R_Free (ctx.k8);
+      R_Free (ctx.k7);
+      R_Free (ctx.k6);
+      R_Free (ctx.k5);    /* reverse order R_Freeing too increase chances */
+      R_Free (ctx.k4);    /* of efficient dynamic memory managing       */
+      R_Free (ctx.k3);
+      R_Free (ctx.k2);
+      R_Free (ctx.k1);
+      R_Free (ctx.yy1);
+      if (ctx.indir)
+        R_Free (ctx.indir);
+      if (ctx.rcont8)
         {
-          R_Free (rcont8);
-          R_Free (rcont7);
-          R_Free (rcont6);
-          R_Free (rcont5);
-          R_Free (rcont4);
-          R_Free (rcont3);
-          R_Free (rcont2);
-          R_Free (rcont1);
+          R_Free (ctx.rcont8);
+          R_Free (ctx.rcont7);
+          R_Free (ctx.rcont6);
+          R_Free (ctx.rcont5);
+          R_Free (ctx.rcont4);
+          R_Free (ctx.rcont3);
+          R_Free (ctx.rcont2);
+          R_Free (ctx.rcont1);
         }
 
       return idid;
     }
 
 } /* dop853 */
-
-
-/* dense output function */
-double contd8 (int ii, double x)
-{
-  int i;
-  double       s, s1;
-
-  i = INT_MAX;
-
-  if (!indir)
-    i = ii;
-  else
-    i = indir[ii];
-
-  if (i == INT_MAX)
-    {
-      Rprintf (_("no dense output available for %uth component"), ii);
-      return 0.0;
-    }
-
-  s = (x - xold) / hout;
-  s1 = 1.0 - s;
-
-  return rcont1[i]+s*(rcont2[i]+s1*(rcont3[i]+s*(rcont4[i]+s1*(rcont5[i]+
-							       s*(rcont6[i]+s1*(rcont7[i]+s*rcont8[i]))))));
-
-} /* contd8 */
diff --git a/src/dop853.h b/src/dop853.h
index cbb50223b..e7c1d19ab 100644
--- a/src/dop853.h
+++ b/src/dop853.h
@@ -180,6 +180,17 @@ nfcnRead    Number of function calls.
 typedef void (*FcnEqDiff)(int *nptr, double x, double *y, double *f);
 typedef void (*SolTrait)(long int nr, double xold, double x, double* y, int *nptr, int* irtrn);
 
+/* Context struct holding all solver state, enabling thread-safe reentrant calls.
+   Previously these were file-scope static variables in dop853.c. */
+typedef struct {
+  long int     nfcn, nstep, naccpt, nrejct;
+  double       hout, xold, xout;
+  int          nrds, *indir;
+  double       *yy1, *k1, *k2, *k3, *k4, *k5, *k6, *k7, *k8, *k9, *k10;
+  double       *rcont1, *rcont2, *rcont3, *rcont4;
+  double       *rcont5, *rcont6, *rcont7, *rcont8;
+} dop853_ctx_t;
+
 
 extern int dop853
  (int *nptr,      /* dimension of the system <= INT_MAX-1*/
@@ -207,15 +218,3 @@ extern int dop853
   int* icont, /* indexes of components for which dense output is required, >= nrdens */
   int licont  /* declared length of icon */
  );
-
-extern double contd8
- (int ii,     /* index of desired component */
-  double x             /* approximation at x */
- );
-
-extern long int nfcnRead (void);   /* encapsulation of statistical data */
-extern long int nstepRead (void);
-extern long int naccptRead (void);
-extern long int nrejctRead (void);
-extern double hRead (void);
-extern double xRead (void);

From 97671af9a39bbac2a62fff3ec8a459c93c678137 Mon Sep 17 00:00:00 2001
From: Bill Denney <wdenney@humanpredictions.com>
Date: Mon, 2 Mar 2026 21:04:50 -0500
Subject: [PATCH 08/13] Fix linCmtB macro collision: use J macro instead of
 lcb.J

The macro `#define J lcb.J` causes `lcb.J` on line 683 to expand to
`lcb.lcb.J`, which fails because linB_t has no member named lcb.
Use the macro name `J` directly instead.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 src/linCmt.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linCmt.cpp b/src/linCmt.cpp
index 743c7d873..0de48b85d 100644
--- a/src/linCmt.cpp
+++ b/src/linCmt.cpp
@@ -680,7 +680,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
       }
     }
   }
-  lc.getJacCp(lcb.J, fx, theta, Jg);
+  lc.getJacCp(J, fx, theta, Jg);
   return lc.adjustF(fx, theta, ind->linCmtHV);
 #undef fx
 #undef J

From 0e1db6a22ca2850876871339377eff5eeae312ab Mon Sep 17 00:00:00 2001
From: Matthew Fidler <matthew.fidler@gmail.com>
Date: Mon, 16 Mar 2026 13:05:23 -0500
Subject: [PATCH 09/13] Update src/linCmt.cpp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 src/linCmt.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/linCmt.cpp b/src/linCmt.cpp
index 0de48b85d..67eacfddf 100644
--- a/src/linCmt.cpp
+++ b/src/linCmt.cpp
@@ -486,7 +486,7 @@ extern "C" double linCmtB(rx_solve *rx, int id,
                           // Oral parameters
                           double ka) {
   // Per-thread linCmtB state (matching linCmtA pattern for thread safety)
-  linB_t &lcb = __linCmtB[omp_get_thread_num()];
+  linB_t &lcb = __linCmtB[rx_get_thread(__linCmtB.size())];
 #define fx        lcb.fx
 #define Jg        lcb.Jg
 #define lc        lcb.lc

From c58e2521931ef2a929f895ae628e5a4f31441261 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 18:06:01 +0000
Subject: [PATCH 10/13] Initial plan


From 59e435692f4378791cdf80a0b8e2279cd7f710bb Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 18:06:19 +0000
Subject: [PATCH 11/13] Initial plan


From 0d4f8a2ae3618080614eb68a541368059dab3342 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 18:10:04 +0000
Subject: [PATCH 12/13] Update dop853.h documentation to match thread-safe
 refactoring

- Update contd8() signature from 2 to 3 arguments (ctx, i, s)
- Replace deleted accessor functions (xRead, hRead, etc.) with
  equivalent dop853_ctx_t struct field documentation
- Update "Remarks about C version" to mention context struct
- Fix typo: 'stff' -> 'stiff'

Co-authored-by: mattfidler <514778+mattfidler@users.noreply.github.com>
---
 src/dop853.h | 30 ++++++++++++++++--------------
 1 file changed, 16 insertions(+), 14 deletions(-)

diff --git a/src/dop853.h b/src/dop853.h
index e7c1d19ab..ac09f7e66 100644
--- a/src/dop853.h
+++ b/src/dop853.h
@@ -23,8 +23,8 @@ Version of Mai 2, 1994.
 Remarks about the C version : this version allocates memory by itself, the
 iwork array (among the initial FORTRAN parameters) has been splitted into
 independant initial parameters, the statistical variables and last step size
-and x have been encapsulated in the module and are now accessible through
-dedicated functions; the variable names have been kept to maintain a kind
+and x have been encapsulated in a dop853_ctx_t context struct for thread-safe
+reentrant usage; the variable names have been kept to maintain a kind
 of reading compatibility between the C and FORTRAN codes; adaptation made by
 J.Colinge (COLINGE@DIVSUN.UNIGE.CH).
 
@@ -74,10 +74,11 @@ solout   A pointer to the output function called during integration.
 	 Continuous output : during the calls to solout, a continuous solution
 	 for the interval (xold,x) is available through the function
 
-	   contd8(i,s)
+	   contd8(ctx, i, s)
 
 	 which provides an approximation to the i-th component of the solution
-	 at the point s (s must lie in the interval (xold,x)).
+	 at the point s (s must lie in the interval (xold,x)), where ctx is
+	 the dop853_ctx_t context pointer.
 
 iout     Switch for calling solout :
 	   iout=0 : no call,
@@ -145,7 +146,7 @@ Memory requirements
 OUTPUT PARAMETERS
 -----------------
 
-y       numerical solution at x=xRead() (see below).
+y       numerical solution at x=ctx->xout (see dop853_ctx_t below).
 
 dopri5 returns the following values
 
@@ -154,21 +155,22 @@ dopri5 returns the following values
 	-1 : input is not consistent,
 	-2 : larger nmax is needed,
 	-3 : step size becomes too small,
-	-4 : the problem is probably stff (interrupted).
+	-4 : the problem is probably stiff (interrupted).
 
 
-Several functions provide access to different values :
+After a successful call to dop853, the following fields of the dop853_ctx_t
+context struct provide solver statistics and output values :
 
-xRead   x value for which the solution has been computed (x=xend after
+ctx->xout    x value for which the solution has been computed (x=xend after
 	successful return).
 
-hRead   Predicted step size of the last accepted step (useful for a subsequent
-	call to dop853).
+ctx->hout    Predicted step size of the last accepted step (useful for a
+	subsequent call to dop853).
 
-nstepRead   Number of used steps.
-naccptRead  Number of accepted steps.
-nrejctRead  Number of rejected steps.
-nfcnRead    Number of function calls.
+ctx->nstep   Number of used steps.
+ctx->naccpt  Number of accepted steps.
+ctx->nrejct  Number of rejected steps.
+ctx->nfcn    Number of function calls.
 
 
 */

From 512b2fc8db53f6f05924a1e9f34d95cbda5020e9 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Mon, 16 Mar 2026 18:11:39 +0000
Subject: [PATCH 13/13] Preserve first error's diagnostic code in op->naTime
 using read-check-write pattern

Co-authored-by: mattfidler <514778+mattfidler@users.noreply.github.com>
---
 inst/include/rxode2parseGetTime.h    | 36 +++++++++++++++++++++-------
 inst/include/rxode2parseHandleEvid.h |  9 +++++--
 2 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/inst/include/rxode2parseGetTime.h b/inst/include/rxode2parseGetTime.h
index c7f402d72..01feba1c2 100644
--- a/inst/include/rxode2parseGetTime.h
+++ b/inst/include/rxode2parseGetTime.h
@@ -39,11 +39,16 @@ static inline double getLag(rx_solving_options_ind *ind, int id, int cmt, double
   double ret = LAG(id, cmt, time);
   if (ISNA(ret)) {
     int newBadSolve = 1;
-    int newNaTime = 1 + 10*cmt;
 #pragma omp atomic write
     op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 1 + 10*cmt;
 #pragma omp atomic write
-    op->naTime = newNaTime;
+      op->naTime = newNaTime;
+    }
   }
   return ret;
 }
@@ -54,11 +59,16 @@ static inline double getRate(rx_solving_options_ind *ind, int id, int cmt, doubl
   double ret = RATE(id, cmt, dose, t);
   if (ISNA(ret)){
     int newBadSolve = 1;
-    int newNaTime = 2 + 10*cmt;
 #pragma omp atomic write
     op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 2 + 10*cmt;
 #pragma omp atomic write
-    op->naTime = newNaTime;
+      op->naTime = newNaTime;
+    }
   }
   return ret;
 }
@@ -70,11 +80,16 @@ static inline double getDur(rx_solving_options_ind *ind, int id, int cmt, double
   double ret = DUR(id, cmt, dose, t);
   if (ISNA(ret)){
     int newBadSolve = 1;
-    int newNaTime = 3 + 10*cmt;
 #pragma omp atomic write
     op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 3 + 10*cmt;
 #pragma omp atomic write
-    op->naTime = newNaTime;
+      op->naTime = newNaTime;
+    }
   }
   return ret;
 }
@@ -347,11 +362,16 @@ static inline double handleInfusionItem(int idx, rx_solve *rx, rx_solving_option
     if (ISNA(f)){
       rx_solving_options *op = &op_global;
       int newBadSolve = 1;
-      int newNaTime = 4 + 10*ind->cmt;
 #pragma omp atomic write
       op->badSolve = newBadSolve;
+      int curNaTime;
+#pragma omp atomic read
+      curNaTime = op->naTime;
+      if (curNaTime == 0) {
+        int newNaTime = 4 + 10*ind->cmt;
 #pragma omp atomic write
-      op->naTime = newNaTime;
+        op->naTime = newNaTime;
+      }
     }
     double durOld = (getAllTimes(ind, ind->idose[infEidx]) -
                      getAllTimes(ind, ind->idose[infBidx]));
diff --git a/inst/include/rxode2parseHandleEvid.h b/inst/include/rxode2parseHandleEvid.h
index e4679a707..312c29a17 100644
--- a/inst/include/rxode2parseHandleEvid.h
+++ b/inst/include/rxode2parseHandleEvid.h
@@ -286,11 +286,16 @@ static inline double getAmt(rx_solving_options_ind *ind, int id, int cmt,
   if (ISNA(ret)){
     rx_solving_options *op = &op_global;
     int newBadSolve = 1;
-    int newNaTime = 5 + 10*cmt;
 #pragma omp atomic write
     op->badSolve = newBadSolve;
+    int curNaTime;
+#pragma omp atomic read
+    curNaTime = op->naTime;
+    if (curNaTime == 0) {
+      int newNaTime = 5 + 10*cmt;
 #pragma omp atomic write
-    op->naTime = newNaTime;
+      op->naTime = newNaTime;
+    }
   }
   return ret;
 }