Fortran - floating-point DO loops

Author

floating-point DO loops

robert.corbett@sun.com

2006-09-12, 7:03 pm

I know that floating-point DO loops have been removed from
the Fortran standard. I am among those who think they
never should have been added to the standard in the first
place. Nonetheless, I know that programmers continue to
use floating-point DO loops. What I would like to know is
what kinds of problems are floating-point DO loops used
to solve.

I know that floating-point DO loops are often used in
numerical integration or quadrature. If floating-point DO
loops are implemented in certain ways, they provide a
convenient way of expressing algorithms for numerical
integration. If they are implemented in some other ways,
they can produce unexpectedly bad results.

What other uses are there for floating-point DO loops?

Bob Corbett

glen herrmannsfeldt

2006-09-12, 7:03 pm

robert.corbett@sun.com wrote:
> I know that floating-point DO loops have been removed from
> the Fortran standard. I am among those who think they
> never should have been added to the standard in the first
> place. Nonetheless, I know that programmers continue to
> use floating-point DO loops. What I would like to know is
> what kinds of problems are floating-point DO loops used
> to solve.

It isn't hard to write a loop with a known loop count,
such as

DO A=0, 100.05, 0.1

The 0.05 such that even with rounding error it will
still stop at the desired point.

Some of my early Fortran (and PL/I) programs were translations
of BASIC programs, which often have floating point loops.
(Many versions of BASIC don't have integer variables.)

When I was in high school I took a programming course after
I knew way more Fortran than they were teaching. One assignment
was specifically designed to use floating point values so one
couldn't use a DO loop. (Before 1977.) I wrote mine to
calculate the appropriate number of iterations for a DO
loop, and then calculate the floating point value from the
DO variable. Because of round-off it would do one fewer
iterations than it was supposed to, so I just added one.
(I probably should have redone the iteration count properly
rounded, but adding one worked well enough.)

Because of cumulative rounding error, it is probably better
to use an integer loop and calculate the appropriate value.

Note that many other languages don't restrict loop variables,
C, C++, and Java with their generalized for() construct, and
PL/I with a real DO that allows floating point (and even COMPLEX
if you don't put a TO clause in.)

-- glen

Arjen Markus

2006-09-13, 4:01 am

glen herrmannsfeldt schreef:

>
> Because of cumulative rounding error, it is probably better
> to use an integer loop and calculate the appropriate value.
>

I have seen reals used as a DO-loop variable only to avoid
that one extra statement (or at least I think that was the
reason - it was the first time ever I saw that reals could
be used like that in FORTRAN 77).

I have never used them myself in that way, at least not
recently enough to remember it, except for one time:
to see what would happen with different compilers :).

> Note that many other languages don't restrict loop variables,
> C, C++, and Java with their generalized for() construct, and
> PL/I with a real DO that allows floating point (and even COMPLEX
> if you don't put a TO clause in.)
>

For-loops in those languages are just shorthands for while():
you combine the initialisation, the condition and the advancing
in one statement. Fortran's do-loops are quite different as the
number of iterations is fixed and you can not use more than
one loop variable.

Regards,

Arjen

David Flower

2006-09-13, 4:01 am

glen herrmannsfeldt

2006-09-13, 4:01 am

Arjen Markus wrote:

(snip)

> For-loops in those languages are just shorthands for while():
> you combine the initialisation, the condition and the advancing
> in one statement. Fortran's do-loops are quite different as the
> number of iterations is fixed and you can not use more than
> one loop variable.

While it is defined with the iteration count fixed at the beginning,
the rules also don't allow modifying the DO variable. It is, then,
allowed to use the DO variable to count the iterations, with undefined
results if it is actually changed. The end value and increment value
do have to be copied to temporary variables if they can possibly change
inside the loop. As far as I know, they don't have to be copied if the
compiler believes that they don't change inside the loop, again with
undefined effects if they actually do change.

-- glen

Brooks Moses

2006-09-13, 4:01 am

robert.corbett@sun.com wrote:
> I know that floating-point DO loops are often used in
> numerical integration or quadrature. If floating-point DO
> loops are implemented in certain ways, they provide a
> convenient way of expressing algorithms for numerical
> integration. If they are implemented in some other ways,
> they can produce unexpectedly bad results.
>
> What other uses are there for floating-point DO loops?

The case where I find that I would most often want them is for setting
up cartesian grids and the like -- for example, I might want to create
an array of gridpoints going from 0.0 to 1.0 in steps of 0.05.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

Ron Shepard

2006-09-13, 7:02 pm

In article <4507BFC6.1030102@cits1.stanford.edu>,
Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:

> The case where I find that I would most often want them is for setting
> up cartesian grids and the like -- for example, I might want to create
> an array of gridpoints going from 0.0 to 1.0 in steps of 0.05.

I think this is one of the situations where floating point DO loops
should *not* be used. If you have an array, you presumably know
ahead of time how many elements are in that array. The DO loop
could result in less elements, that many elements, or more elements.
Unexpected things can happen in the first case, and illegal memory
references can occur in the last case.

The only time a floating point DO loop should be used is when the
exact iteration count is not important and the exact value of the DO
variable itself each iteration is not important. There are
situations like this, but because the construct can be abused in so
many other cases, it is probably a good thing that it is now dropped
from the language.

$.02 -Ron Shepard

Ron Shepard

2006-09-13, 7:02 pm

In article < RPSdnVVMDLQXKJrYnZ2dnUVZ_smdnZ2d@comcast
.com>,
glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:

> While it is defined with the iteration count fixed at the beginning,
> the rules also don't allow modifying the DO variable. It is, then,
> allowed to use the DO variable to count the iterations, with undefined
> results if it is actually changed. The end value and increment value
> do have to be copied to temporary variables if they can possibly change
> inside the loop. As far as I know, they don't have to be copied if the
> compiler believes that they don't change inside the loop, again with
> undefined effects if they actually do change.

I don't think this last sentence is correct. The lower bound, upper
bound, and increment in a DO loop are expressions (not variables),
and the variables used in those expressions are allowed to be
modified within the do loop (unless there are other reasons why they
cannot be modified). By "allowed" I mean that the result of
modifying those variables is defined by the standard.

I don't have the standard handy for a reference, but that is the way
I remember it.

$.02 -Ron Shepard

robert.corbett@sun.com

2006-09-13, 7:02 pm

Ron Shepard wrote:
>
> The only time a floating point DO loop should be used is when the
> exact iteration count is not important and the exact value of the DO
> variable itself each iteration is not important. There are
> situations like this, but because the construct can be abused in so
> many other cases, it is probably a good thing that it is now dropped
> from the language.

In the main, I agree. The problem I have seen is that when
people try to use integer DO loops to range over real values,
they often do a worse job than the compiler would.

Part of the reason I asked how people use floating-point
DO loops is that I am thinking about writing a note explaining
how to emulate floating-point DO loops using integer DO loops.
For the specific case of floating-point DO loops for numerical
integration, I know the properties the translated loops should
have. For other applications, other properties might be
desirable, hence my question.

Bob Corbett

glen herrmannsfeldt

2006-09-13, 7:02 pm

Ron Shepard <ron-shepard@nospam.comcast.net> wrote:
(snip on DO loops)

[color=darkred]
> I don't think this last sentence is correct. The lower bound, upper
> bound, and increment in a DO loop are expressions (not variables),
> and the variables used in those expressions are allowed to be
> modified within the do loop (unless there are other reasons why they
> cannot be modified). By "allowed" I mean that the result of
> modifying those variables is defined by the standard.

If they are expressions, they should be evaluate once and temporaries
used in the loop. If they are variables and there are no statements
to change them in the loop, and no statements that should be able to
indirectly change them, I believe it is legal not to generate
temporaries. If one changes anyway, maybe due to not following
aliasing rules, going outside an array, or from another thread
(in Fortran versions that don't allow multithreading) then the
loop count will be wrong. For a program following all the
requirements of the standard, modifying loop parameter variables
or expressions must give the right results.

As far as allowing expressions, I never liked C loops like:

for(i=0;i<sqrt(a);i++)

which require the sqrt() to be evaluated every iteration.
In cases like that, a temporary variable is especially nice.

-- glen

-- glen

dpb

2006-09-13, 7:02 pm

robert.corbett@sun.com wrote:
....
>
> What other uses are there for floating-point DO loops?

In my former place of reading/maintaining code in days of yore they
were invariably the mistakes you provided examples for! :(

I can see some ways to be helpful in quadrature but I can't recall I've
ever actually written one "in anger".

Brooks Moses

2006-09-13, 7:02 pm

robert.corbett@sun.com wrote:
> Ron Shepard wrote:
>
> In the main, I agree. The problem I have seen is that when
> people try to use integer DO loops to range over real values,
> they often do a worse job than the compiler would.

And I agree as well -- my example was a case where the first naive
implementation that comes to mind is a real-valued DO loop, not a case
where I would actually choose to use one if I could.

> Part of the reason I asked how people use floating-point
> DO loops is that I am thinking about writing a note explaining
> how to emulate floating-point DO loops using integer DO loops.
> For the specific case of floating-point DO loops for numerical
> integration, I know the properties the translated loops should
> have. For other applications, other properties might be
> desirable, hence my question.

I don't know if this is directly useful, but I often find myself doing a
lot of things like:

-------------------------------------------------------------
real :: x1 = 0.0, x2 = 10.0 ! endpoints of interval
real :: dx_max = 0.3 ! maximum step size
real :: dx, x
integer :: i, n_steps

n_steps = ceiling((x2-x1)/dx_max)
dx = (x2-x1)/real(n_steps)

do i = 1, n_steps
x = x1 + i * dx
...
end do
-------------------------------------------------------------

This has a couple of problems that I can see. First, due to round-off,
n_steps might be 1 larger than needed if dx_max evenly divides into the
interval. This could be fairly easily fixed by adding a small epsilon
to the argument of ceiling in the n_steps calculation.

Second, the value of x at i=n_steps may not be exactly equal to x2. I'm
not sure how to fix this except with an explicit end-of-loop clause.

So, any note that contained advice on how to implement this sort of
thing in a way that avoids those problems (or ones that I haven't
noticed, particularly!) would be one that I would find useful.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

Rich Townsend

2006-09-13, 7:02 pm

Brooks Moses wrote:
> robert.corbett@sun.com wrote:
>
> And I agree as well -- my example was a case where the first naive
> implementation that comes to mind is a real-valued DO loop, not a case
> where I would actually choose to use one if I could.
>
>
> I don't know if this is directly useful, but I often find myself doing a
> lot of things like:
>
> -------------------------------------------------------------
> real :: x1 = 0.0, x2 = 10.0 ! endpoints of interval
> real :: dx_max = 0.3 ! maximum step size
> real :: dx, x
> integer :: i, n_steps
>
> n_steps = ceiling((x2-x1)/dx_max)
> dx = (x2-x1)/real(n_steps)
>
> do i = 1, n_steps
> x = x1 + i * dx
> ...
> end do
> -------------------------------------------------------------
>
> This has a couple of problems that I can see. First, due to round-off,
> n_steps might be 1 larger than needed if dx_max evenly divides into the
> interval. This could be fairly easily fixed by adding a small epsilon
> to the argument of ceiling in the n_steps calculation.
>
> Second, the value of x at i=n_steps may not be exactly equal to x2. I'm
> not sure how to fix this except with an explicit end-of-loop clause.
>
> So, any note that contained advice on how to implement this sort of
> thing in a way that avoids those problems (or ones that I haven't
> noticed, particularly!) would be one that I would find useful.

Once you know n_steps and x1, x2, this code is reasonably bomb-proof -- in that
it gets the endpoints spot-on:

do i = 1,n_steps
x = (x1*(n_steps-i) + x2*(i-1))/(n_steps-1)
...
end do

I use this idiom everywhere -- basically, x is recalculated each time as a
linear interpolation between x1 and x2. Not as fast as the x = x1 + i*dx
approach, of course, but immune to the roundoff problem.

cheers,

Rich

Rich Townsend

2006-09-13, 7:02 pm

Rich Townsend wrote:
> Brooks Moses wrote:
>
> Once you know n_steps and x1, x2, this code is reasonably bomb-proof -- in that
> it gets the endpoints spot-on:
>
> do i = 1,n_steps
> x = (x1*(n_steps-i) + x2*(i-1))/(n_steps-1)
> ...
> end do
>
> I use this idiom everywhere -- basically, x is recalculated each time as a
> linear interpolation between x1 and x2. Not as fast as the x = x1 + i*dx
> approach, of course, but immune to the roundoff problem.

Actually, on closer inspection, my claim that the endpoints are 'spot-on' is
wrong -- there is no guarantee that, for instance

x1 == (x1*n_steps)/n_steps

So I'm all ears for a better solution too.

cheers,

Rich

glen herrmannsfeldt

2006-09-13, 7:02 pm

Rich Townsend <rhdt@barvoidtol.udel.edu> wrote:
(snip)
(snip)[color=darkred]
> Actually, on closer inspection, my claim that the endpoints are 'spot-on' is
> wrong -- there is no guarantee that, for instance

> x1 == (x1*n_steps)/n_steps

How about :

x = x1*(float(n_steps-i)/float(n_step-1))+x2*(float(i-1)/float(n_steps-1))

I don't know if it is more or less accurate for intermediate points,
but for the end points it only depends on float(i)/float(i) being
1.0 which seems pretty reliable (for i.NE.0).

-- glen

Steven G. Kargl

2006-09-13, 7:02 pm

In article <ee9k6j$aef$1@scrotar.nss.udel.edu>,
Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:
> Rich Townsend wrote:
>
> Actually, on closer inspection, my claim that the endpoints are 'spot-on' is
> wrong -- there is no guarantee that, for instance
>
> x1 == (x1*n_steps)/n_steps
>
> So I'm all ears for a better solution too.
>

A long time ago when I was more naive about Fortran,
I used real do-loop variables. Now days, I tend to
use

dx = (x2 - x1) / (n - 1)
do i =1, n
x = x1 + (i - 1) * dx
end do

x = x1 is spot on and x = x2 will at most have the
rounding error of a floating point addition and
multiplication. The difference is of course that
you specify the number of steps and let dx be whatever
value is computed. If you specific dx, then you're
stuck with int, nint, ceiling, floor, etc. and x = x2
will most like have a larger error.

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

Richard E Maine

2006-09-13, 7:02 pm

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

> How about :
>
> x = x1*(float(n_steps-i)/float(n_step-1))+x2*(float(i-1)/float(n_steps-1))

Any other matters aside, I'd recommend updating to at least Fortran 77.
:-)

The float intrinsic is a holdover from Fortran 66. While it is still in
the standard (even f2003 - and I recall no proposal to delete it from
f2008 either) for compatibility, you'd have to work darned hard to find
any case where I'd say it made sense to use it today. This isn't even
close to such a case.

--
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

jwm

2006-09-13, 7:02 pm

Rich Townsend wrote:
> Rich Townsend wrote:
>
> Actually, on closer inspection, my claim that the endpoints are 'spot-on' is
> wrong -- there is no guarantee that, for instance
>
> x1 == (x1*n_steps)/n_steps
>

I really don't understand what you guys are talking about (basically,
because I'm used to the F90/95 Standards, and I've just ignored what's
deemed as obsolecent/deleted), so what I'm asking is an explanation of
what's the advantage (besides the number of lines/variables) in using
the following code:
...
real curr
do curr = 0.0, 10.0, 3.0

!block of code

enddo

And what's wrong with using this one:
...
real incr, curr, last
curr = 0.0
incr = 3.0
last = 10.0
do

!same block of code

curr = curr + incr
if (curr > last) exit
enddo

Thanks in advace for the explanation.

Dick Hendrickson

2006-09-13, 7:02 pm

jwm wrote:
> Rich Townsend wrote:
>
>
>
> I really don't understand what you guys are talking about (basically,
> because I'm used to the F90/95 Standards, and I've just ignored what's
> deemed as obsolecent/deleted), so what I'm asking is an explanation of
> what's the advantage (besides the number of lines/variables) in using
> the following code:
> ...
> real curr
> do curr = 0.0, 10.0, 3.0

There's no problem or advantage with this specific loop.
Floating point arithmetic on small integers is exact on
any realistic machine. The advantage, or problem, comes
with using fractions rather than integers. Something like

do curr = 0.0, 1.0, .1

How many iterations? What will the final value of curr
be? Only the FP divider knows for sure ;)

