Add optional solver statistics output (nfev, njev, nlu) to stiff3 (#20)

* Add optional stiff3 statistics counters and test

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/5ea33038-53d1-472f-8bea-c63bb8576909

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

* Refine stats test robustness and initialize optional stats

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/5ea33038-53d1-472f-8bea-c63bb8576909

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

* Address validation feedback on stats implementation

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/5ea33038-53d1-472f-8bea-c63bb8576909

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

* Strengthen stats test with callback count verification

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/5ea33038-53d1-472f-8bea-c63bb8576909

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>

Author: Copilot

README.md

@@ -10,6 +10,7 @@ This repository provides a refactored version with a simplified procedural inter
 
 - Third-order Rosenbrock-type (semi-implicit Runge-Kutta) integrator suitable for stiff ODE systems
 - Adaptive stepsize control with error tolerance per component
+- Optional runtime statistics output: rhs evaluations (`nfev`), Jacobian evaluations (`njev`), and LU decompositions (`nlu`)
 - Requires an exact user-supplied Jacobian
 - Depends on BLAS and LAPACK for linear algebra operations
 - Supports two build systems: [CMake](https://cmake.org/) and [Fortran Package Manager (fpm)](https://github.com/fortran-lang/fpm)
@@ -77,7 +78,7 @@ 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
+  integer :: irtrn, stats(3)
 
 ! initial value
   y = [1.0_wp, 1.0_wp]
@@ -93,7 +94,8 @@ program vanpol
   irtrn = 0
   call out(0,x0,x0,y,0,0.0_wp,irtrn)
 ! integrate system of ODEs
-  call stiff3(n,fun,jac,x0,x1,h0,eps,w,y,solout=out)
+  call stiff3(n,fun,jac,x0,x1,h0,eps,w,y,solout=out,stats=stats)
+  print '(A,3(I0,1X))', 'nfev njev nlu: ', stats
 
 contains
 

src/stiff3_solver.f90

@@ -74,7 +74,7 @@ subroutine output_sub(nr,xold,x,y,iha,qa,irtrn)
 
   !> Semi-implicit Runge-Kutta integrator routine
   !
-  subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
+  subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout,stats)
     integer, intent(in) :: n
       !! Number of equations to be integrated.
     procedure(rhs_sub) :: fun
@@ -95,6 +95,8 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
       !! dependent variables at `x1`.
     procedure(output_sub), optional :: solout
       !! User supplied subprogram for output.
+    integer, intent(out), optional :: stats(3)
+      !! Statistics array with `[nfev, njev, nlu]`.
 
     real(wp), dimension(n) :: yk1, yk2, ya, yold, yold1, f, fold
       !! Workspace for solution vector and right-hand side
@@ -104,6 +106,7 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
       !! Workspace for the pivot array
 
     integer :: icon, iha, i, j, nr, irtrn
+    integer :: nfev, njev, nlu
     real(wp) :: x, xold, h, e, es, q, qa
 
   ! icon = 0 except for last step which ends exactly at x1
@@ -112,6 +115,9 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
     nr = 0
     x = x0
     h = h0
+    nfev = 0
+    njev = 0
+    nlu = 0
 
     outer: do
 
@@ -131,7 +137,9 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
     ! evaluate function and jacobian
 
       call fun(n,y,f)
+      nfev = nfev + 1
       call dfun(n,y,df)
+      njev = njev + 1
 
     ! keep values which are used in half-step integration
 
@@ -145,7 +153,7 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
 
     ! perform full integration step
 
-      call sirk3(n,fun,ip,f,y,yk1,yk2,df,2*h)
+      call sirk3(n,fun,ip,f,y,yk1,yk2,df,2*h,nfev,nlu)
 
       do i = 1, n
         ya(i) = y(i)
@@ -163,13 +171,15 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
       inner: do
         iha = iha + 1
 
-        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h)
+        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu)
         call fun(n,y,f)
+        nfev = nfev + 1
         call dfun(n,y,df)
+        njev = njev + 1
 
         yold1 = y
 
-        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h)
+        call sirk3(n,fun,ip,f,y,yk1,yk2,df,h,nfev,nlu)
 
       ! half step integration finished
       ! compute deviation and compare with tolerance
