Pointer mapping between real and complex arrays

[FedericoPerini]

I do apologize if I bring back a topic that’s already been discussed several times, but I couldn’t find much discussion on the approach I’m trying to use here.

Working on a Schur decomposition function for LAPACK, I want to spare an allocation by splitting a complex output array into two real arrays, without additional internal reallocations/copies (see i.e. how dgees outputs two real arrays for the complex eigenvalues). I came up with this approach (try it at compiler explorer):

   complex(rk), allocatable, target :: arr(:)
   real(rk), pointer :: rreal(:),rimag(:)
   allocate(arr(10))
   rreal => arr%re      ! works in gfortran, wow!
   rimag => arr%im

I think this may be a great and legal way to cast one complex array to two reals, at least it seems consistent to a pointer to a component in a derived type, so I think it may be legal. However:

• works in gfortran and flang but not in ifx (internal compiler error) nor ifort (internal compiler error)
• In both gfortran and flang, there 8 bytes between real values of the pointer array, as I would expect the memory layout of the complex array to be [re,im,re,im,re,im,…]
• gfortran reports is_contiguous==.true. while flang reports is_contiguous==.false.. Does that mean that gfortran stores memory as [re,re,re,…,im,im,im,…]? or is it just a compiler bug?

So I would like to ask any of the compiler and standard experts on here what is your view! Clearly it does not seem ready for production at this point, (too much uncertainty in compiler support), and I will revert to having one slower additional allocation with data copy.

[Euler-37]

I think flang is correct. If use c_loc to point the real pointer to arr firstly,and then use c_f_pointer transfer the type, the result is correct.

program main
    use iso_c_binding
    implicit none
    complex,target::arr(2)
    real,pointer::p(:)
    real,pointer::pre(:),pim(:)
    arr(1)=(1.0,2.0)
    arr(2)=(3.0,4.0)
    call c_f_pointer(c_loc(arr),p,shape(arr)*2)
    pre=>p(::2)
    pim=>p(2::2)
    print*,pre,pim
    print*,is_contiguous(pre),is_contiguous(pre)
end program main

   1.00000000       3.00000000       2.00000000       4.00000000    
 F F

[ivanpribec]

From the Fortran Wringer Tests:

Take care with newer features

• The %RE and %IM syntax references components of complex variables. Some compilers just treat them as if they were expressions equivalent to REAL() and AIMAG() intrinsic function references. The distinction matters when %RE and %IM references are associated with REAL dummy arguments. Exactly one compiler correctly implements a pointer initialization target that applies %RE or %IM to an array slice.

---

The case of real to complex association/mapping with the whole complex array (and not just a subpart like c%re) will hopefully be resolved in F202Y:

• Paper 24-129: allow complex pointer to reals and vice versa (https://j3-fortran.org/doc/year/24/24-129.txt) by @PierU
• Formal Requirements US25 complex/real association (https://j3-fortran.org/doc/year/24/24-173.txt)

The minutes from the J3 meeting in October contain the following:

** Motion 24-173 "Formal Requirements US25 complex/real association"
   [Cohen](Cohen/Johnson)
   1. UC

Does anyone know what UC stands for?

---

For gfortran, there are some corner bugs with the complex c%ref syntax:

• 114023 – complex part%ref of complex named constant array cannot be used in an initialization expression.
• 114022 – ICE with a complex part%ref and nested structure constructor of complex array.

I tried adapting one of the tests:

   ! assume interleaved complex storage
   ! [re(1), im(1), re(2), im(2), ...]
   complex(kind=8),  target :: cmp1(3) = [(1,2),(3,4),(5,6)]
   real(kind=8), pointer :: rr(:) => null(), qq(:) => null()

   rr => cmp1%re
   qq => cmp1%im

   if (any(int(rr) /= [1,3,5])) stop 1
   if (any(int(qq) /= [2,4,6])) stop 2

   if (loc(rr) /= loc(cmp1)) stop 3
   if (loc(qq) /= loc(cmp1(1)%im)) stop 4

   if (is_contiguous(rr)) stop 5    ! stops here (gfortran 14)
   if (is_contiguous(qq)) stop 6

end

Looking at the gfortran dump tree, the stride is wrong:

  static complex(kind=8) cmp1[3] = {__complex__ (1.0e+0, 2.0e+0), __complex__ (3.0e+0, 4.0e+0), __complex__ (5.0e+0, 6.0e+0)};
  static struct array01_real(kind=8) qq = {.data=0B};
  static struct array01_real(kind=8) rr = {.data=0B};

  rr.span = 16;
  rr.dtype = {.elem_len=8, .version=0, .rank=1, .type=3};
  rr.dim[0].lbound = 1;
  rr.dim[0].ubound = 3;
  rr.dim[0].stride = 1;
  rr.data = (void *) &REALPART_EXPR <cmp1[0]>;
  rr.offset = -1;
  rr.dim[0].ubound = rr.dim[0].ubound + (1 - rr.dim[0].lbound);
  rr.offset = rr.offset - (1 - rr.dim[0].lbound) * rr.dim[0].stride;
  rr.dim[0].lbound = 1;

causing the contiguity test to fail later:

  L.8:;
  {
    struct array01_real(kind=8) D.4367;

    D.4367 = rr;
    if (D.4367.dim[0].stride == 1)
      {
        _gfortran_stop_numeric (5, 0);
      }
  }

[cmaapic]

UC - unanimous consent?

[FedericoPerini]

Thanks @ivanpribec @Euler-37, so flang is correct, and gfortran is also OK except for the pointer descriptor, that has a known bug on the stride.

With the C trick that @Euler-37 is suggesting, gfortran also works, probably because it takes the stride from the “manual” real-on-real slicing.

[ivanpribec]

@FedericoPerini

and gfortran is also OK except for the pointer descriptor, that has a known bug on the stride.

I believe a processor could choose to use separate complex storage ([re(1), re(2), ..., im(1), im(2),...] (or even something completely different).

The issue is not that absolute value of is_contiguous is .true. (that would be okay with separate storage of real and complex parts), but that it is not self-consistent with what the addresses show currently.

[FedericoPerini]

I agree with you Ivan, the apparent bug is not in the is_contiguous function, but rather in pointer descriptor creation. Although in practice, I believe all compilers use the (re,im) structure pattern.

[zedthree]

Yep, all compilers I’ve come across use this layout. C and C++ standards require it, in fact, so any compiler suite that includes both C and Fortran almost certainly uses it.

Maybe the committee should formally standardise this layout? C++ standardised two’s complement for similar reasons, as everyone was using it in practice.

[ivanpribec]

@zedthree

in fact, so any compiler suite that includes both C and Fortran almost certainly uses it.

It would certainly need to be used for complex(c_float) and complex(c_double). It does seem complicated to have a compiler with multiple layouts potentially requiring conversions back and forth with an associated runtime penalty.

There were some lengthy discussions on the topic of complex storage before:

• Complex type storage (again)
• Allow complex pointers to real arrays and vice-versa
• Implicit real - complex conversion in fortran

Split storage is not that uncommon in other software products:

• MATLAB versions R2017b and earlier use split storage (MATLAB Support for Interleaved Complex API in MEX Functions)
• the FFT routines in the Apple vDSP library expect split storage: Understanding data packing for Fourier transforms | Apple Developer Documentation

[PierU]

@ivanpribec

I believe a processor could choose to use separate complex storage ([re(1), re(2), ..., im(1), im(2),...] (or even something completely different).

With the current standard text I don’t think so, or at least it would be quite a hassle (and at least for the default complex type).

That means that in this code, a compiler has no other choice than the interleaved storage:

complex c(10)
real r(20)
equivalence (r,c)

[ivanpribec]

But couldn’t that also be a (phase, magnitude) representation?

[PierU]

I think it has been discussed in the discussions you have linked from here, and my conclusion is that although there’s no explicit prohibition of the (phase,magnitude) representation, the various constraints on the complex type make it impossible.

[septc]

   rr => cmp1%re
   qq => cmp1%im

   if (any(int(rr) /= [1,3,5])) stop 1
   if (any(int(qq) /= [2,4,6])) stop 2

   if (loc(rr) /= loc(cmp1)) stop 3
   if (loc(qq) /= loc(cmp1(1)%im)) stop 4

   if (is_contiguous(rr)) stop 5    ! stops here (gfortran 14)

  rr.span = 16;
  rr.dtype = {.elem_len=8, .version=0, .rank=1, .type=3};
  rr.dim[0].lbound = 1;
  rr.dim[0].ubound = 3;
  rr.dim[0].stride = 1;

In the above code, does rr.span = 16 mean that Gfortran represents the array pointer for the real parts as something like a stride-1 array with 8-byte padding between elements? If so, is_contiguous() should be modified as stride == 1 && span == elem_len or similar…?

(Also, if the above interpretation of “span”, “stride”, etc is correct, it seems simpler to just use stride = 2 and span = elem_len = 8 for array pointer, but I don’t know why it is not implemented so… Possibly, treated like components of derived-type arrays like objs(:) % foo?)

[zedthree]

I reviewed complex representations in a wide range of languages, file formats, and applications/libraries for my netCDF complex extension and found almost everyone uses (re, im) (interleaved) layout with Cartesian representation.

Radar/sonar people seem to like (phase, magntidue), and some applications use the split layout for storage on disk, but that’s about it.

[ivanpribec]

@septc

In the above code, does rr.span = 16 mean that Gfortran represents the array pointer for the real parts as something like a stride-1 array with 8-byte padding between elements? If so, is_contiguous() should be modified as stride == 1 && span == elem_len or similar…?

The descriptor contents can be found here: gcc/libgfortran/libgfortran.h at 70999668a1305571d3b5fc57168fcb060a976418 · gcc-mirror/gcc · GitHub

Looking at the commit message that introduced .span it is a property set only for pointer arrays. Skimming through the source, it looks like the array assignment a => c gets handled by gfc_conv_expr_descriptor which internally calls gfc_get_array_span: gcc/gcc/fortran/trans-array.cc at 70999668a1305571d3b5fc57168fcb060a976418 · gcc-mirror/gcc · GitHub

[RonShepard]

@ivanpribec

But couldn’t that also be a (phase, magnitude) representation?

I think this has been discussed before, and the answer is no. The f77 standard has this text:

4.6 Complex Type
A complex datum is a processor approximation to the value of a complex number.
The representation of a complex datum is in the form of an ordered pair of real
data. The first of the pair represents the real part of the complex datum and the
second represents the imaginary part. Each part has the same degree of
approximation as for a real datum. A complex datum has two consecutive
numeric storage units in a storage sequence; the first storage unit is the real part
and the second storage unit is the imaginary part.

Later in the description of function results in Table 5, a complex function result is described as

(6) A complex value is expressed as an ordered pair of reals, (ar,ai), where ar is the real part and ai is the imaginary part.

Then in 17.1.1 where storage sequence association is discussed, there is

A variable or array element of type double precision or complex has a storage
sequence of two numeric storage units. In a complex storage sequence, the real
part has the first storage unit and the imaginary part has the second storage unit.

A little further in 17.1.3, there is a figure showing the storage layout of the declarations

REAL A(4),B
COMPLEX C(2)
DOUBLE PRECISION D
EQUIVALENCE (C(2),A(2),B), (A,D)
...
storage unit   |     1     |     2     |     3     |     4     |      5    |
               |    ----C(1)----       |       ----C(2)----    |
                           |  A(1)     |   A(2)    |   A(3)    |    A(4)   |
                                       |   --B--   |
                           |      -----D-----      |

That last part was typed manually, not cut and paste, but I think it reproduces the f77 figure correctly.

As discussed previously regarding this topic, this text regarding storage sequence association has been shortened in subsequent standards (which I don’t have availble at this moment for comparison), so it is less explicit. Also, f77 had only one complex type+kind, while f90 and later allows multiple kinds of the complex type, so the storage sequence association of these other kinds has never been explicitly laid out like this in such a simple and clear way. Finally, f77 did not allow (or show) storage sequence association between character types and numeric types (they were not even allowed to exist in the same common block), so I think that may still be undefined within the fortran standards.

[PierU]

@RonShepard

Then in 17.1.1 where storage sequence association is discussed, there is

A variable or array element of type double precision or complex has a storage sequence of two numeric storage units. In a complex storage sequence, the real part has the first storage unit and the imaginary part has the second storage unit.

It seems that this constraint is no longer in the current standard.

[PierU]

@FedericoPerini

So I would like to ask any of the compiler and standard experts on here what is your view! Clearly it does not seem ready for production at this point, (too much uncertainty in compiler support), and I will revert to having one slower additional allocation with data copy.

I’m no standard expert, but the standard actually says very little about the “complex part designators”.

To me, among the 4 Fortran compilers on godbolt.org, only flang gets it completely right. gfortran wrongly reports array sections made of complex part designators as contiguous (and always make a temporary copy when they are passed as arguments, even if the dummy argument is assumed shape), ifx messes up everything, and nvfortran doesn’t compile pointers to complex part designators.

[ivanpribec]

@Euler-37

   1.00000000       3.00000000       2.00000000       4.00000000    
 F F

NAG Fortran output:

~/complex> nagfor test_euler-37.f90 
NAG Fortran Compiler Release 7.2(Shin-Urayasu) Build 7203
[NAG Fortran Compiler normal termination]
~/complex> ./a.out
   1.0000000   3.0000000   2.0000000   4.0000000
 F F

Also the tests from above pass:

   if (any(int(rr) /= [1,3,5])) stop 1
   if (any(int(qq) /= [2,4,6])) stop 2

!   if (loc(rr) /= loc(cmp1)) stop 3
!   if (loc(qq) /= loc(cmp1(1)%im)) stop 4

   if (is_contiguous(rr)) stop 5
   if (is_contiguous(qq)) stop 6

(nagfor doesn’t appear to have the non-standard loc intrinsic, so I had to comment those lines)

[cmaapic]

Through C interop you could use loc.