# 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
• 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

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
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
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
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
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: https://godbolt.org/z/a698n1E4b.

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?