Why does temporary variable give the wrong results?

[Bobby_G_]

I'm working with this code (only for example):

Example 1:

DO I=1, TFACE1

    Dtmp2=H(ELFACE2(I))-H(ELFACE1(I))                           
    Dtmp1=DHX(ELFACE2(I))*RPQX(I)+DHY(ELFACE2(I))*RPQY(I)-Dtmp2   

    !Computing at point D 
    IF(Dtmp1*Dtmp2>0.0) THEN

       HRD(I)=H(ELFACE2(I))-RDQPQ(I)*&
              ((Dtmp1**2.0+E)*Dtmp2+(Dtmp2**2.0+E)*Dtmp1)/(Dtmp1**2.0+Dtmp2**2.0+2.0*E)

    ELSE

       HRD(I)=H(ELFACE2(I))

    END IF

    Dtmp1=DHX(ELFACE1(I))*RPQX(I)+DHY(ELFACE1(I))*RPQY(I)-Dtmp2

    !Computing at point D 
    IF(Dtmp1*Dtmp2>0.0) THEN

       HLD(I)=H(ELFACE1(I))+RPDPQ(I)*&
              ((Dtmp1**2.0+E)*Dtmp2+(Dtmp2**2.0+E)*Dtmp1)/(Dtmp1**2.0+Dtmp2**2.0+2.0*E)

    ELSE

       HLD(I)=H(ELFACE1(I))

    END IF


    Dtmp2=ELV(ELFACE2(I))-ELV(ELFACE1(I))                         
    Dtmp1=DELVX(ELFACE2(I))*RPQX(I)+DELVY(ELFACE2(I))*RPQY(I)-Dtmp2

    !Computing at point D 
    IF(Dtmp1*Dtmp2>0.0) THEN

       ELVRD(I)=ELV(ELFACE2(I))-RDQPQ(I)*&
              ((Dtmp1**2.0+E)*Dtmp2+(Dtmp2**2.0+E)*Dtmp1)/(Dtmp1**2.0+Dtmp2**2.0+2.0*E)

    ELSE

       ELVRD(I)=ELV(ELFACE2(I))

    END IF

    Dtmp1=DELVX(ELFACE1(I))*RPQX(I)+DELVY(ELFACE1(I))*RPQY(I)-Dtmp2

    !Computing at point D 
    IF(Dtmp1*Dtmp2>0.0) THEN

       ELVLD(I)=ELV(ELFACE1(I))+RPDPQ(I)*&
              ((Dtmp1**2.0+E)*Dtmp2+(Dtmp2**2.0+E)*Dtmp1)/(Dtmp1**2.0+Dtmp2**2.0+2.0*E)

    ELSE

       ELVLD(I)=ELV(ELFACE1(I))

    END IF


  !Computing at point M
    HR(I)=HRD(I)+DHX(ELFACE2(I))*RDMX(I)+DHY(ELFACE2(I))*RDMY(I)
    HL(I)=HLD(I)+DHX(ELFACE1(I))*RDMX(I)+DHY(ELFACE1(I))*RDMY(I)

    ELVR(I)=ELVRD(I)+DELVX(ELFACE2(I))*RDMX(I)+DELVY(ELFACE2(I))*RDMY(I)
    ELVL(I)=ELVLD(I)+DELVX(ELFACE1(I))*RDMX(I)+DELVY(ELFACE1(I))*RDMY(I)

 END DO

Example 2:
 