@@ -223,6 +233,7 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
         call solout(nr,xold,x,y,iha,qa,irtrn)
         if (irtrn < 0) then
           h0 = h
+          if (present(stats)) stats = [nfev, njev, nlu]
           return
         end if
       end if
@@ -231,6 +242,7 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
 
       if (icon == 1) then
         h0 = h
+        if (present(stats)) stats = [nfev, njev, nlu]
         return
       end if
 
@@ -241,7 +253,7 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout)
 
   !> Single-step semi-implicit integration
   !
-  subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h)
+  subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h,nfev,nlu)
     integer, intent(in) :: n
       !! Size of the system of ODEs
     procedure(rhs_sub) :: fun
@@ -260,6 +272,10 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h)
       !! On output contains the factorized matrix (I - h a J) = LU
     real(wp), intent(in) :: h
       !! Step size of the independent variable
+    integer, intent(inout) :: nfev
+      !! Number of right-hand side evaluations
+    integer, intent(inout) :: nlu
+      !! Number of LU decompositions
 
     integer :: i
 
@@ -281,13 +297,15 @@ subroutine sirk3(n,fun,ipiv,f,y,yk1,yk2,df,h)
     ! perform triangular decomposition and evaluate k1
     !
     call lu(df,ipiv)
+    nlu = nlu + 1
     call back(df,f,ipiv)
 
     do i = 1, n
       yk1(i) = h*f(i)
       yk2(i) = y(i) + 0.75_wp * yk1(i)
     end do
     call fun(n,yk2,f)
+    nfev = nfev + 1
     call back(df,f,ipiv)
 
     !

test/CMakeLists.txt

@@ -1,4 +1,4 @@
-foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants)
+foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants ode_stats)
   add_executable(${test_program} ${test_program}.f90)
   target_link_libraries(${test_program} PRIVATE stiff3)
   add_test(NAME ${test_program} COMMAND ${test_program})

test/ode_stats.f90

@@ -0,0 +1,65 @@
+module ode_stats_problem
+
+  use stiff3_solver, only: wp => stiff3_wp
+
+  implicit none
+
+  real(wp), parameter :: decay_rate = 2.0_wp
+  integer :: fun_calls = 0
+  integer :: jac_calls = 0
+
+contains
+
+  subroutine fun(n,y,f)
+    integer, intent(in) :: n
+    real(wp), intent(in) :: y(n)
+    real(wp), intent(inout) :: f(n)
+    fun_calls = fun_calls + 1
+    f(1) = -decay_rate*y(1)
+  end subroutine
+
+  subroutine jac(n,y,df)
+    integer, intent(in) :: n
+    real(wp), intent(in) :: y(n)
+    real(wp), intent(inout) :: df(n,n)
+    jac_calls = jac_calls + 1
+    df(1,1) = -decay_rate
+  end subroutine
+
+end module
+
+program ode_stats
+
+  use stiff3_solver, only: stiff3, wp => stiff3_wp
+  use ode_stats_problem, only: fun, jac, fun_calls, jac_calls
+
+  implicit none
+
+  integer, parameter :: n = 1
+  real(wp) :: y(n), w(n), x0, x1, h0, eps
+  integer :: stats(3)
+
+  y = [1.0_wp]
+  w = 1.0_wp
+  x0 = 0.0_wp
+  x1 = 1.0_wp
+  h0 = 0.01_wp
+  eps = 1.0e-10_wp
+  fun_calls = 0
+  jac_calls = 0
+
+  call stiff3(n,fun,jac,x0,x1,h0,eps,w,y,stats=stats)
+
+  if (stats(1) /= fun_calls .or. stats(2) /= jac_calls) then
+    print '(A,3(I0,1X),A,2(I0,1X))', &
+      'stats mismatch [nfev njev nlu]: ', stats, &
+      ' callback counts [fun jac]: ', fun_calls, jac_calls
+    error stop 1
+  end if
+
+  if (stats(3) <= 0) then
+    print '(A,3(I0,1X))', 'expected positive nlu in stats [nfev njev nlu], got: ', stats
+    error stop 1
+  end if
+
+end program