Modern Fortran interface to ODEPACK

I made a modern Fortran interface to the LSODA routine in ODEPACK: GitHub - Nicholaswogan/odepack . ODEPACK solves ordinary differential equation initial value problems. LSODA can deal with stiff and non-stiff problems by automatically switching between appropriate methods.

This repository modifies the original ODEPACK source code. The main difference is that I have replaced common blocks with a passed derived type. This makes the code thread safe

Features include

  • easy to use object-oriented interface
  • root-finding
  • thread safe
  • Jacobian can be dense or banded

I am considering replacing the linear algebra matrix routines, which are LINPACK, with calls to more highly optimized LAPACK routines. I’m guessing this would matter for big problems. Any thoughts on the value of this are appreciated.

Any other feedback/issue is also welcome.

14 Likes

3 posts were split to a new topic: Conditionally linking LAPACK with fpm

Thanks a lot for this updated interface, @nicholaswogan. I am big fan of odepack and it’s really nice that we can now employ this set of ODE solvers in a modularized and object-oriented way.
I noticed that you “only” (please don’t take this negatively) made interfaces to LSODA and LSODAR. Are you planning to extend the work to other solvers, e.g. LSODE? For problems that I know for sure are stiff, I prefer calling LSODE with the correct method flag, rather than letting LSODA discover that on its own.
Thanks again for this contribution. Really cool.

3 Likes

No problem! I might put the same effort into LSODE, but not in the near future. I need routines which have root finding, and unfortunately, LSODE doesn’t have root finding.

3 Likes

Are you willing to accept pull requests with the minimum refactorings needed to get rid of the -std=legacy flag? The computational kernel is quite a ball of spaghetti as it stands currently… The first step of course would be to get decent test coverage, and potentially also measure speed, binary size, etc.

3 Likes

Yes, definitely would accept that kind of pull request.

When I remove the -std=legacy flag, everything still compiles, but gives many warnings that look like this.

/Users/nicholas/Documents/Research_local/random/lsoda_numba/lsoda_refactored/src/odepack_sub1.f:123:72:

  123 |  110      PC(I) = PC(I-1) + FNQM1*PC(I)
      |                                                                        1
Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 110 at (1)
1 Like

See this repo: GitHub - jacobwilliams/odepack: Work in Progress to refactor and modernize the ODEPACK Library

For the beginnings of an attempt to actually refactor odepack. Mostly from @urbanjost who used the SPAG tool to get it started. It is very much a work in progress but if somebody wants to contribute we are willing to accept help!

3 Likes

I think that you may find more people willing to help if the refactored code were restricted to, say, F95+TR15581, with as few modules as possible in each project, and explicit interfaces being provided only when necessary.

With those restrictions, the refactored code could be compiled and debugged with Lahey/Fujitsu and Silverfrost, which are my tools of choice for debugging. Unfortunately, these same tools do not handle submodules, the availability of which can make reusing modules easier.

2 Likes

It is one module, if that helps. The individual procedures are in separate INCLUDE files for easier editing, which makes it appear otherwise. A few procedures not completely converted to be able to be in a module are still included into the single file at the end, which are part of the WIP aspects of the project.
INCLUDE was used to split up the source code as a module was desired but submodules, as you mentioned, are not universally supported.

There are a few uses of post-95 usage; The conversion package had options to restrict what the rewrite would write, so you could pick a standard level. Not sure if I added the newer features or if part of that came from the package without looking back through the history.

2 Likes

No I don’t agree. We shouldn’t limit ourselves to 30 year old variants of Fortran. That’s not how we get new users to Fortran. My goal in the modernizations that I do is to try and bring the code up to the present day (i.e., what would it look like if it were written today). The end result is code that we can show people not familiar with Fortran as an example of why Fortran should be used in the first place. Keeping the old FORTRAN 77 spaghetti and adding module wrappers to it is not how we should advocate to new users that code be written today. The algorithms can be liberated from the terrible fixed-form FORTRAN 77 that they have been imprisoned in for all these decades.

6 Likes

If there could be noticeable performance gain by doing this, it will be really awesome!

1 Like

I don’t think it matters for small problem. But it should matter when matrices gets large.

1 Like

I downloaded the package, and then built and ran some of the examples. Impressive piece of work!

There is one remaining anachronistic feature that got carried over from the old F77 code: passing a section of an integer array as an actual argument when the corresponding dummy argument was a double precision array. In a previous attempt to polish the Odepack code, I had stumbled on this very issue, and gave up on finding any simple solution, since some of these work arrays sometimes return status variables that are to be used elsewhere, packed in some delicate order that will cause the package to fail if meddled with.

3 Likes

@nicholaswogan : I wonder if you could modify the lsoda_class interface or the rhs function to accept an optional user-defined data/parameter array. This will allow users to pass necessary parameters to specify their rhs functions without resorting to global variables or having to extend the derived type.

Why not just extend the lsoda_class derived type? Or you can also just specify a function which is contained in another function. Here are the two ways I would do it:

module example
  use odepack_mod, only: lsoda_class, dp
  implicit none

  ! some data to be passed
  type :: MyData
    integer :: n
    real(dp), allocatable :: arr(:)
  end type

  ! customize the lsoda class
  type, extends(lsoda_class) :: lsoda_custom
    type(MyData), pointer :: dat => null()
  end type

