Well for c_float and c_double (and their complex equivalents, c_float_complex and c_float_double) this is not an issue since, section 6.2.5, paragraph 17, page 39, of the draft C23 standard N3096, states:
Each complex type has the same representation and alignment requirements as an array type containing exactly two elements of the corresponding real type; the first element is equal to the real part, and the second element to the imaginary part, of the complex number
The loophole in my eyes, is that we need to piggy-back on C functions and language semantics to do this. The reason it works then is because all Fortran compilers are hosted on top of C-favouring platforms, where the default integer and real kind, âhappenâ to be interoperable with C types.
This plays back exactly to what Wirth says with âstatic typing that cannot be breached.â
use, intrinsic :: iso_c_binding
integer, parameter :: wp = c_double
real(wp), target :: a(6)
complex(c_double), pointer :: b(:)
a = [1,2,3,4,5,6]
call c_f_pointer(c_loc(a),b,[3])
print *, b
end program
If the programmer later decided single precision was sufficient, the above code would not report the type mismatch. But with a true Fortran pointer, it would. Hence c_loc is a loophole because it subverts the type system. (Here the mistake is obvious, but imagine that wp, a and b are declared in different files.)
The Fortran pointer case (not complex to real however):
integer, parameter :: wp = c_float
real(wp), target :: a(6)
real(c_double), pointer :: c(:)
c => a
$ gfortran -Wall test_complex.f90
test_complex.f90:9:0:
9 | c => a
|
Error: Different types in pointer assignment at (1); attempted assignment of REAL(4) to REAL(8)
integer, parameter :: wp = c_double
real(wp), target :: a(6)
complex(wp), pointer :: b(:)
a = [1,2,3,4,5,6]
call c_f_pointer(c_loc(a),b,[3])
print *, b
end program
Then, it is safe against any changes in the kind parameter, and the only implicit assumption is really only the storage equivalence of 2 reals and one complex.
And now try calling say a routine that performs a moving average on the [1., 2., 3. 4.]âŚ
The point is not about accessing the components of complex array, but about viewing the complex array as if it was a real array.
Do you mean that c_f_pointer(c_loc(x),y,shape) is legal whatever the (non-polymorphic) types of x and y? I had once an extensive discussion on Stack Overflow with someone saying the opposite, and he somehow convinced me (and his answers generally demonstrate an accurate knowledge of the standard⌠I do not doubt you have an accurate knowledge either, but Iâm troubled by the disagreement on this point).
There are several things going on in these discussions that are intertwined. In the above call, the complex output matrix C(n,n) is computed fully: the real and imaginary parts of every matrix element. The potential inefficiencies arise when the programmer knows in advance that A(:,:), B(:,:) or C(:,:) has special properties. It might be known that C(:,:) is real (in which case the imaginary parts would not need to be computed), or purely imaginary (the real parts would not need to be computed), or symmetric (only half of the matrix elements need to be computed), or hermitian (only half of the off-diagonal matrix elements need to be computed, and only the real parts of the diagonal elements are needed). There can also be combinations of those, such as purely real and symmetric, etc. In addition to the number of floating point operations involved in these various combinations, there is also the storage issues. Does the programmer really want to store a hermitian matrix in square form, or does he want to store it in a packed form (i.e. half the storage)? Maybe he wants to store it in packed form for most operations, but expand it into square form occasionally for specific operations?
It is true that for some combinations of matrix dimensions all of that can be ignored, and zgemm() can be used blindly with no loss of efficiency. But for other combinations, particularly when the dimensions are large, the optimal solution is usually to exploit the known properties of the input and output matrices. That is when it is useful to work with the %re and %im parts of the input and output arrays separately, and that is when the c_loc()/c_f_pointer() semantics become useful.
Almost. Hereâs what the standard says regarding C_LOC:
Argument. X shall have either the POINTER or TARGET attribute. It shall not be a coindexed object. It shall be a variable with interoperable type and kind type parameters, an assumed-type variable, or a nonpolymorphic variable that has no length type parameter. If it is allocatable, it shall be allocated. If it is a pointer, it shall be associated. If it is an array, it shall be contiguous and have nonzero size. It shall not be a zero-length string.
and then for C_F_POINTER:
CPTR shall be a scalar of type C_PTR. It is an INTENT (IN) argument. Its value shall be
⢠the C address of an interoperable data entity,
⢠the result of a reference to C_LOC with a noninteroperable argument, or
⢠the C address of a storage sequence that is not in use by any other Fortran entity.
The value of CPTR shall not be the C address of a Fortran variable that does not have the TARGET attribute.
If the value of CPTR is the result of a reference to C_LOC with a noninteroperable effective argument X, FPTR shall be a nonpolymorphic pointer with the same type and type parameters as X. In this case, X shall not have been deallocated or have become undefined due to execution of a RETURN or END statement since the reference. If X is scalar, FPTR becomes pointer associated with X; if FPTR is an array, SHAPE shall specify a size of 1. If X is an array and FPTR is scalar, FPTR becomes associated with the first element of X. If both X and FPTR are arrays, SHAPE shall specify a size that is less than or equal to the size of X, and FPTR becomes associated with the first PRODUCT (SHAPE) elements of X (this could be the entirety of X).
This does say that the type of FPTR needs to be the same as that of the argument to C_LOC, which is a catch I had not remembered, not that anyone is checking that.
Regardless, this does seem an unfortunate condition.
There it says that x and y shall have the same type, and my question was about x and y having different types. Your quote is about the âCase (ii)â in the text of the standard, but the âCase (i)â is much trickier to interpret.
Yeah, I was editing my reply when you posted. Case (i) isnât that hard to interpret - it says that the type of the C_LOC argument must be interoperable with the type of the FPTR. So, REAL and COMPLEX are not interoperable with each other.
How is using c_f_pointer to refer to a real array as complex better or different than using EQUIVALENCE to do so? The analog of the code with c_f_pointer
use iso_c_binding
implicit none
integer, parameter :: wp = c_float, nc=2
real(wp), target :: a(2*nc)
complex(wp), pointer :: b(:)
integer :: i
a = [(i,i=1,2*nc)]
call c_f_pointer(c_loc(a),b,[nc])
b = b*2
print*,b
print*,a
end program
is
use iso_c_binding
implicit none
integer, parameter :: wp = c_float, nc=2
real(wp) :: a(2*nc)
complex(wp) :: b(nc)
integer :: i
equivalence (a,b)
a = [(i,i=1,2*nc)]
b = b*2
print*,b
print*,a
end program
But thatâs the point. A complex array is not a real array. And the quotes shared by @ivanpribec above are correct.
To be a bit more specific about my earlier point
This would break backwards compatibility in at least one way
module do_stuff_m
private
public :: do_stuff
interface do_stuff
module procedure do_stuff_real
module procedure do_stuff_complex
end interface
contains
subroutine do_stuff_real(r_)
real, pointer, intent(in) :: r_(:)
print *, r_
end subroutine
subroutine do_stuff_complex(c_)
complex, pointer, intent(in) :: c_(:)
print *, c_
end subroutine
end module
use do_stuff_m
real, target :: r(10)
complex, target :: c(5)
r = 1.
c = (2., 3.)
call do_stuff(r)
call do_stuff(c)
end
If the proposal was added, which of the specific procedures should be called for either case? If youâre going to say that a real array can be viewed as a complex array and vice-versa, what was distinguishable now becomes ambiguous.
You seem to claim that this is something impossible in principle because it is ambiguous. The best way to alleviate your concerns would be to turn to some examples of programming languages where such idea was implemented, is commonly used and found to be of advantage. The first what comes to my mind is the Numpy package in Python, with its concepts of copies and views and the numpy.astype.
And the proposal is actually to be less strict on this point. Technically a complex is made of two reals, and many programmers do use the c_f_pointer trick.
Thatâs a very valid objection. But that has 2 possible solutions:
prioritize the disambiguation by selecting the procedure with the exact match if it exists (but Iâm afraid there would likely be ambiguous cases even with that)
do not authorize the real/complex mismatches in the procedure arguments, and only allow the pointer association between real and complex types (which was the initial point of the proposal, the argument association being just an option in the proposal).
EQUIVALENCE can only be used in fairly limited and specific situations, such as with common arrays and local variables. c_loc() and c_f_pointer() are more general in many ways, allowing for example module arrays, dummy argument arrays, local automatic arrays, allocatable arrays, and so on.
The common aspect in all these cases is the storage sequence association. The fortran standard goes to some effort to define the storage sequence association of real and complex arrays. For example, it explains how a complex variable may be defined by separately assigning its real and imaginary parts. This is in contrast, for example, to what happens when an integer and a real are equivalenced â in this case, when the integer is defined the real becomes undefined and vice versa.
I donât claim itâs impossible, just that it violates the abstraction provided by the language. Being able to say complex means you donât have to say â64 bits, where the first 32 bits are interpreted as the real part of a complex number in the floating point representation of a real variable, and where the last 32 bits are interpreted as the imaginary part of a complex number with the same floating point representationâ. The standard doesnât even say 64 bits, it only says, 2 storage units, the first of which is the real part, and the second is the imaginary part.
But other languages afford different abstractions. Python doesnât even have static typing, so its abstractions and what is useful to define in terms of them are very different. Pythonâs abstractions are everything is a dictionary and mutable, and data sharing is fine. Fortranâs abstractions are static typing and any data sharing is explicit.
My point was that itâs not necessarily a good idea. Just because âmany programmers do itâ doesnât mean its good design, and certainly doesnât mean its good language design.
Iâve heard it said this was a big mistake on the part of C++, that it often leads to lots of developer confusion.
Thatâs basically what my example demonstrated is already possible. You can just only do it one direction. I.e. the complex declaration defines the whole storage sequence, and you can view the real parts or the imaginary parts as real arrays.
Iâd like a bit of pragmatism here. Sometimes a practice is just the result of âoh, I donât want to bother with being more rigorousâ, but here itâs not just that: the alternative is duplicating the data with all the consequences on the performances, which includes not only the time needed to copy the data but also being able to load and process twice less data at once into the available RAM, meaning in turns larger relative overlaps between chunks of data.
I explained above why it does not solve the problem. Having separate pointers to the real and imaginary parts does not enable processing the data as unified real array.
How have we gone with standardising or even allowing âF77 wrappersâ which were so common before ALLOCATE was provided. They are still used today in many legacy codes, with appropriate compiler options, while assuming contiguous memory.
But what if one introduce a new language feature like:
complex :: wz(N)
select complex as real ( wr => wz)
! where here wr is real and has shape (2, N) not necessary contiguous in
! order to accomodate different memory layout
end select
@egio there is a misunderstanding of what the proposal is about. Itâs (mostly) not about accessing the real and imaginary parts of complex numbers: for this %re and %im already do the job. Itâs about representing the same data in different domains, using the same storage.
Typically you have a series of real numbers, for instance a time series in audio processing. You read it in a real :: sr(n) array because this is definitely what the series is at this point. You want to apply some processing steps on it in the time domain (e.g. a dynamic equalization), so you really need to have this real array, not a sr(2,*) array that means nothing here.
But you may also want to apply some processing steps in the Fourier domain (e.g. a spectral balance). You can have a separate complex array sc(n/2), perform a real-to-complex FFT to this array, do your stuff on the complex numbers, and perform the complex-to-real inverse FFT back to the sr array. Fine.
But there are some cases where the sr array can be really huge (and often multidimensional), and the last thing you want is to duplicate it into another array of the same size. in practice itâs very common to have in-place FFT routines, that output the Fourier transform in the original and real sr(:) array. But you know these are complex numbers in it now, and you want to perform complex arithmetic. Hence the proposed sc => sr pointer association, where sc is a complex pointer (with size(sc)==size(sr)/2)
Just for illustration, here is how the same thing can be achieved in Python/NumPy and C++:
# test_complex.py
import numpy as np
r = np.array([1,2,3,4],dtype=np.double)
c = r.view(dtype=np.cdouble)
print("r = ", r)
print("c = ", c)
print("do r and c alias? ", np.shares_memory(r,c))
// test_complex.cpp
//
// gcc -Wall -std=c++20 test_complex.cpp
//
#include <complex>
#include <iostream>
#include <span>
#include <array>
#include <type_traits>
#include <cassert>
template<typename Container,
std::enable_if_t<
std::is_floating_point<typename Container::value_type>::value,bool> = true>
auto as_complex_span(Container &a) {
using T = typename Container::value_type;
assert(a.size() % 2 == 0);
std::span<std::complex<T>> s{
reinterpret_cast<
std::complex<T> *>(a.data()),
a.size()/2
};
return s;
}
int main() {
std::array<float,4> a{1, 2, 3, 4};
// Retrieve complex view of the array
auto c = as_complex_span( a );
for (auto &z: c) {
std::cout << z << ' ';
}
std::cout << '\n';
return 0;
}
(For completeness, the C++ converter function template should also check the container fulfils the contiguous container requirements so it can be converted to a span. That goes beyond my C++ knowledge.)