- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:
- 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.
- 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.
- 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?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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!

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page