# Modern Fortran interface to ODEPACK

**URL:** <https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352>\
**Category:** Announcements\
**Created:** [September 19, 2022, 3:42pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352 "2022-09-19T15:42:47Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [September 19, 2022, 3:42pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/1 "2022-09-19T15:42:47Z")

</div>

I made a modern Fortran interface to the `LSODA` routine in ODEPACK: [GitHub - Nicholaswogan/odepack](https://github.com/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.

---

<div class="post-metadata">

**Author:** ![awvwgk](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/awvwgk/32/154_2.png) [@awvwgk](https://fortran-lang.discourse.group/u/awvwgk)\
**Post date:** [September 28, 2022, 8:39pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/2 "2022-09-28T20:39:43Z")

</div>

3 posts were split to a new topic: [Conditionally linking LAPACK with fpm](https://fortran-lang.discourse.group/t/conditionally-linking-lapack-with-fpm/4425)

---

<div class="post-metadata">

**Author:** ![HugoMVale](https://avatars.discourse-cdn.com/v4/letter/h/e19adc/32.png) [@HugoMVale](https://fortran-lang.discourse.group/u/HugoMVale)\
**Post date:** [October 3, 2022, 11:37am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/3 "2022-10-03T11:37:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 3, 2022, 3:06pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/4 "2022-10-03T15:06:41Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [October 5, 2022, 6:40pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/5 "2022-10-05T18:40:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 5, 2022, 7:29pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/6 "2022-10-05T19:29:52Z")

</div>

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.

```nohighlight
/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)

```

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [October 5, 2022, 10:45pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/7 "2022-10-05T22:45:41Z")

</div>

See this repo: [GitHub - jacobwilliams/odepack: Work in Progress to refactor and modernize the ODEPACK Library](https://github.com/jacobwilliams/odepack)

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!

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [October 5, 2022, 11:23pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/8 "2022-10-05T23:23:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![urbanjost](https://avatars.discourse-cdn.com/v4/letter/u/0ea827/32.png) [@urbanjost](https://fortran-lang.discourse.group/u/urbanjost)\
**Post date:** [October 6, 2022, 2:26am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/9 "2022-10-06T02:26:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [October 6, 2022, 3:13pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/10 "2022-10-06T15:13:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\
**Post date:** [October 7, 2022, 1:53am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/11 "2022-10-07T01:53:56Z")

</div>

> [@nicholaswogan](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 7, 2022, 4:47am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/12 "2022-10-07T04:47:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [October 7, 2022, 2:20pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/13 "2022-10-07T14:20:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![general\_rishkin](https://avatars.discourse-cdn.com/v4/letter/g/a9a28c/32.png) [@general\_rishkin](https://fortran-lang.discourse.group/u/general_rishkin)\
**Post date:** [October 13, 2022, 4:32am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/14 "2022-10-13T04:32:22Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 13, 2022, 5:57am UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/16 "2022-10-13T05:57:36Z")

</div>

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:

```fortran
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

```

---

<div class="post-metadata">

**Author:** ![general\_rishkin](https://avatars.discourse-cdn.com/v4/letter/g/a9a28c/32.png) [@general\_rishkin](https://fortran-lang.discourse.group/u/general_rishkin)\
**Post date:** [October 13, 2022, 2:20pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/17 "2022-10-13T14:20:27Z")

</div>

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(.....)`.

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [October 13, 2022, 2:34pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/18 "2022-10-13T14:34:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 13, 2022, 3:23pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/19 "2022-10-13T15:23:32Z")

</div>

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.

```fortran
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

```

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [October 13, 2022, 5:06pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/20 "2022-10-13T17:06:03Z")

</div>

> [@nicholaswogan](#):
>
> Maybe this is a better interface for the rhs function? Not sure.

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](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).

---

<div class="post-metadata">

**Author:** ![nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@nicholaswogan](https://fortran-lang.discourse.group/u/nicholaswogan)\
**Post date:** [October 13, 2022, 6:10pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352/21 "2022-10-13T18:10:27Z")

</div>

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

[Next page](https://fortran-lang.discourse.group/t/modern-fortran-interface-to-odepack/4352.md?page=2)