Example of calling Fortran from C++ using the enhanced C/Fortran interoperability

cxxfi.hpp

#pragma once

#include <complex>
#include <vector>
#include <array>

#include <type_traits>

#include <span> // C++ 20
#include <iostream>

#include <ISO_Fortran_binding.h> // F2018
// for a peek into the internals from a particular vendor:
// https://github.com/gcc-mirror/gcc/blob/master/libgfortran/ISO_Fortran_binding.h

namespace cfi {

namespace cfi_internal {

/**
 * Function templates to convert a template argument 
 * into an integer type code
 */ 
template<typename T>
constexpr CFI_type_t type();

// Non-interoperable structure (default)
template<typename T>
constexpr CFI_type_t type(){ return CFI_type_other; } 

// Specializations for supported types
template<>
constexpr CFI_type_t type<char>(){ return CFI_type_char; }
template<>
constexpr CFI_type_t type<bool>(){ return CFI_type_Bool; } 
template<>
constexpr CFI_type_t type<float>(){ return CFI_type_float; } 
template<>
constexpr CFI_type_t type<double>(){ return CFI_type_double; } 
template<>
constexpr CFI_type_t type<std::complex<float>>(){ return CFI_type_float_Complex; } 
template<>
constexpr CFI_type_t type<std::complex<double>>(){ return CFI_type_double_Complex; }
template<>
constexpr CFI_type_t type<int>(){ return CFI_type_int; } 
template<>
constexpr CFI_type_t type<long>(){ return CFI_type_long; } 
template<>
constexpr CFI_type_t type<long long>(){ return CFI_type_long_long; } 
template<>
constexpr CFI_type_t type<std::size_t>(){ return CFI_type_size_t; }
template<>
constexpr CFI_type_t type<std::int8_t>(){ return CFI_type_int8_t; }
template<>
constexpr CFI_type_t type<std::int16_t>(){ return CFI_type_int16_t; }
//template<>
//constexpr CFI_type_t type<std::int32_t>(){ return CFI_type_int32_t; }
//template<>
//constexpr CFI_type_t type<std::int64_t>(){ return CFI_type_int64_t; }
template<>
constexpr CFI_type_t type<void*>(){ return CFI_type_cptr; }


} // namespace cfi_internal

/**
 *  Enumerator class for attributes
 */
enum class attr : CFI_attribute_t {
    other       = CFI_attribute_other,
    allocatable = CFI_attribute_allocatable,
    pointer     = CFI_attribute_pointer
};

/**
 *  C++-descriptor class encapsulating Fortran array
 */
template<typename T, int rank_ = 1, attr attr_ = cfi::attr::other,
    bool contiguous = false>
class cdesc_t {
public:

    static_assert(rank_ >= 0, "Rank must be non-negative");
    static_assert(rank_ <= CFI_MAX_RANK, 
        "The maximum allowed rank is 15");

    using value_type = T;
    using size_type = std::size_t;

    using reference = T&;
    using const_reference = const T&;

    using pointer = T*;
    using const_pointer = const T*;

    constexpr CFI_type_t type() const { return cfi_internal::type<T>(); };
    constexpr CFI_rank_t rank() const { return rank_; };

    // Version of ISO_Fortran_binding.h 
    constexpr int version() const { return this->get()->version; }

    // Element length in bytes
    std::size_t elem_len() const { return this->get()->elem_len; }


    // Base constructor for assumed-shape arrays
    cdesc_t(T* ptr, size_t n) :
        ptr_(this->get_descptr())
    {
        static_assert(rank_ == 1,
            "Rank must be equal to 1 to construct descriptor from pointer and length");

        CFI_index_t extents[1] = { static_cast<CFI_index_t>(n) };

        [[maybe_unused]] int status = CFI_establish(
            this->get(),
            ptr,
            static_cast<attribute_type>(attr_),
            this->type(),
            sizeof(T),
            rank_,
            extents
        );

        //std::cout << "Descriptor has been established\n";

    }

    // Constructor for static array
    template<std::size_t N>
    cdesc_t(T (&ref)[N]) : cdesc_t(ref,N) {
        static_assert(attr_ == cfi::attr::other);
        static_assert(rank_ == 1,
            "Rank must be equal to 1 to construct descriptor from static array");
    }

    // Constructor from std::vector
    cdesc_t(std::vector<T> &buffer) : cdesc_t(buffer.data(),buffer.size()) {
        static_assert(attr_ == cfi::attr::other);
        static_assert(rank_ == 1,
            "Rank must be equal to 1 to construct descriptor from std::vector");
    }

