Is there a good Fortran package for SDE, DDE and ODE?

1

Hello all,

I just wanted to ask, has anyone have used, any good Fortran package/libraries or packages or subroutines, that contains

  1. stochastic differential equation (SDE) solver,
  2. delayed differential equation (DDE) solver, and
  3. ordinary differential equations (ODE) solver?

If so how is the performance compared with Julia’s differential equation solver? Is Julia’s differential equation solver is best on the market now?

I just wanted to find good Fortran differential solver so that I can use it in my current Fortran code.
I am not going to completely translate my Fortran code to Julia because based on my experience it seems that, it is not very easy to write a Julia code whose speed is the same as Fortran.

I am basically asking has anyone used any of the three (SDE/ODE/DDE) Fortran subroutines or solvers or libraries, and have any recommendations.

Thank you very much!

2 Likes
3

I am basically asking has anyone used any of the three (SDE/ODE/DDE) Fortran subroutines or solvers or libraries, and have any recommendations if any.

4

I wrote my own for a specific equation:

It’s faster than other codes I tried. But if somebody has a faster solver, definitely let me know.

2 Likes
5

The ones I am aware of for delay differential equations are:

For other types of ODE’s the page of Hairer includes several codes for both stiff and non-stiff problems. I believe many of them are heavily used across various fields.

Another classic resource are the ODE solvers published at netlib. Unfortunately, most of them are F77 style codes which are not thread-safe. An exception is the rksuite which has a Fortran 90 interface and was later modified and included into the NAG library.

The blog post from Christopher Rackauckas (also listed above by @kargl) lists a few more Fortran codes including ODEPACK, Sundials, and FATODE. Many of these are still heavily used today. Some of the solvers by Hairer and from ODEPACK are wrapped in scipy.integrate.

Concerning performance and design, as Christopher states:

DifferentialEquations.jl does offer pretty much everything from the other suite combined, but that’s no accident: our software organization came last and we used these suites as a guiding hand for how to design ours.

Some recent papers like this one found performance of Julia slightly better, but Fortran was comparable. Another paper found Fortran to be faster for a Runge-Kutta method.

As far as stiff problems are concerned, many of the older Fortran packages don’t use LAPACK for the matrix factorizations but their own factorization routines. I believe with some small refactoring performance could be improved. Still it would take a prodigious effort to reach the same level of comprehensiveness that has gone into DifferentialEquations.jl. The author of that package lives and breathes differential equation solvers.

3 Likes
6

Is any of the code for the second paper publicly available? It appears the paper isn’t/

1 Like
7

The code is made available on bitbucket: fortran vs julia.

2 Likes
8

Interestingly, the Julia code isn’t valid Julia.

function RKp6n1(...)
    ...
    local const a21 = 1/2
    ...
end

This is from RK.jl, and isn’t valid Julia syntax. That’s not entirely reassuring that the authors actually know what they’re doing.

3 Likes
9

Also, their code is incredibly inefficient. They wrote their test function in Julia to return a Vector of 2 values, while in Fortran it modifies the inputs. Had the Julia version simply returned 2 values (as a Tuple), they would have removed all of the garbage they were generating, and their function would have been about 10x faster.

2 Likes
11

Just to be on the clear side, I was only quoting part of their conclusion. They do also make some remarks in favor of Julia. But perhaps the rebuttal of that paper can go to a new thread or the Julia Discourse.

To return to the original question, I have personally used the following ODE solvers:

In a distant past I also used LSODE from ODEPACK which gave good results, but the API is rather cumbersome. Another problem with many of the Fortran ODE solvers is they are not thread-safe.

4 Likes
12

I remember Sourcery Institute also had a nice differential equation solver library in modern Fortran. I searched their GitHub pages, but could not find it immediately. Perhaps @rouson could help more on this.

2 Likes
13

You’re probably referring to FOODiE by Zhagi, Curcic et al.:

https://github.com/Fortran-FOSS-Programmers/FOODIE

I don’t think it was ever under Sourcery Institute, but it uses the Abstract Calculus pattern from the book I co-authored and I helped come up with the name.

Damian

3 Likes
14

