Transfer Float64 matrix to double **arg in C

I’ve been trying to transfer a Julia matrix of Float64 to a C function that has a “double **arg” but it does not work as intended so help would be appreciated.

C function I want to call (simplified body just to get idea across):

double emd2(int n_x, double *weight_x, int n_y, double *weight_y, double **cost) {
    double sum = 0.0;
    int i, j;

    for(i = 0; i < n_x; i++) {
        sum = sum + weight_x[i];
    }
    for(i = 0; i < n_y; i++) {
        sum = sum + weight_y[i];
    }
    for(j = 0; j < n_y; j++) {
        for(i = 0; i < n_x; i++) {
            sum = sum + cost[j][i];
        }
    }
    return sum;
}

then I compile this into a lib called libemd and define a low-level function to call it like so:

libemd_emd2(n_x, weight_x, n_y, weight_y, cost) =
    ccall((:emd2,libemd), Cdouble,
          (Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Ptr{Cdouble}}),
          n_x, pointer(weight_x), n_y, pointer(weight_y), pointer(cost))

but when I then try to use it:

w1 = [0.4, 0.2, 0.2, 0.1, 0.1]
w2 = [0.6, 0.2, 0.2]
cost = ones(Float64, 5, 3)
cost_ref = [Ref(cost,i) for i=1:size(cost,1):length(cost)]
libemd_emd2(length(w1), w1, length(w2), w2, cost_ref)

The result is some large number rather than the sums of the values of the matrix and the two arrays. Seems it is reading some random place in memory so I guess I set up the pointers the wrong way or something. I based my solution on this stackoverflow entry: Calling a C function from Julia and passing a 2D array as a pointer of pointers as argument - Stack Overflow

What am I missing here? Thanks for any advice.

This doesn’t really answer your question, but if the C code is something you can edit I would rewrite it to use double *cost as an n_x by n_y column-major array. Just change double **cost to double *cost and change cost[j][i] to cost[j*n_x + i]. Not only will this be faster and easier to call from Julia, but it will probably also be faster in C as well because it avoids the double pointer indirection.

Thanks for the quick reply Steven. I do have access to the C file and it is under an open-source license so I guess I can write a small script to automatically do this conversion, yes. It would simplify and speed up since I will also not need to transpose the matrix (which might be a big one). Only downside is if the original author evolves the code and my script might not be “matched up”… But I’ll see if this is a viable way forward. Thanks!

Note that if you can change the C code, you should change it to a multi-dimensional C array and you don’t need to compute the index yourself.

See either solution at Redirecting to Google Groups or Calling a c-function that returns a 2-dimensional array - #8 by yuyichao

Be careful about column/row ordering…

Somehow this question comes up much more often than I though…

Thanks to both of you. I went with Steven’s proposal since I find it easier to get the column/row ordering right this way. It works now, thanks. If this is a common question maybe some comments about multi-dimensional arrays should be added to the Julia docs for interfacing with C?