DO WHILE and LESS THAN in REAL numbers

[Leos_P_]

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?

[DavidWhite]

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.

[Leos_P_]

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

[Arjen_Markus]

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 :).

[jimdempseyatthecove]

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