What’s the most modern way to extract the real/imaginary parts of a complex number in Fortran? I’m aware of real and aimag, and %re/%im - what’s the most preferred method?
2 Likes
Indeed, with compilers conformant with Fortran 2008 and later revisions, %RE and %IM is the way to go in my “book” too. Note, however, it is processor-dependent whether one gets direct reference to the data, or the address of a copy of it when it comes to the more common scenarios with complex numbers and which is arrays of them.
module m
use, intrinsic :: iso_c_binding, only : c_loc, c_size_t
character(len=*), parameter :: fmtz = "(g0,1x,z0)"
integer(c_size_t), parameter :: mold = 0
contains
impure elemental function cmag( x, y ) result(r)
real, intent(in), target :: x, y
real :: r
print fmtz, "mag: address of x (hex): ", transfer( c_loc(x), mold=mold )
r = sqrt(x**2 + y**2)
end function
end
use m
complex, target :: z(3)
real :: q( size(z) )
z%re = [ 1.0, 2.0, 3.0 ]
z%im = 0.0
print fmtz, "main: address of z(1)%re (hex): ", transfer( c_loc(z(1)%re), mold=mold )
print fmtz, "main: address of z(2)%re (hex): ", transfer( c_loc(z(2)%re), mold=mold )
print fmtz, "main: address of z(3)%re (hex): ", transfer( c_loc(z(3)%re), mold=mold )
q = cmag( z%re, z%im )
print *, q
end
C:\Temp>ifort /standard-semantics i.f90 -o i.exe
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.2.0 Build 20210228_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.Microsoft (R) Incremental Linker Version 14.28.29337.0
Copyright (C) Microsoft Corporation. All rights reserved.-out:i.exe
-subsystem:console
i.objC:\Temp>i.exe
main: address of z(1)%re (hex): 72478FFAE8
main: address of z(2)%re (hex): 72478FFAEC
main: address of z(3)%re (hex): 72478FFAF0
mag: address of x (hex): 72478FFAF4
mag: address of x (hex): 72478FFAF4
mag: address of x (hex): 72478FFAF4
1.000000 2.000000 3.000000C:\Temp>gfortran i.f90 -o gcc-i.exe
C:\Temp>gcc-i.exe
main: address of z(1)%re (hex): 87FDB0
main: address of z(2)%re (hex): 87FDB8
main: address of z(3)%re (hex): 87FDC0
mag: address of x (hex): 87FDB0
mag: address of x (hex): 87FDB8
mag: address of x (hex): 87FDC0
1.00000000 2.00000000 3.00000000C:\Temp>
3 Likes
Thanks to you both. I noticed that you can’t do the following, though: (a+b)%re which is unfortunate. At least, gfortran calls that an unclassifiable statement.
1 Like
Fortran allows the somewhat more verbose ASSOCIATE block with type-inference:
associate ( t => a + b )
! consume t%re as read-only here
end associate
4 Likes
If you want to convert the real or imaginary part of the complex data entity to a different KIND, I would use REAL or AIMAG. Otherwise, %re and %im are the new, shiny approach.
Equally unfortunate, perhaps, is that one cannot expect
("abc" // "def")(3:4)
to yield “cd”. Users moving to a compiled language such as Fortran after extensive usage of interpreted languages will probably encounter and comment on such restrictions in the beginning.
2 Likes
Or:
(A+B+sin(C))(5:10)
where A, B and C are arrays. I will admit that it is not super readable, compared to Python’s (A+B+sin(C))[5:10], was readability the reason it is not allowed?
Indeed. if the function is not elemental, then it would have to compute 1 million elements only to throw away most of them to get (5:10) out. This is exactly how Python would do it also. So it can be wasteful, and that’s probably the reason it is not allowed.
I think there is a bug in gfortran when using %re and %im of complex type components. Consider the following minimal program:
program main
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
type complex_wrap
complex(dp), dimension(:), allocatable :: z
end type
type(complex_wrap) :: w
allocate(w%z(2))
w%z(1) = (1, 2)
w%z(2) = (3, 4)
print*,w%z%re
call print_arr(w%z%re)
contains
subroutine print_arr(x)
real(dp), intent(in), dimension(:) :: x
print*,x
end subroutine
end program main
The print statement in the main program prints “1 , 3” whereas when printing from the subroutine the result is “1, 2”. ifort compiles the program without warnings and prints “1, 3” in both cases.
I’ve submitted a bug report:
Bug 102891
4 Likes
I can confirm the bug with GFortran 11.0.1 20210403 (experimental) on Apple M1 installed using Conda.
The %Im must have the same attributes as a linear mathematical operator.
That means that
a%Im + b%Im
MUST be the same as
(a + b)%Im.
This is not “unfortunate” it is a fundamental mistake.
That reminds me. I have had issues with several compilers, including passing x%im as an argument, using x%im=expression where x%im is on the LHS, … which have all been resolved with the compilers I have had access to except this one…
program testit
implicit none
complex :: a = (1.0, 3.0), b = (2.0, 4.0), x
real :: im
!im(x)=imag(x)
im(x) = x%im
write (*, *) im(a + b)
end program testit
That is, using x%im in a statement function. Does anyone want to weigh in on whether that is standard or not? ifx gets an ICE, and gfortran says it is recursive, which seems incorrect. Definitely a dusty corner. Curious what other compilers report. An error message (if it is an error) is definitely a preferred response instead of an ICE so that needs reported either way; but I do not see why it should be considered recursive so I feel the response from gfortran is misleading at best; but not quite convinced it is standard-conforming anyway.
The Intel forum login has been so unreliable and cumbersome and requires SMA just to report a bug, so I am done with that forum. If anyone wants to report an ICE to Intel here is a minimal reproducer.
program testit
implicit none
complex :: x
real :: im
im(x) = x%im
write (*, *) im((1.0,3.0) + (2.0,4.0))
end program testit
intel_complex_bug.f90(5): error #5623: **Internal compiler error: internal abort** Please report this error along with the circumstances in which it occurred in a Software Problem Report. Note: File and line given may not be explicit cause of this error.
im(x) = x%im
-----------^
compilation aborted for intel_complex_bug.f90 (code 3)
nagfor accepts this just fine and outputs the correct answer. I think it’s valid code.
1 Like
Note that the bug does not happen because of a name collision: the compiler fails whatever the name given to the statement function.
1 Like
Seeing this post I thought renaming the function would likely eliminate the gfortran error where it fails saying it is recursive but it did not:
program testit
implicit none
complex :: x
real :: img
img(x) = x%im
write (*, *) img((1.0,3.0) + (2.0,4.0))
end program testit
gfortran main.f90
main.f90:5:9:
5 | img(x) = x%im
| 1
Error: Statement function at (1) is recursive
So
- nagfor runs
- ifx/ifort ICE
- gfortran produces “recursive” error
so far. Nice to know nagfor works; as two compiler failures made me question if the case was standard although I could not think of a reason
it was not.
Went to report on GNU bugzilla site and found someone already did the dirty work. Thanks!
3 Likes
I’m not sure either. There is this phrase (in 15.6.4 of the f2023 draft), which suggests it might not be:
“[…] the function shall not be a transformational intrinsic, and the result shall be scalar.”
This is referring to a function reference on the rhs of the statement. I’m unsure because I certainly have defined statement functions in the past that violated the transformational intrinsic, and I have never had problems with that. For example, before the modern version of the REAL() intrinsic with its KIND argument, I would routinely do something like this in my codes.
REAL*8 RX
RX(i) = (i) ! integer to working precision real conversion.
The transformational intrinsic was implicit in the RHS in this case, but it was there nonetheless.
In the examples in this discussion, the result is also an array, which seems to be not allowed either. This constraint is C1589, which I think means that the compiler must detect and report violations. So if the programmer wants an elemental function that works with both scalars and arrays, I think he cannot use a statement function, he must use an actual procedure or an internal procedure (with the elemental attribute).
I don’t think this is correct. This is a component reference, not an operator, thus it need not distribute the way an operator might be defined to distribute. In the example,
I would say that the first expression only operates on the imaginary parts of the two arrays, to give the real vector sum. But the second expression involves adding both the real and imaginary parts of the arrays, and then throwing away the real parts and extracting the imaginary parts to give the real vector result. Of course, an optimizing compiler might recognize that the real parts are ignored afterwards and eliminate those operations, but the two statements are not equivalent as far as the language semantics. The other difference between %im and AIMAG() is that the former is modifiable while the latter is not. For example, the former can appear on the LHS of an expression, while the function reference cannot, and the former can be used as an actual argument that is modified by a subprogram, but the latter cannot.
1 Like
The question is then: “is %im a transformational intrinsic”? It is not IMO. Regarding the scalar constraint, the MRE above is just dealing with scalars.
1 Like