>
> !block of code
>
> enddo
>
> And what's wrong with using this one:
> ...
> real incr, curr, last
> curr = 0.0
> incr = 3.0
> last = 10.0
> do
>
> !same block of code
>
> curr = curr + incr
> if (curr > last) exit
> enddo
>
> Thanks in advace for the explanation.

There's nothing "wrong" with it. But, compiler vendors
have spent millions of man hours learning how to do
optimizations on DO loops where the iteration count can
be determined from the DO statement. Writing the loop
this way potentially loses optimzations, especially if the
compiler can't prove that curr, incr, or last are not
changed somewhere else. This form also has an uncertainity
if incr is fractional, say .01. Will the last value be
..99000000000000000000001 or .999999999999999999999999?

Dick Hendrickson
>

glen herrmannsfeldt

2006-09-13, 7:02 pm

Richard E Maine <nospam@see.signature> wrote:

> The float intrinsic is a holdover from Fortran 66. While it is still in
> the standard (even f2003 - and I recall no proposal to delete it from
> f2008 either) for compatibility, you'd have to work darned hard to find
> any case where I'd say it made sense to use it today. This isn't even
> close to such a case.

Yes. In this case, I wanted to emphasize the division of the floating
point values converted from integers. That is, as a description
of the algorithm. For an actual implementation, one could do it
differently.