    // Constructor from std::array
    template<std::size_t N>
    cdesc_t(std::array<T,N> &buffer) : cdesc_t(buffer.data(),N) {
        static_assert(attr_ == cfi::attr::other);
        static_assert(rank_ == 1,
            "Rank must be equal to 1 to construct descriptor from std::array");
    }

    // Constructor from an external descriptor 
    // (typically a dummy argument originating in Fortran)
    cdesc_t(CFI_cdesc_t *desc) : ptr_(desc) {

        // TODO: Runtime checks for matching type, rank, and attributes
        // TODO: Assumed-size should be forbidden (no size available)

    };

    // Implicit cast to C-descriptor pointer
    operator CFI_cdesc_t* () { return this->get(); }

    // Implicit cast to std::span (only for rank-1 arrays, 
    // otherwise use .flatten())
#if __cpp_lib_span
    operator std::span<T> const () {
        static_assert(rank_ == 1,
            "Rank must be equal to 1 to convert implicitly to std::span");
        size_type n = this->get()->dim[0].extent;
        return {this->data(),n};
    }

    // Return a flattened view of the array.
    // TODO: handle assumed-size arrays
    std::span<T> flatten() {
        size_type nelem = 1;
        for (int i = 0; i < rank_; ++i) {
            nelem *= this->get()->dim[i].extent;
        }
        return {this->data(),nelem};
    }
#endif

    // Return pointer to the underlying descriptor
    // (the name get() is inspired by the C++ smart pointer classes)
    constexpr auto get(){ return ptr_; }

    // Array subscript operators
    T& operator[](std::size_t idx) {
        static_assert(rank_ == 1,
            "Rank must be 1 to use array subscript operator");
            return *(data() + idx);
        }
    const T& operator[](std::size_t idx) const {
        static_assert(rank_ == 1,
            "Rank must be 1 to use array subscript operator");
        return *(data() + idx);
    }

    // Iterator support
    T* begin() { return &this->operator[](0); }
    T* end() { return &this->operator[](this->get()->dim[0].extent); }
    const T* cbegin() { return &this->operator[](0); }
    const T* cend() { return &this->operator[](this->get()->dim[0].extent); }

private:

    using attribute_type = typename std::underlying_type<attr>::type;

    // Return base address of the Fortran array
    // cast to the correct type
    inline pointer data() { return static_cast<pointer>(this->get()->base_addr); }

    constexpr auto get_descptr() {
        return reinterpret_cast<CFI_cdesc_t *>(&desc_);
    } 

private:
    // Descriptor containing the actual data
    CFI_CDESC_T(rank_) desc_;
    // Pointer to opaque descriptor object
    CFI_cdesc_t *ptr_{nullptr};
};

} // namespace cfi

dot.f90

!> dot
!>     
!> Dot product of two contiguous arrays
!> The arrays must have the same length
!>
real(c_double) function dot(a,b) bind(c,name='dot')
   use, intrinsic :: iso_c_binding, only: c_double
   implicit none
   real(c_double), intent(in), contiguous :: a(:), b(:)
   dot = dot_product(a,b)
end function

main.cpp

#include <array>
#include <iostream>

#include "cxxfi.hpp"
using namespace cfi;

// real(c_double) function dot(a,b) bind(c,name='dot')
//    use, intrinsic :: iso_c_binding, only: c_double
//    real(c_double), intent(in), contiguous :: a(:), b(:)
//    dot = dot_product(a,b)
// end function
extern "C" double dot(CFI_cdesc_t *a, CFI_cdesc_t *b);

int main() {
    
   std::array<double,3> a{1.,2.,3.};
   std::vector<double> b = {1.,2.,3.};

   { /* fortran operations */

      // cfi::cdesc_t<type,rank> f_obj(cpp_obj);
      // This class is an adaptor from C++ to Fortran with explicit type and rank

      cfi::cdesc_t<double,1> fa(a);
      cfi::cdesc_t<double,1> fb(b);

      // 1^2 + 2^2 + 3^2 = 1 + 4 + 9 = 14
      std::cout << "dot(a,b) = " << dot(fa,fb) << '\n';

      // Range-based for loop  
      for (auto &bitem : fb) {
         bitem += 1.0;
      }

      // 1*2 + 2*3 + 3*4 = 2 + 6 + 12 = 20
      std::cout << "dot(a,b) = " << dot(fa,fb) << '\n';

      for (auto bitem : fb) {
         std::cout << bitem << '\n';
      }
   }

   return 0;
}