That explains why I could not find it on Sourcery’s page. Thanks for your correction.

1 Like
15

Others:

There are certainly a lot of ODE solvers floating around. A good stdlib project would be to modernize, clean up and organize them all into a nice modern API.

6 Likes
16

The FLINT page says

For performance comparison, the latest FLINT code is tested against Julia’s DifferentialEquations package (https://docs.sciml.ai/release-2.0/index.html) and FLINT appears to be faster by at least an order of magnitude as shown in the following screenshot. The Julia test code along with results are provided in the media folder on the FLINT’s GitHub repository GitHub - princemahajan/FLINT: Fortran Library for numerical INTegration of differential equations.

4 Likes
17

We should add FLINT to the packages here: Fortran Packages - Fortran Programming Language and the others too!

2 Likes
18

pull request with some of them here: added 2 interpolation and 4 ode packages by jacobwilliams · Pull Request #306 · fortran-lang/fortran-lang.org ·

3 Likes
19

Thank you very much indeed guys!

My purpose is to see if we can find corresponding Fortran packages/subroutines/libraries, which can be combined and have the same functions as in the Julia’s DifferentialEquations.jl

For now I see some ODE, DDE recommendations, and may I ask, do you have any recommendations for stochastic differential equation (SDE) solver?

1 Like
20

Note that this benchmark is really flawed. I believe the results they see come from the fact that their Julia code sucks in ways their Fortran doesn’t. For example, the Julia version defines

function JacobiConstant(odesol)
    
    X = odesol[1:3,:]
    Xdot = odesol[4:6,:]

    ~,n = size(X)
    
    # distances
    R1 = copy(X)
    R2 = copy(X)
    R1[1,:] = R1[1,:] .+ mu
    R2[1,:] = R2[1,:] .- 1 .+ mu
    r1 = [norm(R1[:,i]) for i = 1:n]
    r2 = [norm(R2[:,i]) for i = 1:n]
    
    # radius and velocity nagnitudes
    R =  [sqrt(sum((X.^2)[:,i])) for i = 1:n]
    V = [sqrt(sum((Xdot.^2)[:,i])) for i = 1:n]

    # pseudo-potential
    Omega = 1/2*(R.^2) .+ (1-mu)./r1 .+ mu./r2
    
    # Jacobi's Constant
    C = (V.^2)/2 .- Omega
    
    return C
end;

While in Fortran, they define

function JacobiC(mu, X)
        real(wp), intent(in) :: mu
        real(wp), dimension(:,:), intent(in) :: X
        real(wp), dimension(size(X,2)) :: JacobiC
        real(wp), dimension(size(X,2)) :: Omega
        real(WP), dimension(size(X,2)) :: r1, r2
        
        r1 = sqrt((X(1,:) + mu)**2 + X(2,:)**2 + X(3,:)**2)
        r2 = sqrt((X(1,:) - 1.0 + mu)**2 + X(2,:)**2+X(3,:)**2)
        Omega = 1.0/2.0*sum(X(1:3,:)**2,1) + (1.0-mu)/r1 + mu/r2
        ! Jacobian's Constant
        JacobiC = sum(X(4:6,:)**2,1)/2 - Omega
    end function JacobiC

The Julia version here is incredibly inefficient because when you do something like [sqrt(sum((X.^2)[:,i])) for i = 1:n], you are computing X.^2 n times only to throw away all but 1 column, and allocating n matrices and n vectors to throw those away by summing over them. Just using @views appropriately, and doing less stupid stuff should easily regain the performance defecit here. Since this is the core of the benchmark, it’s kind of hard for me to believe that the authors are trying to make a fair comparison, when the Julia code is so obviously stupidly written (for this function, they would have gotten way better performance by just following the same logic that their fortran version used, instead of randomly computing a ton of info, allocating a ton of arrays, and then throwing them away).

21

Maybe it is something about Julia that prevents people from writing efficient code. I doubt so many people and smart scientists are so “stupid” to write inefficient “stupid stuff” so easily in a language as you describe.

1 Like
22

They litterally could have written the exact same thing, except with (X[:,i].^2)) instead of (X.^2)[:,i]) and it would have been multiple times faster.