Add required idid status flag to stiff3 APIs

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/cfc168cc-31b3-4e0c-9237-d65abc89893e

Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>

Author: Copilot

README.md

@@ -15,6 +15,7 @@ This repository provides a heavily refactored version of the original implementa
 
 - Third-order Rosenbrock-type (semi-implicit Runge-Kutta) integrator suitable for stiff ODE systems
 - Adaptive stepsize control with error tolerance per component
+- Required integer exit flag (`idid`) reporting success and early-exit conditions
 - Optional maximum absolute half-step size (`hmax`) to cap step growth
 - Optional runtime statistics output: successful steps (`nacc`), rejected steps (`nrej`), rhs evaluations (`nfev`), Jacobian evaluations (`njev`), LU decompositions (`nlu`), and linear-system solves (`nsol`)
 - Optional explicit workspace interface via `stiff3(..., rwork, iwork, ...)` for caller-managed memory
@@ -103,7 +104,7 @@ program vanpol
   integer, parameter :: n = 2
   real(wp), parameter :: mu = 10.0_wp
   real(wp) :: y(n), w(n), x0, x, x1, h0, eps, hmax
-  integer :: stats(6)
+  integer :: idid, stats(6)
 
 ! initial value
   y = [1.0_wp, 1.0_wp]
@@ -119,8 +120,9 @@ program vanpol
   x = x0
   x1 = 100.0_wp
 ! integrate system of ODEs
-  call stiff3(n,fun,x,y,x1,jac,h0,eps,w,solout=out,stats=stats,hmax=hmax)
-! on successful exit, x == x1
+  call stiff3(n,fun,x,y,x1,jac,h0,eps,w,idid,solout=out,stats=stats,hmax=hmax)
+  if (idid /= 0) error stop 'stiff3 failed'
+! on successful exit, idid == 0 and x == x1
   print '(A)', 'Solver statistics'
   print '(A,I0)', '  nacc: ', stats(1)
   print '(A,I0)', '  nrej: ', stats(2)
@@ -180,7 +182,17 @@ The two tolerance parameters `eps` and `w` together control how accurately the s
 
 A reasonable first choice is `eps = 1.0e-4` with all `w(i) = 1.0`, which requests roughly four significant digits from every component.
 
-For interoperability scenarios where the caller manages memory allocation (e.g. language bindings), `stiff3` also provides an overload with explicit work arrays: `rwork` of size `n*(7 + 2*n)` and `iwork` of size `n`.
+For interoperability scenarios where the caller manages memory allocation (e.g. language bindings), `stiff3` also provides an overload with explicit work arrays: `rwork` of size `n*(7 + 2*n)` and `iwork` of size `n`. The required `idid` argument appears before the optional arguments in both overloads:
+
+- `call stiff3(n, fun, x, y, x1, jac, h0, eps, w, idid, [solout=...], [stats=...], [hmax=...])`
+- `call stiff3(n, fun, x, y, x1, jac, h0, eps, w, rwork, iwork, idid, [solout=...], [stats=...], [hmax=...])`
+
+The integer exit flag `idid` reports solver status:
+
+- `0` — successful completion at `x1`
+- `-1` — LU factorization failed because the Jacobian system matrix was singular
+- `-2` — integration interrupted by the user `solout` callback (`irtrn < 0`)
+- `-3` — step-size underflow occurred during bisection (`abs(h) <= spacing(x_current)`)
 
 When using this explicit-workspace overload together with `solout`, accepted-step dense output is available through the generic interface `stiff3_interp`:
 

example/e5_chemical.f90

@@ -65,13 +65,15 @@ subroutine integrate_case(eps, y, stats, err, elapsed)
     real(wp), intent(out), optional :: err, elapsed
 
     real(wp) :: h0, t0, t1, x
+    integer :: idid
 
     y = y0
     x = x0
     h0 = h0_initial
     call cpu_time(t0)
-    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, stats=stats)
+    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, idid, stats=stats)
     call cpu_time(t1)
+    if (idid /= 0) error stop 'stiff3 failed in e5_chemical example'
     if (present(err)) err = max_relative_error(y, yref)
     if (present(elapsed)) elapsed = max(t1 - t0, 0.0_wp)
   end subroutine

example/hires.f90

@@ -68,13 +68,15 @@ subroutine integrate_case(eps, y, stats, err, cpu_time_seconds)
 
     real(wp), parameter :: min_cpu_time = 1.0e-12_wp
     real(wp) :: h0, cpu_start, cpu_end, x
+    integer :: idid
 
     y = y0
     x = x0
     h0 = h0_initial
     call cpu_time(cpu_start)
-    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, stats=stats)
+    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, idid, stats=stats)
     call cpu_time(cpu_end)
+    if (idid /= 0) error stop 'stiff3 failed in hires example'
     if (present(err)) err = max_relative_error(y, yref)
     if (present(cpu_time_seconds)) cpu_time_seconds = max(min_cpu_time, cpu_end - cpu_start)
   end subroutine

example/lorenz.f90

@@ -7,6 +7,7 @@ program lorenz
   integer, parameter :: n = 3
   real(wp), parameter :: sigma = 10.0_wp, rho = 28.0_wp, beta = 8.0_wp/3.0_wp
   real(wp) :: y(n), w(n), x0, x1, h0, eps
+  integer :: idid
 
 ! initial value
   y = [1.0_wp, 1.0_wp, 1.0_wp]
