Fortran test file compilation request

[mvrht__u3425923]

Hello,

I am in a process of searching for bug which I suspect depends on Fortran compiler version (gfortran vs Intel Fortran Compiler). However, I do not have Intel Fortran Compiler. Maybe someone could be so kind to compile small attached test file with Intel compiler and to post here its outputs?

I will appreciate your help, thanks in advance!

[David_Billinghurst]

I built this using Intel Fortran 2019 u4 under Windows 10 in 64-bit debug mode, with no optimization, full compiler warnings and full run-time checks.

Compiler output:

Compiling with Intel(R) Visual Fortran Compiler 19.0.4.245 [Intel(R) 64]...
test.for
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(22): warning #6717: This name has not been given an explicit type.   [PI]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(23): warning #6717: This name has not been given an explicit type.   [EL]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(24): warning #6717: This name has not been given an explicit type.   
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(24): warning #6717: This name has not been given an explicit type.   
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(25): warning #6717: This name has not been given an explicit type.   [A0]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(28): warning #6717: This name has not been given an explicit type.   [XT]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(30): warning #6717: This name has not been given an explicit type.   [CE1]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(34): warning #6717: This name has not been given an explicit type.   [CR]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(35): warning #6717: This name has not been given an explicit type.   
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(54): warning #6717: This name has not been given an explicit type.   [ZC]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(55): warning #6717: This name has not been given an explicit type.   [ZD]
C:\Users\david\source\repos\mvrht_01\mvrht_01\test.for(56): warning #6717: This name has not been given an explicit type.   [ZDC]

Output:

 (-1.89511781635594,-3.14159265358979)
 (-1.89511781635594,-3.14159265358979)

For interested viewers the code is

        PROGRAM TEST
        complex(8) :: POS, POS_RES, NEG, NEG_RES
        POS=CMPLX(-1, 0)
        call E1Z(POS, POS_RES)
        PRINT *, POS_RES
        NEG=CMPLX(-1, -0.0)
        call E1Z(NEG, NEG_RES)
        PRINT *, NEG_RES
        END


        SUBROUTINE E1Z(Z,CE1)
C
C       ====================================================
C       Purpose: Compute complex exponential integral E1(z)
C       Input :  z   --- Argument of E1(z)
C       Output:  CE1 --- E1(z)
C       ====================================================
C
        IMPLICIT COMPLEX*16 (C,Z)
        IMPLICIT DOUBLE PRECISION (A,D-H,O-Y)
        PI=3.141592653589793D0
        EL=0.5772156649015328D0
        X=DBLE(Z)
        A0=CDABS(Z)
C       Continued fraction converges slowly near negative real axis,
C       so use power series in a wedge around it until radius 40.0
        XT=-2*DABS(DIMAG(Z))
        IF (A0.EQ.0.0D0) THEN
           CE1=(1.0D+300,0.0D0)
        ELSE IF (A0.LE.5.0.OR.X.LT.XT.AND.A0.LT.40.0) THEN
C          Power series
           CE1=(1.0D0,0.0D0)
           CR=(1.0D0,0.0D0)
           DO 10 K=1,500
              CR=-CR*K*Z/(K+1.0D0)**2
              CE1=CE1+CR
              IF (CDABS(CR).LE.CDABS(CE1)*1.0D-15) GO TO 15
10         CONTINUE
15         CONTINUE
           IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0) THEN
C     Careful on the branch cut -- use the sign of the imaginary part
C     to get the right sign on the factor if pi.
              CE1=-EL-CDLOG(-Z)+Z*CE1-DSIGN(PI,DIMAG(Z))*(0.0D0,1.0D0)
           ELSE
              CE1=-EL-CDLOG(Z)+Z*CE1
           ENDIF
        ELSE
C          Continued fraction https://dlmf.nist.gov/6.9
C
C                           1     1     1     2     2     3     3
C          E1 = exp(-z) * ----- ----- ----- ----- ----- ----- ----- ...
C                         Z +   1 +   Z +   1 +   Z +   1 +   Z +
           ZC=0D0
           ZD=1/Z
           ZDC=1*ZD
           ZC=ZC + ZDC
           DO 20 K=1,500
              ZD=1/(ZD*K + 1)
              ZDC=(1*ZD - 1)*ZDC
              ZC=ZC + ZDC

              ZD=1/(ZD*K + Z)
              ZDC=(Z*ZD - 1)*ZDC
              ZC=ZC + ZDC

              IF (CDABS(ZDC).LE.CDABS(ZC)*1.0D-15.AND.K.GT.20) GO TO 25
20         CONTINUE
25         CE1=CDEXP(-Z)*ZC
           IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0) CE1=CE1-PI*(0.0D0,1.0D0)
        ENDIF
        RETURN
        END

