Calling others languages in

source available at https://github.com/aitzkora/nahasketa/tree/master/julia/interfaces
- with CMake files to build shared objects

disclaimer :
- we address only the Julia calling other code part
- only tested under Linux (more precisely Funtoo)

Outline

calling C code

calling Fortran code

calling Python

Calling C code

the basic julia syntax to interface with C code is

ccall((:funcName, library), returnType, (argType1, argType2, ...), (argVal1, argVal2, ...))
ccall(funcPtr, returnType, (argType1, argType2, ...), (argVal1, argVal2, ...))

Remarks

library is a formally a string :
- you could use "./mylib.so"
- but ⚠ you could not use string(pwd(),"/mylib.so") ⚠
to use a library which is not in ., add the path to LD_LIBRARY_PATH before launching Julia
using dlopen and dlsym one could directly use the function pointer call

Basic example

we want to call a scalar function taking one int and returning one int
we write the C function and compile it



In [1]:

    
io = open("C/skel.c","w")
write(io, "int ajoute2(int x) { return x+2; }")
close(io)
run(`gcc -o C/ajoute2.so --shared C/skel.c`);

and we call it in Julia



In [2]:

    
w = ccall((:ajoute2, "ajoute2.so"), Int32, (Int32,), 12)
run(`rm C/ajoute2.so C/skel.c`)
println("w = $w")

Example with a matrix

We can call the following C code

#include <inttypes.h>
double mysum(int64_t m, int64_t n, double *x)
{
 double s = 0.;
 int i,j;
 for(i = 0; i < m; ++i)
     for(j = 0; j < n; ++j)
          s += x[i * n + j];
 return s;
}

with the julia code

a = reshape(collect(1.:10.), 5, 2)
s = ccall((:mysum, "libmysum.so"), Float64, (Int64, Int64, Ptr{Float64}), size(a, 1), size(a, 2), a)
println("s = $s")



In [7]:

    
include("C/test_mysum.jl")

Notes :

we use Ptr{Float64} to pass a matrix to the function
using standards types like int64_t is not mandatory, but enforces compatibility trough more portability
we could also use Cdouble and Cint in Julia code

Calling a Julia Function in the C code called by Julia

using a function pointer in C code

#include <inttypes.h>

double apply_f(double (*fun)(double), int64_t x_n , double *x, double * y)
{
  int64_t i;
  for(i = 0; i < x_n; ++i)
      i[y] = (*fun)(i[x]);
}

and the @cfuntion julia macro as in

a = collect(1.:10.)
fun = @cfunction(x->x*x, Float64, (Float64,))
b = zeros(10)
s = ccall((:apply_f, "libjuliafun.so"), Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}), fun, 10, a, b)
println("b = $b")

one could call julia function into the C code



In [3]:

    
include("C/test_apply_f.jl")









    



b = [1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0]

Basic C++ support

we can use C++ code, as long as we provided a extern "C" interface for the entry function

#include <cinttypes>
#include <vector>
#include <numeric>

extern "C" {
double mysum2(int64_t n, double *x);
};

double mysum2(int64_t n, double *x)
{
    double s = 0;
    std::vector<double> t(x, x + n);
    for(auto & z : t) s+= z;
    return s;
}

and we could call with the julia code

a = collect(1.:10.)
s = ccall((:mysum2, "libmysum2.so"), Float64, (Int64, Ptr{Float64}), size(a, 1), a)
println("s = $s")



In [4]:

    
include("C/test_mysum2.jl")

Persistance across multiple calls

#define ALLOC_MAX 4
double *table[ALLOC_MAX];
int64_t compteur = 0;

int64_t my_alloc(int64_t size) {
    double * p;
    printf("comp = %d\n", compteur);
    if (compteur <  ALLOC_MAX )
    {
        p = (double*) (calloc(size, sizeof(double)));
        if (p)
        {
            table[compteur] = p;
            return compteur++;
        }
        else
            return -1; // bad calloc
    }
    else
      return -2; // allocation table is full
}

void my_free(int64_t num)
{
     if (table[num]) { free(table[num]); table[num] = NULL; }
}

int64_t set_values(int64_t size, int64_t num, double * val)
{
    if (table[num])
        memcpy((void*)(table[num]), (void*)(val), sizeof(double)*size);
    else
        return -1;
}

int64_t get_values(int64_t size, int64_t num, double * val)
{
    if (table[num])
        memcpy((void*)(val), (void*)(table[num]), sizeof(double)*size);
    else
        return -1;
}

using Libdl: dlopen, dlsym
lib_alloc = dlopen("libmyalloc")
all_ptr, fre_ptr, set_ptr, get_ptr = map(x->dlsym(lib_alloc, x), [:my_alloc, :my_free, :set_values, :get_values])