(Also, that I don't really like the reuse of the REAL function from
its traditional and mathematical use as the real part of a complex value.)

It could have been worse, I could have used DFLOAT!

-- glen

jwm

2006-09-13, 7:02 pm

Dick Hendrickson wrote:
> jwm wrote:
>
> There's no problem or advantage with this specific loop.
> Floating point arithmetic on small integers is exact on
> any realistic machine. The advantage, or problem, comes
> with using fractions rather than integers. Something like
>
> do curr = 0.0, 1.0, .1
>
> How many iterations? What will the final value of curr
> be? Only the FP divider knows for sure ;)

According to this, I should avoid using:

do curr = 0.0, 1.0, 0.001

and rather trying:

do icurr = 0, 1000
curr = 0.001 * real(icurr)
enddo

OK, this one shouldn't be used, since it's similar to what was stated
in the previous posts. But, because of what you said about floating
point arithmetic on "small integers", will this one work just fine?:

do fcurr = 0.0, 1000.0, 1.0
curr = 0.001 * fcurr
enddo

Dick Hendrickson

2006-09-13, 7:02 pm

jwm wrote:
> Dick Hendrickson wrote:
>
>
>
> According to this, I should avoid using:
>
> do curr = 0.0, 1.0, 0.001
>
> and rather trying:
>
> do icurr = 0, 1000
> curr = 0.001 * real(icurr)
> enddo
>
> OK, this one shouldn't be used, since it's similar to what was stated
> in the previous posts. But, because of what you said about floating
> point arithmetic on "small integers", will this one work just fine?:
>
> do fcurr = 0.0, 1000.0, 1.0
> curr = 0.001 * fcurr
> enddo
>
It should. The key is that integers less than 2**23
(depends on the machine and precision, etc.) or so
are exactly representable in floating point formats
and thus +, - and * should be exact. But, most common
decimal fractions, like .1 or .05, are not exactly
representable in floating poing and thus arithmetic
on them will not be exact. Note that divides of
floating point integers often produces results that
are not exactly representable. The serious round off
problems occur when you repeatedly add up inexact
values. The loop
curr = -0.001
do icurr = 0, 1000
curr = curr + 0.001
enddo
is clearly (well, to me anyhow ;) ) inferior to the
others because you are repeatedly adding an inexact
number. In reality, this loop is equivalent to
do curr = 0.0, 1.0., 0.001
because the compiler will undoubtable implement it
as
curr = curr + 0.001
as the incrementation step.

