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

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.


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

1 Like

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


Interestingly, the Julia code isn’t valid Julia.

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

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.


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.


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.


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.


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

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.



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

1 Like


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.


The FLINT page says

For performance comparison, the latest FLINT code is tested against Julia’s DifferentialEquations package ( 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<.


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


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


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

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

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

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

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.

@oscardssmith I think if the FLINT author really compare FLINT with Julia’s differentialequation.jl package, then the result is what it is.

I can be wrong, but it is difficult to believe the people who wrote, test, and maintain the Julia’s differentialequation.jl package did not aware of the stuff you mentioned. If those can really help they will do so already. Perhaps you could also make your comments and suggestions to Julia Forum and to see how they respond.

I can be absolutely wrong. But based on my very shallow experience with Julia, no offense, I feel that with or without the Julia’s optimization tricks you suggested, the speed of the Julia code will not really change much. 10-30% faster is possible, but like you said easily 10x faster is unlikely. Because as long as they follow the performance tips, Performance Tips · The Julia Language
recent versions of Julia should be intelligent enough to optimize those `stupid’ part of the code you pointed out automatically. If Julia cannot handle those optimization automatically and requires the user to write code with that level of awareness and carefulness you suggested, it will be really funny and seems lost its purpose because not everyone is software engineer level.

Like even in Fortran, based on my shallow experience, with modern compilers and hardware, even if I wrote things like A( i, : ) instead of column major A(:, i ) in the loop, 10-20% slower is possible (actually sometimes the difference is even smaller and even negligible). But the code will highly unlikely slow down by 10 times. In Julia I believe things are similar, the compiler will do some optimization work automatically and the modern hardware is fast enough. Besides, it is difficult to always maintain column major everywhere in the code. Sometime you just have to have some A(i, : ) no matter what.

This isn’t about code in DifferentialEquations.jl this is about the benchmarking code that the FLINT authors used to benchmark Julia against their FLINT. They wrote their own model separately in fortran and in julia, so their inefficiency in their code for the Julia model isn’t something DifferentialEquations is able to optimize for.

Thank you! @oscardssmith
Eh, which jl test file in FLINT are you referring to, it seems I only find one jl file,

It seems the Julia script they wrote was not that bad.
Perhaps you could quickly do a test on you computer, with their original Julia script, and with your optimized script, to see if there is a 10x speedup?