Incomprehensible (at least for me) behaviour in a do loop

[DomBernard]

I've extracted from a larger code the following code where the DO loop is the cause of a bug:

PROGRAM MODIF_GEOM3D
IMPLICIT NONE
INTEGER :: NX,NY,NZ,I
REAL(KIND=4) :: IVF1,IVF2
! REAL(KIND=8) :: IVF1,IVF2
! INTEGER :: IVF1,IVF2
NX = 400**3
DO I=1,NX
IVF2 = IVF1
IVF1 = I
IF(IVF1==IVF2) THEN
WRITE(6,*) I,IVF1,IVF2
STOP
ENDIF
ENDDO
WRITE(6,*) 'IVF1 = ',IVF1
END PROGRAM MODIF_GEOM3D

After compiling with gfortran on 3 different computers I've obtained the following:

[uxmalrog_DB/Geom/src] bernard% gfortran Simple_Test.FOR
[uxmalrog_DB/Geom/src] bernard% ./a.out
16777217 16777216. 16777216.

If IVF1 and IVF2 are declared as REAL 8 or INTEGER the result is now:

[uxmalrog_DB/Geom/src] bernard% gfortran Simple_Test.FOR
[uxmalrog_DB/Geom/src] bernard% ./a.out
IVF1 = 64000000.000000000

This is of course the logical one.

Any explanation?

 

[xwb]

IVF1 is not initialized

 

[DomBernard]

Completely right.
I've deleted so many lines to look for the core of the problem that the initialisation line disappeared.

PROGRAM MODIF_GEOM3D
IMPLICIT NONE
INTEGER :: NX,NY,NZ,I
REAL(KIND=4) :: IVF1,IVF2
! REAL(KIND=8) :: IVF1,IVF2
! INTEGER :: IVF1,IVF2
NX = 400**3

IVF1 = 0

DO I=1,NX
IVF2 = IVF1
IVF1 = I
IF(IVF1==IVF2) THEN
WRITE(6,*) I,IVF1,IVF2
STOP
ENDIF
ENDDO
WRITE(6,*) 'IVF1 = ',IVF1
END PROGRAM MODIF_GEOM3D

Result is unchanged.

Thanks for the basic and relevant remark.

 

[xwb]

Real numbers have a precision of approximately 6/7 digits, double precision about 15. Since 400**3 has 8 digits, and it is being stored in a number with 6 digit precision, the program stops as soon as the 6th or 7th digit stops changing.

 

[Beilstwh]

This is working as designed. Real and double is an approximation of the number and the larger the number, I less certain it is. Use INTEGER or "DOUBLE PRECISION" to get the correct values. While DOUBLE PRECISION is also a floating number, it uses 8 bits and is much more likely to support larger numbers.

PROGRAM MODIF_GEOM3D
IMPLICIT NONE
INTEGER :: NX,NY,NZ,I
! REAL(KIND=4) :: IVF1,IVF2
DOUBLE PRECISION :: IVF1,IVF2
! REAL(KIND=8) :: IVF1,IVF2
! INTEGER :: IVF1,IVF2
NX = 400**3

IVF1 = 0

DO I=1,NX
IVF2 = IVF1
IVF1 = I
IF(IVF1==IVF2) THEN
WRITE(6,*) I,IVF1,IVF2
STOP
ENDIF
ENDDO
WRITE(6,*) 'IVF1 = ',IVF1
END PROGRAM MODIF_GEOM3D

Bill
Lead Application Developer
New York State, USA

 

[DomBernard]

Thanks for the suggestions.
I want to say first that the objective is not to run this piece of code.
This is only a test code reproducing a behaviour that generates a problem in a larger one.
I'm testing a lot of hypothesis around this problem and the number of correct digits has been one.
Here is a variant of the code:

PROGRAM MODIF_GEOM3D
IMPLICIT NONE
INTEGER :: NX,NY,NZ,I,ICONT
REAL(KIND=4) :: IVF1,IVF2
! REAL(KIND=8) :: IVF1,IVF2
! INTEGER :: IVF1,IVF2
NX = 400**3
ICONT = 0
IVF1 = 0
DO I=1,NX
IVF2 = IVF1
IVF1 = I
IF(ICONT /= 0) WRITE(6,*) I,IVF1,IVF2
IF(IVF1==IVF2) THEN
WRITE(6,'(1x,I10,3x,F20.7,3x,F20.7)') I,IVF1,IVF2
READ(5,*) ICONT
IF(ICONT==0) STOP
ENDIF
ENDDO
WRITE(6,*) 'IVF1 = ',IVF1
END PROGRAM MODIF_GEOM3D

with the corresponding results:

[bernard@hpc2 Evolution_Diffusion]$ Test.Test
16777217 16777216.0000000 16777216.0000000
1
16777218 16777218. 16777216.
16777219 16777220. 16777218.
16777220 16777220. 16777220.
16777220 16777220.0000000 16777220.0000000
1
16777221 16777220. 16777220.
16777221 16777220.0000000 16777220.0000000
1
16777222 16777222. 16777220.
16777223 16777224. 16777222.
16777224 16777224. 16777224.
16777224 16777224.0000000 16777224.0000000
1
16777225 16777224. 16777224.
16777225 16777224.0000000 16777224.0000000
1
16777226 16777226. 16777224.
16777227 16777228. 16777226.
16777228 16777228. 16777228.
16777228 16777228.0000000 16777228.0000000
0

My impression is that it is not a problem of digit.
Many thanks for helping.

 

[DomBernard]

After a few minutes it appears to me that the previous test was not relevant.
In fact, REAL 4 are coded on 32bits, 1 for the sign, 8 for the exponent and 23 for the fraction.
The maximum value given by my code is 16777216 = 2**24 => this explain the results.

Thanks for helping.

 

[Beilstwh]

Exactly. Your need to use DOUBLE PRECISION if you are going to have fractions in your code or INTEGER if you are only using whole numbers. If you are using a large INTEGER try using

integer(kind=int64) :: my_variable

which will allow you to use up to
9223372036854775807

Bill
Lead Application Developer
New York State, USA

 

[DomBernard]

To close this question I want to add that it is not exactly a problem of range (you can use numbers larger than 400**3 with real 32b), but a problem of discontinuity in the number representation. By step of one I reached the limit of 2**24 and adding one was not enough to change the exponent so the representation in 32b remained unchanged.
For compting Integer is the best choice because the granularity of the number representation is constant and equal to one. For large number the use of double precision Integer is the solution.
Thanks again.