Does Fortran have inherent limitations on numerica

Question

While working on a simple programming exercise, I produced a while loop (DO loop in Fortran) that was meant to exit when a real variable had reached a precise value.

I noticed that due to the precision being used, the equality was never met and the loop became infinite. This is, of course, not unheard of and one is advised that, rather than comparing two numbers for equality, it is best see if the absolute difference between two numbers is less than a set threshold.

What I found disappointing was how low I had to set this threshold, even with variables at double precision, for my loop to exit properly. Furthermore, when I rewrote a "distilled" version of this loop in Perl, I had no problems with numerical accuracy and the loop exited fine.

Since the code to produce the problem is so small, in both Perl and Fortran, I'd like to reproduce it here in case I am glossing over an important detail:

Fortran Code

PROGRAM precision_test
IMPLICIT NONE

! Data Dictionary
INTEGER :: count = 0 ! Number of times the loop has iterated
REAL(KIND=8) :: velocity
REAL(KIND=8), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0

velocity = 0.5 * MACH_2_METERS_PER_SEC ! Initial Velocity
DO
        WRITE (*, 300) velocity
        300 FORMAT (F20.8)
        IF (count == 50) EXIT
        IF (velocity == 5.0 * MACH_2_METERS_PER_SEC) EXIT
!       IF (abs(velocity - (5.0 * MACH_2_METERS_PER_SEC)) < 1E-4) EXIT
        velocity = velocity + 0.1 * MACH_2_METERS_PER_SEC
        count = count + 1 
END DO

END PROGRAM precision_test

Perl Code

#! /usr/bin/perl -w
use strict;

my $mach_2_meters_per_sec = 340.0;

my $velocity = 0.5 * $mach_2_meters_per_sec;

while (1) {
        printf "%20.8f\n", $velocity;   
        exit if ($velocity == 5.0 * $mach_2_meters_per_sec);
        $velocity = $velocity + 0.1 * $mach_2_meters_per_sec;
}

The commented-out line in Fortran is what I would need to use for the loop to exit normally. Notice that the threshold is set to 1E-4, which I feel is quite pathetic.

The names of the variables come from the self-study-based programming exercise I was performing and don't have any relevance.

The intent is that the loop stops when the velocity variable reaches 1700.

Here are the truncated outputs:

Perl Output

    170.00000000
    204.00000000
    238.00000000
    272.00000000
    306.00000000
    340.00000000

...

   1564.00000000
   1598.00000000
   1632.00000000
   1666.00000000
   1700.00000000

Fortran Output

    170.00000000
    204.00000051
    238.00000101
    272.00000152
    306.00000203
    340.00000253

...

   1564.00002077
   1598.00002128
   1632.00002179
   1666.00002229
   1700.00002280

What good is Fortran's speed and ease of parallelization if its accuracy stinks? Reminds me of the three ways to do things:

1. The Right Way

2. The Wrong Way

3. The Max Power Way

"Isn't that just the wrong way?"

"Yeah! But faster!"

All kidding aside, I must be doing something wrong.

Does Fortran have inherent limitations on numerical accuracy compared to other languages, or am I (quite likely) the one at fault?

My compiler is gfortran (gcc version 4.1.2), Perl v5.12.1, on a Dual Core AMD Opteron @ 1 GHZ.

Answer 1

As already answered, "plain" floating point constants in Fortran will default to the default real type, which will likely be single-precision. This is an almost classic mistake.

Also, using "kind=8" is not portable -- it will give you double precision with gfortran, but not with some other compilers. The safe, portable way to specify precisions for both variables and constants in Fortran >= 90 is to use the intrinsic functions, and request the precision that you need. Then specify "kinds" on the constants where precision is important. A convenient method is to define your own symbols. For example:

integer, parameter :: DR_K = selected_real_kind (14)

REAL(DR_K), PARAMETER :: MACH_2_METERS_PER_SEC = 340.0_DR_K

real (DR_K) :: mass, velocity, energy

energy = 0.5_DR_K * mass * velocity**2

This can also be important for integers, e.g., if large values are needed. For related questions for integers, see Fortran: integer*4 vs integer(4) vs integer(kind=4) and Long ints in Fortran

Answer 2

Your assignment is accidentally converting the value to single precision and then back to double.

Try making your 0.1 * be 0.1D0 * and you should see your problem fixed.