contains

  subroutine test_1()
    type(lsoda_custom) :: ls
    integer :: neq, itask, istate
    real(dp) :: y(2), t, tout, rtol, atol(1)
    type(MyData), target :: dat

    ! custom data
    dat%n = 3
    allocate(dat%arr(dat%n))
    dat%arr(:) = 1.0_dp

    neq = 2
    call ls%initialize(rhs_1, neq, istate=istate)
    if (istate < 0) then
      error stop '"test" failed'
    endif
    ls%dat => dat

    y(:) = [5.0_dp, 0.8_dp]
    t = 0.0_dp
    tout = 10.0_dp
    rtol = 1.0e-8_dp
    atol = 1.0e-8_dp
    itask = 1
    istate = 1
    call ls%integrate(y, t, tout, rtol, atol, itask, istate)
    if (istate < 0) then
      error stop '"test" failed'
    endif

    print*,y

  end subroutine

  subroutine rhs_1(self, neq, t, y, ydot, ierr)
    class(lsoda_class), intent(inout) :: self
    integer, intent(in) :: neq
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(neq)
    real(dp), intent(out) :: ydot(neq)
    integer, intent(out) :: ierr

    type(MyData), pointer :: dat

    select type (self)
    class is (lsoda_custom)
      dat => self%dat
    end select

    ydot(1) = y(1)-y(1)*y(2)*sum(dat%arr)
    ydot(2) = y(1)*y(2)-y(2)
    ierr = 0
  end subroutine

  subroutine test_2()
    type(lsoda_class) :: ls
    integer :: neq, itask, istate
    real(dp) :: y(2), t, tout, rtol, atol(1)
    type(MyData) :: dat

    ! custom data
    dat%n = 3
    allocate(dat%arr(dat%n))
    dat%arr(:) = 1.0_dp

    neq = 2
    call ls%initialize(rhs_tmp, neq, istate=istate)
    if (istate < 0) then
      error stop '"test" failed'
    endif

    y(:) = [5.0_dp, 0.8_dp]
    t = 0.0_dp
    tout = 10.0_dp
    rtol = 1.0e-8_dp
    atol = 1.0e-8_dp
    itask = 1
    istate = 1
    call ls%integrate(y, t, tout, rtol, atol, itask, istate)
    if (istate < 0) then
      error stop '"test" failed'
    endif

    print*,y

  contains
    subroutine rhs_tmp(self_, neq_, t_, y_, ydot_, ierr_)
      class(lsoda_class), intent(inout) :: self_
      integer, intent(in) :: neq_
      real(dp), intent(in) :: t_
      real(dp), intent(in) :: y_(neq_)
      real(dp), intent(out) :: ydot_(neq_)
      integer, intent(out) :: ierr_
      call rhs_2(neq_, t_, y_, ydot_, ierr_, dat)
    end subroutine
  end subroutine

  subroutine rhs_2(neq, t, y, ydot, ierr, dat)
    integer, intent(in) :: neq
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(neq)
    real(dp), intent(out) :: ydot(neq)
    integer, intent(out) :: ierr
    type(MyData), intent(inout) :: dat
    ydot(1) = y(1)-y(1)*y(2)*sum(dat%arr)
    ydot(2) = y(1)*y(2)-y(2)
    ierr = 0
  end subroutine

end module

program main
  use example, only: test_1, test_2
  implicit none
  call test_1()
  call test_2()
end program

Extending the lsoda_class or using a sub-function are definitely options. I just believe that these put more burden on the user of the library (plus I am not a fan of inheritance in any language unless absolutely necessary).

It would be nice to be able to use the functionality just as is done in MATLAB, Scipy and Julia, i.e., just call odeint(.....).

When the ODE problem being solved involves parameters, there are some delicate questions that need to be answered first. Does the “solver” (in this case, one of the ODEPACK solvers such as LSODAR) need to know that there are parameters to be used? If not, the parameters p in \dot{y} = f(t, y; p) do not need to be passed to the solver routines in their argument lists. The user may place the parameters in a module (or COMMON blocks, in F77 and earlier) for sharing between the driver and the problem function routines (f, the Jacobian routine, etc.), or using polymorphic variable types, as discussed earlier in this thread.

The group led by Hindmarsh at LLNL also published the SUNDIALS package, which I recollect does include the capability for performing sensitivity analysis w.r.t. parameters. Other solvers address parameter optimization, as well.

In summary, the solution may have to be adjusted to fit the circumstances, and some planning will be needed.

I agree to some extent. I think “the right way” to do something should be very obvious. This will result in better code produced by the community.

Maybe this is a better interface for the rhs function? Not sure.

subroutine rhs(t, y, ydot, user_data, ierr)
  real(dp), intent(in) :: t
  real(dp), intent(in) :: y(:)
  real(dp), intent(out) :: ydot(size(y))
  class(*), intent(inout), optional :: user_data
  integer, intent(out) :: ierr
end subroutine

Since you need to use a select type block to restore the type of user_data, this adds a jump in the function evaluation to perform the type check: Compiler Explorer.

Furthermore the compilers add a small number of instructions to check for contiguity. I’d also point out that in your interface y could be non-contiguous, while ydot must be contiguous since it is assumed-size, which seems kind of inconsistent. For large functions (n >> 3) such overheads probably don’t really matter (assuming you are passing a contiguous dummy variable so no copy in/out is performed).

I’ve seen you post about the most ideal callback interface somewhere. What do you think is best in situations like this?