DO WHILE and LESS THAN in REAL numbers

DO WHILE and LESS THAN in REAL numbers

I have two codes to reproduce the situation:

Code1

PROGRAM test
	REAL(KIND=8) dt, t, tmax
	INTEGER i, imax
	PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10)
	t = 0
	i = 0
	DO WHILE (t .LT. tmax)
		PRINT *, "before", t, i
		t = t + dt
		i = i + 1
		PRINT *, "after", t, i
	END DO
	
END PROGRAM test

Code 2

PROGRAM test
	REAL(KIND=8) dt, t, tmax
	INTEGER i, imax
	PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10)
	t = 0
	i = 0
	DO WHILE (i .LT. imax)
		PRINT *, "before", t, i
		t = t + dt
		i = i + 1
		PRINT *, "after", t, i
	END DO
	
END PROGRAM test

If I look at the output, the last 3 terms code 1 prints:

after   10.2000000000000               51
before   10.2000000000000               51
after   10.4000000000000               52

and the code 2

after   1.80000000000000                9
before   1.80000000000000                9
after   2.00000000000000               10

The code 2 is doing what it is supposed to, code 1 is doing one more loop it is not supposed to do. I believe the reason is the machine representation of numbers. But what do I do to get around it, to get a real "less than" condition?

 

Thread Topic: 

Question
5 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

You could determine the loop count before you start and use an integer counter.

E.G. nloops = NINT((tmax-tmin)/dt)

assuming tmin is set.

Testing reals for equality is always going to lead to grief unless you allow for the possibility of round off error.

Thank you for that. But then, there is an interesting (and for me embarassing). How do you determine the loop cound (I think  you should have INT( (tmax-tmin)/dt ) +1 in your post, but anyway, this is valid only for the loop of the type where I wanted DO WHILE (t .LT. tmax)  but this no longer works for the loop type DO WHILE (t .LE. tmax) and there is a difference from the above number of loops by 1 and this depends on dt by I have not been able to figure out how. It is always one less or the same as INT( (tmax-tmin)/dt ) +1

To clarify, I want the same number of loops as would have happened in the case t .LE. tmax in the case we were able to do precise arithmentics (no round off errors etc).

The phenomenon you are seeing is precisely the reason why real numbers are tricky to use as a loop variable. One improvement:

Use t = tmin + dt * (i-1)

instead of accumulating the variable t. Accumulation means that also the rounding error is accumulated.

If you seek a lower-than loop:

n = int( (tmax - tmin - 0.5d0*dt) / dt)

ought to do it

For a lower-equal:

n = int( (tmax - tmin + 0.5d0*dt) / dt)

Caveat: I have not tested these formulae in any formal sense, so I may be wrong :).

To avoid the accumulation of roundoff errors...
you avoid accumulation altogether.

PROGRAM test
 REAL(KIND=8) dt, t, tmax
 INTEGER i, imax
 PARAMETER (tmax = 10.2D0, dt = 0.2D0, imax = 10)
 t = 0
 i = 0
 DO WHILE (i*dt .LT. tmax)
  PRINT *, "before", t, i
  t = t + dt
  i = i + 1
  PRINT *, "after", t, i
 END DO
 
END PROGRAM test

Jim Dempsey

Leave a Comment

Please sign in to add a comment. Not a member? Join today