function alloc(n::Int)
    ret = ccall(all_ptr, Int64, (Int64, ), n)
    if (ret == -2) println("allocation table is full") end
    if (ret == -1) println("calloc could not allocate another bloc") end
    return ret
end

function free(num::Int)
    ccall(fre_ptr, Int64, (Int64, ), num)
end

function set_values(num::Int, vals::Array{Float64, 1})
    ret = ccall(set_ptr, Int64, (Int64, Int64, Ptr{Float64}), size(vals, 1), num, vals)
    return ret
end

function get_values(num::Int, vals::Array{Float64, 1})
    ret = ccall(get_ptr, Int64, (Int64, Int64, Ptr{Float64}), size(vals, 1), num, vals)
    return ret
end

we now could use the C code to allocate some double memory blocks

h1 = alloc(10); h2 = alloc(5)
a = collect(1.:10.)
set_values(h1, a)
b = zeros(3)
ret2 = get_values(h1, b)
println("b = $b")
free(h1)
h3 = alloc(1); h4 = alloc(1); h5 = alloc(1)



In [5]:

    
include("C/test_myalloc.jl")









    



b = [1.0, 2.0, 3.0]
allocation table is full






    Out[5]:





-2

Using a reference for a function modifying some arguments

if a function modifying an argument which is not a pointer like
```
int increment(int * x) { (*x)++; }
```

you have to use a julia reference such as in

a = Ref{Cint}(12)
ccall((:increment, "libincrement.so"), Cvoid, (Ref{Cint},), a)
println(a[])

Calling Fortran code

the principle is the same for C, modulo mangling of the name → use bind to mangle correctly the name

module heat
    use iso_c_binding, only: c_int32_t, c_double
    public :: kernel
contains
    subroutine kernel(m, n, u_in,  u_out, error) bind( C, name="heatKernel" )

        implicit none
        integer(c_int32_t), intent(in) :: m, n
        real(c_double), dimension( 1:m, 1:n ), intent(in) :: u_in
        real(c_double), dimension( 1:m, 1:n ), intent(out) :: u_out
        real(c_double), intent(out) :: error

        integer(c_int32_t) :: i, j

        error = 0.d0
        u_out(2:m-1,2:n-1) = 4.d0 * u_in(2:m-1, 2:n-1) &
                                  - u_in(1:m-2, 2:n-1) &
                                  - u_in(3:m, 2:n-1)   &
                                  - u_in(2:m-1,1:n-2)  &
                                  - u_in(2:m-1,3:n)
        error  =  sum((u_out(:,:) - u_in(:,:))**2)

    end subroutine kernel
end module heat

could be mapped on a julia function with

function heat(src::Array{Float64,2}, dest::Array{Float64,2})
    @assert size(dest) == size(src)
    (size_x, size_y) = size(dest)
    size_x = Int32(size_x) # could overflows
    size_y = Int32(size_y)
    err = Ref{Float64}(0.) # ✏️ VERY IMPORTANT ✏️ 
    ccall((:heatKernel, "./libheatKernel.so"), Cvoid, (Ref{Int32}, Ref{Int32}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}), size_x, size_y, src, dest, err)
    return err[]
end

Remarks

since fortran pass value by address, all values are mapped on reference
err is modify in the subroutine, we must declare as a Ref{Float64}
arrays could be passed as Ptr{Float64}
we also use here standards like c_double to enforce portability

Example with MPI

sometimes, we want to call a parallel code using MPI

subroutine solve(n_x, n_y, p_x, p_y, snapshot_step, snapshot_size, iter_max, solution) bind( C, name="solve" )
  use iso_c_binding, only: c_int32_t, c_double
  use communications
  use heat
  use mpi
  implicit none
  integer(c_int32_t), intent(in) :: n_x, n_y ! sizes of the global matrix
  integer(c_int32_t), intent(in) :: p_x, p_y ! nb of processes (in each dimensions)
  integer(c_int32_t), intent(in) :: snapshot_step, snapshot_size
  integer(c_int32_t), intent(inout) ::  iter_max ! max number of iterations
  real(c_double), intent(inout) :: solution(1:n_x, 1:n_y, 1:snapshot_size)

  ! (...)
  call set_bounds( coords, p_x, p_y, u_in)
  call set_bounds( coords, p_x, p_y, u_out)

  do i=1, iter_max

      call stencil_4( h_x, h_y, d_t, u_in, u_out, error_loc )
      call MPI_ALLREDUCE( error_loc, error, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr )
      error = sqrt( error )

      if (mod( i, snapshot_step ) == 0)  then
          !(...)
          call gather_solution( sol_space, n_x, n_y, u_in, ndims, comm2D, is_master )
          if (is_master) solution(:, :, i / snapshot_step) = sol_space(:, :)
          !(...)
      end if

      call ghosts_swap( comm2D, type_row, neighbour, u_in )
      u_in = u_out
      if (error <= prec) exit
  end do

  if (i < iter_max ) iter_max = i ! to know if we do all loops
  !(...) 

  end subroutine solve

