Sachin
Sachin

Reputation: 48

Properly passing multidimensional C array to fortran such that size function (in fortran) gets the correct size in each dimension?

I have a fortran subroutine like this which expects a 3d array and prints the size along each dimension:

subroutine printarray(arr)
real*8 :: arr(:,:,:)

print*, "Dimensions of arr -> ", size(arr(:,1,1)), size(arr(1,:,1)), size(arr(1,1,:))

end subroutine printarray

Now I want to pass a 3d array from a C program to this subroutine.

I attempted it in many ways, none of them worked.

Attempt 1:

#include <stdio.h>
#include <stdlib.h>

void printarray_(double *);

int main(int argc, char *argv[])
{
  double *ap;

  ap = (double *) malloc(sizeof(double) * 10 * 10 * 2);
  printarray_(ap);

  if(ap != NULL)
    free(ap);
  return 0;
}

Output is : 1 1 1

One issue here is I haven't explicitly specified the size of each dimension in my C code itself, so I cannot expect Fortran to figure it out. So I tried it differently:

Attempt 2:

#include <stdio.h>
#include <stdlib.h>

void printarray_(double ***);

int main(int argc, char *argv[])
{
  int i;
  int j;
  double ***a;
  a = (double ***) malloc(sizeof(double **) * 200);
  for(i = 0;i < 10;i++)
    {
      a[i] = (double **) malloc(sizeof(double *) * 10);
      for(j = 0;j < 10;j++)
      {
        a[i][j] = (double *) malloc(sizeof(double) * 2);
      }
    }

  printarray_(a);

  for(i = 0; i < 10;i++)
    {
      for(j = 0;j < 10;j++)
      {
        if(a[i][j] != NULL)
          free(a[i][j]);
      }
      if(a[i] != NULL)
        free(a[i]);
    }
  if(a != NULL)
    free(a);
  return 0;
}

Here, the output is : 289 289 1

What is the correct way of passing multidimensional arrays to Fortran subroutines from C programs, such that the size function in fortran picks up correct size of each dimension?

Upvotes: 1

Views: 708

Answers (2)

user4490638
user4490638

Reputation:

If your Fortran compiler supports ISO/IEC TS 29113:2012 , you can use the ISO_Fortran_binding.h header for interoperability of Fortran arrays with C. I am not sure which compilers support the TS; unfortunately, gfortran does not do so, in this respect. So, assumed-shape arrays like arr(:,:,:) might not be an option for you.

Failing that, you should use Fortran's Iso C binding feature.

Note that Fortran arrays consist of consecutive storage locations. For 1D-arrays, this is not an issue for C interop. For arrays with two or more dimensions, this creates problems because the normal C way of accessing arrays through pointers to pointers to ... has no good equivalent in Fortran (and is inefficient anyway).

The closest you can get in C to Fortran arrays is through the use of variable-length arrays. Unfortunately, these have been relegated to optional in the last version of the C standard, so C support for multidimensional arrays has fallen back behind Fortran 66 again. However, I am assuming that you have a compiler which supports VLAs (gcc does) and you do not hit the limitations that sometimes plague the implementaions that, for example, VLAs are only allocated on the stack and that the stack is too small for your problem size...

So, you should be able to declare

  int m,n,k;
  // Obtain a value for m,n,k
  {
     double arr[10][10][2];
     // Do something with it
     printarray(a,k,m,n)  // Fortran arrays have the opposite ordering of indices versus C
  }

... and have on the Fortran side

subroutine printarray(arr,k,m,n) bind(C,name="printarray")
  use iso_c_binding
  integer(c_int), value :: k, m, n
  real(c_double), dimension(k,m,n) :: arr

Upvotes: 1

user3528438
user3528438

Reputation: 2817

You just can not do it in C, because arrays in C decay to pointers during a function call.

There are a few ways to do it: Function parameter as array with declared size

But best and most common practice is still to pass in dimensions as separate inputs, and also pass in the address of of the beginning element and let subroutine figure out the rest itself.

And your second attempt is dangerous because dimensions are not contiguous but your function might expect/assume it to be so.

Upvotes: 1

Related Questions