INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Jobs

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

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

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

(OP)
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:

[uxmal:Prog_DB/Geom/src] bernard% gfortran Simple_Test.FOR
[uxmal:Prog_DB/Geom/src] bernard% ./a.out
16777217 16777216. 16777216.

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

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

This is of course the logical one.

Any explanation?


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

IVF1 is not initialized

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Resources

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close