[David_Billinghurst]

From your two inputs you appear to have issues with signed zeroes.  What do you see with this?

program signed_zero_01
  ! test signed zero support
  implicit none
  integer, parameter :: k = SELECTED_REAL_KIND(13) ! usually IEEE double precision
  write(*,*) sign(1.0_k, 0.0_k), sign(1.0_k,-0.0_k)
end

With gfortran 7.4.0 using default compiler flags on cygwin/Windows 10

   1.0000000000000000       -1.0000000000000000

intel fortran 2019u4 on Windows 10 with default compiler flags

   1.00000000000000        1.00000000000000

The behaviour can depend on your system math library, compiler version and compiler flags.  For intel fortran look at the assume minus0 compiler option.

[David_Billinghurst]

Perhaps it is time to explicitly test for signed zeroes with the ieee_arithmetic module.

program signed_zero_02
  ! test signed zero support
  use, intrinsic :: ieee_arithmetic
  implicit none
  integer, parameter :: k = SELECTED_REAL_KIND(13) ! usually IEEE double precision
  write(*,*) IEEE_CLASS( 0.0_k)==IEEE_POSITIVE_ZERO
  write(*,*) IEEE_CLASS(-0.0_k)==IEEE_POSITIVE_ZERO  
  write(*,*) IEEE_CLASS( 0.0_k)==IEEE_NEGATIVE_ZERO
  write(*,*) IEEE_CLASS(-0.0_k)==IEEE_NEGATIVE_ZERO
end

Both gfortran and ifort with the standard compiler flags give

 T
 F
 F
 T

I will have a play this evening.

[mvrht__u3425923]

Hi David,

thank you very much for helping!

David Billinghurst wrote:

From your two inputs you appear to have issues with signed zeroes.  What do you see with this?

program signed_zero_01
  ! test signed zero support
  implicit none
  integer, parameter :: k = SELECTED_REAL_KIND(13) ! usually IEEE double precision
  write(*,*) sign(1.0_k, 0.0_k), sign(1.0_k,-0.0_k)
end

With gfortran 7.4.0 under Linux it prints as yours

1.0000000000000000       -1.0000000000000000

Actually, I am not Fortran programer at all and

        PROGRAM TEST
        complex(8) :: POS, POS_RES, NEG, NEG_RES
        POS=CMPLX(-1, 0)
        call E1Z(POS, POS_RES)
        PRINT *, POS_RES
        NEG=CMPLX(-1, -0.0)
        call E1Z(NEG, NEG_RES)
        PRINT *, NEG_RES
        END

is my first Fortran program :)

SUBROUTINE E1Z(Z,CE1) is copied directly from SciPy project. Standard SciPy test suite has check

assert sc.exp1(complex(-1, 0)).imag == (-sc.exp1(complex(-1, -0.0)).imag)

which under the hood calls provided Fortran routine. SciPy compiled in Linux with MKL using gfortran successfully passes this test; SciPy taken from pip in Windows with OpenBlas and compiled with Mingw also passes the test; SciPy taken from https://www.lfd.uci.edu/~gohlke/pythonlibs/ compiled in windows with MKL using iFort failing this test. I was searching for the reason. So, thanks to your investigation, I guess now I have these options:

1. From man pages https://software.intel.com/en-us/fortran-compiler-developer-guide-and-reference-assume I can now understand that it is intended behavior for ifort (The default is -assume nominus0). Though it is old (Fortran 90/77) semantics and newer ones distinguishes the difference between -0.0 and +0.0. So, I guess, it is not a bug from ifort perspective.
2. I can try to ask SciPy team to add this flag to their build system's scripts, but considering that it is not generic flag for all Fortran compilers and recommended compiler is opensource gfortran (for Windows also) they might reject it.
3. And the last option is to ask the maintainer of http://www.lfd.uci.edu/~gohlke/pythonlibs to provide this flag explicitly while building SciPy in Windows with ifort.

David, am I right or am I missing something else?

[David_Billinghurst]

mvrht, u3425923 wrote:

David, am I right or am I missing something else?

A more interesting first Fortran program the the usual "Hello, world".  Well done.

The standard conforming way would be to use the IEEE_CLASS function from the ieee_arithmetic module to explicitly test for -0.0 (as above).  This was introduced in Fortran 2003 and is supported by recent gfortran and ifort.  This is a small change.  Leave it with me and I will have a look over the weekend.

[David_Billinghurst · Accepted solution]

I had a look at this over lunch.  Here is what I came up with.  The only changes to the subroutes E1Z are to use the ieee_arithmetic module and immediately below the "Careful on the branch cut" comment.  

