Expand ODE solver stats and align solout grid-point callback semantics (#36)

Copilot · ivan-pi · web-flow · commit 0b299ecc9491 · 2026-05-11T16:56:08.000+02:00
* Expand solver stats to include steps and back-substitutions Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/562c29e2-5602-4cc7-965b-83224e87db8c Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Refine stats test expectation for linear-system solve counter Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/562c29e2-5602-4cc7-965b-83224e87db8c Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Adjust solout initial callback and reorder solver stats Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/8e92236a-2f09-487c-9169-cb3dc8700e1a Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Refine stats test readability and README stats output format Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/ccb23e84-4807-42ce-a5c1-5b2811603ce3 Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Normalize README stats label spacing Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/ccb23e84-4807-42ce-a5c1-5b2811603ce3 Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> --------- Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>
diff --git a/README.md b/README.md
@@ -16,7 +16,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
 - Optional maximum absolute half-step size (`hmax`) to cap step growth
-- Optional runtime statistics output: rhs evaluations (`nfev`), Jacobian evaluations (`njev`), and LU decompositions (`nlu`)
+- 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
 - Dense output helper `stiff3_interp` for accepted steps in `solout`
 - Requires an exact user-supplied Jacobian
@@ -95,7 +95,7 @@ program vanpol
   integer, parameter :: n = 2
   real(wp), parameter :: mu = 10.0_wp
   real(wp) :: y(n), w(n), x0, x1, h0, eps, hmax
-  integer :: irtrn, stats(3)
+  integer :: stats(6)
 
 ! initial value
   y = [1.0_wp, 1.0_wp]
@@ -109,12 +109,15 @@ program vanpol
 ! time interval
   x0 = 0.0_wp
   x1 = 100.0_wp
-! output initial condition
-  irtrn = 0
-  call out(0,x0,x0,y,0,0.0_wp,irtrn)
 ! integrate system of ODEs
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out,stats=stats,hmax=hmax)
-  print '(A,3(I0,1X))', 'nfev njev nlu: ', stats
+  print '(A)', 'Solver statistics'
+  print '(A,I0)', '  nacc: ', stats(1)
+  print '(A,I0)', '  nrej: ', stats(2)
+  print '(A,I0)', '  nfev: ', stats(3)
+  print '(A,I0)', '  njev: ', stats(4)
+  print '(A,I0)', '  nlu: ', stats(5)
+  print '(A,I0)', '  nsol: ', stats(6)
 
 contains
 
diff --git a/example/lorenz.f90 b/example/lorenz.f90
@@ -7,7 +7,6 @@ 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 :: irtrn
 
 ! initial value
   y = [1.0_wp, 1.0_wp, 1.0_wp]
@@ -19,9 +18,6 @@ program lorenz
 ! time interval
   x0 = 0.0_wp
   x1 = 40.0_wp
-! output initial condition
-  irtrn = 0
-  call out(0,x0,x0,y,0,0.0_wp,irtrn)
 ! integrate system of ODEs
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
 
diff --git a/example/pendant_drop.f90 b/example/pendant_drop.f90
@@ -7,6 +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)
 
 ! dimensionless pressure offset and density difference
   pL = 3.3888_wp
@@ -27,10 +28,10 @@ program pendant_drop
   eps = 1.0e-6_wp
   w = 1.0_wp
 
-! output initial condition at s = 0
-  write(*,'(4(E18.12,2X),I4,2X,G0)') 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0, 0.0_wp
 ! integrate Young-Laplace system
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out,stats=stats)
+  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)
 
 contains
 
diff --git a/example/predator_prey.f90 b/example/predator_prey.f90
@@ -8,7 +8,6 @@ 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 :: irtrn
 
 ! initial value
   y = [10.0_wp, 5.0_wp]
@@ -20,9 +19,6 @@ program predator_prey
 ! time interval
   x0 = 0.0_wp
   x1 = 40.0_wp
-! output initial condition
-  irtrn = 0
-  call out(0,x0,x0,y,0,0.0_wp,irtrn)
 ! integrate system of ODEs
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
 
diff --git a/example/robertson.f90 b/example/robertson.f90
@@ -49,7 +49,6 @@ program main
   real(wp) :: y(n), w(n)
   real(wp) :: h0, eps, x0, x1
 
-  integer :: irtrn
 
 ! initial value
   y = [1.0_wp, 0.0_wp, 0.0_wp]
@@ -66,8 +65,6 @@ program main
   x0 = 0.0_wp
   x1 = 10.0_wp
 
-  irtrn = 0
-  call output(0,x0,x0,y,0,0.0_wp,irtrn)
   call stiff3(n,fun,x0,y,x1,dfun,h0,eps,w,solout=output)
 
 contains
diff --git a/example/three_equation_system.f90 b/example/three_equation_system.f90
@@ -8,7 +8,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 :: irtrn, nsteps
+  integer :: stats(6)
 
 ! initial value
   y = [1.0_wp, 1.0_wp, 0.0_wp]
