Fortran/C++ compatibility library

[ivan-pi]

With the introduction of Fortran 2018 enhancements to interoperability with C it is possible to write C (and C++) code which interoperates with

• optional dummy arguments
• assumed-length character dummy arguments
• assumed-shape arrays
• allocatable dummy arguments
• pointer dummy arguments

To handle such objects the standard introduced the concept of C descriptors which are defined in the header "ISO_Fortran_binding.h" header provided by the Fortran processor.

While very powerful, using these low-level tools to access Fortran objects feels clunky and error-prone. It would be nice to provide some higher level tools in the form of templated C++ wrapper classes. A similar concept is followed by the popular and fast-growing Rcpp library for R and C++ integration.

Here is a minimal example of such a wrapper class (kudos to @LKedward for helping me overcome a problem with C descriptor lifetimes in the Discourse thread) :

// fortran_array.cpp

#include "ISO_Fortran_binding.h"
#include <iostream>

template<typename T>
struct FortranArray {
  typedef T real_t;

  const void *base_addr;
  const CFI_dim_t *dim;

    const CFI_cdesc_t *obj; // A const reference to rank 2 Fortran array (can be discontiguous with regular strides).

    FortranArray(const CFI_cdesc_t *obj_) : : base_addr(obj_->base_addr), dim(obj_->dim) {}

    // Returns the i'th component of the j'th column in the array:
    inline T operator()(const size_t i, const size_t j) const {
        char *elt = (char *) base_addr;
        elt += (i*dim[0].sm + j*dim[1].sm);
        return *(T *) elt;
    }

    size_t size(const size_t d) {
        return obj->dim[d].extent;
    }

};

extern "C" {

void print_array(CFI_cdesc_t *array) {

    // View into a Fortran array
    FortranArray<double> a(array);

    for (size_t i = 0; i < a.size(0); i++) {
        for (size_t j = 0; j < a.size(1); j++) {
            std::cout << a(i,j) << " ";
        }
        std::cout << "\n";
    }
}

}

Given the popularity of C++ and it's strong presence in scientific computing, I think these tools would go a long way to increase adoption of Fortran as a viable co-language. Ideally, one would implement an STL-compatible container interface, allowing to use C++ algorithms directly on Fortran arrays.

A few questions:

• Do you see benefit in tighter integration between Fortran and C++?
• Do you think we can find support from the C++ community?
• Should this go into a repository separate from stdlib?
• Does it make sense to port the Fortran run-time library routines? (e.g. for a Fortran array living in C++, do we just call std::sin(array(i)) or a Fortran version fortran::sin(array(i)))

[ivan-pi]

As a small further example, I have added an attribute parameter to the FortranArray class:

namespace cfi {
  typedef CFI_attribute_t attribute_t;
  enum attribute {
    allocatable = CFI_attribute_allocatable,
    pointer = CFI_attribute_pointer,
    other = CFI_attribute_other
  };
}

template <typename T, cfi::attribute_t attribute = cfi::attribute::other>
struct FortranArray {

  // ...

  // Constructor
  FortranArray(const CFI_cdesc_t *obj) :
      base_addr(obj->base_addr), dim(obj->dim) {
        // Make sure true attribute and declared attribute match
        assert(obj->attribute == attribute); 
      }
}

using cfi::attribute::allocatable;

void process_allocatable_array(const CFI_cdesc_t *array) {

    FortranArray<double, allocatable> a(array); // If attributes don't match, assertion fails.

}

E.g. trying to call process_allocatable_array with a non-allocatable arrays leads to the error:

main_foo: foo.cpp:31: FortranArray<T, attribute>::FortranArray(const CFI_cdesc_t*) [with T = double; signed char attribute = 1; CFI_cdesc_t = CFI_cdesc_t]: Assertion `obj->attribute == attribute' failed.

Program received signal SIGABRT: Process abort signal.

The error message could be improved with some preprocessor syntax.

Similar run-time checks of the shadow C++ class could be done for the rank and type of the array. Type-checking can be done like this:

        // Check type
        if (std::is_same<real_t, double>::value) {
          assert(obj->type == CFI_type_double);
        } else if (std::is_same<real_t, float>::value) {
          assert(obj->type == CFI_type_float);
        } else {
          std::abort();
        }

[ivan-pi]

So I just added the following iterator code:

  struct Iterator
  {
    using iterator_category = std::forward_iterator_tag;
    using difference_type   = std::ptrdiff_t;
    using value_type        = real_t;
    using pointer           = real_t *;
    using reference         = real_t &;

    Iterator(pointer ptr) : m_ptr(ptr) {}

    reference operator*() const { return *m_ptr; }
    pointer operator->() { return m_ptr; }
    Iterator& operator++() { m_ptr++; return *this; }
    Iterator operator++(int) { Iterator tmp = *this; ++(*this); return tmp; }
    friend bool operator== (const Iterator& a, const Iterator& b) { return a.m_ptr == b.m_ptr; };
    friend bool operator!= (const Iterator& a, const Iterator& b) { return a.m_ptr != b.m_ptr; };

  private:
    pointer m_ptr;
  };

  Iterator begin() { return (address(0)); }
  Iterator end()   { return (address(this->size())); }

  inline real_t *address(const size_t i) {
    real_t *val = const_cast<real_t *>(
        static_cast<const real_t*>(this->base_addr));
    return &val[i];
  }

into the FortranArray class (no guarantees this is legitimate C++...; I just followed some tutorial).

A subroutine to multiply by a scalar can be written as:

void multiply_by_two(CFI_cdesc_t *array) {
  Vector f_array(array);
  for (auto &element: f_array) {
    element *= 2.0;
  }
}

For the demo program:

  A = [1,2,3,4,5]

  print *, A
  call multiply_by_two(A)
  print *, A

I get the output:

$ ./main_foo 
   1.0000000000000000        2.0000000000000000        3.0000000000000000        4.0000000000000000        5.0000000000000000     
   2.0000000000000000        4.0000000000000000        6.0000000000000000        8.0000000000000000        10.000000000000000   

[LKedward]

Given the popularity of C++ and it's strong presence in scientific computing, I think these tools would go a long way to increase adoption of Fortran as a viable co-language. Ideally, one would implement an STL-compatible container interface, allowing to use C++ algorithms directly on Fortran arrays.

I agree completely. This was my thinking when you posted on the Discourse; this will definitely be useful to have in a library. I'm tempted to say it should be its own library separate to stdlib.

[ivan-pi]

I'm interested in whether this would also allow straightforward usage of the oneAPI parallel STL directly on the memory of Fortran objects.

[ivan-pi]

I just seized GTC 2021 as an opportunity to ask @brycelelbach about using C++ standard parallel algorithms on Fortran array objects using the new NVC++ and NVFORTRAN compilers. His answer:

Yes, that should be possible. If it's dynamically allocated memory and you have automatic unified memory enabled.

Great way to get the best of both worlds.