Wrong result of logical expression evaluation

[Josué_B_]

Hello,

I'm facing a very bizarre problem: a logical expression of the form (A .and. B) .or. (C .and. D) has to be evaluated in a do while construct. But if gives the wrong result. I mean: when printing the expressions one by one, they are correct. When printing the left part, and the right part, they are also correct. When the full expression is used, its wrong.

A, B, C, D are : T, F, T, F

(A .and. B) gives F

(C .and. D) gives F

BUT: (A .and. B) .or. (C .and. D) gives T !!!!

Here is a copy of the code:

100   IF ( (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
     &     .OR. (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO) ) THEN

      print*,'*** start of time loop**'
      print*,'T, Tend',T,Tend
      print*,'test result',(Direction.GT.0).AND.((T-Tend)+Roundoff.LE.
     & ZERO) .OR. (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO) 

! (list of statements, which increase T towards Tend)

      print'(a,2f25.18)','T, Tend',T,Tend
      print*,'Roundoff =',Roundoff
      AA = (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
      BB = (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO)

      print*,'test result',
     & (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
     & .OR.
     & (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO)

      print*,'test with previous evalution',AA.OR.BB
      print*,'each part A,B (for A .OR. B test)',AA,BB
      print*,""
      print*,'subtest p1',
     & (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
      print*,'subtest p2',
     & (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO)

      print*,'split 1',(Direction.GT.0),((T-Tend)+Roundoff.LE.ZERO)
      print*,'split 2',(Direction.LT.0),((Tend-T)+Roundoff.LE.ZERO)

      print*,'*** END of time loop**'

      GOTO 100

      END IF

and here is the output:

 *** start of time loop**
 T, Tend   19.3310000000000        20.0000000000000     
 test result T
T, Tend    20.000000000000000000    20.000000000000000000
 Roundoff =  2.220446049250313E-016
 test result T
 test with previous evalution T
 each part A,B (for A .OR. B test) T F
 
 subtest p1 F
 subtest p2 F
 split 1 T F
 split 2 F T
 *** END of time loop**
 *** start of time loop**
 T, Tend   20.0000000000000        20.0000000000000     
 test result T
forrtl: error (73): floating divide by zero

Look: when evaluation first the logical expressions, and storing them in logical variables AA and BB, they are said to be T (which is wrong!) and F (which is correct).

Now, when evaluated inside the print*, command (see "subtest p1/2"), they are F and F (which is now correct).

How is that possible ?

Note that "Roundoff" is defined as epsilon(1.d0).

The ifort version used is 9.0.

Many thanks for your help.

[mecej4]

I do not have a complete picture of what your program does, but there is one part of it that is flawed. The expression (T-Tend)+Roundoff.LE.ZERO will probably not work as you seem to expect. Machine epsilon is defined relative to 1.0, not 20 or some such number much larger than 1.0.

The following program should illustrate the point.

program tst
implicit none
double precision :: T,Tend,dT,expr
!
Tend=20d0
dT=20d0*epsilon(1d0)
do
   T=Tend-dT
   expr=(Tend-T)+epsilon(1d0)
   write(*,10)T,dT,expr
   if(expr.le.0d0)exit
   dT=0.5d0*dT
end do
10 format(1p,2D22.15,4x,D11.3)
end program

[Josué_B_]

Thanks a lot for your answer.

I'm affraid I don't understand everything:

- firstly, I'm concerned by the output of your test program, which gives (on my machine):

 2.000000000000000D+01 4.440892098500626D-15      3.553D-15
 2.000000000000000D+01 2.220446049250313D-15      3.553D-15
 2.000000000000000D+01 1.110223024625157D-15      0.000D+00

Why is the third value smaller than the second one?

In your example, expr = (Tend-T)+epsilon(1d0) and T=Tend-dT. Thus expr could be rewritten as dT+epsilon(1d0).

So expr (the third printed value) should be larger than dT (the second printed value), right?

- secondly, I don't understand why you say I'm comparing epsilon with 20. In fact, epsilon is rather compared with (Tend-T), which tends towards 0. Thus, I would rather say that epsilon is compared to something close to zero. BUT, as far as I understood, the machine dependant value returned by the intrinsic function "epsilon" was a value which could be distinguished from any real, when summed or substracted to this real (basically, R + epsilon /= R whatever R). Did I misunderstood the objective of this epsilon function ?

- last, there is still a remaining point which appears unbelievable: in my first message, why are the outputs of print*,(...) AA (line 22) and print*,'split 1' (line 29) different, while they refer to the same logical instruction!?!? The first one is evaluated as true, the second one is evaluated as false.

Many thanks

[mecej4]

Josué B. wrote:

... "epsilon" was a value which could be distinguished from any real, when summed or subtracted to this real (basically, R + epsilon /= R whatever R).

That is incorrect. One definition of machine epsilon is "the smallest number that, when added to 1.0, yields a result different from 1.0". Because of finite precision, while 1 + ε differs from 1, 2 + ε is equal to 2. To firm up these notions, you may wish to write some test programs, and look up the related intrinsic functions NEAREST() and SPACING(). See also the Wikipedia article https://en.wikipedia.org/wiki/Machine_epsilon .

Regarding your question concerning logical expressions, it would be best if you presented a complete program that involves the question. I could then address the question in the specific context, instead of speculating.

[jimdempseyatthecove]

      AA = (Direction.GT.0).AND.((T-Tend)+Roundoff.LE.ZERO)
      BB = (Direction.LT.0).AND.((Tend-T)+Roundoff.LE.ZERO)

In looking at the above two statements, when the absolute value of the difference of T and Tend is less than the precision capability of the floating point format holding Roundoff, then

((T-Tend)+Roundoff) == .((Tend-T)+Roundoff) == Roundoff

The above would be computationally true even though T != Tend

IOW the fractional difference cannot be held within the floating point precision of Roundoff.

Jim Dempsey

[Josué_B_]

mecej4 wrote:
That is incorrect. One definition of machine epsilon is "the smallest number that, when added to 1.0, yields a result different from 1.0". Because of finite precision, while 1 + ε differs from 1, 2 + ε is equal to 2. To firm up these notions, you may wish to write some test programs, and look up the related intrinsic functions NEAREST() and SPACING(). See also the Wikipedia article https://en.wikipedia.org/wiki/Machine_epsilon .

Thanks, I will have a look to these functions

mecej4 wrote:
Regarding your question concerning logical expressions, it would be best if you presented a complete program that involves the question. I could then address the question in the specific context, instead of speculating.

Please see the attached file, 'rosenbrock.f'. It is part of KPP-2.2.3 distribution (KPP stands for Kinetic PreProcessus).

I have discovered something using your test program. First, I modified it cause there was a small mistake in it : when defining T as T=Tend-dT, this means that T is smaller than Tend and approach Tend with lower values (this is indeed my case). However, to keep consistent with the formulation of "expr", it should be written expr = (T-Tend)+epsilon, and exit the do loop for values .ge.0. The modified code that I use is as follows:

program tst
implicit none
double precision :: T,Tend,dT,expr
!
Tend=20d0
dT=20d0*epsilon(1d0)
do
   T=Tend-dT
   expr=(T-Tend)+epsilon(1d0)
   write(*,10)T,dT,expr
   if(expr.ge.0d0)exit
   dT=0.5d0*dT
end do
10 format(1p,2D22.15,4x,D11.3)
end program

Now, I can see different results depending on the compiler used.

When using ifort:

 2.000000000000000D+01 4.440892098500626D-15     -3.553D-15
 2.000000000000000D+01 2.220446049250313D-15     -3.553D-15
 2.000000000000000D+01 1.110223024625157D-15      0.000D+00

When using gfortran:

 2.000000000000000D+01 4.440892098500626D-15     -3.331D-15
 2.000000000000000D+01 2.220446049250313D-15     -3.331D-15
 2.000000000000000D+01 1.110223024625157D-15      2.220D-16

This comes as another proof that something is going wrong with ifort: the 3rd value of the 3rd line should be close to epsilon, which is 2.22d-16 on my machine.

(To be accurate: in the modified program tst, I should have used a ".gt." test instead of ".ge", the first one being obviously the negation of ".le.". But because of the strange behaviour of the program when compiled with ifort, this leads to infinite do loop (and that's my problem, even if formulated in another way in the original program "rosenbrock.f")

[Josué_B_]

jimdempseyatthecove wrote:

In looking at the above two statements, when the absolute value of the difference of T and Tend is less than the precision capability of the floating point format holding Roundoff, then

((T-Tend)+Roundoff) == .((Tend-T)+Roundoff) == Roundoff

The above would be computationally true even though T != Tend

IOW the fractional difference cannot be held within the floating point precision of Roundoff.

Thank you for your answer.

I agree. But my problem is that the result of the calculation (which should be ~Roundoff), is not! See my previous answer to mecej4, and the output comparing the same program output, when compiled either with ifort and with gfortran.

To go back to the original file (rosenbrock.f), the problem is that after some integration steps, T actually becomes equal to Tend. BUT the result of the logical expression (T-Tend)+Roundoff.LE.ZERO stays true, while it should be false.

[mecej4]

There are some INCLUDE statements in rosenbrock.f. The included files are needed in order to compile the program.

[Josué_B_]

Hello,

I made progress regarding this problem. In the test program suggested by mecej4, I simply split the expr calculation in 2 steps:

! original
! expr=(T-Tend)+epsilon(1d0)
!
expr = (T-Tend)
expr = expr + epsilon(1d0)

Now, the results with ifort are the same than those I get with gfortran, and the one I expect:

 2.000000000000000D+01 4.440892098500626D-15     -3.331D-15
 2.000000000000000D+01 2.220446049250313D-15     -3.331D-15
 2.000000000000000D+01 1.110223024625157D-15      2.220D-16

And thanks to the explanations of mecej4 about the spacing function, I can understand the reason of the bug:

in the original expression (expr=(T-Tend)+epsilon(1d0)), the program doesn't care of the parenthesis, and computes first Tend+epsilon(1d0). Since Tend is greater than 1, epsilon(1d0) is smaller than spacing(Tend), and the result of the calculation is Tend+epsilon(1d0) = Tend.

Then, T-Tend gives 0 when dT becomes smaller than spacing(Tend).

My next (and maybe the last) question will be: is there a way to force ifort to care of the parenthesis evenif they are not mandatory?

Many thanks

[Kevin_D_Intel]

Yes, see -fprotect-parens

[Josué_B_]

Thanks for your answer, Kevin D.

I think that this option is not available with the old version of ifort I use.

I fount another option: -fltconsistency which looks similar, and allowed the model to run for a few integration steps (it was crashing at the first step previously) but crashes after a while.

So I rephrased the initial instruction:

DO WHILE (ABS(Tend-T) .LE. Roundoff)

Like this I'm sure that Tend-T is evaluated separately as I expect, without needing a further compilation option. And the model works well like this.

Cheers

[TimP]

I hope you aren't using such an old version of ifort that /assume:protect_parens isn't supported.

/standard-semantics is documented as the option for standards compliance (including observance of parentheses).  It does appear to have a few side effects which aren't fully documented, varying somewhat among compiler updates.  Without that option, several of these behaviors are supposed to be compatible with predecessor compilers prior to f95, including violation of parentheses which seems to have been intended to help optimization for x87/ia32 mode.
 

[Josué_B_]

Tim P. wrote:

I hope you aren't using such an old version of ifort that /assume:protect_parens isn't supported.

I'm afraid I do use a version of ifort which doesn't support this option. The version I use is ifort 9.0

I tried the -fprotect-parens option suggested by Kevin D., and also the -fp-model precise which was mentioned in the manual page refered to by Kevin D.: https://software.intel.com/en-us/node/579289

None of these options are working with ifort 9.0

[Kevin_D_Intel]

ifort 9.0 is about 10+ years old now and lots of new features/options/bug fix content has gone into the compiler since then. You should consider updating. You can keep your current ifort 9.0 installed and install a separate instance of our latest PSXE 2016 assuming you have a valid license to permit that.

I don't know whether this applies to your situation, https://software.intel.com/en-us/qualify-for-free-software