fixes for flang bug on arithmetic if flang-compiler/flang#844 #179

src/NWints/hondo/hnd_rt123.F

@@ -14,6 +14,7 @@ subroutine hnd_rt123
       equivalence (w9(5),ww5)
       double precision r12, pie4, r22, w22, r13, r23, w23, r33, w33
       double precision y, f1, f2, t1, t2, t3, a1, a2, g
+      integer*4 n4
       data r12,pie4/2.75255128608411d-01, 7.85398163397448d-01/
       data r22,w22/ 2.72474487139158d+00, 9.17517095361369d-02/
       data r13/     1.90163509193487d-01/
@@ -23,7 +24,12 @@ subroutine hnd_rt123
       if(yy.gt.1.0d+00) go to 30
       if(yy.gt.3.0d-07) go to 20
 c     yy is approximately zero.         nroots=1,2, or 3
+#if 1
+      n4=nroots-2
+      if(n4) 11,12,13
+#else
       if(nroots-2) 11,12,13
+#endif
    11 rt1= 0.5d+00 -yy/5.0d+00
       ww1= 1.0d+00 -yy/3.0d+00
       return
@@ -203,7 +209,12 @@ subroutine hnd_rt123
      2+2.4645596956002d-01)/yy-4.9984072848436d-01)/yy
      3-3.1501078774085d-06)*g + dsqrt(pie4/yy)
       f1=(ww1-g)/(yy+yy)
+#if 1
+      n4=nroots-2
+      if(n4) 51,52,53
+#else
       if(nroots-2) 51,52,53
+#endif
    51 rt1=f1/(ww1-f1)
       return
    52 y=yy-7.5d+00
@@ -247,7 +258,12 @@ subroutine hnd_rt123
    60 ww1=       (((-1.8784686463512d-01/yy+2.2991849164985d-01)/yy
      1-4.9893752514047d-01)/yy-2.1916512131607d-05)*g + dsqrt(pie4/yy)
       f1=(ww1-g)/(yy+yy)
+#if 1 
+      n4=nroots-2
+      if(n4) 61,62,63
+#else
       if(nroots-2) 61,62,63