we can call this code using MPI package → using MPI will do the job

import MPI

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    commSize = MPI.Comm_size(comm)
    commRank = MPI.Comm_rank(comm)
    nArg = size(ARGS, 1)
    # (..) some checks about nargs 
    nX, nY, pX, pY, iterMax, snapshotStep = map(x->parse(Int32, x), ARGS)
    # (...) some checks than pX*pY == commSize
    @assert snapshotStep < iterMax
    snapshotSize = max(iterMax ÷ snapshotStep, 1)
    iter = Ref{Int32}(iterMax) # very important !!
    solution = zeros(nX, nY, snapshotSize)

    ccall((:solve, "./libheat_solve.so"), Cvoid,
    (Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32},  Ref{Int32},  Ref{Int32}, Ref{Int32}, Ptr{Float64}),
     nX, nY, pX, pY, snapshotStep, snapshotSize,  iter, solution
    )
    iter = iter[] # very important !!
    if (commRank == 0)
        indexToDisplay = (iter == iterMax ? snapshotSize : (iter ÷ snapshotStep))
        show(IOContext(stdout, :limit=>true), "text/plain", solution[:, : , indexToDisplay]);
    end
    MPI.Finalize()
end
main()

after we can call the code in a shell doing

mpirun -np 4 julia heat_solve.jl 6 6 2 2 30 10

Calling Python from Julia

we could calling Python using the package PyCall : e.g. PyPlot package binding method uses PyCall

we could import a python symbol trough @pyimport and julia will convert the result



In [2]:

    
using PyCall
@pyimport numpy
a = transpose(numpy.arange((5)))









    Out[2]:





1×5 LinearAlgebra.Transpose{Int64,Array{Int64,1}}:
 0  1  2  3  4

using directly python code is possible



In [3]:

    
py"""2 in range(10)"""









    Out[3]:





true

to more complex objects than primitives types, we use ["name_field"] for methods or members



In [4]:

    
so = pyimport("scipy.optimize")
filter(x->occursin("newton", x),convert(Array{String}, so["__all__"]))









    Out[4]:





2-element Array{String,1}:
 "newton"       
 "newton_krylov"



In [6]:

    
so["newton"](x->log(x) - 1 , 1)









    Out[6]:





2.718281828459031

Example with VTK

vtk has C++ bindings (hard to bind due to C++ classes) but also Python bindings 😅
➡️ let's write a simple mesh Reader using PyCall

struct Mesh
    points::Array{Float64,2}
    cells::Array{Array{Int}}
end
using PyCall
@pyimport vtk
@pyimport vtk.util as vu
@pyimport vtk.util.numpy_support as vuns
function readMeshFromFileVtk(fileName::String)
    reader=vtk.vtkUnstructuredGridReader()
    reader["SetFileName"](fileName)
    reader["Update"]()
    out = reader[:GetOutput]()
    pts = out[:GetPoints]()
    pointsData = pts[:GetData]()
    points = convert(Array{Float64,2}, vuns.vtk_to_numpy(pointsData))
    nbPoints = size(points, 1)
    clls = out[:GetCells]()
    cllData = clls[:GetData]()
    cellsRaw = vuns.vtk_to_numpy(cllData)
    i = Int32(1)
    nbCells = out[:GetNumberOfCells]()
    cellIdx = 1
    cells = fill(Int32[], nbCells) # beware 32 bits cast
    while i < size(cellsRaw,1)
         currentCellSize = cellsRaw[i]
         cells[cellIdx] = cellsRaw[i+1:i+currentCellSize]
         i += currentCellSize+1
         cellIdx += 1
    end
    return Mesh(points, cells)
end

⚠ you need a VTK lib compiled against the same python executable used by PyCall ⚠

Conclusion

making interfaces remains simple in Julia : no glue langage, types directly mapped

do not forget use Ref{T} when Julia variables are modify by the routine

the calling Julia code must do some checks before passing parameters

Interfacing with is possible : see package RCall

References
- https://github.com/JuliaPy/PyCall.jl (Well documented)
- https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/

Calling others languages in

Outline

calling C code

calling Fortran code

calling Python

Calling C code

Remarks

Basic example

Example with a matrix

Notes :

Calling a Julia Function in the C code called by Julia

Basic C++ support

Persistance across multiple calls

Using a reference for a function modifying some arguments

Calling Fortran code

Remarks

Example with MPI

Calling Python from Julia

Example with VTK

⚠ you need a VTK lib compiled against the same python executable used by PyCall ⚠

Conclusion