CODEWITHSUNDEEP
fortrando-loops
Roger
Roger

Reputation: 21

Fortran loop not working

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?

Upvotes: 2

Views: 418

Answers (1)

Kai Guther
Kai Guther

Reputation: 793

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.

Upvotes: 3

Related Questions