It’s certainly possible. Previously, I’ve proposed a Fortran/C++ compatibility library (Fortran/C++ compatibility library · Issue #325 · fortran-lang/stdlib · GitHub) but the issue didn’t get much attention. Presumably there is low community interest.
At least for vectors of intrinsic kinds, it’s fairly simple to have iterator access when working from C++. This also allows you to use STL routines which rely on iterator access.
For tensors (arrays of rank 2 and higher) the situation is more precarious. If you consider how many C++ matrix libraries are out there: Eigen, Armadillo, Blaze, Boost UBLAS,… , you may infer there are (too) many different ways to achieve the same thing with various trade-offs with respect to expressivity and performance. The amount of effort which has gone into these libraries is much higher than what a single programmer could meaningfully spend without backing from an institution or unless persistent enough to develop them over a long time period.
I would note that Eigen and also Armadillo have facilities to interface with a raw C buffer. This can be used also for Fortran arrays. In both cases of C/Fortran calling you need to be careful with who allocates the memory. Here’s an example of using Eigen from Fortran:
/* multiply_by_two.cc */
#include<Eigen/Dense>
#ifdef __cplusplus
extern "C" {
#endif
// Multiplies a C or Fortran array by the scalar 2
//
// Inputs:
// nx, ny - dimensions of the array
// a - pointer to a float array
//
void multiply_by_two(int nx, int ny, float* a) {
Eigen::Map<Eigen::Matrix<float,Eigen::Dynamic,Eigen::Dynamic>> a_m(a,nx,ny);
a_m *= 2;
}
#ifdef __cplusplus
}
#endif
! call_eigen.f90
!
program call_eigen
use iso_c_binding, only: c_float, c_int
implicit none
interface
subroutine multiply_by_two(nx,ny,a) bind(c,name="multiply_by_two")
import c_float, c_int
integer(c_int), intent(in), value :: nx, ny
real(c_float), intent(inout) :: a(nx, ny)
end subroutine
end interface
real(c_float) :: a(3,2)
a = reshape([real(c_float) :: 1, 2, 3, 4, 5, 6], [3, 2])
call multiply_by_two(3,2,a)
print *, a
end program
$ g++ -c multiply_by_two.cc -I /usr/include/eigen3/
$ gfortran call_eigen.f90 multiply_by_two.o -lstdc++
$ ./a.out
2.00000000 4.00000000 6.00000000 8.00000000 10.0000000 12.0000000
If interested I also have an example of defining a sparse matrix in CSR format in Fortran, and then using one of Eigen’s iterative solvers to solve for a given RHS.
With the enhanced interoperability features, the Fortran interface could be
subroutine multiply_by_two(a) bind(c,name="multiply_by_two")
import c_float
real(c_float), intent(inout), contiguous :: a(:,:)
end subroutine
and on the C++ side:
void multiply_by_two(CFI_cdesc_t* a) {
const void * a_ptr = a->base_addr;
const CFI_index_t nx = a->dim[0].extent;
const CFI_index_t ny = a->dim[1].extent;
Eigen::Map<Eigen::Matrix<float,Eigen::Dynamic,Eigen::Dynamic>> a_m(a,nx,ny);
a_m *= 2;
}
Ideally, you would write a small helper template function, which would return the Eigen::Map instance, something like:
auto a_m = asEigenMat<float>(a); // returns Eigen::Map instance