The simple rule is that floating point Do loops are
like cliches: avoid them like the plague!

Dick Hendrickson

robert.corbett@sun.com

2006-09-13, 7:02 pm

Dick Hendrickson wrote:

> In reality, this loop is equivalent to
> do curr = 0.0, 1.0., 0.001
> because the compiler will undoubtable implement it
> as
> curr = curr + 0.001
> as the incrementation step.

I am not sure what you are saying. You appear to be
saying that compilers undoubtedly implement the
increment step of a floating-point DO loop by summation.
Some compilers do it that way, but most of the
compilers to which I have access do not. As you note,
that implementation of floating-point DO loops is
unsuitable in many situations.

Bob Corbett

Brooks Moses

2006-09-13, 10:05 pm

jwm wrote:
> According to this, I should avoid using:
>
> do curr = 0.0, 1.0, 0.001
>
> and rather trying:
>
> do icurr = 0, 1000
> curr = 0.001 * real(icurr)
> enddo

Yes.

> OK, this one shouldn't be used, since it's similar to what was stated
> in the previous posts. But, because of what you said about floating
> point arithmetic on "small integers", will this one work just fine?:
>
> do fcurr = 0.0, 1000.0, 1.0
> curr = 0.001 * fcurr
> enddo

In theory, that exact form will work without trouble -- so long as your
compiler happens to ignore the fact that it's not legal Fortran code
according to the Fortran 95 standard.

