Fortran - Re: Random position on sphere generator -- numerics question

Author

Re: Random position on sphere generator -- numerics question

Ron Shepard

2006-09-02, 3:59 am

In article <1157136993.969477.56210@e3g2000cwe.googlegroups.com>,
mvukovic@nycap.rr.com wrote:

> Ron Shepard wrote:
>
> You are right. My post had a mistake in it. Here is the actual
> fortran code (one, two and half are pre-defined constants):
> DO
> call random_number( rn_angles )
> !print *,'RNs:',rn_angles
> x1=two*(rn_angles(1)-half)
> x2=two*(rn_angles(2)-half)
> ! print *,x1,x2,(x1**2+x2**2)
> if ((x1**2+x2**2).le. 1) exit
> END DO
> sq_arg=one-x1**2-x2**2
> cos1=two*x1*sqrt((sq_arg)) ! should use (abs(sq_arg))
> cos2=two*x2*sqrt((sq_arg)) ! to prevent sqrt(neg.num)
> cos3=one-two*(x1**2+x2**2)
> rn_point_on_sphere=(/cos1,cos2,cos3/)
>
> However, due to round-offs, the argument of the square root (on the
> cos1,2 lines) would sometimes be negative. So, my way of dealing with
> that was to just take the sqrt of the absolute value.

Ok, I think I see what might be the problem. In the test inside the
loop, the comparison is being computed as something equivalent to

one-(x1**2+x2**2) >= 0.0

Then, outside the loop, it is probably being computed as

(one-x1**2)-x2**2

In exact arithmetic, the two ways of computing that value would be
the same, but in finite floating point, the two values may be
different. I'm not sure, but I think the first way of computing the
value is probably better because there is only one subtraction
rather than two, and it is the subtractions of nearly equal values
that causes the roundoff problems.

What if you change the test inside the DO loop to

sq_arg=one-(x1**2+x2**2)
if (sq_arg .gt. 0.0) exit

This way you will test on exactly the argument that is used in the
sqrt() call. Furthermore, sq_arg is computed only once this way (it
is effectively computed twice in your code). Also, sqrt(sq_arg)
should be computed only once (it is computed twice in your code).

Another observation is that in the
http://mathworld.wolfram.com/SpherePointPicking.html link that you
gave earlier, it says that you should reject points for which
x1**2+x2**2==1.0. In your original code, you are accepting such
points; the above modification rejects those points, meaning that
rn_point_on_sphere(:) will never be (/0.0,0.0,-1.0/).

Another possibility has to do with mixed-precision arithmetic. Are
all of the variables and parameter constants (x1, x2, one, two,
half, cos1, cos2, cos3, rn_angles, rn_point_on_sphere) defined to be
the same TYPE and KIND? If they aren't, then that is a possible
source of these kinds of errors. You might be computing a slightly
different result in one place than the other because of different
ways of computing the expressions involving mixed-precision
variables. For this reason, I have cross-posted this to
comp.lang.fortran in case someone there can recognize what might be
the problem (your original code is written in fortran). In that
newsgroup, you will probably be asked to post the entire code,
including declarations.

$.02 -Ron Shepard

Ron Shepard

2006-09-06, 4:00 am

In article <1157483940.899987.96740@m79g2000cwm.googlegroups.com>,
mvukovic@nycap.rr.com wrote:

>
> I was aware of the multiple sqrt evaluations, but I thought it best to
> leave it to the compiler to optimize away instead of introducing
> another variable.

In general, it is easier for the compiler to optimize away a memory
store than a function reference. If you have questions about what a
particular compiler does, or doesn't do, you can examine the
assembler listing.

$.02 -Ron Shepard

glen herrmannsfeldt

2006-09-06, 4:00 am

Ron Shepard wrote:

> In article <1157483940.899987.96740@m79g2000cwm.googlegroups.com>,
> mvukovic@nycap.rr.com wrote:

(snip, someone wrote)
[color=darkred]
[color=darkred]
> In general, it is easier for the compiler to optimize away a memory
> store than a function reference. If you have questions about what a
> particular compiler does, or doesn't do, you can examine the
> assembler listing.

Common subexpression elimination has been part of optimizing Fortran
compilers for a long time, including intrinsic functions. Exactly how
good any particular compilers is at it is always a question. If the
compiler can figure it out, it has been suggested that it is better to
let the compiler do it. In this case, though, I agree that a temporary
variable makes sense to be sure that the expression is always evaluated
consistently.

-- glen

Ron Shepard

2006-09-06, 7:01 pm

In article <1157555141.009493.134930@b28g2000cwb.googlegroups.com>,
mvukovic@nycap.rr.com wrote:

> Are there directives in fortran (not the compiler option switches) to
> disable optimization of some parts of the code? Like saying: "leave
> this variable alone, don't mess with it"

It is not clear exactly what you mean by "mess with it." Some
compilers do have these kinds of directives, but they are not
standardized. There is also the VOLATILE keyword in f2003 (and it
is an extension in some earlier compilers) which inhibits register
optimization of the associated variable (i.e. it forces memory
address reads or writes whenever it is referenced).

Also, fortran has always been required to respect parentheses in
expressions. This eliminates the need for many of these kinds of
directives. Many languages are not required to respect parentheses
in arithmetic expressions. For example, C was not required to do so
until either the 89 or the 99 revision of the language (I forget
which).

$.02 -Ron Shepard

Richard Maine

2006-09-06, 7:01 pm

Ron Shepard <ron-shepard@NOSPAM.comcast.net> wrote:

> There is also the VOLATILE keyword in f2003 (and it
> is an extension in some earlier compilers) which inhibits register
> optimization of the associated variable (i.e. it forces memory
> address reads or writes whenever it is referenced).

