Passing COMPLEX to REAL (or vice versa)

[Thomas Koenig]

I just investigated a few warnings issued by gfortran with
LAPACK, and stumbled across an interesting case.

gfortran -flto complained about (severely simplified and de-f77'd)

subroutine foo
complex (kind=8) :: a(5)
call bar(a,10)
end

and in another file

subroutine bar(a,n)
real(kind=8) :: a(n)
end

I thought this was legal, but apparently not in F18 (at least):

15.5.2.4 Ordinary dummy variables

2 The dummy argument shall be type compatible with the actual argument.

and 7.3.2.3

5 A nonpolymorphic entity is type compatible only with entities of the
same declared type.

Although probably a common idiom (or it would probably not have
made it into LAPACK), my guess would be it was never legal.

[JCampbell]

I would expect that this was never "legal" Fortran but is typical usage for a F77 wrapper. I expect in this case it was an easy way of accessing the real and imaginary components of complex(8) ( complex*16 ) :: A(5) by equivalencing to real(8) ( real*8 ) :: A(10)

Isn't this transformation typical of F77 LAPACK ?

[FortranFan]

On Saturday, April 20, 2019 at 8:19:20 AM UTC-4, Thomas Koenig wrote:

> ..

> Although probably a common idiom (or it would probably not have
> made it into LAPACK), my guess would be it was never legal.

For whatever it's worth, even ANSI X3.9-1966 document toward FORTRAN 66 had statements such as:

"The actual arguments, which constitute the argument list, must agree in order, number, and type with the corresponding dummy arguments in the defining program unit"

and that document makes clear "real, double precision, complex" are all *different types*.

[dpb]

On 4/20/2019 7:19 AM, Thomas Koenig wrote:
> I just investigated a few warnings issued by gfortran with
> LAPACK, and stumbled across an interesting case.
>
> gfortran -flto complained about (severely simplified and de-f77'd)
>
> subroutine foo
> complex (kind=8) :: a(5)
> call bar(a,10)
> end
>
> and in another file
>
> subroutine bar(a,n)
> real(kind=8) :: a(n)
> end
>
> I thought this was legal, but apparently not in F18 (at least):

...

> Although probably a common idiom (or it would probably not have
> made it into LAPACK), my guess would be it was never legal.

Indeed, it is not "legal" but rampant and widespread use for type
"spoofing" in order to retrieve the real/imaginary components directly.

It is required to work correctly by association rules for COMMON and
size matching between REAL and COMPLEX variables as far as memory is
concerned even though it isn't "pure" from a computer science
programming viewpoint and compliant to Standard Fortran.

It's one of those cases where convenience and history far overrides
"political correctness" in my view--I'd turn those warnings off if
could, selectively, if it were me or certainly ignore them rather than
worry about "fixing" LAPACK or any other similar libraries/user code.

There are far bigger fish to fry than to worry about some things...and
eliminating all warnings is one of those other things. :)

--

[robin....@gmail.com]

On Saturday, April 20, 2019 at 10:19:20 PM UTC+10, Thomas Koenig wrote:

> I just investigated a few warnings issued by gfortran with
> LAPACK, and stumbled across an interesting case.
>
> gfortran -flto complained about (severely simplified and de-f77'd)
>
> subroutine foo
> complex (kind=8) :: a(5)

de-F77'd ? Double precision complex was not around in F77,
except perhaps as an exension.

However, you might try
double precision a(*)
in BAR.

[FortranFan]

> ..

> subroutine foo
> complex (kind=8) :: a(5)
> call bar(a,10)
> end
>

> ..

> subroutine bar(a,n)
> real(kind=8) :: a(n)
> end
>
> I thought this was legal, but apparently not in F18 (at least):

> ..


If I'm not mistaken, the following which uses "RE" component of a COMPLEX data type is conformant starting with Fortran 2008, in case any modern Fortran coders are interested:

--- begin example ---
module m
implicit none
integer, parameter :: RK = selected_real_kind( p=12 )
contains
subroutine foo( c )
complex(kind=RK), intent(inout) :: c(:)
call bar( c%re )
end subroutine
subroutine bar( r )
real(kind=RK), intent(inout) :: r(:)
print *, "In bar: r = ", r
end subroutine
end module
use m
complex(kind=RK) :: a(2)
a = [ (1.0_rk,0.0_rk), (2.0_rk,0.0_rk) ]
call foo( a )
end
--- end example ---

Upon execution with gfortran:
In bar: r = 1.0000000000000000 2.0000000000000000


Though a temporary is likely with the above given the storage association rules in the standard for which a coder pays a price which can be contrasted with the earlier non-standard argument association usage as employed (presumably) in LAPACK, etc. where the coder has to carefully handle the dummy argument array of twice the size and make sure to consume the real and imaginary parts of said array appropriately.

[FortranFan]

On Saturday, April 20, 2019 at 8:19:20 AM UTC-4, Thomas Koenig wrote:

> ..

> subroutine foo
> complex (kind=8) :: a(5)
> call bar(a,10)
> end

> ..

> subroutine bar(a,n)
> real(kind=8) :: a(n)
> end
>

> .. probably a common idiom (or it would probably not have
> made it into LAPACK) ..


One would guess codes such as LAPACK, if they employed such idioms, either used ASSUMED SIZE array e.g., a(*) in the dummy argument declaration in 'bar' or automatic arrays of twice the size e.g., a(2*n).

And if a subprogram such as 'bar' only intended to work with REAL part of the complex data, then it must have had logic to explicitly skip over the imaginary parts in the array 'a'.

[Ron Shepard]

On 4/20/19 10:45 AM, FortranFan wrote:
> On Saturday, April 20, 2019 at 8:19:20 AM UTC-4, Thomas Koenig wrote:
>
>> ..
>> subroutine foo
>> complex (kind=8) :: a(5)
>> call bar(a,10)
>> end
>> ..
>> subroutine bar(a,n)
>> real(kind=8) :: a(n)
>> end
>>
>> .. probably a common idiom (or it would probably not have
>> made it into LAPACK) ..
>
>
> One would guess codes such as LAPACK, if they employed such idioms, either used ASSUMED SIZE array e.g., a(*) in the dummy argument declaration in 'bar' or automatic arrays of twice the size e.g., a(2*n).

No, they worked mostly as the poster described. The actual array was
declared as 5 and the dummy argument is declared as 10. It was up to the
programmer to put in the factor of 2.

They also worked the other way around, with a REAL actual argument and a
COMPLEX dummy argument. Again, it was up to the programmer to get the
factor of 1/2 right in the calling sequence.

>
> And if a subprogram such as 'bar' only intended to work with REAL part of the complex data, then it must have had logic to explicitly skip over the imaginary parts in the array 'a'.
>

The storage sequence rules required the odd indexes to be the real part
and the even indexes to be the complex part.

The reason for all of this type punning is that many compilers had
notoriously poor implementations of complex arithmetic, both with
protections against unnecessary overflow and underflow and with
efficiency. Some operations, such as division of a complex by a real
were especially inefficient. So it was often the case that the
programmer could take the real and imaginary parts of the complex number
and do the work himself much better than the compiler could do it.

Another reason was that it was not possible to reference just the real
or imaginary part of a complex number in f77, or prior to f2003, for
that matter. So if you wanted to do, for example, a simple assignment of
just the imaginary part, leaving the real part unchanged, that was not
possible. Instead, you had to assign the whole complex value,
overwriting the real part unnecessarily with its original value. This
was a design flaw in the language, it was not just a matter of compiler
efficiency. This was finally corrected in f2003, but of course that was
decades after millions of lines of code had been written in various
linear algebra and signal processing libraries that had been based on
the storage sequence trick.

And as others have pointed out, f77 only had single precision COMPLEX
type, it did not define the behavior of anything corresponding to double
precision or other extended precision types. Although double precision
complex, or the alternative complex*16, were common extensions, they
were not standard. So that was yet another reason for writing complex
arithmetic libraries in terms of real variables and arrays.

$.02 -Ron Shepard

[FortranFan]

On Saturday, April 20, 2019 at 12:36:20 PM UTC-4, Ron Shepard wrote:

> ..

> No, they worked mostly as the poster described. The actual array was
> declared as 5 and the dummy argument is declared as 10. It was up to the
> programmer to put in the factor of 2.
>
> They also worked the other way around, with a REAL actual argument and a
> COMPLEX dummy argument. Again, it was up to the programmer to get the
> factor of 1/2 right in the calling sequence.

> ..


Ok, my limited experience with old dialect FORTRAN involving COMPLEX and REAL data types is similar to that described at this webpage with Fast Fourier Transform subroutine from Numerical Recipes:
http://hosting.astro.cornell.edu/~akgun/Fortran/Numerique/fourier.f

where one notices "'data' is a complex array of length nn or, equivalently, a real array of length 2*nn":
subroutine four1(data,nn,isign)
..
double precision data(2*nn)
..
n = 2*nn
..
do i = 1,n,2
..

meaning the details of the array length are handled in the 'library' code rather in the caller. Note also the DO loop that increments by 2.

Anyways, one would think most use cases involving complex/real data are better handled now using generic interfaces and perhaps INCLUDE files for common code, though the standards body will be remiss if they fail to offer a good generics-based solution in the next revision, 202X.

[Ron Shepard]

On 4/20/19 12:43 PM, FortranFan wrote:
> where one notices "'data' is a complex array of length nn or, equivalently, a real array of length 2*nn":
> subroutine four1(data,nn,isign)
> ..
> double precision data(2*nn)
> ..
> n = 2*nn
> ..
> do i = 1,n,2
> ..
>
> meaning the details of the array length are handled in the 'library' code rather in the caller. Note also the DO loop that increments by 2.

In this case, one might have in the calling program

double precision :: data(n) ! n is even
...
call four1( data, n/2, isign )

So no matter which convention is adopted, the programmer must account
correctly for the factor of two. This was a perpetual point of confusion
and source of error when mixing code using different libraries or from
different programmers.

$.02 -Ron Shepard

[ga...@u.washington.edu]

On Saturday, April 20, 2019 at 9:36:20 AM UTC-7, Ron Shepard wrote:
> On 4/20/19 10:45 AM, FortranFan wrote:

(snip on mixing REAL and COMPLEX actual/dummy arguments)

> > One would guess codes such as LAPACK, if they employed such idioms,
> either used ASSUMED SIZE array e.g., a(*) in the dummy argument
> declaration in 'bar' or automatic arrays of twice the size e.g., a(2*n).

> No, they worked mostly as the poster described. The actual array was
> declared as 5 and the dummy argument is declared as 10. It was up to the
> programmer to put in the factor of 2.

(snip)

> The reason for all of this type punning is that many compilers had
> notoriously poor implementations of complex arithmetic, both with
> protections against unnecessary overflow and underflow and with
> efficiency. Some operations, such as division of a complex by a real
> were especially inefficient. So it was often the case that the
> programmer could take the real and imaginary parts of the complex number
> and do the work himself much better than the compiler could do it.

This is especially true for FFT routines, which tend to process
the real and imaginary parts separately.

> Another reason was that it was not possible to reference just the real
> or imaginary part of a complex number in f77, or prior to f2003, for
> that matter. So if you wanted to do, for example, a simple assignment of
> just the imaginary part, leaving the real part unchanged, that was not
> possible. Instead, you had to assign the whole complex value,
> overwriting the real part unnecessarily with its original value. This
> was a design flaw in the language, it was not just a matter of compiler
> efficiency. This was finally corrected in f2003, but of course that was
> decades after millions of lines of code had been written in various
> linear algebra and signal processing libraries that had been based on
> the storage sequence trick.

Reminds me that PL/I has an interesting way to separately
reference the real and imaginary parts:

DCL I FIXED BIN(31) COMPLEX;
I=0;
DO IMAG(I) = 1 TO 100;
PRINT SKIP LIST(I, SQRT(I));
END;

Note that in this case IMAG is both an intrinsic
function and pseudo variable. You can also extract both
at the same time with:

DCL ((X,Y) REAL, Z COMPLEX) FLOAT BIN(53);
Z = SQRT(2I);
COMPLEX(X,Y) = Z;

Fortran now has many features from PL/I, but not yet this one.

As for the original question, you can EQUIVALENCE to the
appropriate type, or TRANSFER it. I have no idea if the TRANSFER
can be done in an actual argument without a temporary copy.

[robin....@gmail.com]

The intrinsic functions AIMAG and REAL gave access to the
imaginary and real parts of a complex number -- since FORTRAN IV
days (1960s) at least.

[James Van Buskirk]

"FortranFan" wrote in message
news:5335209f-39f1-42c2...@googlegroups.com...

> Ok, my limited experience with old dialect FORTRAN involving
> COMPLEX and REAL data types is similar to that described at
> this webpage with Fast Fourier Transform subroutine from
> Numerical Recipes:
> http://hosting.astro.cornell.edu/~akgun/Fortran/Numerique/fourier.f

> where one notices "'data' is a complex array of length nn or,
> equivalently, a real array of length 2*nn":
> subroutine four1(data,nn,isign)
> ..
> double precision data(2*nn)
> ..
> n = 2*nn
> ..
> do i = 1,n,2
> ..

Where a real-half complex in-place DFT is performed there is
necessarily some type punning required.

[ga...@u.washington.edu]

On Sunday, April 21, 2019 at 9:01:45 PM UTC-7, James Van Buskirk wrote:
> "FortranFan" wrote in message
> news:5335209f-39f1-42c2...@googlegroups.com...

> > Ok, my limited experience with old dialect FORTRAN involving
> > COMPLEX and REAL data types is similar to that described at
> > this webpage with Fast Fourier Transform subroutine from
> > Numerical Recipes:
> > http://hosting.astro.cornell.edu/~akgun/Fortran/Numerique/fourier.f

(snip)

> Where a real-half complex in-place DFT is performed there is
> necessarily some type punning required.

I suppose the routine could accept a COMPLEX array of
half the length, with two real values in each element.

That complicates, and makes ugly, the calling routine, though.

But FFT routines themselves are much easier to write using
real arrays.