It is, however, an unwise thing to do. You don't gain anything by using
a real number there, and aside from the fact that it makes your code
nonstandard for no good reason, there's also the temptation to change
that stepsize of 1.0 to something that isn't an integer -- particularly
since in practice it will most likely be coming from a variable
somewhere else.

Besides which, if you have a number that represents integer data (the
number of iterations) and needs to only be set to integer values for the
program to work right, then it should be declared integer. That's half
the reason why the type system exists.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

robert.corbett@sun.com

2006-09-13, 10:05 pm

Brooks Moses wrote:
> jwm wrote:
>
> Yes.
>
>
> In theory, that exact form will work without trouble -- so long as your
> compiler happens to ignore the fact that it's not legal Fortran code
> according to the Fortran 95 standard.
>
> It is, however, an unwise thing to do. You don't gain anything by using
> a real number there, and aside from the fact that it makes your code
> nonstandard for no good reason, there's also the temptation to change
> that stepsize of 1.0 to something that isn't an integer -- particularly
> since in practice it will most likely be coming from a variable
> somewhere else.
>
> Besides which, if you have a number that represents integer data (the
> number of iterations) and needs to only be set to integer values for the
> program to work right, then it should be declared integer. That's half
> the reason why the type system exists.

The translation of the floating-point DO loop using a floating-point
DO loop will be more efficient on many machines because of the
cost of doing the integer-to-floating-point conversion before
multiplying.

