Calling Fortran functions which return KIND=10 complex values from C++

Question

I have a bunch of old F77 source code (usually compiled on a x86_64 with gfortran -std=legacy). It contains quite a few functions in form:

      double complex function f(x, y, i)
      double precision x, y
      integer i
      f = cmplx(x, y) * i
      return
      end

I need to call these functions from some C++ code (usually compiled on a x86_64 with g++).

1. It works with the default Fortran KIND=8:

   extern "C" { std::complex<double> f_(double *x, double *y, int *i); }
   

2. It works when I enforce the default Fortran KIND=4 using the -freal-8-real-4 option:

   extern "C" { std::complex<float> f_(float *x, float *y, int *i); }
   

3. It works when I enforce the default Fortran KIND=16 using the -freal-8-real-16 option (and in C++ #include <quadmath.h>):

   extern "C" { __complex128 f_(__float128 *x, __float128 *y, int *i); }
   

   To my surprise, in this case, it also seems to work with (the returned value is in *z):

   extern "C" { void f_(__complex128 *z, __float128 *x, __float128 *y, int *i); }
   

   Which of these two above prototypes is the (more?) proper one?

4. My problem is that I cannot get it working with my desired default Fortran KIND=10 using the -freal-8-real-10 option. Inside of Fortran, the kind, precision, range and sizeof return values which directly correspond to the C++ long double. So, I tried:

   extern "C" { std::complex<long double> f_(long double *x, long double *y, int *i); }
   extern "C" { void f_(std::complex<long double> *z, long double *x, long double *y, int *i); }
   extern "C" { void f_(long double *x, long double *y, int *i, std::complex<long double> *z); }
   

   But I cannot get it working at all.

   Maybe I need to add some special flags to gfortran and / or g++ calls in order to let C++ retrieve Fortran KIND=10 complex values? Note: I don't think I can use -ff2c.

Update (2020.08.04): I have been able to fool the C++ compiler so that it seems to generate proper code for any Fortran KIND=4,8,10. The trick is to use the ISO C99 _Complex in C++ (note: this trick is required only for KIND=10, but it actually works for KIND=4,8, too):

#include <complex.h>
#define C99KIND long double /* it can be "float", "double" or "long double" */
extern "C" { C99KIND _Complex f_(C99KIND *x, C99KIND *y, int *i); }

Note that, in C++, you cannot use e.g. long double complex but fortunately long double _Complex is still fine.

The usability of the ISO C99 _Complex in C++ is rather limited. For example, with -std=c++11 (or newer) even the most basic creal* and cimag* functions disappear.

So, the best idea is to immediately copy the returned value into some standard C++ templated complex variable, e.g using something like (note: f_ returns C99KIND _Complex):

std::complex<C99KIND> z = f_(&x, &y, &i);

Answer 1

If you want to do the job correctly, then learn how to actually use Fortran's kind type parameters and do a proper port of your Fortran code to REAL(10). Yes, I know 10 is not portable; we are however, discussing a particular Fortran processor.

Consider,

function f(x, y, i) result(r) bind(c, name='f')
  use iso_c_binding, only : ep => c_long_double
  implicit none
  complex(ep) r
  real(ep), intent(in), value :: x, y
  integer, intent(in), value :: i
  r = cmplx(x, y, ep) * i
end function f

and, since I don't do C++ but you should be able to update the C to your needs

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

long double complex f(long double, long double, int);

int
main(void)
{
  int i;
  long double x, y;
  long double complex z;

  i = 42;
  x = 1;
  y = 2;
  printf("%.10Lf %.10Lf\n", x, y);

  z = f(x, y, i);
  x = creall(z);
  y = cimagl(z);
  printf("%.10Lf %.10Lf\n", x, y);
  return 0;
}

% gfortran -c a.f90
% gcc -o z b.c a.o -lm
% ./z
1.0000000000 2.0000000000
42.0000000000 84.0000000000

OP states that he cannot be bothered with a proper port of his/her Fortran code, and so therefore, must use a compiler option to magically (and yes it is magic) to do the port. Here's an example

% cat a.f90
double complex function f(x, y, i)
  implicit none
  double precision :: x, y
  integer i
  f = cmplx(x, y) * i   ! cmplx my not do what you expect
end function f

% cat b.c
#include <complex.h>
#include <stdio.h>
#include <stdlib.h>

long double complex f_(long double *, long double *, int *);

int
main(void)
{
  int i;
  long double x, y;
  long double complex z;

  i = 42;
  x = 1;
  y = 2;
  printf("%.10Lf %.10Lf\n", x, y);

  z = f_(&x, &y, &i);
  x = creall(z);
  y = cimagl(z);
  printf("%.10Lf %.10Lf\n", x, y);
  return 0;
}

% gfcx -c -freal-8-real-10 a.f90
% gcc -o z b.c a.o -lm
% ./z
1.0000000000 2.0000000000
42.0000000000 84.0000000000

Might as well go for the trifecta. Here's the C++ code to go with the file a.f90 above in the second example.

#include <iostream>
#include <complex>
#include <cmath>

extern "C" { std::complex<long double> f_(long double *, long double *, int *); }

int main()
{
  std::complex<long double> z;
  long double x, y;
  int i;

  i = 42;
  x = 1;
  y = 2;
  z = f_(&x, &y, &i);

  std::cout << z << '\n';
 }