exponent function gives different result for the same input

[Roderik]

I'm currently busy refactoring a script i got from somebody and a really strange thing happens. At a certain point in the code, a simple calculation is made. In order to check if the result and input is the same, I wrote the following in both versions:

open(310,file='testFile.txt',status='new')
write(310,*)'Begin output phi'
do px=1,Nx
  write(310,*) phix(px)
end do
do px=1,Nx
  do py=1,Ny
    Hx(px,py)=2.0D0*exp(-1.0D0*phix(px))
  end do
end do
write(310,*)'Begin output Hx'
do px=1,Nx
  write(310,*) Hx(px,1)
end do
close(310)

If I look at the output files for both runs, the value of phix is exactly the same for all px. The value of Hx is also exactly the same for all px, except for px = 6. In the original version, Hx(6,1) = 0.4861756328364343, and in the new version, Hx(6,1) = 0.4861756328364341.

If I change the calculation of Hx to

  Hx(px,py)=2.0D0*exp(-1.1D0*phix(px))

the samething occurs, but for px = 7, with Hx(7,1) = 0.4220554448918039 in the original version and Hx(7,1) = 0.4220554448918036 in the new version.

Does somebody know what could be the reason of this?

 

[ArkM]

You can't garantee that all phix elements values are the same. It's possible to get changed LSB of mantissa with refactored code (explicit or imlicit through optimization changes of calculations order, for example).

Don't forget: floating-point mantissa is a binary (not decimal) fraction. Printed values are decimal rounded fractions of true number values. Now you may get small variations (in LSB, as we can see) of result values, especially after exp function (as usually, it's implemented as polynomial approximation of math exp function).

It seems no reasons for any disturbance. There is a good chance to remember that floating point data are real number approximations only...

 

[Roderik]

Thanks for your comment! I figured out what is the reason of my problem, but I do not know how to solve it. This post can be somewhat long, but it takes me sometime to explain what the problem is, sorry for that.

Basically, the first version of the program was written something like this:

PROGRAM omega
	IMPLICIT NONE	
	INTEGER,PARAMETER :: Nx=10,Ny=15
        DOUBLE PRECISION,DIMENSION(1:Nx,1:Ny) :: Hx
        ! And the rest of the program
END

My goal for now is to alter this program, such that Nx and Ny can be set in an input file. The previous lines are replaced by

PROGRAM omega
	use parameters
        IMPLICIT NONE	

        call loadParameters
        allocate( Hx(1:Nx,1:Ny) )
        ! And the rest of the program
END

with the parameters module defined as

module parameters
	IMPLICIT NONE
	
	integer, public :: Nx, Ny
        DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:) :: Hx

        contains

        subroutine loadParameters
		IMPLICIT NONE
                namelist/Collocation/Nx,Ny
                open( unit=10, file='parameters.ini', status='old', action='read', iostat=fileError )
                read(10, nml=Collocation,iostat=fileError)
	end subroutine loadParameters
end module parameters

This gives no errors, and Hx has the same size as before. However, during the calculations done in the program, some digits change just as described in my first post. The problem is however that for some reason the problem becomes much larger and after a lot of calculations 1.23456789 in the original model is -2.3217865 in the new model (these numbers are just to give an example, but it really is in this order of magnitude).

As a final check, I gave the Nx and Ny in the parameters module a completely different name and changed the original program to

PROGRAM omega
	use parameters
        IMPLICIT NONE	
        INTEGER,PARAMETER :: Nx=10,Ny=15

        call loadParameters
        allocate( Hx(1:Nx,1:Ny) )
        ! And the rest of the program
END

Now, my results are completely the same as in the original model.

Long story short: By changing Nx and Ny from being PARAMETER to an input parameter, all calculations are effected a little bit, but after lots of calculations, this is a real problem!

Does somebody know what the problem could be?

 

[gullipe]

Hi Roderik

What happens if you change the function exp(...) to dexp (...) and declare the vector phix as double precision ?

 

[Roderik]

I tried that, but the same thing happens. Turning back to using INTEGER,PARAMETER's instead of parameters from a file solved the trick, but I do want these parameters in a file!

 

[Roderik]

I've managed to isolate the problem to a very small piece of code.

The original version:

program omega
  implicit none
  integer,parameter :: Nx=10
  double precision :: pi
  integer px

  pi = 4.0D0*DATAN(1.0D0)

  px = Nx - 1
  write(*,*)cos(pi*(6-1)/(px))
  write(*,*)cos(pi*(6-1)/(Nx-1))
end

This produces:
-0.1736481776669303
-0.1736481776669301

while this new version (the only difference is that PARAMETER is ommited)

program omega
  implicit none
  integer :: Nx=10
  double precision :: pi
  integer px

  pi = 4.0D0*DATAN(1.0D0)

  px = Nx - 1
  write(*,*)cos(pi*(6-1)/(px))
  write(*,*)cos(pi*(6-1)/(Nx-1))
end

produces -0.1736481776669303 twice.

Can someone explain me why these results are different?

 

[mikrom]

I tried your first example (with integer,parameter :: Nx=10) with 2 compilers (g95 and gfortran) and got the same results for both computations.
See here:

$ g95 omega.f90 -o omega  

$ omega
 -0.1736481776669303
 -0.1736481776669303

$ gfortran omega.f90 -o omega

$ omega
 -0.17364817766693030     
 -0.17364817766693030

What compiler are you using?

 

[mikrom]

The second example (with integer :: Nx=10) gives with g95 and gforran exact the same results.

 

[gullipe]

Hi both

Even this F77 program (compiled with Microsoft Fortran optimizing compiler, version 5.10):

program omega

implicit none
integer Nx
double precision pi
integer px
parameter (Nx=10)

pi = 4.0D0*DATAN(1.0D0)
px = Nx - 1
write(*,*)cos(pi*(6-1)/(px))
write(*,*)cos(pi*(6-1)/(Nx-1))
end

gave the result (same in both cases):
-1.736481776669303E-001
-1.736481776669303E-001

 

[xwb]

Which compiler/OS are you using?

Are you building it with the math library or the alternative math library. MS used to have an alternative math library for machines that didn't have an FPU.

It could be a fault in the FPU. Sun and MS sometimes have fixes for these hardware problems in their service packs.