Bob Corbett

glen herrmannsfeldt

2006-09-13, 10:05 pm

robert.corbett@sun.com wrote:

(snip)

(snip)
[color=darkred]
> The translation of the floating-point DO loop using a floating-point
> DO loop will be more efficient on many machines because of the
> cost of doing the integer-to-floating-point conversion before
> multiplying.

The floating point compare, and getting the result out to do
the conditional branch, may be slower than an integer compare
and branch. I might guess that the floating point loop is
still faster, but it might not be.

For x87, there is a complicated system to transfer the floating
point flags to AX, then load them into the flags register for
the conditional branch, the last step using an instruction left
over from the 8080. On the other hand, from CX to the x87 goes
through memory, unless they recently added a more direct path.

-- glen

Tim Prince

2006-09-14, 8:01 am

glen herrmannsfeldt wrote:
> robert.corbett@sun.com wrote:
>
> (snip)
>
>
> (snip)
>
>
> The floating point compare, and getting the result out to do
> the conditional branch, may be slower than an integer compare
> and branch. I might guess that the floating point loop is
> still faster, but it might not be.
>
> For x87, there is a complicated system to transfer the floating
> point flags to AX, then load them into the flags register for
> the conditional branch, the last step using an instruction left
> over from the 8080. On the other hand, from CX to the x87 goes
> through memory, unless they recently added a more direct path.
>
> -- glen
>
All recent compilers for IA (last 3 years or so) know to use SSE
instructions for these translations, even for the CPUs like pentium-m
which generally prefer x87. I believe a similar time has elapsed since
the last CPU which supported x87 but not SSE2 was made. It would be
difficult to find meaningful cases where the real DO index could improve
performance.

My understanding of the x87 standard is that the DO loop involves a
calculation of an exact number of trips before entering the loop, even
when the counter is real. In principle, then, the number of trips would
not change simply because of changes in compiler options, unless one
chooses an option which promotes real to double or the like.

glen herrmannsfeldt

2006-09-14, 10:00 pm

Tim Prince <timothyprince@sbcglobal.net> wrote:
(snip about floating point DO loops)

> All recent compilers for IA (last 3 years or so) know to use SSE
> instructions for these translations, even for the CPUs like pentium-m
> which generally prefer x87. I believe a similar time has elapsed since
> the last CPU which supported x87 but not SSE2 was made. It would be
> difficult to find meaningful cases where the real DO index could improve
> performance.