DO I=1, TFACE1

    FC1=ELFACE1(I)   
    FC2=ELFACE2(I)

    DHC=H(FC2)-H(FC1)   
    DELVC=ELV(FC2)-ELV(FC1)     

    DHUPWR=DHX(FC2)*RPQX(I)+DHY(FC2)*RPQY(I)-DHC 
    DEUPWR=DELVX(FC2)*RPQX(I)+DELVY(FC2)*RPQY(I)-DELVC  

    DHUPWL=DHX(FC1)*RPQX(I)+DHY(FC1)*RPQY(I)-DHC          
    DEUPWL=DELVX(FC1)*RPQX(I)+DELVY(FC1)*RPQY(I)-DELVC    

   !Computing the function
    AR=DHUPWR
    AL=DHUPWL
    B =DHC

      IF(AR*B>0.0) THEN

         LVAHR=((AR**2+E)*B+(B**2+E)*AR)/(AR**2+B**2+2*E)

      ELSE

         LVAHR=0.0

      END IF

      IF(AL*B>0.0) THEN

         LVAHL=((AL**2+E)*B+(B**2+E)*AL)/(AL**2+B**2+2*E)

      ELSE

         LVAHL=0.0

      END IF


    AR=DEUPWR
    AL=DEUPWL
    B =DELVC

      IF(AR*B>0.0) THEN

         LVAER=((AR**2+E)*B+(B**2+E)*AR)/(AR**2+B**2+2*E)

      ELSE

         LVAER=0.0

      END IF


      IF(AL*B>0.0) THEN

         LVAEL=((AL**2+E)*B+(B**2+E)*AL)/(AL**2+B**2+2*E)

      ELSE

         LVAEL=0.0

      END IF


  !Computing at point D        
    HRD(I)=H(FC2)-RDQPQ(I)*LVAHR
    HLD(I)=H(FC1)+RPDPQ(I)*LVAHL

    ELVRD(I)=ELV(FC2)-RDQPQ(I)*LVAER
    ELVLD(I)=ELV(FC1)+RPDPQ(I)*LVAEL

  !Computing at point M
    HR(I)=HRD(I)+DHX(FC2)*RDMX(I)+DHY(FC2)*RDMY(I)
    HL(I)=HLD(I)+DHX(FC1)*RDMX(I)+DHY(FC1)*RDMY(I)

    ELVR(I)=ELVRD(I)+DELVX(FC2)*RDMX(I)+DELVY(FC2)*RDMY(I)
    ELVL(I)=ELVLD(I)+DELVX(FC1)*RDMX(I)+DELVY(FC1)*RDMY(I)

 END DO

FYI, these are only some parts of our code. I run the program with -O3 flag, and somehow I found that Example 1 gave the infinity results, while Example 2 gave the correct one, even I think (to my knowledge) at first glance they are same.

Does anyone know why it is different?

PS: if it's only just because I wrote any parts wrongly, please also let me know.

[mecej4]

You are asking a question that is tough to answer, because many important details are missing.

You show two code fragments, each about 80 lines long, both of which contain calls to external functions whose properties you have not informed us about: for example, are the functions pure? Since you did not show any declarations, are we to assume that all variables and functions are implicitly typed? Some variables appear to have been renamed between the two versions, and that imposes an extra burden on the person doing a visual comparison of the two versions. 

You said that one version gave infinite 'results'. Which of the numerous variables in the code are the 'results'? Do you obtain infinite values always, regardless of the values of those variables that should have been defined prior to executing the code fragment?

Compiler version, etc.?

Having a complete test case (that can be compiled and run) would facilitate understanding the problem.

[Bobby_G_]

@mecej4: Thanks for the answer. I'm sorry if the code is too bad to see. The final variable results are not written on these codes. I can't post here all parts of codes, but the important thing is that if I used code on example 1 ( keeping all other parts of the codes same when compiling using example 2), the final variable results were infinity, while example 2 gave the correct results.

Example 1: The variable declarations are:

INTEGER::NEM,TFACE1,TFACE3W

REAL::ERROR,Dtmp1,Dtmp2

INTEGER, DIMENSION(TFACE3W)::ELFACE1,ELFACE2

INTEGER, DIMENSION(TFACE1)::LMIN,RMIN,INDEKS5

REAL, DIMENSION(NEM)::DHX,DHY,DELVX,DELVY,DUX,DUY,DVX,DVY,H,ELV

REAL, DIMENSION(TFACE3W)::RPQX,RPQY,RDMX,RDMY,RPDPQ,RDQPQ,HR,HL,ELVR,ELVL,UR,UL,VR,VL

REAL, DIMENSION(TFACE1)::HRD,HLD,ELVRD,ELVLD,URD,ULD,VRD,VLD

Example 2: The variable declarations are:

INTEGER::NEM,TFACE1,TFACE3W,FC1,FC2

