Yes, I realise that; thank you. Your first comment led me to question what I had interpreted the line to state and have solved that one (and similar).
I can’t get the first one to clear, though. I tried this to no avail:
if (z00 == 0d0 .or. z01 == 0d0 .or. z11 == 0d0 .or. z10 == 0d0) then
z00 et al are declared as double precision and they all derive their individual values via another variable, also declared as double precision. However, this latter one (zkap()) appears in a read, where the data type being read in might well be defined as a real?
This is there as a reminder that exact equality of floating point values is very rarely true, and you should make sure you understand that if you’re going to do it. For example, what do you think the following program will print? T or F?
program floating_point_math
print *, 0.1 + 0.2 == 0.3
end program
Due to binary representation, and assuming these are 8-bit real values, it’d come out as false.
My issue is that I am dealing with a large program that I did not originally write. It was written in F77 and I have, slowly, moved it across to F90. However, there are places where I now find the logic runs the risk of failing, yet the code works - for very large sequences (50000 steps, taking a day or two to process).
What those kinds of comparisons generally end up doing is serving as a proxy for “Did some other branch in the code get taken”. I.e.
program roundabout
implicit none
real :: x
real, parameter :: magic_value = 42.0
x = magic_value
if (the_stars_align) then
x = something_else()
end if
if (x == magic_value) then
print *, "The stars didn't align"
else
print *, "The stars aligned"
end if
end program
because there’s virtually zero chance that something_else will return the magic_value. Most projects just end up defaulting to the convention that the magic value happens to be zero, because they didn’t really think about it (or at least in that way). Now if you see something like
program has_a_bug
if (sum(things) == 100.0) then
...
end if
end program
It’s pretty much guaranteed to be a bug, because again, the chances of that summation ever being exactly 100 are pretty much nil.
This is the recommended approach, rather than adding D exponents. If you do this with a parameter constant, then when you later decide to use one of the other supported KINDs, it would just involve changing one parameter value, rather than numerous D exponents scattered throughout the code.
Just that you can define something like dp and use if as the suffix for all the constants, and then
if you wanted to change to another kind all you have to do is change the definition of dp (of course that is oversimplifying as your computations could call a function that only takes certain kinds, etc …). As @kargl noted comparing to zero is generally OK and so on and that the compiler is being “too helpful”. The second example IF below would probably remove your messages, but (probably) would have performance hit compared to your original if called intensively, but it shows everything done in terms of the kind “dp” which is generally a good idea.
integer,parameter :: dp=kind(0.0d0)
real(kind=dp) :: z00, z01, z02, z03, tol=2*epsilon(0.0_dp)
z01=0.1_dp
z02=0.2_dp
z03=0.3_dp
z00=z01+z02
if (z00 == z03 .or. any([z01,z02] == 0.0_dp) )then
write(*,*)'OK',z00,z03,z00-z03,epsilon(0.0_dp)
else
write(*,*)'BAD',z00,z03,z00-z03,epsilon(0.0_dp)
endif
! use of functions could be slower, but should not produce warnings
if (abs(z00 - z03) < tol .or. any(abs([z01,z02]) < tol) )then
write(*,*)'OK NOW',z00,z03,z00-z03,tol
else
write(*,*)'BAD NOW',z00,z03,z00-z03,tol
endif
end
not saying this is a great solution, it is just for demonstration purposes.
Thank you for the clarity. Curious that you have defined dp as an integer, yet then appear to define it in terms of a double?
However, more relevant is that I think I need to find out what the code is trying to achieve and remove these lines altogether. The model already takes considerable time to run and I am trying to make it more efficient. Hahahaha!
All kinds are integer values. The KIND() function returns the integer value representing the kind possessed by the argument. The value of the argument is irrelevant; it is just returning the kind for that type of argument. You may see them by predefined names like INT8 or REAL64 but those are just integer constants.
In floating point computation, the problem you have with comparing reals is how small does the value need to become to be considered zero ?
You might consider replacing the test by using a logical function IS_ZERO. This could be used as:
if ( is_zero (z00) .or. is_zero (z01) .or. is_zero (z11) .or. is_zero (z10) ) then
where:
logical function IS_ZERO ( value )
real(8) :: value
real(8) :: typical_value = 1.0
is_zero = abs ( value ) <= tiny ( typical_value )
end function IS_ZERO
The advantage of this test is that it removes all the compiler warnings.
The problem remains : what is a typical value ?
In reality, most often, we have initialised a real variable ( as Z00 = 0.) and then later we simply want to test if it has been changed. Zero is a significant number for any real kind, so the z00 == 0. is a safe test.
Perhaps a change could be to define “real(8) :: zero = 0”
Z00 = zero
… ( computation)
then use “if ( Z00 == zero ) …”, Although this would produce more compiler warnings !!
Be aware that most real numbers are not exact.
! Try the following in gfortran
real(4) :: one_4 = 0.1_4
real(8) :: one_8 = 0.1_8
real(8) :: one_x
one_x = one_4
if ( one_4 /= one_8 ) write (*,*) "one_4 is not one_8", one_4, one_8, one_x
if ( one_4 == one_x ) write (*,*) "one_4 is one_x", one_4, one_x, one_8
if ( one_4 == one_8 ) write (*,*) "unexpected one_4 is one_8"
if ( one_4 /= one_x ) write (*,*) "unexpected one_4 is not one_x"
end
Then change to the following and see the result. ( this would not have happened in most F77 compilers, eg LF77 )
real(4) :: one_4 = 0.1
real(8) :: one_8 = 0.1