Add stiff3 workspace overloads with shared core implementation (#25)

Copilot · ivan-pi · web-flow · commit d63052363135 · 2026-05-08T15:18:39.000+02:00
* Add stiff3 workspace overload and core routine Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/a3bd8265-d95a-447c-96ff-c8de05acbfbd Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Improve ode_workspace mismatch diagnostics Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/a3bd8265-d95a-447c-96ff-c8de05acbfbd Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Apply review feedback for workspace argument mapping Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/0bfd1f1c-6066-4bb0-9df3-d8078e1d0164 Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Polish workspace argument indexing in stiff3_work Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/0bfd1f1c-6066-4bb0-9df3-d8078e1d0164 Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com> * Simplify workspace offset expressions in stiff3_work Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/0bfd1f1c-6066-4bb0-9df3-d8078e1d0164 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
@@ -12,6 +12,7 @@ This repository provides a refactored version with a simplified procedural inter
 - 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 explicit workspace interface via `stiff3(..., rwork, iwork, ...)` for caller-managed memory
 - 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)
@@ -151,6 +152,8 @@ 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`.
+
 You can also optionally pass `hmax` to limit the absolute half-step size used by the adaptive controller. If `hmax` is omitted or set to `0`, the default is `abs(x1-x0)`. If provided and positive, the solver uses `min(hmax, abs(x1-x0))`. Negative values are rejected.
 
 **How the tolerance is applied:** after each step the solver estimates the local error in each component and checks:
diff --git a/src/stiff3_solver.f90 b/src/stiff3_solver.f90
@@ -16,7 +16,7 @@ module stiff3_solver
   implicit none
   private
 
-  public :: stiff3, stiff3_wp
+  public :: stiff3, stiff3_auto, stiff3_work, stiff3_wp
   public :: rhs_sub, jacobian_sub, output_sub
 
   !> Kind parameter for working precision of stiff3 reals
@@ -66,6 +66,11 @@ subroutine output_sub(nr,xold,x,y,iha,qa,irtrn)
 
   end interface
 
+  interface stiff3
+    module procedure stiff3_auto
+    module procedure stiff3_work
+  end interface
+
 
 contains
 
@@ -74,7 +79,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,stats,hmax)
+  subroutine stiff3_auto(n,fun,dfun,x0,x1,h0,eps,w,y,solout,stats,hmax)
     integer, intent(in) :: n
       !! Number of equations to be integrated.
     procedure(rhs_sub) :: fun
@@ -108,6 +113,81 @@ subroutine stiff3(n,fun,dfun,x0,x1,h0,eps,w,y,solout,stats,hmax)
     integer :: ip(n)
       !! Workspace for the pivot array
 
+    call stiff3_core(n,fun,dfun,x0,x1,h0,eps,w,y, &
+                     yk1,yk2,ya,yold,yold1,f,fold,df,dfold,ip, &
+                     solout,stats,hmax)
+
+  end subroutine
+
+
+  !> Semi-implicit Runge-Kutta integrator routine with explicit workspace
+  !
+  subroutine stiff3_work(n,fun,dfun,x0,x1,h0,eps,w,y,rwork,iwork,solout,stats,hmax)
+    integer, intent(in) :: n
+      !! Number of equations to be integrated.
+    procedure(rhs_sub) :: fun
+      !! User supplied subprogram for function evaluation.
+    procedure(jacobian_sub) :: dfun
+      !! User supplied subprogram for evaluation of the Jacobian.
+    real(wp), intent(in) :: x0, x1
+      !! Limits of the independent variable between which the differential
+      !! equation is solved.
+    real(wp), intent(inout) :: h0
+      !! Suggested initial half-step length. On exit `h0` contains suggested
+      !! value of half-step length for continued integration beyond `x1`,
+      !! or the suggested step size at an interrupted callback return.
+    real(wp), intent(in) :: eps, w(n)
+      !! Tolerance parameters.
+    real(wp), intent(inout) :: y(n)
+      !! Vector of dependent variables at `x0`. On exit `y` is the vector of
+      !! dependent variables at `x1`.
+    real(wp), intent(inout) :: rwork(n*(7 + 2*n))
+      !! Real workspace of size `n*(7 + 2*n)`.
+    integer, intent(inout) :: iwork(n)
+      !! 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]`.
+    real(wp), intent(in), optional :: hmax
+      !! Maximum absolute half-step size. If absent or zero, defaults to
+      !! `abs(x1 - x0)`.
+
+    call stiff3_core(n,fun,dfun,x0,x1,h0,eps,w,y, &
+                     yk1   = rwork(1), &
+                     yk2   = rwork(n+1), &
+                     ya    = rwork(2*n+1), &
+                     yold  = rwork(3*n+1), &
+                     yold1 = rwork(4*n+1), &
+                     f     = rwork(5*n+1), &
+                     fold  = rwork(6*n+1), &
+                     df    = rwork(7*n+1), &
+                     dfold = rwork(7*n+n*n+1), &
+                     ip    = iwork(1), &
+                     solout=solout,stats=stats,hmax=hmax)
+
+  end subroutine
+
+
+  !> Semi-implicit Runge-Kutta integrator routine core implementation
+  !
+  subroutine stiff3_core(n,fun,dfun,x0,x1,h0,eps,w,y, &
+                         yk1,yk2,ya,yold,yold1,f,fold,df,dfold,ip, &
+                         solout,stats,hmax)
+    integer, intent(in) :: n
+    procedure(rhs_sub) :: fun
+    procedure(jacobian_sub) :: dfun
+    real(wp), intent(in) :: x0, x1
+    real(wp), intent(inout) :: h0
+    real(wp), intent(in) :: eps, w(n)
+    real(wp), intent(inout) :: y(n)
+    real(wp), intent(inout) :: yk1(n), yk2(n), ya(n), yold(n), yold1(n), f(n), fold(n)
+    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)
+    real(wp), intent(in), optional :: hmax
+
     integer :: icon, iha, i, j, nr, irtrn
     integer :: nfev, njev, nlu
     real(wp) :: x, xold, h, e, es, q, qa, hmax_used
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -1,4 +1,4 @@
-foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants ode_stats ode_hmax)
+foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants ode_stats ode_hmax ode_workspace)
   add_executable(${test_program} ${test_program}.f90)
   target_link_libraries(${test_program} PRIVATE stiff3)
   add_test(NAME ${test_program} COMMAND ${test_program})
