[J3] [EXTERNAL] [BULK] SCALE(x,i) with complex x
Steven G. Kargl
kargl at uw.edu
Mon Nov 6 20:09:53 UTC 2023
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>
> 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
>
--
Steve
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/x-pkcs7-signature
Size: 3950 bytes
Desc: not available
URL: <https://mailman.j3-fortran.org/pipermail/j3/attachments/20231106/2f0ce5dc/attachment.bin>
More information about the J3
mailing list