Allow complex pointers to real arrays and vice-versa

Language enhancement
9

Would it be more “legal” to implement the type-casting in C or C++?

// complex_cast.cpp

#include <complex>
#include <ISO_Fortran_binding.h>

extern "C"
void c_f_complex_1(CFI_cdesc_t *c, CFI_cdesc_t *d) {

// Fortran -- complex(c_double)
// --> C   -- double _Complex
// --> C++ -- std::complex<double>

    auto cptr = static_cast<std::complex<double> *>(c->base_addr);
    auto dptr = reinterpret_cast<double *>(cptr);
    CFI_index_t dshape[1] = { c->dim[0].extent * 2 };

    CFI_establish(d,
        dptr,
        CFI_attribute_pointer,
        CFI_type_double,
        0 /* ignorer */,
        (CFI_rank_t) 1,
        dshape);
}

extern "C"
int c_f_complex_associated(CFI_cdesc_t *c, CFI_cdesc_t *d) {
    return c->base_addr == d->base_addr;
}

In C++, the layout of the complex type is constrained to allow array-oriented access using a reinterpret_cast: std::complex - cppreference.com, and for compatibility with the C layout.

In C the following applies (see Arithmetic types - cppreference.com):

Each complex type has the same object representation and alignment requirements as an array of two elements of the corresponding real type (float for float complex, double for double complex, long double for long double complex). The first element of the array holds the real part, and the second element of the array holds the imaginary component. [emphasis added]

Judging by this description, the Fortran type complex(c_double), which is C interoperable, can not be implemented as a (magnitude, phase) pair, but must hold floating point values of the real and imaginary components. So as long as kind(c_double) == kind(1.0d0) is true, we should be able to get away with it.

For sake of completeness, here is the Fortran program, demonstrating the C++ function above:

! test.f90
program test
use, intrinsic :: iso_c_binding

interface c_f_complex
    subroutine c_f_complex_1(c,d) bind(c)
        import c_double
        complex(c_double), pointer, intent(in) :: c(:)
        real(c_double), pointer, intent(out) :: d(:)
    end subroutine
end interface

complex(c_double), target :: a(4)
real(c_double), pointer :: b(:) => null()

character(*), parameter :: cfmt='(a,*(1x,f0.0))'

call c_f_complex(a, b)
print *, associated(b)
!print *, associated(b,a) ! Not allowed, pointer and target must have same type
print *, c_f_complex_associated(b,a)

b = [1,2,3,4,5,6,7,8]

write(*,cfmt) 'a', a

contains

    ! Check association status between a real pointer and complex target
    logical function c_f_complex_associated(pointer,target) result(stat)
        real(c_double), pointer, intent(in) :: pointer(:)
        complex(c_double), pointer, intent(in) :: target(:)
        interface
            function c_f_complex_associated_(c,d) &
                    bind(c,name="c_f_complex_associated")
                import c_double, c_int
                complex(c_double), pointer, intent(in) :: c(:)
                real(c_double), pointer, intent(in) :: d(:)
                integer(c_int) :: c_f_complex_associated_
            end function
        end interface
        ! Convert the C int result to a Fortran logical scalar
        stat = c_f_complex_associated_(target,pointer) == 1_c_int
    end function

end program

(Caveat: I didn’t include any status checks against misuse of attributes, but the strong Fortran typing should take care of this in principle.)