I made examples of wrapping jl code from Fortran & C; how generalize?

I am interested in calling Julia from Fortran (eg F90), and since I couldn’t find any, I built an example of that which goes through a C wrapper layer; see here:

It adapts code from the manual (Ch 30: Embedding) plus a SO post which explains how to call general J code. F is responsible for preallocating output arrays, in the classical F style. Multithreading works.

It works but feels too clunky, involving text commands and a lot of jl_* calls particular to the function’s signature. The thought of wrapping multiple Julia functions this way is daunting.

  1. Does anyone know if there are cleaner, more general such wrappers or examples? In particular, C/F calling J code.
  2. In the case of multiple calls to J, I imagine the F driver code should call jl_init() and load any J modules once at start, and jl_atexit_hook(0) once at finish, so that only one session is created. Is that right?
  3. I don’t understand where JULIA_DEFINE_FAST_TLS() should be called if the C code is merely a wrapper layer to F.

Anyway, I hope these are useful and others will improve them. Thanks, Alex

5 Likes

Hi Alex, welcome to Julia!

What I usually recommend here, rather than doing a zillion jl_* calls to pack/unpack arguments, is to just use @cfunction on the Julia side to get a compiled C function pointer for whatever signature you want on the C/Fortran side. Basically, you want to do all of the conversions on the Julia side where it is easy, and then on the C side you have an ordinary function pointer that you can call like any ordinary C function.

(I think you can use C function pointers directly from Fortran as well using C_F_PROCPOINTER?)

I posted an example of this on stackoverflow recently in answer to someone who wanted to call the Julia SpecialFunctions.polygamma function from C++.

I feel like this should be added to the Julia manual: explain @cfunction in Embedding Julia manual · Issue #38932 · JuliaLang/julia · GitHub

7 Likes

If it is of your interest I have some examples calling Fortran. They are just minimal interfaces to get started. Also, there are examples interfacing C, C++, Rust, and Python.

1 Like

Yes.

I think this works with some linker magic (it just defines a couple of functions with special compiler-attribute directives), so you only need to stick it once in any of your linked files. So yes, stick it in one of your C files, even if it is just being used as a wrapper layer from Fortran.

1 Like

Dear @stevengj, thanks for the help.

While able to get cfunction working with scalar I/O, I cannot for even one array. I have tried putting explicit types on the Julia side. No luck. Am I messing up the cfunction? I have added it to the repo from my original post, under ccallj/testcfuncfail.c. It’s a minimally complete C example which passes an array and prints the max value:

#include <julia.h>
typedef void (*test2_ptr)(double*, int);
test2_ptr test2 = NULL;
int main()
{
  jl_init();
  jl_eval_string("function test2(x,n) println(maximum(x[1:n])) end");
  jl_eval_string("x=rand(10); test2(x,6)");  // test it from jl side 
  jl_value_t *ret = jl_eval_string("@cfunction(test2, Cvoid, (Ref{Cdouble}, Cint))");
  test2 = (test2_ptr) jl_unbox_voidpointer(ret);     // convert to C func ptr
  double a[4] = {1.0,2.0,7.0,3.0};             // input data
  test2(a,4);
  jl_atexit_hook(0);
  return 0;
}

Compilation links to -ljulia as usual (see repo). Execution gives:

0.8654571660294361
fatal: error thrown and no exception handler available.
MethodError(f=typeof(Base.getindex)(), args=(1, Base.UnitRange{Int64}(start=1, stop=4)), world=0x0000000000006ca2)
jl_method_error_bare at /buildworker/worker/package_linux64/build/src/gf.c:1798
jl_method_error at /buildworker/worker/package_linux64/build/src/gf.c:1816
jl_lookup_generic_ at /buildworker/worker/package_linux64/build/src/gf.c:2373 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:2394
test2 at ./none:1
unknown function (ip: 0x7fbaf0073b73)
main at ./testcfuncfail (unknown line)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
_start at ./testcfuncfail (unknown line)

Maybe the size of x is not understood on the jl side? Maybe I have the wrong paradigm x[1:n] for passing an array with length controlled on the C side? But I can’t find any example codes with arrays.
Being able to pass an array from C to Julia is pretty essential for me (and if cfunction is the best way to wrap). Could someone kindly help me fix this example?

Thanks so much, Alex

Ok, letting the student (ie, me) play around, I seem to have figured it out! The Julia function has to repackage the incoming pointer into an array. Here is a complete working example:

#include <julia.h>
typedef void (*myfun_ptr)(double*, int);     // ptr type for func's C signature
myfun_ptr myfun = NULL;                      // instantiate a ptr
int main()
{
  jl_init();
  jl_eval_string("function myfun(xptr,n) x=unsafe_wrap(Array,xptr,(n,); own=false); println(maximum(x)); end");
  jl_eval_string("x=rand(10); myfun(pointer(x),6)");    // test it in jl. I don't understand why Ref(x) and Ptr(x) fail here :(
  jl_value_t *ret = jl_eval_string("@cfunction(myfun, Cvoid, (Ptr{Cdouble},Cint))");
  myfun = (myfun_ptr) jl_unbox_voidpointer(ret);    // convert to C func ptr
  double a[4] = {1.0,2.0,7.0,3.0};                  // input array
  myfun(a,4);                                       // call it
  printf("%s: should have printed 7.0\n",__FILE__); // chk by eye
  jl_atexit_hook(0);
  return 0;
}

Assuming the above is cfunctionarr.c, compilation with

export JULIA_DIR=/usr/local/julia-1.5.2/
gcc-9 -Wall -g -o cfunctionarr -fPIC cfunctionarr.c -I$JULIA_DIR/include/julia -L$JULIA_DIR/lib -Wl,-rpath=$JULIA_DIR/lib -Wl,-rpath=$JULIA_DIR/lib/julia -ljulia
./cfunctionarr

gives expected output:

0.9694413805077366
7.0
cfunctionarr.c: should have printed 7.0

I am committing tutorial versions of this to my repo at the top of this thread. Hope they are useful to those trying to get the needed language interoperability that we’ve been used to between C<->Fortran for decades.

  1. Please let me know if there’s something bad about the above approach.
  2. Is it needed to insert return nothing in the Julia function? I recall getting an error about wrong return type, but can’t reproduce it.
  3. I also have the noob question about why pointer is needed not Ptr or Ref in the Julia-side initial test?

Best, Alex

2 Likes

The return type of the function needs to match what you declare in the @cfunction. Note that if you don’t have an explicit return statement, then the return value of a function is the value of the last statement in the function.

Here, your @cfunction is declaring the return type as Cvoid, so the function needs to return nothing in order to match this. You need an explicit return nothing if the last statement has some other type. In this example, your last statement is a call to println, which already returns nothing, so a return nothing statement would be redundant.

Ptr is the return type of pointer(array), but that doesn’t mean it has to be name of the function you use to get the pointer to the data in the array. (length(array) returns an Int, but that doesn’t mean we should call the function Int(array).)

When to use the the type name as a function (i.e. a “constructor”) is somewhat of a matter of taste, but the constructor idiom is usually employed when the constructed object is in some sense a synonym for the argument(s). pointer(array), in contrast, loses information like the dimensions of the array.

As for Ref, that is a completely different type — it is a supertype for objects like RefValue that wrap another Julia value, and are convertable to a Ptr (e.g. in a ccall), but are not identical to a bare pointer object. (An important property of a Ref object, unlike a bare pointer, is that it is “rooted” — the Julia garbage collector knows about the underlying object and will not free it out from under you while you are using the Ref.)

3 Likes