The driver program compares points on either side of the branch cut with nearby points.  It passes with ifort and gfortran using the default compiler flags

      program test_e1z
      ! Check subroutine e1z near the negative x axis over [-40,0)
      ! Compare values for
      !   z = (x,+0) and (x,+eps)
      !   z = (x,-0) and (x,-eps)
      ! for eps = 8*epsilon(double) ~ 1.0e-16 for IEEE-754 double precision
      ! The relative difference (last column) between the function values should be ~ 1.0e-15
      implicit none
      integer, parameter :: dp=8
      real(kind=dp) :: x, y, eps, delta, reldiff, xmin
      complex(kind=dp) :: za, zb, ra, rb
      integer side
      ! check along negative real axis over (0,40]
      eps = 4 * epsilon(x) ! look really close the x-axis
      delta = 0.1_dp
      xmin = -40.0_dp
      do side = -1, 1, 2
        write(*,*) 'side = ', side
        x = 0.0_dp
        do
          x = x-delta
          !za = cmplx(x,side*0.0_8) ! oops, -1*0.0_8 /= -0.0_8
          if (side==1) then; y=0.0_dp; else; y=-0.0_dp; end if
          za = cmplx(x,y)
          call e1z(za,ra)
          zb = cmplx(x,side*eps)
          call e1z(zb,rb)
          reldiff = abs(ra-rb)/max(abs(ra),abs(rb))
          write(*,'(2(x,f7.2),5(x,es14.7))') real(za), imag(za), 
     &      real(ra), imag(ra), real(rb), imag(rb), reldiff
        if (x < xmin+delta/2) exit
        end do
        write(*,*)
      end do
      END


        SUBROUTINE E1Z(Z,CE1)
C
C       ====================================================
C       Purpose: Compute complex exponential integral E1(z)
C       Input :  z   --- Argument of E1(z)
C       Output:  CE1 --- E1(z)
C       ====================================================
C
        use, intrinsic :: ieee_arithmetic
        IMPLICIT COMPLEX*16 (C,Z)
        IMPLICIT DOUBLE PRECISION (A,D-H,O-Y)
        PI=3.141592653589793D0
        EL=0.5772156649015328D0
        X=DBLE(Z)
        A0=CDABS(Z)
C       Continued fraction converges slowly near negative real axis,
C       so use power series in a wedge around it until radius 40.0
        XT=-2*DABS(DIMAG(Z))
        IF (A0.EQ.0.0D0) THEN
           CE1=(1.0D+300,0.0D0)
        ELSE IF (A0.LE.5.0.OR.X.LT.XT.AND.A0.LT.40.0) THEN
C          Power series
           CE1=(1.0D0,0.0D0)
           CR=(1.0D0,0.0D0)
           DO 10 K=1,500
              CR=-CR*K*Z/(K+1.0D0)**2
              CE1=CE1+CR
              IF (CDABS(CR).LE.CDABS(CE1)*1.0D-15) GO TO 15
10         CONTINUE
15         CONTINUE
           IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0) THEN
C     Careful on the branch cut -- use the sign of the imaginary part
C     to get the right sign on the factor if pi.
              if (IEEE_CLASS(imag(z))==IEEE_POSITIVE_ZERO) then
                CE1=-EL-CDLOG(-Z)+Z*CE1-PI*(0.0D0,1.0D0)
              else
                CE1=-EL-CDLOG(-Z)+Z*CE1+PI*(0.0D0,1.0D0)
              end if
           ELSE
              CE1=-EL-CDLOG(Z)+Z*CE1
           ENDIF
        ELSE
C          Continued fraction https://dlmf.nist.gov/6.9
C
C                           1     1     1     2     2     3     3
C          E1 = exp(-z) * ----- ----- ----- ----- ----- ----- ----- ...
C                         Z +   1 +   Z +   1 +   Z +   1 +   Z +
           ZC=0D0
           ZD=1/Z
           ZDC=1*ZD
           ZC=ZC + ZDC
           DO 20 K=1,500
              ZD=1/(ZD*K + 1)
              ZDC=(1*ZD - 1)*ZDC
              ZC=ZC + ZDC

              ZD=1/(ZD*K + Z)
              ZDC=(Z*ZD - 1)*ZDC
              ZC=ZC + ZDC

              IF (CDABS(ZDC).LE.CDABS(ZC)*1.0D-15.AND.K.GT.20) GO TO 25
20         CONTINUE
25         CE1=CDEXP(-Z)*ZC
           IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0) CE1=CE1-PI*(0.0D0,1.0D0)
        ENDIF
        RETURN
        END

[mvrht__u3425923]

Hi David,

yes, with your proposed changes it works! I will try to propose this fix to SciPy project

Thank you very much, you are extremally helpful!