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

Mon Nov 6 20:19:20 UTC 2023

Steve,

Not sure which case you were asking about, so I did both:

   print*, c, cmplx(2.**2,0.)*c

NAG:  (1.0000000,Inf) (4.0000000,Inf)
Ifort:  (1.000000,Infinity) (4.000000,Infinity)

   print*, c, cmplx(2.**2,0.) * c

NAG: (1.0000000,-1.4012985E-45) (4.0000000,-5.6051939E-45)
ifort: (1.000000,-1.1754942E-38) (4.000000,0.0000000E+00)

From: "Steven G. Kargl" <kargl at uw.edu>
Reply-To: <kargl at uw.edu>
Date: Monday, November 6, 2023 at 3:09 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

Tom,

With your first example, one can get the NAG and Intel result
with gfortran if one uses the -Ofast or -ffast-math option.
This essentially allows gfortran to use a distributive property
2.**2*C -> (4.*C%RE, 4.*C%IM)

What do NAG and Intel do if you change to (2.**2, 0.)*C?
I suspect it wont' change anything, as both compiler likely
detect the 0.

--
steve

On Mon, Nov 06, 2023 at 07:50:51PM +0000, Clune, Thomas L. (GSFC-6101) via J3 wrote:
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<mailto:kargl at uw.edu>>
Reply-To: <kargl at uw.edu<mailto:kargl at uw.edu>>
Date: Monday, November 6, 2023 at 1:19 PM
To: j3 <j3 at mailman.j3-fortran.org<mailto:j3 at mailman.j3-fortran.org>>
Cc: "Clune, Thomas L. (GSFC-6101)" <thomas.l.clune at nasa.gov<mailto: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

--
Steve

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