[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