Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
25 Views

reinterpreting complex as reals

Consider the following trick,

module xxx
  implicit none
  integer, parameter :: stdprec=kind(1.d0)
contains
  function conv2real(carr)
    use, intrinsic :: iso_c_binding
    ! returns real pointer to a complex array
    complex(stdprec), contiguous, intent(inout), target  :: carr(:)
    real(stdprec),contiguous,pointer :: conv2real(:)
    call c_f_pointer(c_loc(carr),conv2real,[size(carr)*2])
  end function conv2real
end module xxx

The function conv2real takes a complex array as argument and , with the help of c binding module, returns a real pointer to it. 

First question: is Fortrand standard OK with this technique? I think it should be.

Second question: Is this violating ansi aliasing? I think not, since complex numbers are just pairs of reals, no weird reinterpreting is going on. 

 

Now, consider the following test case

module xxx
  implicit none
  integer, parameter :: stdprec=kind(1.d0)
contains
  function conv2real(carr)
    use, intrinsic :: iso_c_binding
    ! returns real pointer to a complex array
    complex(stdprec), contiguous, intent(inout), target  :: carr(:)
    real(stdprec),contiguous,pointer :: conv2real(:)
    call c_f_pointer(c_loc(carr),conv2real,[size(carr)*2])
  end function conv2real
end module xxx

program zorglub
  use xxx
  implicit none
  complex(stdprec), contiguous, pointer :: x(:)
  real(stdprec), contiguous, pointer :: rx(:)
  type ccc
    complex(stdprec), contiguous, pointer :: p(:)
    real(stdprec), contiguous, pointer :: rp(:)
  end type ccc
  type(ccc) :: o

  !sequence of steps: allocate, point real to complex, then initialise complex
  allocate(x(100))
  rx=>conv2real(x)
  x(1)=1.23_stdprec
  print *, x(1)
  print *, rx(1) !real part of x(1) == rx(1) CORRECT

  
  !sequence of steps: allocate, point real to complex, then initialise complex
  allocate(o%p(100))
  o%rp=>conv2real(o%p)
  o%p(1)=1.23_stdprec
  print *, 'p', o%p(1)
  print *, 'rp', o%rp(1) !real part of o%p(1) /= o%rp(1) !!! WRONG (BUG?)

  !do some clean up and reuse
  deallocate(o%p)
  nullify(o%p)
  nullify(o%rp)
  
  !sequence of steps: allocate, then initialise complex, then point real to complex
  allocate(o%p(100))
  o%p(1)=1.23_stdprec
  o%rp=>conv2real(o%p)
  print *, 'p', o%p(1)
  print *, 'rp', o%rp(1) !real part of o%p(1) == o%rp(1) !!! CORRECT

end zorglub

I am first allocating the complex x pointer , then pointing rx to it. Afterwards, I initialise x(1) and print both x and rx. Both yield the same , expected result.

However, when I make a derived type containing a complex and a real array pointer, repeat the previous procedure, I get different values for o%p(1) and o%rp(1) . I think this is a bug.

The last bit of the test is a printout of o%p and o%rp, but this time, the first value of o%p is initialised _before_ the o%rp is made to point to o%p . So, the steps are reversed in respect to the preceding examples. Now again the pointer and its target have the same value at the first index. 

Lastly, when I replace 'target' in function conv2real with 'pointer', everything seem to work correctly.

 

The ifort compiler version is 16.0.2

0 Kudos
1 Reply
Highlighted
25 Views

The first bit of code doesn't violate the standard as the dummy argument has the TARGET attribute.

Your full program (which contains a minor syntax error) gives me identical results for all cases in 17.0.1, but not in 16.0.2 (which is what I happen to have on this system.)

0 Kudos