Complex pointer to a real array and vice-versa
PierUgit opened this issue · 23 comments
In many codes that use FFTs and the Fourier domain one need to alternatively work on the same array as real or as complex numbers. Different tricks have to used, none being really satisfactory, e.g.:
real, allocatable :: r(:)
! ...
! some computations on r(:)...
call rfft(r) ! real to complex in-place transform
call compute_complex(r,size(r)) ! this routine MUST NOT have any interface !!!
call irfft(r) ! complex to real in-place inverse transform
! ...
subroutine compute_complex(c,n)
complex :: c(n/2)
! some computations on c(:)
end subroutineInstead, the standard should provide a way to have a complex pointer to a real array:
real, allocatable, target :: r(:)
complex, pointer :: c(:)
! ...
! some computations on r(:)...
call rfft(r) ! real to complex in-place transform
c => cmplx_pointer(r)
! some computations on c(:)...
call irfft(r) ! complex to real in-place inverse transform
! ...And vice-versa
complex, allocatable, target :: c(:)
real, pointer :: r(:)
! ...
r => real_pointer(c)
! some computations on r(:)...
call rfft(r) ! real to complex in-place transform
! some computations on c(:)...
call irfft(r) ! complex to real in-place inverse transform
! ...Of course some constraints would be needed:
- the size of the real array in the first dimension shall be even
- the real or complex array shall have a unit stride along the first dimension (forcing the whole array to be contiguous could be desirable, so that the compiler could catch the error at compile-time)
c(:)%reshall be equivalent tor(1::2)andc(:)%imshall be equivalent tor(2::2), which enforces the real component to be stored first (de-facto standard in all compilers that I know).
One thing that might be relevant here: often times in FFT it is possible to compute the transformation faster if all the real parts are stored first, and then all the imaginary parts. In other words, if c(:)%re is equivalent to r(1:n/2) and c(:)%im is equivalent to r(n/2+1:n). It would be cool if Fortran compilers allowed this alternative layout of complex arrays, perhaps with some attribute such as complex(dp), layout(interleaved), allocatable :: c(:) and layout(block). It seems the language itself almost allows it, except some corner cases.
@certik This goes beyond and looks more complex than my proposal, which is actually very simple.
An alternative syntax has been proposed on the discourse group:
c => complex :: r(:)
Your proposal is only simple if you assume that Fortran prescribes an ABI/layout. In most cases however, the Fortran standard does not prescribe a particular layout.
I don't know, but I think that the memory layout for the complex type is already significantly constrained by the existing rules.
About the layout(block) that you are considering, it would mean that an elementary type could be stored non contiguously. I have to admit that I can't find anything that disallow such a layout, but it really looks weird to me...
- the real or complex array shall have a unit stride along the first dimension (forcing the whole array to be contiguous could be desirable, so that the compiler could catch the error at compile-time)
One difficulty here is if the pointed object is itself a pointer: the compiler has then no mean to check the contiguity at compile time.
EDIT not a problem actually, I forgot that pointers could have the contiguous attribute...
I think that we should extend Fortran to allow to select a memory layout of arrays. Just like in C++ the array support (std::mdspan as well as Kokkos) allows you to select the layout. Now coming to your proposal, we should allow the complex to real array casting, and the order that you get would then depend on the layout. Right now, I think the standard prescribes somewhere (I can't find the exact place) that the complex array layout has to be interleaved.
All I am saying is to keep the door open to allow to select array layouts and for all the features to still work.
Putting aside my proposal, I have some doubt about the feasibility of a "blocked layout" for the complex type:
complex, layout(block), target :: c(n) ! n > 1
complex, pointer :: s ! s is a scalar
s => c(1)With this layout the real part of c(1) would be for instance at address p and the imaginary part at address p+4*n: why not, but it would be also the same for the s scalar. Again it looks weird to me to have a scalar object that is stored non contiguously...
I think this last case is probably why the interleaved layout is currently used. You can imagine passing c(1) to a function that only expects a complex scalar. In the blocked layout it would require a copy.
I am wondering if there could be an alignment issue here...
allocate( r(1000) )
c => r(2:999)
Compilers are supposed to handle the above case anyway, because they are also supposed to handle the similar case:
equivalence (r(2),c(1))
But can it have performance issues? Maybe aligning complex numbers on odd 4-bytes boundaries (that is on 4*(2*k+1) addresses) is a performance killer... Just wondering, I have no clue about that.
If it was, maybe only r => c should be allowed, and not c => r
This initial paper has passed at the june 2024 meeting (https://j3-fortran.org/doc/meeting/233).
Now, what's next?
@PierUgit I think the committee then requires an edits paper, or possibly a few other papers before that (requirements, specification, etc.). @everythingfunctional, can you please summarize what the next steps for @PierUgit should be?
The following aspects need to be addressed. These could be separate papers, but it's also not uncommon to combine 1 and 2, or 3 and 4.
- Use Cases: Outline some examples where the feature would be required, or at least would be highly beneficial, to solving a particular problem
- Requirements: Based on the use cases, define what is needed for a solution to the problem.
- Specifications: Define the characteristics/aspects of a feature that satisfies the previously outlined requirements.
- Syntax: Define a syntax and its semantics that are congruent with the previously defined specifications
- Edits: Define the exact changes to the standard, with page and line number, that would incorporate the feature as outlined by items 3 and 4.
You'll need to coordinate submitting these papers with Malcolm Cohen, the head of the Data Subgroup. He is the one who will bring the paper(s) forward to the committee for vote. Ideally you would request to be a member of the committee (you can be added as an alternate to someone else if your organisation is unwilling to pay the membership fee), and attend the meetings where your papers are discussed.
Thanks @everythingfunctional! @PierUgit go ahead and become a member of the committee as an alternate, and then you can start with the papers 1.-4., and as Brad said, you can combine 1 and 2, as well as 3 and 4. You can submit them first as PRs and ask others for feedback.
Could you point me towards typical papers illustrating the differents points? "Uses cases" and "syntax" is clear to me, but "requirements" and "specifications" not so much...
@PierUgit after you understand it well, can you also please send a PR against this repository and document this process with examples in the README? That would be very helpful for others.
@PierUgit the below are examples of requirements and specifications papers. The gist is requirements are what problems the feature needs to solve or what capabilities it needs to have. Specifications are how the feature satisfies the requirements without necessarily getting into specific syntax.
You can always look at more examples from any of the papers at https://j3-fortran.org/doc/year, but those are the best examples I could think of recently.
Here are the first drafts of two papers, one with uses cases and the requirements, one with the specifications and the syntax. I have somehow some difficulties to fully separate the specifications from the requirements.
Those are looking reasonable, but make sure you stick to the 76 character line length limit.
Sure, I'll do in the final version.
By the way, in the minutes of the last meeting, what does mean:
** Motion 24-129 "allow complex pointer to reals and vice versa"
[Hugonnet](Cohen/Shafran)
This motion is to move the concept, not any of the proposed solutions
1. UC
\JoR
1. No work items
\Interp
When adding things to the worklist we don't necessarily want to imply that we are choosing a particular design for the feature.
I have finalized the files, and they can possibly be uploaded (after appropriate header edits)