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

Wed Nov 8 23:44:58 UTC 2023

Hi Steven,

The Fortran standard does not require that compilers calculate real*complex multiplication in the way gfortran does it. It does not say anything about the algorithm for complex multiplication at all.

And it is not the distributive law - there are no parentheses here. It's just the mathematical meaning of real*complex. Gfortran is the only compiler I've ever seen to insist on using the naïve complex multiply algorithm for real*complex. I think someone just misunderstood the standard here.

Similarly for complex*complex when one or more of the imaginary parts is known to be zero.

The standard does not even require complex multiply to be done using the naïve algorithm when the imaginary part is non-zero. It is really is completely silent! (The naïve algorithm is pretty decent when both parts are nonzero, so I'd expect it to be used then.)

I note that there are also multiple algorithms for complex division, and it is not uncommon for different compilers to use different algorithms there.

I note that there is no IEEE arithmetic standard for complex operations. One might speculate that if there were, it would say that real*complex would not calculate zero times the complex and add it, i.e. that INF*complex would not produce spurious NaNs.

> It's not that much slower.  But it can be more robust.

It can be a lot slower if subnormals are supported, as that means branching, and that means it's harder to simd-vectorise it in software. But mostly I agree that in the cases where one would be using SCALE, one is probably not in a vectorised inner loop anyway.

BTW, it's not that I think the suggestion has zero value, but it doesn't seem to me that it is likely to be widely useful.

And we can do it already. In the common case (scaling by a constant), multiply is easy to write.

In the (I expect) rare case where one is scaling by an unknown factor that might exceed the maximum/minimum exponent, one may write
    Cmplx(Scale(c%re,n),Scale(c%im,n),Kind(c))

Is that more verbose? Sure. But would it occur so often, and is it so intrusive, that Scale(c,n) should be added to the standard? (There are very few zero-cost features.) That is the question, and one that people may differ on.

I don't see a lot of use of even the SCALE that we have, in Fortran programs that cross my desk. Most of the time, people are scaling by less than maxexponent (scaling by 2**digits or 2**(emax/2) are what I see more commonly), and in those cases, multiplication is easy to understand. So, I find the case for extending it to COMPLEX not very convincing.

Finally, I think the gfortran guys should get the free performance boost by calculating the mathematical real*complex. There are no user parentheses to violate here, so it falls clearly within the "mathematical equivalence" rule.

Cheers,
-- 
..............Malcolm Cohen, NAG Oxford/Tokyo.

-----Original Message-----
From: J3 <j3-bounces at mailman.j3-fortran.org> On Behalf Of Steven G. Kargl via J3
Sent: Tuesday, November 7, 2023 3:33 PM
To: Malcolm Cohen <malcolm at nag-j.co.jp>
Cc: Steven G. Kargl <kargl at uw.edu>; 'General J3 interest list' <j3 at mailman.j3-fortran.org>
Subject: Re: [J3] [BULK] SCALE(x,i) with complex x

On Tue, Nov 07, 2023 at 12:39:56PM +0900, Malcolm Cohen wrote:
> >2.**n * (1., 1.) will raise FE_INEXACT.
> 
> Not so. Every intermediate calculation there is exact.
> 

Whoops.  I should have been more careful.

> call test(3)

  call test(-128) // nagfor results?

> end
> subroutine test(n)
> use ieee_arithmetic
> complex c,d
> logical inex
> c = (1,1)
> d = 2.0**n*c
> call ieee_get_flag(ieee_inexact,inex)
> print *,inex,d
> end
> 
> This prints "f (8.000,8.000)"  (with variable numbers of zeroes) on gfortran, intel fortran, nag fortran, ...
> 
> >2.**n * (1., inf) = (nan,inf) will raise FE_INVALID.
> 
> Highly unlikely unless you use a pessimising compiler. Wow, gfortran is
> a pessimising compiler here?! I am very surprised.

"Pessimising" might be a bit strong.  gfortran implements real*complex
multiplication according the Fortran standard and correctly.  Type
conversion gives 2.**n*(x,y) -> (2.**n, 0.) * (x,y).  Short-circuiting
2.**n*(x,y) to (2.**n*x, 2.**n*y) by other compilers seems to suggest
those compilers are appealing to 10.1.5.2.4 where one is dismissing
the possibility that either x or y is exceptional.  Or, those 
compilers are using a distributive property of ordinary math.  Seems
to be QoI.

BTW, you get the NAG and Intel results with gfortran if one
uses the -ffast-math option.

> SCALE is slower than multiplication (unless of course, the compiler
> can see all the values and turn it into a multiply, in which case
> it's merely more opaque). Hardware multiplication is highly optimised.
> We are no longer in the days when cpus sometimes did not have hardware
> multiplication.

It's not that much slower.  But it can be more robust.
What does NAG and Intel do with 2**n*(x,y) if n > maxexponent(x)
and x and y are sufficiently less than 1.?

2.**(maxexponent(t0) + 1) * (0.125, 0.125)

I suppose forcing a user to do 2.**maxexponent(t0) * (0.125, 0.125) * 2
is acceptable.

> 
> Cheers,
> -- 
> ..............Malcolm Cohen, NAG Oxford/Tokyo.
> 
> -----Original Message-----
> From: J3 <j3-bounces at mailman.j3-fortran.org> On Behalf Of Steven G. Kargl via J3
> Sent: Tuesday, November 7, 2023 4:44 AM
> To: Long, Bill F <william.long at hpe.com>
> Cc: Steven G. Kargl <kargl at uw.edu>; General J3 interest list <j3 at mailman.j3-fortran.org>
> Subject: Re: [J3] [BULK] SCALE(x,i) with complex x
> 
> 2.**n * (1., 1.) will raise FE_INEXACT.
> 2.**n * (1., inf) = (nan,inf) will raise FE_INVALID.
> 
> SCALE((1.,1.), 2) is exact, .i.e, no exception raised.
> SCALE((1.,inf), 2) = (4., inf) will raise FE_OVERFLOW.
> 
> SCALE((0.125, 0.125), maxexponent(t0) + 1) = (8.507059173E+37,8.507059173E+37)
> 
> 2.**(maxexponent(t0) + 1) * (0.125, 0.125) does not compile.
> 
> My actual use case is computation of bessel functions
> of a complex argument via downward recursion.  Rescaling
> the recursion sequence allows one to extend the length
> of the sequence while avoiding overflow.
> 
> And, yes, I realize that this would be for F202Y.
> 
> -- 
> steve
> 
> On Mon, Nov 06, 2023 at 04:21:11PM +0000, Long, Bill F wrote:
> > In F90 there was a category of intrinsics "Floating point manipulation functions" that included EXPONENT, FRACTION, NEAREST, RRSPACING, SCALE, SET_EXPONENT, and SPAING.
> > For most of these only REAL arguments make sense. For SCALE, there is also a sensible result definition for a COMPLEX argument.  But given the context, I suspect that was not
> > in the scope of the conversation.  This could be added to the requirements list for F202Y.  (It is too late in the process to add this to F2023.)  But the question for F202Y would be "Why bother?".
> > 
> > Cheers,
> > Bill
> > 
> > ________________________________
> > From: J3 <j3-bounces at mailman.j3-fortran.org> on behalf of Clune, Thomas L. (GSFC-6101) via J3 <j3 at mailman.j3-fortran.org>
> > Sent: Monday, November 6, 2023 6:55 AM
> > To: kargl at uw.edu <kargl at uw.edu>; General J3 interest list <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
> > 
> > 
> > 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.
> > 
> > 
> > 
> > Cheers,
> > 
> > 
> > 
> >   *   Tom
> > 
> > 
> > 
> > From: J3 <j3-bounces at mailman.j3-fortran.org> on behalf of j3 <j3 at mailman.j3-fortran.org>
> > Reply-To: "kargl at uw.edu" <kargl at uw.edu>, j3 <j3 at mailman.j3-fortran.org>
> > Date: Saturday, November 4, 2023 at 11:29 PM
> > To: j3 <j3 at mailman.j3-fortran.org>
> > Cc: "Steven G. Kargl" <kargl at uw.edu>
> > Subject: [EXTERNAL] [BULK] [J3] SCALE(x,i) with complex x
> > 
> > 
> > 
> > CAUTION: This email originated from outside of NASA.  Please take care when clicking links or opening attachments.  Use the "Report Message" button to report suspicious messages to the NASA SOC.
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> 
> -- 
> Steve
> 

-- 
Steve