For "forces" substitute some more wafflish word - perhaps "encourages".
The details of VOLATILE are a bit ill-defined in the standard. Partly
that is because VOLATILE is only useful in conjunction with things
outside of the standard (such as shared memory). The standard's
description of VOLATILE (which I don't think I'll drag out and quote
right now, although it isn't very long) uses, as I recall, forms of the
word "should". In standard-speak, that word means, well, pretty much
what it does in English, except that ISO makes the distinction between
"should" and "shall" explicit and strong, while we sometimes blur over
it in English. In the standard, "should" indicates something that is
recommended, but not required.

Yes, you have the right idea. I'm just being nitpicky about the word
"forces".

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

glen herrmannsfeldt

2006-09-06, 7:01 pm

In comp.lang.fortran Craig Powers <enigma@hal-pc.org> wrote:
(I wrote)

[color=darkred]
> What does VOLATILE have to do with Duff's Device? :

:

The true use for Duff's device is writing to a memory mapped
I/O device. That is, all byte are written to the same memory
location. Without VOLATILE, the store, along with the whole loop,
might get optimized away.

I don't believe jumps into Fortran DO loops are allowed.
I am not sure about WHILE loops.

-- glen

Craig Powers

2006-09-06, 7:01 pm

glen herrmannsfeldt wrote:
> In comp.lang.fortran Craig Powers <enigma@hal-pc.org> wrote:
> (I wrote)
>
>
>
> The true use for Duff's device is writing to a memory mapped
> I/O device. That is, all byte are written to the same memory
> location. Without VOLATILE, the store, along with the whole loop,
> might get optimized away.

Note my reply to myself. :)

While that was the original purpose that led to the development of the
device, it appears to me to be equally usable when the destination is a
conventional variable.

glen herrmannsfeldt

2006-09-06, 7:01 pm

In comp.lang.fortran Craig Powers <enigma@hal-pc.org> wrote:
(snip regarding Duff's Device)

> Note my reply to myself. :)

> While that was the original purpose that led to the development of the
> device, it appears to me to be equally usable when the destination is a
> conventional variable.

I saw your reply after writing mine. As Duff says, the reason is
that C supplies memcpy() which is supposed to be an efficient
implementation routine to copy data from one place to another. Otherwise,
it could be useful for that. I suppose if one was decrementing such
that it reversed the bytes in memory, that would be a useful form
not duplicating memcpy() and not needing volatile.

-- glen

Richard E Maine

2006-09-06, 7:01 pm

glen herrmannsfeldt <gah@seniti.ugcs.caltech.edu> wrote:

> I don't believe jumps into Fortran DO loops are allowed.

They aren't. (Someone will probably point out that there were cases in
f66 where it was allowed. That changed in f77. It is no longer allowed
in any cases, and compilers often enforce that.)

> I am not sure about WHILE loops.

Hint... There is no such thing as a WHILE loop. It is just a special
case of a DO loop. Note that the spelling is not

WHILE (condition)

Instead, it is

DO WHILE (condition)

Therefore, the prohibition against jumping into a DO loop applies.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain| experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Dave Thompson

2006-10-15, 10:02 pm

On Thu, 21 Sep 2006 02:19:46 GMT, "James Giles"
<jamesgiles@worldnet.att.net> wrote:

> Dave Thompson wrote:
> ...
>
> I can't find justification for that statement in the C standard. Can
> you post a section number (or page/line number pair) for the part
> that states this?
>
> In fact, I've talked to C implementors that state that C's requirements
> on expressions are much more restrictive than Fortran's. Not only
> does C have to observe the integrity of parenthesis, but also the
> specified left- or right-associativity of operators. They do *not*
> have a provision that permits "mathematically equivalent" expressions
> to be substituted for the ones actually present like Fortran does.
> That's what they say, and that's what I find in the standard document.
>
> The "as-if" rule applies. That is, if the compiler can tell that the
> result of evaluating in some other order will be the same as if it
> were evaluated in the specified order, it's free to do so. But, since
> few operations are computationally associative (or any of the other
> really convenient properties of "real" mathematics) such optimizations
> must be rare. Given the possibility of overflow, even integer addition
> cannot be reordered.
>
There is (I'm pretty sure) no specific statement; it would be under
the as-if rule plus the escape clause that before or in the absence of
the C99 option for compliance to 60559, much of the 'quality' of
floating point is at the choice of the implementor/tation. I thought I
recalled a DR, but a quick scan doesn't find any. The other
possibility would be discussion on comp.std.c, but (1) that would be
very hard to find, and (2) on reflection I may have misremembered
(more frequent) discussion about sequence points (and most operands
not necessarily being done in any 'proper' (serializable) order).

> Of course, "traditional" (pre-standard) C had much looser rules.
> What you say was true back then. Boy I really didn't like that.

- David.Thompson1 at worldnet.att.net

glen herrmannsfeldt

2006-10-17, 4:00 am

Richard E Maine wrote:

> glen herrmannsfeldt <gah@seniti.ugcs.caltech.edu> wrote:

[color=darkred]
[color=darkred]
[color=darkred]
> Eh? relatively new as in dating back at least to the first Fortran
> standard 40 years ago (f66)? Prior to that, of course, things are more
> vague, but I still think that the defacto-standards defined it.

OK, yes, I looked it up now. I must have been using C too long.

It isn't quite fair, though. The 704, 709, and 7090 are all
sign-magnitude where it is natural to do it that way.
(Divide the magnitude, XOR the sign.)

C is designed to allow what the hardware allows, though most
hardware does what Fortran requires, probably so it can run Fortran.
(Or copied after a processor designed to run Fortran.)

Blauuw & Brooks discuss the question, but don't indicate that any
processors produce only non-negative remainders.

-- glen

Sponsored Links