@@ -21,12 +21,9 @@ program three_equation_system
   x0 = 0.0_wp
   x1 = 50.0_wp
 
-  nsteps = 0
-  irtrn = 0
-  call out(0,x0,x0,y,0,0.0_wp,irtrn)
-  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
-
-  print '(A,I0)', 'Accepted steps required: ', nsteps
+  call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,stats=stats)
+  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)
 
 contains
 
@@ -58,16 +55,4 @@ subroutine jac(n,y,df)
     df(3,3) = -k2*y(1) - k3*y(2)
   end subroutine
 
-  subroutine out(nr,told,t,y,ih,qa,irtrn)
-    integer, intent(in) :: nr
-    real(wp), intent(in) :: told
-    real(wp), intent(in) :: t
-    real(wp), intent(in) :: y(:)
-    integer, intent(in) :: ih
-    real(wp), intent(in) :: qa
-    integer, intent(inout) :: irtrn
-
-    nsteps = nr
-  end subroutine
-
 end program
diff --git a/example/vanpol.f90 b/example/vanpol.f90
@@ -6,7 +6,6 @@ program vanpol
   integer, parameter :: n = 2
   real(wp), parameter :: mu = 10.0_wp
   real(wp) :: y(n), w(n), x0, x1, h0, eps
-  integer :: irtrn
 
 ! initial value
   y = [1.0_wp, 1.0_wp]
@@ -18,9 +17,6 @@ program vanpol
 ! time interval
   x0 = 0.0_wp
   x1 = 100.0_wp
-! output initial condition
-  irtrn = 0
-  call out(0,x0,x0,y,0,0.0_wp,irtrn)
 ! integrate system of ODEs
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,solout=out)
 
diff --git a/src/stiff3_solver.f90 b/src/stiff3_solver.f90
@@ -48,7 +48,7 @@ subroutine jacobian_sub(n,y,df)
     subroutine output_sub(nr,xold,x,y,iha,qa,irtrn)
       import wp
       integer, intent(in) :: nr
-        !! Number of successful steps that have been taken
+        !! Number of the current grid point, starting at 1 for the initial value
       real(wp), intent(in) :: xold
         !! Previous value of the independent variable
       real(wp), intent(in) :: x
@@ -105,8 +105,8 @@ subroutine stiff3_auto(n,fun,x,y,xend,jac,h0,eps,w,solout,stats,hmax)
       !! dependent variables at `xend`.
     procedure(output_sub), optional :: solout
       !! User supplied subprogram for output.
-    integer, intent(out), optional :: stats(3)
-      !! Statistics array with `[nfev, njev, nlu]`.
+    integer, intent(out), optional :: stats(6)
+      !! Statistics array with `[nacc, nrej, nfev, njev, nlu, nsol]`.
     real(wp), intent(in), optional :: hmax
       !! Maximum absolute half-step size. If absent or zero, defaults to
       !! `abs(xend - x)`.
@@ -153,8 +153,8 @@ subroutine stiff3_work(n,fun,x,y,xend,jac,h0,eps,w,rwork,iwork,solout,stats,hmax
       !! Integer workspace of size `n`.
     procedure(output_sub), optional :: solout
       !! User supplied subprogram for output.
-    integer, intent(out), optional :: stats(3)
-      !! Statistics array with `[nfev, njev, nlu]`.
+    integer, intent(out), optional :: stats(6)
+      !! Statistics array with `[nacc, nrej, nfev, njev, nlu, nsol]`.
     real(wp), intent(in), optional :: hmax
       !! Maximum absolute half-step size. If absent or zero, defaults to
       !! `abs(xend - x)`.
@@ -191,18 +191,17 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
     real(wp), intent(inout) :: df(n,n), dfold(n,n)
     integer, intent(inout) :: ip(n)
     procedure(output_sub), optional :: solout
-    integer, intent(out), optional :: stats(3)
+    integer, intent(out), optional :: stats(6)
     real(wp), intent(in), optional :: hmax
 
-    integer :: icon, iha, i, j, nr, irtrn
-    integer :: nfev, njev, nlu
+    integer :: icon, iha, i, j, irtrn
+    integer :: nfev, njev, nlu, nacc, nrej, nsol
     real(wp) :: x_current, xold, h, e, es, q, qa, hmax_used
     logical :: have_f
 
   ! icon = 0 except for last step which ends exactly at x1
     icon = 0
 
-    nr = 0
     x_current = x
     if (present(hmax)) then
       if (hmax < 0.0_wp) error stop 'stiff3: hmax must be a non-negative real value'
@@ -218,8 +217,21 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
     nfev = 0
     njev = 0
     nlu = 0
+    nacc = 0
+    nrej = 0
+    nsol = 0
     have_f = .false.
 
+    if (present(solout)) then
+      irtrn = 0
+      call solout(1,x_current,x_current,y,0,0.0_wp,irtrn)
+      if (irtrn < 0) then
+        h0 = h
+        if (present(stats)) stats = [nacc, nrej, nfev, njev, nlu, nsol]
+        return
+      end if
+    end if
+
     outer: do
 
     ! last step - or first step longer than interval
@@ -259,7 +271,7 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
 
     ! perform full integration step
 
-      call sirk3(n,fun,ip,f,y,yk1,yk2,df,2*h,nfev,nlu)
+      call sirk3(n,fun,ip,f,y,yk1,yk2,df,2*h,nfev,nlu,nsol)
 
       do i = 1, n
         ya(i) = y(i)
@@ -277,15 +289,15 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
       inner: do
         iha = iha + 1
 
-        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu)
+        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu,nsol)
         call fun(n,y,f)
         nfev = nfev + 1
         call jac(n,y,df)
         njev = njev + 1
 
         yold1 = y
 
