2

Fortran returns an error message saying floating point overflow when i run this

program single
implicit none
integer :: n,k
real, dimension(255) :: x
n = 255
x(1) = 1/3.0
x(2) = 1/12.0
do k = 3,n
  x(k) = 2.25*x(k-1) - 0.5*x(k-2)
end do
print *,x
end program

After some testing, I think it is due to the variable index in the do loop but I don't quite understand what the problem is. Can anyone help me?

Roger
  • 21
  • 1

1 Answers1

3

By printing the values of x during the loop, you can get some insight to what is going on. The first few values are

2.08333284E-02
5.20832092E-03
1.30205788E-03
3.25469766E-04
8.12780345E-05
2.01406947E-05
4.67754580E-06
4.54130713E-07
-1.31697880E-06
-3.19026776E-06
-6.51961273E-06
-1.30739954E-05

Note how the sign changes and the value keeps increasing monotonically from there. You might have noticed that for negative values of x, the sequence you defined is growing exponentially, so once you go negative, the exponential growth quickly gets you and leads to a floating point overflow, as x just grows and grows (after 100 iterations its already beyond 1E+21).

However, let's see what happens if we increase the floating point precision, by specifying

Integer, Parameter :: wp = Selected_real_kind( 13, 70 )
real(kind=wp), dimension(255) :: x

If we compile now and run the very same program, we will notice that x(k) goes to 0 for increasing k, but does not change sign within 255 iterations. So in fact, the overflow you see is caused by insufficient accuracy of the floating point representation which leads to an erroneous sign change, and subsequent exponential growth.

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Kai Guther
  • 783
  • 3
  • 11
  • 5
    Please don't teach `(kind=8)` to beginners, its non-standard and may not be portable. Other than that, good answer. – Ross Jul 31 '18 at 22:48
  • An interesting difference I have found is 32-bit and 64-bit compiler results vary for the same selected kind. Not sure if it is the compiler I am using (FTN95) or different round-off between 32-bit and 64-bit instructions. 64-bit FTN95 and gFortran produce the same results. I also used "x(k) = (4.5*x(k-1) - x(k-2))/2" to see if this changed the round-off, but no change. Changing to 8-byte reals removed the -ve result, but 32-bit /64-bit round-off difference still evident. – johncampbell Aug 05 '18 at 13:35