In connection with my effort to implement Kahan summation in a manner
that is robust against the treatments of optimizing compilers I just
posted a request for interpretation to the J3 committee web server.
The request is attached FYI.  --Bas Braams

To: J3                                           08-208
From: Bas Braams
Subject: Interpretation, precise f.p. semantics of the REAL intrinsic
Date: 2008 June 15

REQUEST FOR INTERPRETATION:

Must the intrinsic function REAL with KIND parameter wp return a value
that is a REAL (KIND=wp) floating point number?

RATIONALE FOR THE QUESTION:

Computer hardware may use a wider floating-point format for registers
than for memory; e.g., 80 bits for registers and 64 bits for memory
for the case of standard double precision floating point numbers.
Some algorithms require a high level of control over floating point
semantics.  If the intrinsic function REAL with KIND parameter wp is
guaranteed to return a REAL (KIND=wp) result then a programmer can use
this to force intermediate results into main memory format, never mind
that the optimizing compiler may have placed the intermediate in a
register.

I am interested in a J3 interpretation of this matter, especially a
loud and clear affirmative interpretation, because it appears that
some present Fortran compilers optimize away my explicit use of the
REAL intrinsic.  The context is code for compensated summation (Kahan
summation).  I appreciate that parentheses are inviolable courtesy of
the Fortran standard, but in order to have code that cannot be broken
by an optimizing compiler I seem to need also a language mechanism to
force intermediate results into main memory format.

Bas Braams
Chemistry Department and
Emerson Center for Scientific Computation
Emory University
Atlanta, GA

In article <f910aaa6-6099-4edd-8954-9c56aefd0482@w7g2000hsa.googlegroups.com>,
	"braams@mathcs.emory.edu" <braams@mathcs.emory.edu> writes:
> In connection with my effort to implement Kahan summation in a manner
> that is robust against the treatments of optimizing compilers I just
> posted a request for interpretation to the J3 committee web server.
> The request is attached FYI.  --Bas Braams
> 
> To: J3                                           08-208
> From: Bas Braams
> Subject: Interpretation, precise f.p. semantics of the REAL intrinsic
> Date: 2008 June 15
> 
> REQUEST FOR INTERPRETATION:
> 
> Must the intrinsic function REAL with KIND parameter wp return a value
> that is a REAL (KIND=wp) floating point number?

I think the Fortran 95 standard already answers this question:


13.14 Specifications of the intrinsic procedures

This section contains detailed specifications of the generic intrinsic
procedures in alphabetical order.

The types and type parameters of intrinsic procedure arguments and
function results are determined by these specifications.  A program
is prohibited from invoking an intrinsic procedure under circumstances
where a value to be returned in a subroutine argument or function result
is outside the range of values representable by objects of the specified
type and type parameters.

AFAIK, the standard does not address machine level details such as
the width of FP rregisters.

-- 
Steve
http://troutmask.apl.washington.edu/~kargl/

braams@mathcs.emory.edu <braams@mathcs.emory.edu> wrote:

> In connection with my effort to implement Kahan summation in a manner
> that is robust against the treatments of optimizing compilers I just
> posted a request for interpretation to the J3 committee web server.
....
> Must the intrinsic function REAL with KIND parameter wp return a value
> that is a REAL (KIND=wp) floating point number?

Did you check the existing interprations first? This sounds like almost
a word-for-word copy of one that I recall Van Snyder submitting a few
years ago.