Still, many processors keep the floating point and fixed point
logic on separate parts of the chip. That is partly the cause
of slow conversions between fixed and float, and also tends to
slow down floating point conditional branches. If you really
want to know, you have to study a specific processor, and
the answer will only be true for that processor.

-- glen

robin

2006-09-17, 8:05 am

robert.corbett@sun.com wrote in message <1158106931.914881.89770@p79g2000cwp.googlegroups.com>...
>I know that floating-point DO loops have been removed from
>the Fortran standard. I am among those who think they
>never should have been added to the standard in the first
>place. Nonetheless, I know that programmers continue to
>use floating-point DO loops. What I would like to know is
>what kinds of problems are floating-point DO loops used
>to solve.
>
>I know that floating-point DO loops are often used in
>numerical integration or quadrature. If floating-point DO
>loops are implemented in certain ways, they provide a
>convenient way of expressing algorithms for numerical
>integration. If they are implemented in some other ways,
>they can produce unexpectedly bad results.
>
>What other uses are there for floating-point DO loops?

Any loop where the precise number of iterations is
unimportant;
Any loop executed a great number of times where
the repeated conversion from integer to real (which may include
multiplication and/or division) becomes significant.

robin

2006-09-17, 8:05 am

glen herrmannsfeldt wrote in message ...
>robert.corbett@sun.com wrote:
>
>It isn't hard to write a loop with a known loop count,
>such as
>
> DO A=0, 100.05, 0.1
>
>The 0.05 such that even with rounding error it will
>still stop at the desired point.
>
>Some of my early Fortran (and PL/I) programs were translations
>of BASIC programs, which often have floating point loops.
>(Many versions of BASIC don't have integer variables.)
>
>When I was in high school I took a programming course after
>I knew way more Fortran than they were teaching. One assignment
>was specifically designed to use floating point values so one
>couldn't use a DO loop. (Before 1977.) I wrote mine to
>calculate the appropriate number of iterations for a DO
>loop, and then calculate the floating point value from the
>DO variable. Because of round-off it would do one fewer
>iterations than it was supposed to, so I just added one.
>(I probably should have redone the iteration count properly
>rounded, but adding one worked well enough.)
>
>Because of cumulative rounding error, it is probably better
>to use an integer loop and calculate the appropriate value.
>
>Note that many other languages don't restrict loop variables,
>C, C++, and Java with their generalized for() construct, and
>PL/I with a real DO that allows floating point (and even COMPLEX
>if you don't put a TO clause in.)

And since you mention PL/I, it should be pointed out that the problem
doesn't arise in PL/I because the loop can be written--
DECLARE A FIXED DECIMAL (5,2);
DO A=0 TO 100.05 BY 0.1
...
END;
in which A is a decimal variable, and all values (including
the constants 100.05 and 0.1) are held exactly.

robin

2006-09-17, 8:05 am

Dick Hendrickson wrote in message ...
>jwm wrote:
>There's no problem or advantage with this specific loop.
>Floating point arithmetic on small integers is exact on
>any realistic machine. The advantage, or problem, comes
>with using fractions rather than integers. Something like
>
> do curr = 0.0, 1.0, .1
>
>How many iterations? What will the final value of curr
>be? Only the FP divider knows for sure ;)

>
>There's nothing "wrong" with it. But, compiler vendors
>have spent millions of man hours learning how to do
>optimizations on DO loops where the iteration count can
>be determined from the DO statement. Writing the loop
>this way potentially loses optimzations, especially if the
>compiler can't prove that curr, incr, or last are not
>changed somewhere else. This form also has an uncertainity
>if incr is fractional, say .01. Will the last value be
>.99000000000000000000001 or .999999999999999999999999?

Probably neither, as the final value is likely to be
out by as much as n_steps in the least-significant digit.

robin

2006-09-17, 8:05 am

glen herrmannsfeldt wrote in message ...
>Ron Shepard <ron-shepard@nospam.comcast.net> wrote:
>(snip on DO loops)
>
>
>
>If they are expressions, they should be evaluate once and temporaries
>used in the loop. If they are variables and there are no statements
>to change them in the loop, and no statements that should be able to
>indirectly change them, I believe it is legal not to generate
>temporaries. If one changes anyway, maybe due to not following
>aliasing rules, going outside an array, or from another thread
>(in Fortran versions that don't allow multithreading) then the
>loop count will be wrong.

No it won't and can't.
In Fortran, the loop count is determined BEFORE loop entry,
and changing the value of any variables used for the
initial value. the final value, and the increment,
cannot change the number of iterations.

> For a program following all the
>requirements of the standard, modifying loop parameter variables
>or expressions must give the right results.

Irrelevant C material deleted.

glen herrmannsfeldt

2006-09-18, 8:00 am

robin wrote:

(snip)

> And since you mention PL/I, it should be pointed out that the problem
> doesn't arise in PL/I because the loop can be written--
> DECLARE A FIXED DECIMAL (5,2);
> DO A=0 TO 100.05 BY 0.1
> ...
> END;
> in which A is a decimal variable, and all values (including
> the constants 100.05 and 0.1) are held exactly.

That is why I specifically said floating point, which
FIXED DEC(5,2) is not.

-- glen

robin

2006-09-19, 10:01 pm

glen herrmannsfeldt wrote in message ...
>robin wrote:
>
>
>That is why I specifically said floating point, which
> FIXED DEC(5,2) is not.

Which is why I specifcally mentioned that the problem doesn't arise
in PL/I.

In PL/I, floating-point may be DECIMAL, using decimal hardware
if available, and in that case, it doesn't arise either.

Brooks Moses

2006-09-20, 4:01 am

robin wrote:
> glen herrmannsfeldt wrote in message ...
>
> Which is why I specifcally mentioned that the problem doesn't arise
> in PL/I.
>
> In PL/I, floating-point may be DECIMAL, using decimal hardware
> if available, and in that case, it doesn't arise either.

DO A = 0 to PI by PI/100.

That has problems in decimal, no?

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

glen herrmannsfeldt

2006-09-20, 7:01 pm

Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
> robin wrote:
(snip on floating point DO loops in any language)

[color=darkred]
> DO A = 0 to PI by PI/100.

> That has problems in decimal, no?

I once thought of a floating point format with one bit that meant
to multiply the value by pi. It would be convenient for radian
arguments to trigonometric functions, for one.

Decimal helps in some common cases, but not all even for rational
decimal constants.

-- glen

Kevin G. Rhoads

2006-09-21, 7:04 pm

>In PL/I, floating-point may be DECIMAL, using decimal hardware
>if available, and in that case, it doesn't arise either.

Implementations of Fortran can, and have, supported decimal arithmetic,
either fixed or floating point.

Microsoft Fortran v3.3, for example, had the option of using BCD decimal
coding available on x86's in the RTL, at least.

And many of IBM's Fortran compilers for decimal machines (e.g., IBM 1620)
used decimal arithmetic. (Although IIRC you could load the ADD and MULTIPLY
tables in the 1620 I with values corresponding to octal arithmetic and
thus fake up Decimal Coded Octal. With the model II 1620 the ADD was
in hardware and that trick became less useful ...)

robin

2006-09-28, 7:00 pm

Kevin G. Rhoads wrote in message <4512CD1C.CF5204D1@alum.mit.edu>...
>
>Implementations of Fortran can, and have, supported decimal arithmetic,
>either fixed or floating point.

But it isn't part of the language, and where decimal _was_
provided, it was because none other hardware was available.

>Microsoft Fortran v3.3, for example, had the option of using BCD decimal
>coding available on x86's in the RTL, at least.

Again, only because most PCs of that era didn't have a
numeric co-processor.

>And many of IBM's Fortran compilers for decimal machines (e.g., IBM 1620)
>used decimal arithmetic.

Like I said, only because no other (floating binary) hardware was available.

robin

2006-09-28, 7:00 pm

Brooks Moses wrote in message <45110D9C.6010200@cits1.stanford.edu>...
>robin wrote:
>
>DO A = 0 to PI by PI/100.
>
>That has problems in decimal, no?

That is irrelevant. The question is in regard to
values with two decimal places (and, of course,
similar types of problems).

But if you _did_ want to do that, you'd write:

DO Angle = 0 TO 180 BY 1.8;
A = Angle*Pi/180;
...
END;
where Angle is defined as FIXED DECIMAL (3,1).
A itself would be floating-point, probably.

Alternatively, the increment (viz. 180/pi) could be stored
with 30 decimal places, so that the effects of truncation
errors were very small.

Sponsored Links