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>

Author: Copilot

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`), LU decompositions (`nlu`), successful steps (`nacc`), rejected steps (`nrej`), and linear-system solves (`nsol`)
+- 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(6)
+  integer :: stats(6)
 
 ! initial value
   y = [1.0_wp, 1.0_wp]
@@ -109,13 +109,10 @@ 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(1), stats(2), stats(3)
-  print '(A,3(I0,1X))', 'accepted rejected nsol: ', stats(4), stats(5), stats(6)
+  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
 

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)
 

example/pendant_drop.f90

@@ -28,12 +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,stats=stats)
-  print '(A,3(I0,1X))', 'nfev njev nlu: ', stats(1), stats(2), stats(3)
-  print '(A,3(I0,1X))', 'accepted rejected nsol: ', stats(4), stats(5), stats(6)
+  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
 

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)
 

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

example/three_equation_system.f90

@@ -22,8 +22,8 @@ program three_equation_system
   x1 = 50.0_wp
 
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,stats=stats)
-  print '(A,3(I0,1X))', 'nfev njev nlu: ', stats(1), stats(2), stats(3)
-  print '(A,3(I0,1X))', 'accepted rejected nsol: ', stats(4), stats(5), stats(6)
+  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
 

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)
 

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
@@ -106,7 +106,7 @@ subroutine stiff3_auto(n,fun,x,y,xend,jac,h0,eps,w,solout,stats,hmax)
     procedure(output_sub), optional :: solout
       !! User supplied subprogram for output.
     integer, intent(out), optional :: stats(6)
-      !! Statistics array with `[nfev, njev, nlu, nacc, nrej, nsol]`.
+      !! 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)`.
@@ -154,7 +154,7 @@ subroutine stiff3_work(n,fun,x,y,xend,jac,h0,eps,w,rwork,iwork,solout,stats,hmax
     procedure(output_sub), optional :: solout
       !! User supplied subprogram for output.
     integer, intent(out), optional :: stats(6)
-      !! Statistics array with `[nfev, njev, nlu, nacc, nrej, nsol]`.
+      !! 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)`.
@@ -194,15 +194,14 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
     integer, intent(out), optional :: stats(6)
     real(wp), intent(in), optional :: hmax
 
-    integer :: icon, iha, i, j, nr, irtrn
+    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'
@@ -223,6 +222,16 @@ subroutine stiff3_core(n,fun,x,y,xend,jac,h0,eps,w, &
     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
@@ -342,14 +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, nacc, nrej, nsol]
+          if (present(stats)) stats = [nacc, nrej, nfev, njev, nlu, nsol]
           return
         end if
       end if
@@ -358,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, nacc, nrej, nsol]
+        if (present(stats)) stats = [nacc, nrej, nfev, njev, nlu, nsol]
         return
       end if
 

test/ode_hmax.f90

@@ -27,7 +27,7 @@ program ode_hmax
 
   if (any(stats_zero /= stats_default)) then
     print '(A,6(I0,1X),A,6(I0,1X))', &
-      'expected hmax=0 stats [nfev njev nlu nacc nrej nsol] to match default. got ', stats_zero, &
+      'expected hmax=0 stats [nacc nrej nfev njev nlu nsol] to match default. got ', stats_zero, &
       ' expected ', stats_default
     error stop 1
   end if

test/ode_stats.f90

@@ -50,25 +50,25 @@ program ode_stats
 
   call stiff3(n,fun,x0,y,x1,jac,h0,eps,w,stats=stats)
 
-  if (stats(1) /= fun_calls .or. stats(2) /= jac_calls) then
+  if (stats(3) /= fun_calls .or. stats(4) /= jac_calls) then
     print '(A,6(I0,1X),A,2(I0,1X))', &
-      'stats mismatch [nfev njev nlu nacc nrej nsol]: ', stats, &
+      'stats mismatch [nacc nrej nfev njev nlu nsol]: ', stats, &
       ' callback counts [fun jac]: ', fun_calls, jac_calls
     error stop 1
   end if
 
-  if (stats(3) <= 0) then
-    print '(A,6(I0,1X))', 'expected positive nlu in stats [nfev njev nlu nacc nrej nsol], got: ', stats
+  if (stats(5) <= 0) then
+    print '(A,6(I0,1X))', 'expected positive nlu in stats [nacc nrej nfev njev nlu nsol], got: ', stats
     error stop 1
   end if
 
-  if (stats(4) <= 0 .or. stats(5) < 0) then
-    print '(A,6(I0,1X))', 'expected nonnegative step counters [nfev njev nlu nacc nrej nsol], got: ', stats
+  if (stats(1) <= 0 .or. stats(2) < 0) then
+    print '(A,6(I0,1X))', 'expected nonnegative step counters [nacc nrej nfev njev nlu nsol], got: ', stats
     error stop 1
   end if
 
-  if (stats(6) < stats(3)) then
-    print '(A,6(I0,1X))', 'expected nsol>=nlu in stats [nfev njev nlu nacc nrej nsol], got: ', stats
+  if (stats(6) < stats(5)) then
+    print '(A,6(I0,1X))', 'expected nsol>=nlu in stats [nacc nrej nfev njev nlu nsol], got: ', stats
     error stop 1
   end if