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);
    }
}