REAL::DHC,DELVC,B,E,DHUPWL,DEUPWL,AL,LVAHL,LVAEL,BR,BL,DHUPWR,DEUPWR,AR,LVAHR,LVAER

INTEGER, DIMENSION(TFACE3W)::ELFACE1,ELFACE2

INTEGER, DIMENSION(TFACE1)::LMIN,RMIN,INDEKS5

REAL, DIMENSION(NEM)::DHX,DHY,DELVX,DELVY,DUX,DUY,DVX,DVY,H,ELV

REAL, DIMENSION(TFACE3W)::RPQX,RPQY,RDMX,RDMY,RPDPQ,RDQPQ,HR,HL,ELVR,ELVL,UR,UL,VR,VL

REAL, DIMENSION(TFACE1)::HRD,HLD,ELVRD,ELVLD,URD,ULD,VRD,VLD

In this case the parameters are:

NEM     = 8800
TFACE1  = 12500
TFACE3W = 14000

ERROR   = 10E-16
E       = 10E-16

[jimdempseyatthecove]

I suggest you go back to the code that works, then incrementally introduce your "improvements". Run test at each step. This will inform you when the error was introduced.

Jim Dempsey

[Bobby_G_]

@Jim: Thanks for your suggestion. I've just done that. I did step by step and so when I changed (for example) this one:

new=old * ((AR**2+ERROR)*B+(B**2+ERROR)*AR)/(AR**2+B**2+2*ERROR)

with this one:

LV= ((AR**2+ERROR)*B+(B**2+ERROR)*AR)/(AR**2+B**2+2*ERROR)
new=old * LV

then it worked well.

Do you have any recommendations? Thanks in advance.

[jimdempseyatthecove]

The this seems to indicate that variable LV is used after that statement (iow it is not set in first example).

Also note that if LV is local but SAVE, its prior value may be required upon next call.

Jim Dempsey

[Bobby_G_]

I wrote again some parts of a code and somehow it looks like the following (for simplification)

INTEGER::J,K,TFACE1,TFACE3W,NEM
REAL::Dtmp1,Dtmp2
INTEGER, DIMENSION(TFACE3W)::ELFACE1,ELFACE2
REAL, DIMENSION(TFACE3W)::RDQPQ
REAL, DIMENSION(TFACE1)::HRD
REAL, DIMENSION(NEM)::H

ERROR=10E-16

 DO I=1, TFACE1

    J=ELFACE2(I)
    K=ELFACE1(I)

    Dtmp1=H(J)-H(K)
    Dtmp2=DHX(J)*RPQX(I)+DHY(J)*RPQY(I)-DHC

      IF(Dtmp1*Dtmp2>0.0) THEN

         LV=((Dtmp1**2+ERROR)*Dtmp2+(Dtmp2**2+ERROR)*Dtmp1)/(Dtmp1**2+Dtmp2**2+2*ERROR) 
         HRD(I)=H(J)-RDQPQ(I)*LV

      ELSE

         HRD(I)=H(J)

      END IF

 END DO

The problem is,

1. When I defined the variable LV as REAL(KIND=4) or REAL(KIND=8) or REAL(KIND=16), the program gave the infinity results.
2. But, when I didn't define the variable LV at all, the program gave the correct result. I assume it was because of ERROR = 10E-16, which is closely to machine accuracy, so the computer just wasn't be able to calculate it. But, why did it give me the correct results when I didn't define the variable at all?

3. I have also tried the direct multiply like this:

   HRD(I)=H(J)-RDQPQ(I)*((Dtmp1**2+ERROR)*Dtmp2+(Dtmp2**2+ERROR)*Dtmp1)/(Dtmp1**2+Dtmp2**2+2*ERROR)

But also it gave the infinity results. Well, for this case, it could be concluded since I declared variable Dtmp1 and Dtmp2 as REAL before.

[jimdempseyatthecove]

When you do not define LV
... .AND. ...
you are not using IMPLICIT NONE
... THEN ...
LV is of implicit type INTEGER

This will truncate the REAL expression to integer.

When you switched to the locally defined REAL, there is no truncation. So the value differs (when fractional component present in expresson).

Your former convergence code was likely wrong (unless you really need the truncation).

Jim Dempsey