diff --git a/test/ode_workspace.f90 b/test/ode_workspace.f90
@@ -0,0 +1,62 @@
+program ode_workspace
+
+  use stiff3_solver, only: stiff3, wp => stiff3_wp
+
+  implicit none
+
+  integer, parameter :: n = 2
+  real(wp) :: y_auto(n), y_work(n), w(n), x0, x1, h0, eps
+  real(wp), allocatable :: rwork(:)
+  integer, allocatable :: iwork(:)
+  integer :: stats_auto(3), stats_work(3)
+
+  y_auto = [1.0_wp, -0.5_wp]
+  y_work = y_auto
+  w = 1.0_wp
+  x0 = 0.0_wp
+  x1 = 1.0_wp
+  h0 = 0.02_wp
+  eps = 1.0e-8_wp
+
+  allocate(rwork(n*(7 + 2*n)))
+  allocate(iwork(n))
+
+  call stiff3(n,fun,jac,x0,x1,h0,eps,w,y_auto,stats=stats_auto)
+
+  h0 = 0.02_wp
+  call stiff3(n,fun,jac,x0,x1,h0,eps,w,y_work,rwork,iwork,stats=stats_work)
+
+  if (any(y_auto /= y_work)) then
+    print '(A,2(1X,ES12.4),A,2(1X,ES12.4))', &
+      'solutions differ. auto=', y_auto, ' work=', y_work
+    error stop 1
+  end if
+
+  if (any(stats_auto /= stats_work)) then
+    print '(A,3(I0,1X),A,3(I0,1X))', &
+      'stats mismatch auto [nfev njev nlu]: ', stats_auto, &
+      ' work: ', stats_work
+    error stop 1
+  end if
+
+contains
+
+  subroutine fun(n,y,f)
+    integer, intent(in) :: n
+    real(wp), intent(in) :: y(n)
+    real(wp), intent(inout) :: f(n)
+    f(1) = -2.0_wp*y(1) + y(2)
+    f(2) = -3.0_wp*y(2)
+  end subroutine
+
+  subroutine jac(n,y,df)
+    integer, intent(in) :: n
+    real(wp), intent(in) :: y(n)
+    real(wp), intent(inout) :: df(n,n)
+    df(1,1) = -2.0_wp
+    df(1,2) = 1.0_wp
+    df(2,1) = 0.0_wp
+    df(2,2) = -3.0_wp
+  end subroutine
+
+end program

Original file line number	Diff line number	Diff line change
`@@ -1,4 +1,4 @@`
`1`		`-foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants ode_stats ode_hmax)`
	`1`	`+foreach(test_program ode_linear_decay ode_order ode_logistic ode_stiff_linear ode_robertson_invariants ode_stats ode_hmax ode_workspace)`
`2`	`2`	`add_executable(${test_program} ${test_program}.f90)`
`3`	`3`	`target_link_libraries(${test_program} PRIVATE stiff3)`
`4`	`4`	`add_test(NAME ${test_program} COMMAND ${test_program})`