Send vector of doubles in a C-struct to Fortran Derived Type

21

What about this? Based on what cmaapic said, I am not certain about your claim. I learned a lot from what you said so far otherwise

22

OK, I give up, you clearly don’t want to understand.

1 Like
23

The C-descriptor won’t help because your misalignment is at the struct level. The reason you have undefined behavior is because your types are not interoperable.

Let me start by re-stating that, the struct you’ve shown here,

struct blah //define structure
{
    int  n;
    int* array;
    double* array2;
    std::vector<int> myvec;
    //std::vector<double> myvec2;
    double* D;
};

is NOT interoperable with Fortran, just like @PIerU has said. While the scalar and pointer components are interoperable, the std::vector members are not. Unlike C, where arrays are just “chunks” of memory, the std::vector objects are class instances containing additional private information. Here is a skeleton of what std::vector might contain in practice (vendors can choose to implement std::vector differently):

template <class T, class A = std::allocator<T> >
class vector {
public:
    // public member functions
private:
    T*                    data_;
    typename A::size_type capacity_;
    typename A::size_type size_;
    A                     allocator_;
};

It is these extra members, which cause the misalignment, among other things. In your struct blah, the members myvec and myvec2 are not pointers but values, so trying to interoperate with these as type(c_ptr) objects is futile.

As both @PierU and @cmaapic have shown, you can pass the data component of a std::vector to Fortran, either by using .data() member function which returns a pointer to the underlying storage, or by using the address of the first element directly as in &x[0] (less desirable, IMO). This is known as pass-by-reference.

If you have multiple arrays/vectors of the same length, you can pack them together in one single struct, similar to the one shown by @PierU, or pass the pointers individually, depending how you’d like to have things organised. If the method is not public, you could just pass the pointers directly, without a temporary struct, but the function prototype will be lengthier.

That said, C++ has a few advanced concepts which can make such passing effortless for the consumer, including user-defined conversions:

// main.cpp
#include <vector>

extern "C" {
  // A C- and Fortran-compatible struct
  struct f_multivec {
      int n;
      int *ivec;
      double *dvec;
  };
  // Fortran routine we will invoke from C++
  void send2fortran(struct f_multivec s); // ATTENTION - pass by value!
}

// Our C++ class with ownership of the data
struct Multivec {

   // Constructor of the multi-vector
   Multivec(int n_) : 
        n(n_), ivec(n,0), dvec(n,-1.0) {}

   // Public members
   int n;
   std::vector<int> ivec;
   std::vector<double> dvec;

   // User-defined conversion (implicit)
   operator f_multivec() { 
    return {n,ivec.data(),dvec.data()} ; 
  }
};

int main() {
  auto s = Multivec(3);
  send2fortran(s);  // <-- sneaky implicit conversion
}

! send2fortran.f90

!> void send2fortran(struct multivec s); // pass by value !
subroutine send2fortran(s) bind(c)
  
  use, intrinsic :: iso_c_binding
  implicit none
  type, bind(c) :: multivec
    integer(c_int) :: n
    type(c_ptr) :: ivec, dvec
  end type

  type(multivec), intent(in), value :: s

  integer(c_int), pointer :: ivec(:) => null() 
  real(c_double), pointer :: dvec(:) => null()

  call c_f_pointer(s%ivec,ivec,[s%n])
  call c_f_pointer(s%dvec,dvec,[s%n])

  print *, "ivec = ", ivec
  print *, "dvec = ", dvec

end subroutine

$ gfortran-12 -Wall -c send2fortran.f90 
$ g++-12 -std=c++11 -Wall main.cpp send2fortran.o -lgfortran
$ ./a.out
 ivec =            0           0           0
 dvec =   -1.0000000000000000       -1.0000000000000000       -1.0000000000000000

While implicit conversions can make interfacing very simple in C++, they are hard to spot, which can make them surprising.Item C.164 in the C++ Core Guidelines suggests to avoid them entirely, or make them explicit:

  // Explicit conversion
  explicit operator f_multivec() { 
    return {n,ivec.data(),dvec.data()} ; 
  }

If we now try to compile this, we get a conversion error:

$ g++-12 -std=c++11 -Wall main.cpp send2fortran.o -lgfortran
main.cpp: In function 'int main()':
main.cpp:35:17: error: could not convert 's' from 'Multivec' to 'f_multivec'
   35 |    send2fortran(s);  // <-- sneaky implicit conversion
      |                 ^
      |                 |
      |                 Multivec

To make the code work we must convert explicitly:

int main() {
  auto s = Multivec(3);
  send2fortran(f_multivec(s));  // explicit conversion
}

The downside of this code is that the explicit conversion might give us a false feeling that s won’t be mutated. Even if we attempted to secure things using const modifiers, since the procedure is implemented in Fortran, the modifiers can be by-passed without any warning.

If the Fortran routine is called seldomly, you could simply construct the struct in plain-sight:

  int n = 10;
  std::vector<int> ivec(n,42);
  std::vector<double> dvec(n,-99.0);
  
  /* "Fortran" scope */ 
  { 
    send2fortran( 
      {n, ivec.data(), dvec.data()}
    );
  }

  int *iarray = new int[n];
  double *darray = new double[n];
  
  // ... set values ...

  /* "Fortran" scope */ 
  { 
    send2fortran( 
      {n, iarray, darray}
    );
  }

  // (!) don't forget to free memory
  delete [] iarray;
  delete [] darray;

This makes it obvious pointers to the data are being passed.

3 Likes