Dear forum members,
In general what is the recommended method to pass multi dimensional Fortran arrays to C functions?
Thank you in advance
It's still ugly, but I continue to use the old f2c scheme, a literal implementation of assumed size, which depends on the Fortran array being a sequence type. Technically, in order to use iso_c_binding, I suppose the value passed over to C should be taken by c_loc, and of course the data type would be one of those from iso_c_binding. The sizes would need to be passed (as is usual with assumed size). On the C side, the address is calculated literally by using the Fortran subscripts and sizes to find the offset from the passed address.
As a consequence, the C compiler considers all elements of the array to be aliased with each other, even though the pointer is defined with restrict so that it doesn't alias other arrays (unless some trick such as icc #pragma vectorlength is invoked).
Keep in mind that C "arrays" aren't like Fortran arrays, especially when it comes to more than one dimension. There is no standard C syntax corresponding to a Fortran multidimensional array - something like ARRAY in C is an array of pointers to arrays, which is not directly representable in Fortran. This is why only one-dimensional arrays are considered "interoperable".
First of all I would like to thank you for the valuable comments. As far as I know C arrays are stored row major order in the memory, on the contrary Fortran array are stored column major order. Therefore, one has to transpose them. I figured out that one can reshape the array into 1D and then pass it to C. From this point on one can use the C pointer arithmetic and do things like “*(array+offset)” in order to get the needed element. However I need some more experiments in order to check the consistency of what I did.
I have a follow up question which is related to this array stuff. If you do not mind I would like to pose it here. I am running a Monte-Carlo code and as my particles fly over the domain I update some arrays. In general I change the values of those arrays. The problem is that my arrays are 4D.
I do things like arr1(I,j,k,1:5)=arr2(1:5) I guess that the code is doing some serious loops over the memory in order to go to the right position. If you do it once it is fine but if you do it millions of times well that adds up time. For instance in order to test some MPI broadcast of arrays I put a simple construction such as
if( any(array /= Value)) then
This resulted in a dramatic increase in computation time. More than 5 times, quite a depressing fact after all, my array had 500 integer elements.
This is actually why I was trying to use C functions because I was thinking that if I pass just references things might be faster. I guess this is an ugly solution, and not a wise thing to do. So I will really appreciate if someone can give me any idea. Is it at all possible to fine better solution and speed up things? I am using Intel Fortran compiler under Linux.
Thank you in advance
If your solution for performance issues related to cache locality is to make a transposed copy of the array, that works at least as well in Fortran. You hit the locality issue while performing the transpose, but it may be worth doing if it saves you hitting it many times.
You seem to be going off in several directions with this, so I don't know which one you'd like to pursue, within reasonable limitations on topicality.
Generally speaking, reaching the same level of performance in C as in Fortran, for applications where Fortran is well suited, involves significant additional effort and reduction in readability. If you need to go further than to consider the requirement for dealing with ansi-alias/strict-aliasing options and restrict qualifiers, besides the issues we already touched on about the different representations of multi-dimensioned arrays, there's a lot of contradictory material already available for you to research.
Well the code I am using is Fortran 77, which calls C functions in order to allocate arrays. It is an archaic piece of software, and this is how I stumbled upon all these peculiar problems. Well I think I the best approach would be to rewrite everything in Fortran-90 (use dynamic memory allocation and substitute the commons with modules).
Yes I really went out of the scope of my original question so I would like to apologize for that. I just wanted to get some professional opinion on my problems.
At this point I would like to apologize for the inconvenience of broadening the original topic scope.
Thanks for the help.
>>I do things like arr1(I,j,k,1:5)=arr2(1:5)
Can you change your subscript order such that the preponderance of the equates/manipulations of sequential elements occur on the left most dimension?
Changing the preponderance of the equate operation to the left most dimension is a really good idea. I can redo the indexing as I wish. Well actually my array has 3 dimensions resembling the coordinates (x,y,z) and the last dimension is material composition percentage. I would have to do some alternations in my code, but will definitely try the idea.
Thanks for the idea.