[J3] [EXTERNAL] [BULK] SCALE(x,i) with complex x

[Clune, Thomas L. (GSFC-6101)]

Hi Steve,


The following program returns  (4.0000000,Inf) with NAG and Intel.    GFortran returns (NaN,Inf) which suggests it does not optimize away the 0 multiplies.

program main
   use IEEE_ARITHMETIC
   implicit none

   complex :: c

   c = cmplx(1., IEEE_VALUE(1., IEEE_POSITIVE_INF))
   print*, 2.**2 * c

end program main


The following program returns (4.0000000,-5.6051939E-45) with NAG,  (4.000000,0.0000000E+00) with ifort, and (4.00000000,-2.350988702E-38) with gfortran.

program main
   use IEEE_ARITHMETIC
   implicit none

   complex :: c

   c = cmplx(1., IEEE_VALUE(1., IEEE_NEGATIVE_SUBNORMAL))
   print*, c, 2.**2 * c

end program main

I’ll continue to steer clear of any applications that rely critically on such subtle numerical features.  ☺


  *   Tom

From: "Steven G. Kargl" <kargl at uw.edu>
Reply-To: <kargl at uw.edu>
Date: Monday, November 6, 2023 at 1:19 PM
To: j3 <j3 at mailman.j3-fortran.org>
Cc: "Clune, Thomas L. (GSFC-6101)" <thomas.l.clune at nasa.gov>
Subject: Re: [J3] [EXTERNAL] [BULK] SCALE(x,i) with complex x

On Mon, Nov 06, 2023 at 12:55:20PM +0000, Clune, Thomas L. (GSFC-6101) via J3 wrote:
Hi Steve,

Extending to complex would not seem to pose any challenge.    As
to why (well before my time), my guess would be that someone was
attempting to emulate a specific bit of functionality in another
layer (IEEE?) that did not address the complex case.

Of course back in the day, SCALE() could provide some significant
performance advantage via bit-twiddling vs full multiply.     These
days, though, unless an arg and result are known to be very small
and cache resident, the memory access will dominate.   Given that
I think a variant of your 3rd version is actually preferable in
most implementations:

PROGRAM FOO
    COMPLEX A, B
    A = (1,1)
    B = (2**2) * A   ! parens added for clarity for reader
    PRINT *, B
END PROGRAM FOO

I find this implementation much more obvious than SCALE() as one
needs to 1st remember which argument is which.   Of course, if
I used SCALE() frequently, I might feel differently.

Not opposed to the extension – just explaining why I there may
be mixed views on priority.

Thanks for the thoughtful reply.   Your guess aligns with what I
also suspect.  Fortran 90 was a massive update to Fortran 77,
and introduced the idea of a model number for the REAL type.  SCALE
was likely motivated by the definition of a model number.

While the code you show is indeed clear to a programmer, it
involves the evaluation of (2**2) (1 multiply), a conversion
of that result to COMPLX(4.,0.), and then the complex multiplication
CMPLX(4.*A%RE-0.*A%IM, 4.*A%IM+0.*A%RE) (4 multiplies, 2 additions).

A Fortran processor might optimize out the multiplications with
0, but this assumes that the result of that multiplication is
not exceptional.  F2018, 10.1.5.2.4, gives enough wiggle room
to allow this.  An example to consider (2.**2)*(1,Inf) = (NaN,Inf)
while SCALE((1,inf),2) would be (4., Inf).

A more important case might be a point in the 4th quadrant of the
complex (1,sign(0.,-1.)) = (1, -0.) i.e., a point infinitesimally
close to the real axis.  (2**2)*(1., -0.) yields (4., 0.) a point
that is now in the first quadrant.  SCALE((1.,-0.), 2) = (4.,-0.).

--
Steve

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.j3-fortran.org/pipermail/j3/attachments/20231106/c9b14956/attachment-0001.htm>