I don't think I'll jump into the middle of that debate. (Added later -
Looks like I did after all, but anyway, I won't argue about it). I seem
to recall the one on Van's interp being fairly heated. But I'll note a
few things.

1. Part of what strikes me as so simillar to Van's request is that you
haven't really phrased the question in a way that a direct answer will
be very useful. Of course it returns a real(kind=wp) floating point
number. It is hard to imagine how the standard could be much more
explicit about that than it already is; and it is hard to imagine how an
interp answer could be any more explicit either. I know that doesn't
answer your real question, which is much harder, but that is the
trivially obvious answer to the question literally as you asked it. Your
rationale is a bit more to the point, but it is "dangerous" to ask
questions that have such trivially obvious, but unhelpful answers. The
"danger" is that the committee will avoid the trickier issue by just
answering exactly the question that was asked.

I've seen things like this before, often where the question was probably
intended to be a leading one.

2.Your real question is much more subtle and more controversial. I don't
actually recall the answer to Van's similar question. But issues to keep
in mind are

  a. The Fortran standard doesn't guarantee that much of anything having
to do with floating point arithmetic give exact results. Of course,
that's partly because some things don't have exact results, but even
when they do, the Fortran standard doesn't guarantee that you will get
them. In my opinion, it skips some very major steps to ask about fine
points like this one when the Fortran standard allows 2.0+2.0 to return
the result 5.0. IEEE makes things more "interesting", but you didn't
mention IEEE kinds specifically.

  b. There is the whole "mathematical equivalence" thing. That's what
you are alluding to when you mention parenthesis, as the parenthesis
rule is an exception to the allowance of allowing mathematical
equivalence. Unfortunately, the difference between single and double
precision is, in the terms used by the standard, a computational
difference rather than a mathematical one.

 c. Finally, there is the matter of confusing specification and
implementation. I think that is the place where your question really
fails to get at what you actually mean. That the REAL intrinsic returns
a value of a given kind is part of its specification. But that does't
translate into "the computer must store the bits in that form." It
sounds to me as though you are assuming that translation. I think that's
why the obvious answer "yes, it returns a value of that kind" doesn't
mean "yes, the bits to represent such a value have to be computed".

I might be (probably am) showing my own prejudices in the above points.
As I said, I forget exactly how the previous interp got answered. But
I'm a little afraid that the answer might have been less than the
perfectly clear definitive ruling you might be hoping for. I sure recall
that some of the draft answers weren't; I just don't remember whether
that finally got cleared up in the final answer. I'm half thinking that
I recall a followon interp by people who didn't like the answer to the
first one, trying to rephrase the question to get the "right" answer.
But I could well be confusing floor conversations with formal interps. I
suppose I ought to go look up the interp(s), but I'm feeling on the lazy
side at the moment. Perhaps someone else will. And perhaps it does give
a perfectly clear answer.

Unfortunately, it is sometimes difficult get perfectly clear answers out
of a committee when there is strong irreconcilable disagreement. The
perfectly clear answers tend to cause objections from one side or the
other. The answers that aren't so clear can sometimes slip through
better, particularly if the committee is "tired" of a debate that seems
to be going nowhere. I've seen it happen.

In this case, there is a strong conflict between the ability to optimize
code well and the ability for the user to state exactly what the
implementation must do, down to the last bit. Both things are strongly
desired. There have been several proposals in this area, but I don't
recall that they have done very well. I doubt you are going to be able
to adequately satisfy both needs with just an interpretation of the
existing standard and syntax. People have sure tried, but I don't recall
much luck. Some of the proposals have tried things like adding new
syntax to indicate expressions that should be evaluated exactly as is
without otherwise allowed optimizations.

-- 
Richard Maine                    | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle           |  -- Mark Twain

braams@mathcs.emory.edu wrote:
(snip)

> RATIONALE FOR THE QUESTION:

> Computer hardware may use a wider floating-point format for registers
> than for memory; e.g., 80 bits for registers and 64 bits for memory
> for the case of standard double precision floating point numbers.
> Some algorithms require a high level of control over floating point
> semantics.  If the intrinsic function REAL with KIND parameter wp is
> guaranteed to return a REAL (KIND=wp) result then a programmer can use
> this to force intermediate results into main memory format, never mind
> that the optimizing compiler may have placed the intermediate in a
> register.

For x87, I believe the only way to force memory format is to
actually write to memory.  (Hopefully to cache.)  That is slow
enough that compilers try to avoid doing it.

> I am interested in a J3 interpretation of this matter, especially a
> loud and clear affirmative interpretation, because it appears that
> some present Fortran compilers optimize away my explicit use of the
> REAL intrinsic.  The context is code for compensated summation (Kahan
> summation).  I appreciate that parentheses are inviolable courtesy of
> the Fortran standard, but in order to have code that cannot be broken
> by an optimizing compiler I seem to need also a language mechanism to
> force intermediate results into main memory format.

I wonder if storing into a VOLATILE or ASYNCHRONOUS variable
would be less likely to be optimized away.

-- glen

glen herrmannsfeldt wrote:
> braams@mathcs.emory.edu wrote:
> (snip)
>
>> RATIONALE FOR THE QUESTION:
>
>> Computer hardware may use a wider floating-point format for registers
>> than for memory; e.g., 80 bits for registers and 64 bits for memory
>> for the case of standard double precision floating point numbers.
>> Some algorithms require a high level of control over floating point
>> semantics.  If the intrinsic function REAL with KIND parameter wp is
>> guaranteed to return a REAL (KIND=wp) result then a programmer can
>> use this to force intermediate results into main memory format,
>> never mind that the optimizing compiler may have placed the
>> intermediate in a register.
>
> For x87, I believe the only way to force memory format is to
> actually write to memory.  (Hopefully to cache.)  That is slow
> enough that compilers try to avoid doing it.

Not quite true.  Writing to memory and retrieving again rounds
to the desired precision *and* restricts the values to the exponent
range of the non-extended format.  However, if (as is the case
with Kahan summation) you only need to restrict the precision,
setting the precision bit for the X87 instructions costs only a
few clocks.  And even that small cost is only twice: once before
the summation operation and once to restore the defaults after
the summation.  Back when they still published instruction
timing I think the cost was less than 10 clocks each time.

However, there is no option in most Fortran compilers to set
that precision mode flag. :-(

-- 
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies."   --  C. A. R. Hoare

"Simplicity is prerequisite for reliability"  -- E. W. Dijkstra 

glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

>
> For x87, I believe the only way to force memory format is to
> actually write to memory.  (Hopefully to cache.)  That is slow
> enough that compilers try to avoid doing it.

No, you can instruct the co-processor to use 64-bit double precision
with this (GNU/Linux platform):

#include <fpu_control.h>

void set_double_(void)
{
  unsigned short cw;

  _FPU_GETCW(cw);
  cw &= ~0x300;
  cw |= _FPU_DOUBLE;
  _FPU_SETCW(cw);
}

Then

  external set_double
  call set_double()

from Fortran.

Chip

-- 
Charles M. "Chip" Coldwell
"Turn on, log in, tune out"
GPG Key ID: 852E052F
GPG Key Fingerprint: 77E5 2B51 4907 F08A 7E92  DE80 AFA9 9A8F 852E 052F

On 16 jun, 01:40, "bra...@mathcs.emory.edu" <bra...@mathcs.emory.edu>
wrote:
> In connection with my effort to implement Kahan summation in a manner
> that is robust against the treatments of optimizing compilers I just
> posted a request for interpretation to the J3 committee web server.
> The request is attached FYI. =A0--Bas Braams
>
> To: J3 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
 =A0 =A0 =A0 =A0 08-208
> From: Bas Braams
> Subject: Interpretation, precise f.p. semantics of the REAL intrinsic
> Date: 2008 June 15
>
> REQUEST FOR INTERPRETATION:
>
> Must the intrinsic function REAL with KIND parameter wp return a value
> that is a REAL (KIND=3Dwp) floating point number?
>
> RATIONALE FOR THE QUESTION:
>
> Computer hardware may use a wider floating-point format for registers
> than for memory; e.g., 80 bits for registers and 64 bits for memory
> for the case of standard double precision floating point numbers.
> Some algorithms require a high level of control over floating point
> semantics. =A0If the intrinsic function REAL with KIND parameter wp is
> guaranteed to return a REAL (KIND=3Dwp) result then a programmer can use
> this to force intermediate results into main memory format, never mind
> that the optimizing compiler may have placed the intermediate in a
> register.
>
> I am interested in a J3 interpretation of this matter, especially a
> loud and clear affirmative interpretation, because it appears that
> some present Fortran compilers optimize away my explicit use of the
> REAL intrinsic. =A0The context is code for compensated summation (Kahan
> summation). =A0I appreciate that parentheses are inviolable courtesy of
> the Fortran standard, but in order to have code that cannot be broken
> by an optimizing compiler I seem to need also a language mechanism to
> force intermediate results into main memory format.
>
> Bas Braams
> Chemistry Department and
> Emerson Center for Scientific Computation
> Emory University
> Atlanta, GA

Probably a naive solution, but why not:

function realf( x )
    use precision  ! To import "wp"
    real :: x
    real(wp) :: realf
    realf =3D x
end function realf

(and a similar double precision version).

It might be enough if you put them in a module to prevent
the extra precision from surfacing. But if that does not
work, you can keep them outside of any module, so that
a smart compiler will not optimise them out.

Regards,

Arjen

Charles Coldwell wrote:
> glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

>>For x87, I believe the only way to force memory format is to
>>actually write to memory.  (Hopefully to cache.)  That is slow
>>enough that compilers try to avoid doing it.

> No, you can instruct the co-processor to use 64-bit double precision

You can, but:

It doesn't reduce the exponent range to the double memory format.

It doesn't apply to all operations.

Also, even without those two, does it give the same
result in all cases?  (Maybe, but can you be sure?)

-- glen

Richard,

Thank you for your detailed comments.  Right now I am going to skip
over
much of what you wrote and address just two things.  Recall my
question:

> Must the intrinsic function REAL with KIND parameter wp return a value
> that is a REAL (KIND=wp) floating point number?

I tried to be concise.  You wrote

> Of course it returns a real(kind=wp) floating point number. It is hard
> to imagine how the standard could be much more explicit about that than
> it already is; and it is hard to imagine how an interp answer could be
> any more explicit either.

Let's be a bit more careful then and formulate the question exactly
how
it is meant.  (Except, I'm afraid I'm using the Adams et al. Fortran
95
Handbook as my reference, not the full 2003 Standard.)  Associated
with any
Type is its set of values (Handbook, Ch. 4.2), and I would like to
understand that set of values to be the set of, so to speak,
"permanent"
values that an entity of the given type can have.  Or, operationally
speaking, the set of values that is associated with how the type is
stored in memory.  I take it for granted that implementations may use
wider register arithmetic and still be in conformance with the
standard.
Therefore, it is apparently correct to interpret the standard so as to
allow intermediate values of a given type to belong to a super-set of
the
set of possible permanent values.  I do not object to this, not at
all.

Taking for granted then that the "set of values" of a type is
understood
to be the set of possible permanent values (operationally: the set of
values associated with storage into main memory) then, for the case of
the REAL intrinsic with KIND parameter wp, I ask that the result value
is
guaranteed to lie in the set of values of the REAL (KIND=wp) type.

> In this case, there is a strong conflict between the ability to optimize
> code well and the ability for the user to state exactly what the
> implementation must do, down to the last bit. Both things are strongly
> desired.

There is no conflict at all in this case.  In fact, I think that my
(implicit) proposal offers the user a beautifully simple device,
similar
to the use of explicit parentheses, to force a certain precise
floating point computation only where it matters.  Normally, in code
in
which some intermediate expression E is of type REAL (kind=wp), the
programmer will not be tempted to apply the intrinsic function and
write REAL(E,KIND=wp); it would look odd.  However, just like the
brackets in the expression x+(y-z), it doesn't look odd anymore if,
courtesy of an authoritative interpretation of the Standard, the
expression REAL(E,KIND=wp) means, no matter what is the level of
optimization, the expression E must now be projected onto the set of
permanent values associated with the REAL (KIND=wp) type.

Bas Braams

glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

> braams@mathcs.emory.edu wrote:
> (snip)
>
>> RATIONALE FOR THE QUESTION:
>
>> Computer hardware may use a wider floating-point format for registers
>> than for memory; e.g., 80 bits for registers and 64 bits for memory
>> for the case of standard double precision floating point numbers.
>> Some algorithms require a high level of control over floating point
>> semantics.  If the intrinsic function REAL with KIND parameter wp is
>> guaranteed to return a REAL (KIND=wp) result then a programmer can use
>> this to force intermediate results into main memory format, never mind
>> that the optimizing compiler may have placed the intermediate in a
>> register.
>
> For x87, I believe the only way to force memory format is to
> actually write to memory.  (Hopefully to cache.)  That is slow
> enough that compilers try to avoid doing it.

For gcc/gfortran see 


`-ffloat-store'
     Do not store floating point variables in registers, and inhibit
     other options that might change whether a floating point value is
     taken from a register or memory.

     This option prevents undesirable excess precision on machines such
     as the 68000 where the floating registers (of the 68881) keep more
     precision than a `double' is supposed to have.  Similarly for the
     x86 architecture.  For most programs, the excess precision does
     only good, but a few programs rely on the precise definition of
     IEEE floating point.  Use `-ffloat-store' for such programs, after
     modifying them to store all pertinent intermediate computations
     into variables.


-- 
Charles M. "Chip" Coldwell
"Turn on, log in, tune out"
GPG Key ID: 852E052F
GPG Key Fingerprint: 77E5 2B51 4907 F08A 7E92  DE80 AFA9 9A8F 852E 052F