How do I pass a 2D array in C++ to a Fortran subroutine?

Question

I am writing a small C++ program which passes a 2-D array (of complex numbers) to a Fortran subroutine and receives it back filled with values. I wrote a version which passes and receives a 1-D array, and this works well. The 2-D version doesn't work (my true goal is to write a 4-D version with large dimensions - so these arrays have to be dynamically allocated).

I will post both my working code and non-working code, but first note that I was forced to use structs (simple ones, just containing two doubles) because Fortran seems to interpret these in exactly the same way as its own native complex numbers. This is why my 1-D version works. For the same reason, I don't think this is a 'complex number' issue.

Here is my working code. Passing 1-D array of complex numbers to Fortran subroutine:

The Fortran subroutine:

subroutine carray(A)
complex*16 A(2)

A(1) = cmplx(3,7)
A(2) = cmplx(9,5)

return
end

The C++ code:

include <iostream>
include <complex>
using namespace std;

struct cpx{double r, i;};

extern"C"
{
   void carray_(struct cpx* A);
}

int main()
{
   struct cpx* A;
   A = new struct cpx [2];

   carray_(A);

   complex<double>* P;
   P = new complex<double> [2];

   for(x = 0; x < 2; x++)
   {
      real(P[x] = A[x].r;
      imag(P[x] = A[x].i;
   }

   cout << real(P[0]) << "\t" << imag(P[0]) << "\n";
   cout << real(P[1]) << "\t" << imag(P[1]) << "\n";

   return 0;
}

Compiling with the following commands works without complaint:

gfortran -c CarrayF.f
g++ -c CarrayC.cc
g++ -o Carray CarrayC.o CarrayF.o -lgfortran

So as long as I treat the (native) Fortran complex number like a struct of two doubles, I can put them into the (non-native) C++ complex type. The Fortran subroutine seems to be perfectly happy receiving a pointer where it expects an array. So far so good.

Here is my non-working attempt to pass a 2D array:

The Fortran code:

subroutine carray(A)
complex*16 A(2,2)

A(1,1) = cmplx(3,7)
A(1,2) = cmplx(9,5)
A(2,1) = cmplx(2,3)
A(2,2) = cmplx(4,9)

return
end

The C++ code:

include <iostream>
include <complex>
using namespace std;

struct cpx{double r, i;};

extern"C"
{
   void carray_(struct cpx** A);
}

int main()
{
   struct cpx** A;
   A = new struct cpx* [2];
   for(int x = 0; x < 2; x++)
   {
      A[x] = new struct cpx [2];
   }

   carray_(A);

   complex<double>** P;
   P = new complex<double>* [2];
   for(int x = 0; x < 2; x++)
   {
      P[x] = new complex<double> [2];
   }

   for(x = 0; x < 2; x++)
   {
      for(int y = 0; y < 2; y++)
      {
         real(P[x][y] = A[x][y].r;
         imag(P[x][y] = A[x][y].i;
      }
   }

   cout << real(P[0][0]) << "\t" << imag(P[0][0]) << "\n";
   cout << real(P[0][1]) << "\t" << imag(P[0][1]) << "\n";
   cout << real(P[1][0]) << "\t" << imag(P[1][0]) << "\n";
   cout << real(P[1][1]) << "\t" << imag(P[1][1]) << "\n";

   return 0;
}

This actually compiles without complaint (same compilation procedure as for the 1-D version), but running the executable produces an immediate segmentation fault. Because of the headache of using two languages together, the debugger is not being helpful.

Have I made a trivial error somewhere? I don't seem to be exceeding any array bounds. The Fortran subroutine is happy to receive a pointer, but obviously it doesn't understand what to do with a pointer to a pointer. Ordinarily, Fortran would simply deal in terms of array names, even for multidimensional arrays, but I need to understand how Fortran deals with 2D arrays.

Answer 1

The way to do this in this era is to use the ISO C Binding. Officially this is part of Fortran 2003, but is is widely supported by Fortran compilers, including gfortran. It even includes interoperability of complex types! See the list of types in the gfortran manual under the chapter "Intrinsic Modules", section "ISO_C_Binding". Plus there are several examples in the chapter "Mixed-Language Programming". On the C++ side, use extern C for the routine to be called in order to use C calling conventions without name-mangling. On the Fortran side, use the ISO C Binding. Then you will have compatibility at the language level rather than having to make assumptions about how the compilers work.

Answer 2

Multidimensional arrays in Fortran are stored as flat column-major arrays in memory. They are not pointers to pointers—they're really just 1D arrays where you have to do some extra arithmetic to compute array indices.

The correct way to pass the array from C++ is like this:

extern "C" void carray_(struct cpx *A);
...
struct cpx A[2*2];
carray_(A);

complex<double> P[2][2];
for(int i = 0; i < 2; i++)
{
    for(int j = 0; j < 2; j++)
    {
        // note order of indicies: A is column-major
        P[i][j] = std::complex<double>(A[j*2 + i].r, A[j*2 + i].i);
    }
}

Answer 3

New advice

A two-dimensional array in c++ is not--I repeat not--the same as a pointer to a pointer{*}!

When you do

struct cpx** A;

you are setting up to construct a so called "ragged array" for which there is not a fortran equivalent. You want something like

struct cpx *A[2][2] = new struct cpx[2][2];

which is a pointer to a two-dimensional array with rows 2 long.

{*} Yes, you can access a pointer-to-pointer structure using the two-dimensional array notation, but they are laid out differently in memory. ::grumble:: People who tell other people that arrays and pointers are the same thing in c need to meet the Big Foam Clue Bat.

• A two dimensional array is a single allocation of <row-dimension>*<column-dimension>*sizeof(<type>). It's name decays to a pointer-to-type.
• A ragged array is 1+ allocation. The one is <column-dimension<*sizeof(<type>*) and the other are <row-dimension>*sizeof(<type>).

Old, correct, but non-applicable advice

The thing to be aware of here is that c++ and fortran believe arrays are stored differently in memory.

C++ thinks that the memory position after a[1][1] is a[1][2], while fortran believes that is is a[2][1]. The distinction is between row-major (c, c++, etc) and column major (fortran, matlab, a few others).

Note that this is separate from fortran indexing arrays from 1 by default.

Answer 4

I will skip the unnecessary complication of the complex type. It is just a trivial change to plug it back. I will also change the dimensions to 3x4, because 2x2 makes things confusing.

Static array

If the array is static, you just swap the dimensions [3][4] to (4,3). See how the modern Fortran binding with C is used:

subroutine carray(array) bind(C, name="carray")
  use iso_c_binding, only: c_double
  implicit none
  real(c_double), intent(in) :: array(4,3)
  integer :: j
  
  do j = 1, 3
    print '(*(f8.2))', array(:,j)
  end do
      
end subroutine

extern "C"
{
   void carray(double arr[3][4]);
}

int main()
{
   double array[3][4];
   
   for(int i = 0; i < 3; i++)
   {
      for(int j = 0; j < 4; j++)
      {
         array[i][j] = i + j/100.;
      }
   }

   carray(array);

   return 0;
}

> g++ -Wall carray.f90 main.cc -lgfortran
> ./a.out 
    0.00    0.01    0.02    0.03
    1.00    1.01    1.02    1.03
    2.00    2.01    2.02    2.03

Dynamic array

You really want to use a contiguous 2D array and not a jagged array in this case. Just a one big blob of memory, not rows scattered over many places of memory. If you want to index it using [][], see the link just given. I will show the simple case of a 1D array indexed in two dimensions manually (or using a macro). Again, the extension is just mechanical and almost trivial:

subroutine carray(array, n, m) bind(C, name="carray")
  use iso_c_binding, only: c_double, c_int
  implicit none
  integer(c_int), value :: n, m
  real(c_double), intent(in) :: array(m, n)
  integer :: j
  
  do j = 1, n
    print '(*(f8.2))', array(:,j)
  end do
      
end subroutine

extern "C"
{
   void carray(double* arr, int n, int m);
}

int main()
{
   int n = 3;
   int m = 4;
   double* array = new double[n*m];
   
   for(int i = 0; i < n; i++)
   {
      for(int j = 0; j < m; j++)
      {
         array[j+i*m] = i + j/100.;
      }
   }

   carray(array, n, m);

   return 0;
}

> g++ -Wall carray.f90 main.cc -lgfortran
> ./a.out 
    0.00    0.01    0.02    0.03
    1.00    1.01    1.02    1.03
    2.00    2.01    2.02    2.03

Please note the value attribute in Fortran. If you do not want it or cannot put it there, you must pass addresses (int* and &n) in C or C++ instead.