Allow complex pointers to real arrays and vice-versa

Language enhancement
1

Following recurrent discussions on the topic, I have prepared a draft proposal to allow complex pointers to real arrays and vice-versa:

Any feedback is appreciated.

2

You could actually emulate this feature, it seems, via:

call c_f_pointer( c_loc(c), r, [2] )

and similarly vice versa.

3

This is indeed what many people do (including me), but this is not valid Fortran. The standard explicitely says that in the above intruction r and c shall have the same type.

4

Consider this program:

program noncontig
   use, intrinsic :: iso_c_binding, only: c_loc, c_f_pointer
   implicit none
   real, target :: t(8)=[1,2,3,4,5,6,7,8]
   real, pointer :: a(:,:,:), rp(:,:,:)
   complex, pointer :: cp(:,:)
   character(*), parameter :: cfmt='(a,*(1x,f0.0))'
   !cp(1:1,1:4) => a    ! contiguous
   write(*,cfmt) ' t', t
   a(1:2,1:2,1:2) => t
   write(*,cfmt) ' a', a
   call c_f_pointer( c_loc(a), cp, [1,4] )
   write(*,cfmt) 'cp', cp

   !cp(1:1,1:2) => a(:,1:1,1:2) ! noncontiguous
   write(*,cfmt) ' a', a(1:2,1:1,1:2)
   rp => a(1:2,1:1,1:2)
   write(*,cfmt) 'rp', rp
   cp => cp(1:1,1:4:2)
   write(*,cfmt) 'cp', cp
end program noncontig

$ a.out
 t 1. 2. 3. 4. 5. 6. 7. 8.
 a 1. 2. 3. 4. 5. 6. 7. 8.
cp 1. 2. 3. 4. 5. 6. 7. 8.
 a 1. 2. 5. 6.
rp 1. 2. 5. 6.
cp 1. 2. 5. 6.

Note that all of the real pointer assignments are legal and are straightforward within the confines of the language. It is only the complex assignments that run afoul of the restrictions on the programmer.

However, through storage sequence association, the real and complex array elements are required by the standard to match up. Therefore it seems like an artificial restriction on the programmer that he is not allowed to take advantage of that storage sequence association.

5

Clear enough :innocent: So this enhancement would endorse the practice that is already there.

6

Indeed, this proposal essentially consists in standardizing the current practice.

Actually, to be consistent, passing a real actual argument to a dummy complex argument and vice versa should also be allowed, I think.

7

Actually, to be consistent, passing a real actual argument to a dummy complex argument and vice versa should also be allowed, I think.

I agree with that. We’ve used the following code for years:

do j=1,n
  do i=1,j-1
    c(i,j)=c(i,j)+zdotc(l,a(1,i),1,b(1,j),1)
  end do
  c(j,j)=c(j,j)+ddot(2*l,a(1,j),1,b(1,j),1)
end do

to set up a Hermitian matrix. The diagonal components are real so ddot can be used instead of zdot. The compiler obviously cannot know this optimization, but also has to be kept ignorant of the standards violation.

This would still work even if the real and imaginary parts are not interleaved but stored in blocks in the same contiguous array. However, it wouldn’t work if the compiler stored complex numbers as (magnitude, phase) pairs.

8

In the previous discussion on that topic, it turned out that such a storage was actually not possible (one of the reasons is that r => c%re would not be possible, then).

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.)

10

Probably, yes. But obsiouvsly a legal 100% Fortran solution is desirable.

1 Like
11

It is unfortunate that it is still not possible to write this operation in an equivalent standard-conforming way.

c(j,j)%re = c(j,j)%re + dot_product(a(:,j)%re,b(:,j)%re) + dot_product(a(:,j)%im,b(:,j)%im)

is close, but there are two function calls instead of one. Also, I remember a previous discussion about %re and %im where some compilers generate temp arrays for the argument association, which would be problematic for efficiency.