+#endif
    61 rt1=f1/(ww1-f1)
       return
    62 rt1=      ((((-1.01041157064226d-05*yy+1.19483054115173d-03)*yy
@@ -289,7 +305,12 @@ subroutine hnd_rt123
       ww1=        (( 1.9623264149430d-01/yy-4.9695241464490d-01)/yy
      1-6.0156581186481d-05)*g + dsqrt(pie4/yy)
       f1=(ww1-g)/(yy+yy)
+#if 1
+      n4=nroots-2
+      if(n4) 71,72,73
+#else
       if(nroots-2) 71,72,73
+#endif
    71 rt1=f1/(ww1-f1)
       return
    72 rt1=      ((((-1.14906395546354d-06*yy+1.76003409708332d-04)*yy
@@ -331,7 +352,12 @@ subroutine hnd_rt123
       go to 301
 c     yy = 33.0 to infinity             nroots=1,2, or 3
    90 ww1= dsqrt(pie4/yy)
+#if 1
+      n4=nroots-2
+      if(n4) 91,92,93
+#else
       if(nroots-2) 91,92,93
+#endif
    91 rt1=0.5d+00/(yy-0.5d+00)
       return
    92 if(yy.gt.40.0d+00) go to 102

src/NWints/ints_sp/sinfo.F

@@ -52,6 +52,7 @@ subroutine sinfo(
       double precision p2, temp, p3333, sinp, p3131, csai, cpai
 c
       integer ijshlo, ijato, klshlo, klato
+      integer*4 jtype4
       save ijshlo, ijato, klshlo, klato ! STATIC variables
 c     
       data ijshlo, ijato, klshlo, klato /4*0/
@@ -75,6 +76,7 @@ subroutine sinfo(
 *      klatom = kat*1000  + lat
 *      ijshel = ish*10000 + jsh
 *      klshel = ksh*10000 + lsh
+      jtype4=jtype
       ijatom = 1
       klatom = 1
       ijshel = 1
@@ -323,10 +325,20 @@ subroutine sinfo(
                if (qqq-sixty) 18000,18000,10000
 10000          ismlp(ind) = 2
                dp00p(ind) = dzero
+#if 1
+               jtype4=jtype-3
+               if (jtype4) 36000,36000,12000
+#else
                if (jtype-3) 36000,36000,12000
+#endif
 12000          dp01p(ind) = dzero
                conp(ind) = dzero
+#if 1
+               jtype4=jtype-5
+               if (jtype4) 14000,14000,16000
+#else
                if (jtype-5) 14000,14000,16000
+#endif
 14000          bpp(ind) = bpp(ind)*gab
                go to 36000
 16000          dp10p(ind) = dzero
@@ -343,9 +355,19 @@ subroutine sinfo(
 26000          ismlp(ind) = 2
 28000          qqqq = pito52*qq
                dp00p(ind) = qqqq*csai*csb(j)
+#if 1
+               jtype4=jtype-3
+               if (jtype4) 36000,36000,30000
+#else
                if (jtype-3) 36000,36000,30000
+#endif
 30000          dp01p(ind) = qqqq*csai*cpb(j)
+#if 1
+               jtype4=jtype-5
+               if (jtype4) 32000,32000,34000
+#else
                if (jtype-5) 32000,32000,34000
+#endif
 32000          conp(ind) = dp01p(ind)*eab
                dp00p(ind) = dp00p(ind)*gab/dp01p(ind)
                bpp(ind) = bpp(ind)*gab

src/NWints/texas/big_service.F

@@ -90,6 +90,7 @@ subroutine mxmtr(a,b,c,m,n,p)
 c     subroutine sdiag2v(m,n,a,d,x)
       subroutine sdiag2 (m,n,a,d,x)
       implicit real*8 (a-h,o-z)
+      integer*4 f4,h4,l4
 ckwol parameter (mxdim=401)
 ckwol parameter (mxdim=601)
       parameter (mxdim=2000)
@@ -164,18 +165,33 @@ subroutine sdiag2 (m,n,a,d,x)
          l=i-2
          h=0.0
          g=x(i,i-1)
+#if 1
+         l4=l
+         if (l4) 140,140,20
+#else
          if (l) 140,140,20
+#endif
    20    do 30 k=1,l
    30    h=h+x(i,k)**2
          s=h+g*g
          if (s.ge.tol) go to 40
          h=0.0
          go to 140
+#if 1
+   40    h4=h
+         if (h4) 140,140,50
+#else
    40    if (h) 140,140,50
+#endif
    50    l=l+1
          f=g
          g=dsqrt(s)
+#if 1
+         f4=f
+         if (f4) 70,70,60
+#else
          if (f) 70,70,60
+#endif
    60    g=-g
    70    h=s-f*g
          x(i,i-1)=f-g
@@ -410,6 +426,7 @@ subroutine sdiag2sc (m,n,a,d,x)
 c
       dimension a(m,m), d(m), x(m,m)
       dimension e(160)
+      integer*4 f4,h4,l4
 c
 c     correct adjustment for ieee floating point numbers (64 bits)
 c
@@ -433,18 +450,33 @@ subroutine sdiag2sc (m,n,a,d,x)
          l=i-2
          h=0.0d0
          g=x(i,i-1)
+#if 1
+         l4=l
+         if (l4) 140,140,20
+#else
          if (l) 140,140,20
+#endif
    20    do 30 k=1,l
    30    h=h+x(i,k)**2
          s=h+g*g
          if (s.ge.tol) go to 40
          h=0.0d0
          go to 140
+#if 1
+   40    h4=h
+         if (h4) 140,140,50
+#else
    40    if (h) 140,140,50
+#endif
    50    l=l+1
          f=g
          g=dsqrt(s)
+#if 1
+         f4=f
+         if (f4) 70,70,60
+#else
          if (f) 70,70,60
+#endif
    60    g=-g
    70    h=s-f*g
          x(i,i-1)=f-g