How to build a wrapper around a C-style array using Julia C API?

General Usage
, ,
1

Hi, I need to use the Julia C API to build a wrapper around a C-style array to use this array in Julia (the array is actually a Python NumPy array). The thing is that the array is a matrix (2D) and Julia uses Fortran-style arrays.

I have the following basic code:

            jl_value_t *arr_type =
                jl_apply_array_type((jl_value_t *)jl_float64_type, c_array->nd);
            jl_value_t *dims =
                build_julia_tuple_from_size_t_array(c_array->dimensions, c_array->nd);
            bool own_buffer = false;
            cur_julia_arg = (jl_value_t *)jl_ptr_to_array(arr_type, c_array->data,
                                                          (jl_value_t *)dims, own_buffer);

where c_array is a structure like this:

// This structure closely follows PyArray_Object that describes NumPy arrays.
typedef struct {
    // Number of dimensions in the array.
    int nd;
    // Size of each axis, i = 0, .., nd-1.
    intptr_t *dimensions;
    // Pointer to actual data.
    double *data;
    // Flags that describe the array: things like row-major order, etc.
    int flags;
} CArrayF64;

What I do not understand, is how to adjust the strides in the Julia wrapper array, as they must be in the opposite order to what they are using the above code.

An auxiliary question is, does it make sense at all? Do Julia libraries understand that arrays can be C- or Fortran-like?
I have done the opposite (wrapping a Julia array using Python C API), and because NumPy/SciPy use BLAS/LAPACK under the hood, they handle both styles of arrays transparently. But I am not sure, if this is the case for Julia.

Thank you in advance!

1 Like
2

PyCall.jl and PythonCall.jl already implement new AbstractArray types that can wrap around row-major numpy arrays, why not use them?

Otherwise, you can just use transpose(unsafe_wrap(Array, pointer_to_row_major_array, dims..., own=false)) to wrap a lazy transposition around it. A big advantage here is that many of the linear-algebra routines have specialized methods for Transpose array wrappers, allowing them to perform well on row-major arrays without resorting to generic fallbacks.

(I would urge you to write as much of your interface code in Julia as possible … it’s really a lot easier, not to mention better documented, than working directly with the low-level C API. If you need to call something in Julia from C, you can just make a @cfunction from a Julia routine to get a C-callable function pointer.)

1 Like
3

Thank you for the answer!

I am working on this software: GitHub - MaRDI4NFDI/open-interfaces: Measure 2.2 “Open Interfaces” of the Mathematical Research Data Initiative (MaRDI) project
and do not quite understand how to do what I do using Julia directly. Essentially, I need to be able to import Julia modules programmatically (because the module names are known only in runtime as strings) and be able to invoke function in those modules while also being able to convert a generic list of void pointers with ids to corresponding Julia types, something like this:

call_interface_impl(
    impl:String,
    method:String,
    args:List{Ptr{Cvoid}},
    argsTypeIds:List{Int}
)

and I can understand how to do it with embedded Julia, but do not quite understand how would I do it in Julia directly.

ADDED Forgot to answer the first question.

PyCall.jl and PythonCall.jl already implement new AbstractArray types that can wrap around row-major numpy arrays, why not use them?

The idea is that the code should be agnostic to the fact that the array is from Python, because it can be from another language, e.g., C.

4

eval can import Julia modules programmatically. And you can convert pointers (with metadata indicating the size and type of the data) to Julia types by ordinary Julia code (in the same way as you would read data from a file, or convert data from Python), e.g. using unsafe_wrap as shown above.

All this will be a lot easier in Julia than using the C API.

5

Why would you need to do this? If the wrapper is just forwarding to the indexing or whatever other operations the c_array supports, then it only needs to process inputs and outputs between Julia and C, no need for separate or derivative copies of metadata. AbstractArray code often won’t care what the memory layout is, it just indexes on the semantic level. Efficient indexing orders can be (not always possible, not always implemented) abstracted away with eachindex, other times the user has to be aware of the concrete type they’re using e.g. order of axes in nested loops.

For example, while A::Array has a strictly column-major layout, that also means that axes-reversing wrapper types like PermutedDimsArray(A, reverse(1:length(axes(A)))) have a strictly row-major layout. Such a row-major array doesn’t separately store strides, nor do we reverse the dimensions to match the column-major parent. Permuting dimensions is more general than reversing axes so PermutedDimsArray doesn’t leverage efficient indexing in eachindex, but the point is row-major array types do naturally exist in Julia (transpose is the more common one strictly for matrices). ArrayLayouts.jl is an effort to map layouts to array types for dispatch to optimized linear algebra routines.

The only possible place where I could imagine swapping axes in a wrapper concerns how matrices tend to be “stacked” in higher dimensional arrays. While axes are conventionally ordered by strides, we also expect rows to precede columns in contiguous matrices, hence the strangely almost reversed axes between row-major (…, depth2, depth, rows, columns) and column-major (rows, columns, depth, depth2, …). That adds a headache when loading stacked matrices in the wrong layout because simply reversing the axes doesn’t get it right for the matrices’ axes. Your wrapper could automatically permute dimensions for convenience if you have a narrow context, but otherwise I’d leave it up to the users to just learn their types and what the ordered axes mean.

1 Like