@@ -19,7 +20,8 @@ program lorenz
   x0 = 0.0_wp
   x1 = 40.0_wp
 ! integrate system of ODEs
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,idid,solout=out)
+  if (idid /= 0) error stop 'stiff3 failed in lorenz example'
 
 contains
 

example/oregonator.f90

@@ -70,7 +70,7 @@ subroutine integrate_case(eps, y, stats, cpu_time_seconds, solution_filename)
     real(wp), intent(out) :: cpu_time_seconds
     character(len=*), intent(in), optional :: solution_filename
 
-    integer :: ios
+    integer :: ios, idid
     real(wp), parameter :: min_cpu_time = 1.0e-12_wp
     real(wp) :: cpu_start, cpu_end, h0, x
     logical :: write_solution
@@ -87,8 +87,9 @@ subroutine integrate_case(eps, y, stats, cpu_time_seconds, solution_filename)
     end if
 
     call cpu_time(cpu_start)
-    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, solout=output_solution, stats=stats)
+    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, idid, solout=output_solution, stats=stats)
     call cpu_time(cpu_end)
+    if (idid /= 0) error stop 'stiff3 failed in oregonator example'
     cpu_time_seconds = max(min_cpu_time, cpu_end - cpu_start)
 
     if (write_solution) then

example/pendant_drop.f90

@@ -7,7 +7,7 @@ program pendant_drop
   integer, parameter :: n = 3
   real(wp) :: y(n), w(n), x0, x1, h0, eps
   real(wp) :: pL, drho, ds
-  integer :: stats(6)
+  integer :: stats(6), idid
 
 ! dimensionless pressure offset and density difference
   pL = 3.3888_wp
@@ -29,7 +29,8 @@ program pendant_drop
   w = 1.0_wp
 
 ! integrate Young-Laplace system
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out,stats=stats)
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,idid,solout=out,stats=stats)
+  if (idid /= 0) error stop 'stiff3 failed in pendant_drop example'
   print '(A,3(I0,1X))', 'accepted rejected nfev: ', stats(1), stats(2), stats(3)
   print '(A,3(I0,1X))', 'njev nlu nsol: ', stats(4), stats(5), stats(6)
 

example/predator_prey.f90

@@ -8,6 +8,7 @@ program predator_prey
   real(wp), parameter :: alpha = 1.1_wp, beta = 0.4_wp
   real(wp), parameter :: delta = 0.1_wp, gamma = 0.4_wp
   real(wp) :: y(n), w(n), x0, x1, h0, eps
+  integer :: idid
 
 ! initial value
   y = [10.0_wp, 5.0_wp]
@@ -20,7 +21,8 @@ program predator_prey
   x0 = 0.0_wp
   x1 = 40.0_wp
 ! integrate system of ODEs
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,idid,solout=out)
+  if (idid /= 0) error stop 'stiff3 failed in predator_prey example'
 
 contains
 

example/ring_modulator.f90

@@ -76,13 +76,15 @@ subroutine integrate_case(eps, y, stats, err, cpu_time_seconds)
 
     real(wp), parameter :: min_cpu_time = 1.0e-12_wp
     real(wp) :: h0, cpu_start, cpu_end, x
+    integer :: idid
 
     y = y0
     x = x0
     h0 = initial_step_scale*eps
     call cpu_time(cpu_start)
-    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, stats=stats)
+    call stiff3(n, fun, x, y, x1, jac, h0, eps, w, idid, stats=stats)
     call cpu_time(cpu_end)
+    if (idid /= 0) error stop 'stiff3 failed in ring_modulator example'
     if (present(err)) err = max_relative_error(y(1:nphys), yref)
     if (present(cpu_time_seconds)) cpu_time_seconds = max(min_cpu_time, cpu_end - cpu_start)
   end subroutine

example/robertson.f90

@@ -48,6 +48,7 @@ program main
   integer, parameter :: n = 3
   real(wp) :: y(n), w(n)
   real(wp) :: h0, eps, x0, x1
+  integer :: idid
 
 
 ! initial value
@@ -65,7 +66,8 @@ program main
   x0 = 0.0_wp
   x1 = 10.0_wp
 
-  call stiff3(n,fun,x0,y,x1,dfun,h0,eps,w,solout=output)
+  call stiff3(n,fun,x0,y,x1,dfun,h0,eps,w,idid,solout=output)
+  if (idid /= 0) error stop 'stiff3 failed in robertson example'
 
 contains
 

example/three_equation_system.f90

@@ -9,7 +9,7 @@ program three_equation_system
   real(wp), parameter :: k2 = 1000.0_wp
   real(wp), parameter :: k3 = 2500.0_wp
   real(wp) :: y(n), w(n), x0, x1, h0, eps
-  integer :: stats(6)
+  integer :: stats(6), idid
 
 ! initial value
   y = [1.0_wp, 1.0_wp, 0.0_wp]
@@ -22,7 +22,8 @@ program three_equation_system
   x0 = 0.0_wp
   x1 = 50.0_wp
 
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,stats=stats)
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,idid,stats=stats)
+  if (idid /= 0) error stop 'stiff3 failed in three_equation_system example'
   print '(A,3(I0,1X))', 'accepted rejected nfev: ', stats(1), stats(2), stats(3)
   print '(A,3(I0,1X))', 'njev nlu nsol: ', stats(4), stats(5), stats(6)