Hey everyone,
I am fairly new to Julia and am interested to do machine learning with it. So far I learnt how to code in Julia and how to use the Flux library. Yet, My goal is to use trained neural networks in a fortran code. I saw on the Julia page that it is possible to embed Julia code in C/C++, but I can’t manage to find the way to do it in fortran. Does anyone have an idea?
A cleaner way is to bind to ONNXRuntime from Fortran so you don’t need Julia runtime for inference, after some searches, it appears that this doesn’t exist..
It is much easier to interface the other way around. It is conceivable as a workflow to have a main julia code that calls your fortran routines?
In general Julia will be much friendlier for file manipulation, visualization, and the higher-level stuff, and usually one could let some lower level tasks to Fortran if there are nice routines available for one’s purpose.
Thank you very much for your answers, so I agree with the fact of using Julia for higher level stuff, yet I am working on a fortran code that already exists and we want to add some machine learning to replace some operations in the code. In order to use Julia as the main code to use machine learning algorithms, I would have to rewrite most of the existing optimized fortran code, which isn’t really what we want, hence my question of embedding Julia code in fortran.
Since there is no better answer, unfortunately, I would just point that probably you would need “only” a wrapper for the main Fortran program. You could convert the main fortran program into a subroutine, and call it from Julia. Or create a new “main” program in Julia and call all the subroutines of Fortran from it. Of course how hard is that depends a lot on the the complexity of the program.
The alternative, easy one, is just to use each code (Julia and Fortran) as just independent program to execute, and interface them by writing and reading stuff from and to files. Of course that’s only viable if the cost of this slow interface is not important.
One more option would be not embedding Julia, but calling Julia. Meaning you have a Julia server and some ways of communication between your Fortran and Julia programs (e.g. ZMQ).
What are your requirements for throughput and latency?
In Fortran, you call a C function (Interoperability with C (The GNU Fortran Compiler) ). So maybe you can call the jl_init function (and others) by declaring a interface with the bind c attribute in Fortran. The problem would be C macros and potentially C data structures not suppoted by the Fortran’s C binding. So maybe you need some small C wrapper around your Julia code, that you can call from Fortran.
I have never tried this myself, but I would be interested to hear about your approach.
There’s one other option, call some other language from Fortran, I suggest to Python.
I googled a bit, found this, not sure if it’s the best, seems good, see its blog post:
I don’t mean this necessarily as a substitute for Julia and Flux (though would be helpful for e.g. Tensorflow or PyTorch), but rather to then call from Python to Julia, with Pythoncall.jl (it’s also to call from Python, or use pyjulia). Learning to call to or from Python (for Julia and likely any for language) is very useful, likely to be helpful to know, also independently of what you’re trying to do.
This might be the easiest intermediate language, and then you don’t need to code anything for interop, just call a bit differently. And since calling to or from Python is in-process (and I suppose all three), it’s fast, most likely.
There are some issues with coding your own interop infrastructure, that I’m not to up-to-speed on. I just know you get rid of them for some minor inconvenience. If some of the code is multi-threaded, that might be an issue, but likely will be with and without this solution.
There are some other options for the middle language. I wouldn’t recommend e.g. Java (not in-process), Rust should work. It’s possible to call to Julia from it and from C++ (library available to make easier from C++ and another from Rust). Both would be an option to call from, but it’s likely easier to call to Python than them.
It’s not mandatory interacting with the C macros or data structures to embed Julia. The approach in GitHub - GunnarFarneback/DynamicallyLoadedEmbedding.jl: Embed Julia with dynamical loading of libjulia at runtime. uses neither of those and minimizes the number of libjulia functions called.
If I remember correctly, @sloede has a working prototype for this.
Thanks for pinging me, @ranocha!
Yes, we do indeed have a working prototype for calling Julia from Fortran. Since Fortran is not nice (some may call it horrible) when playing with others, I found that for me it was much easier to write a small shim library in C that itself is called from Fortran, and then have the C library interact with Julia.
You can fin examples for pure C-to-Julia usage and Fortran-to-C-to-Julia usage in this repository:
Please feel free to ping me here (or better yet, in the repository using an issue) if there are questions/problems.
For reference, here is some example code that I just updated for julia 1.13.
Previously (julia 1.11 and 1.12), I could call julia directly from Fortran without using C code but something changed in the declaration of global variables like jl_main_module, that you can no longer access them via from Fortran as:
type(c_ptr), bind(C, name="jl_main_module") :: jl_main_module
But it is possible to get these pointers using a small C wrapper importing all declaration from julia.h.
You can but I definitly won’t recommand for begginers :
The neural network is compiled into a native library, with a custom ReLU activation passed as a Fortran function pointer.
One small hurdle: KernelAbstractions 0.10 (PoCL backend) caused unresolved GPU symbols during linking. Downgrading to 0.9 fixed it. And if you wonder you won’t be able to use gpu’s yet.
the files :
main.f90 main.jl Manifest.toml Project.toml
environement:
[587475ba] Flux v0.16.15
⌃ [63c18a36] KernelAbstractions v0.9.45 ⚲
using Flux
struct CActivation <: Function
ptr::Ptr{Cvoid}
end
@inline (f::CActivation)(x::Float64) =
ccall(f.ptr, Cdouble, (Cdouble,), x)
function build_nn(ni::Int64, no::Int64, nh::Int64, act::Ptr{Cvoid})
Flux.f64(Flux.Chain(
Flux.Dense(ni => nh, CActivation(act)),
Flux.Dense(nh => no)
))
end
const NNType = typeof(build_nn(1, 1, 1, C_NULL))
const NETWORKS = Dict{Int64, NNType}()
const NEXT_ID = Ref{Int64}(0)
Base.@ccallable function make_nn(
ni::Int64, no::Int64, nh::Int64,
act::Ptr{Cvoid}
)::Int64
(ni > 0 && no > 0 && nh > 0 && act != C_NULL) ||
return 0
nn = build_nn(ni, no, nh, act)
id = (NEXT_ID[] += 1)
NETWORKS[id] = nn
return id
end
Base.@ccallable function run_nn(
id::Int64,
px::Ptr{Float64},
py::Ptr{Float64},
batchsize::Int64
)::Cint
haskey(NETWORKS, id) || return -1
(px != C_NULL && py != C_NULL && batchsize >= 0) ||
return -1
nn = NETWORKS[id]
ni = size(nn.layers[1].weight, 2)
no = size(nn.layers[2].weight, 1)
x = unsafe_wrap(Matrix{Float64}, px, (ni, batchsize))
y = unsafe_wrap(Matrix{Float64}, py, (no, batchsize))
copyto!(y, nn(x))
return 0
end
Base.@ccallable function free_nn(id::Int64)::Cvoid
delete!(NETWORKS, id)
return nothing
end
fortran :
program test_nn
use iso_c_binding
implicit none
interface
function make_nn(ni, no, nh, act) &
bind(C, name="make_nn") result(id)
import :: c_int64_t, c_funptr
integer(c_int64_t), value :: ni, no, nh
type(c_funptr), value :: act
integer(c_int64_t) :: id
end function
function run_nn(id, px, py, batchsize) &
bind(C, name="run_nn") result(status)
import :: c_int64_t, c_ptr, c_int
integer(c_int64_t), value :: id, batchsize
type(c_ptr), value :: px, py
integer(c_int) :: status
end function
subroutine free_nn(id) bind(C, name="free_nn")
import :: c_int64_t
integer(c_int64_t), value :: id
end subroutine
end interface
integer(c_int64_t), parameter :: ni = 2, no = 1
integer(c_int64_t), parameter :: nh = 4, batch = 3
real(c_double), target :: x(ni,batch), y(no,batch)
integer(c_int64_t) :: id
integer(c_int) :: status
x(:,1) = [1.0d0, 2.0d0]
x(:,2) = [3.0d0, 4.0d0]
x(:,3) = [5.0d0, 6.0d0]
id = make_nn(ni, no, nh, c_funloc(relu))
if (id == 0) error stop "make_nn failed"
status = run_nn(id, c_loc(x), c_loc(y), batch)
if (status /= 0) error stop "run_nn failed"
print *, "Inputs:"
print *, x
print *, "Outputs:"
print *, y
call free_nn(id)
contains
function relu(v) bind(C) result(res)
real(c_double), value :: v
real(c_double) :: res
res = max(v, 0.0d0)
end function
end program test_nn
juliac --output-lib libcarre --project=. main.jl
then
gfortran main.f90 -L. -lcarre -o test_nn
finally
./test_nn
Inputs:
1.0000000000000000 2.0000000000000000 3.0000000000000000 4.0000000000000000 5.0000000000000000 6.0000000000000000
Outputs:
-1.0558196649046820 -1.9435691153369969 -2.8313185657693118
edit : oh my thats a 4 years revive of the thread lol.
version Lux.jl:
using Lux, Random
struct CActivation <: Function
ptr::Ptr{Cvoid}
end
@inline (f::CActivation)(x::Float64) =
ccall(f.ptr, Cdouble, (Cdouble,), x)
const PS = @NamedTuple{
weight::Matrix{Float64},
bias::Vector{Float64}
}
struct Network{L1,L2}
l1::L1
l2::L2
p1::PS
p2::PS
ni::Int64
no::Int64
end
const L1 = typeof(Lux.Dense(1 => 1, CActivation(C_NULL)))
const L2 = typeof(Lux.Dense(1 => 1))
const NNType = Network{L1,L2}
const NETWORKS = Dict{Int64,NNType}()
const NEXT_ID = Ref{Int64}(0)
function build_nn(ni::Int64, no::Int64, nh::Int64,
act::Ptr{Cvoid})
rng = Random.Xoshiro(0)
l1 = Lux.Dense(ni => nh, CActivation(act))
l2 = Lux.Dense(nh => no)
p1 = (
weight = randn(rng, Float64, nh, ni) ./ sqrt(ni),
bias = zeros(Float64, nh)
)
p2 = (
weight = randn(rng, Float64, no, nh) ./ sqrt(nh),
bias = zeros(Float64, no)
)
return Network(l1, l2, p1, p2, ni, no)
end
Base.@ccallable function make_nn(
ni::Int64, no::Int64, nh::Int64,
act::Ptr{Cvoid}
)::Int64
(ni > 0 && no > 0 && nh > 0 && act != C_NULL) ||
return 0
nn = build_nn(ni, no, nh, act)
id = (NEXT_ID[] += 1)
NETWORKS[id] = nn
return id
end
Base.@ccallable function run_nn(
id::Int64,
px::Ptr{Float64},
py::Ptr{Float64},
batchsize::Int64
)::Cint
haskey(NETWORKS, id) || return -1
(px != C_NULL && py != C_NULL && batchsize > 0) ||
return -1
nn = NETWORKS[id]
x = unsafe_wrap(Matrix{Float64}, px, (nn.ni, batchsize))
y = unsafe_wrap(Matrix{Float64}, py, (nn.no, batchsize))
h, _ = nn.l1(x, nn.p1, NamedTuple())
result, _ = nn.l2(h, nn.p2, NamedTuple())
copyto!(y, result)
return 0
end
Base.@ccallable function free_nn(id::Int64)::Cvoid
delete!(NETWORKS, id)
return nothing
end
with trimming:
% juliac --output-lib libcarre --project=. main.jl --trim
✓ Compiling...
% gfortran main.f90 -L. -lcarre -o test_nn
% ./test_nn
Inputs:
1.0000000000000000 2.0000000000000000 3.0000000000000000 4.0000000000000000 5.0000000000000000 6.0000000000000000
Outputs:
0.11682082443574934 0.26708492893181329 0.41734903342787721