-        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu)
+        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu,nsol)
 
       ! half step integration finished
       ! compute deviation and compare with tolerance
@@ -314,6 +326,7 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
 
         h = h/2.0_wp
         icon = 0
+        nrej = nrej + 1
 
       end do inner
 
@@ -338,13 +351,13 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
 
     ! perform output if appropriate
 
-      nr = nr + 1
+      nacc = nacc + 1
       if (present(solout)) then
         irtrn = 0
-        call solout(nr,xold,x_current,y,iha,qa,irtrn)
+        call solout(nacc+1,xold,x_current,y,iha,qa,irtrn)
         if (irtrn < 0) then
           h0 = h
-          if (present(stats)) stats = [nfev, njev, nlu]
+          if (present(stats)) stats = [nacc, nrej, nfev, njev, nlu, nsol]
           return
         end if
       end if
@@ -353,7 +366,7 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
 
       if (icon == 1) then
         h0 = h
-        if (present(stats)) stats = [nfev, njev, nlu]
+        if (present(stats)) stats = [nacc, nrej, nfev, njev, nlu, nsol]
         return
       end if
 
@@ -441,7 +454,7 @@ subroutine stiff3_interp_vector(xold,x,y,rwork,xeval,yeval)
 
   !> Single-step semi-implicit integration
   !
-  subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
+  subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu,nsol)
     integer, intent(in) :: n
       !! Size of the system of ODEs
     procedure(rhs_sub) :: fun
@@ -464,6 +477,8 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
       !! Number of right-hand side evaluations
     integer, intent(inout) :: nlu
       !! Number of LU decompositions
+    integer, intent(inout) :: nsol
+      !! Number of linear-system solves (back substitutions)
 
     integer :: i
 
@@ -487,6 +502,7 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
     call lu(df,ipiv)
     nlu = nlu + 1
     call back(df,f,ipiv)
+    nsol = nsol + 1
 
     do i = 1, n
       yk1(i) = h*f(i)
@@ -495,6 +511,7 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
     call fun(n,yk2,f)
     nfev = nfev + 1
     call back(df,f,ipiv)
+    nsol = nsol + 1
 
     !
     ! evaluate k2
@@ -510,6 +527,7 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
     ! for convenience stored in yk2
     !
     call back(df,yk2,ipiv)
+    nsol = nsol + 1
     do i = 1, n
       y(i) = y(i) + yk2(i)
     end do
diff --git a/test/ode_hmax.f90 b/test/ode_hmax.f90
@@ -10,7 +10,7 @@ program ode_hmax
   real(wp), parameter :: h0_start = 0.5_wp
   real(wp) :: y_default(n), y_zero(n), y_capped(n), w(n), x0, x1, h0, eps_test
   real(wp) :: max_h_seen
-  integer :: stats_default(3), stats_zero(3), stats_capped(3)
+  integer :: stats_default(6), stats_zero(6), stats_capped(6)
 
   w = 1.0_wp
   x0 = 0.0_wp
@@ -26,8 +26,8 @@ program ode_hmax
   call stiff3(n,fun,x0,y_zero,x1,jac,h0,eps_test,w,stats=stats_zero,hmax=0.0_wp)
 
   if (any(stats_zero /= stats_default)) then
-    print '(A,3(I0,1X),A,3(I0,1X))', &
-      'expected hmax=0 stats [nfev njev nlu] to match default. got ', stats_zero, &
+    print '(A,6(I0,1X),A,6(I0,1X))', &
+      'expected hmax=0 stats [nacc nrej nfev njev nlu nsol] to match default. got ', stats_zero, &
       ' expected ', stats_default
     error stop 1
   end if
diff --git a/test/ode_stats.f90 b/test/ode_stats.f90
diff --git a/test/ode_workspace.f90 b/test/ode_workspace.f90
