
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
ljfarrugia
Joined: 06 Jun 2016 Posts: 37

Posted: Mon Oct 07, 2024 11:47 am Post subject: Confusing difference 


The trivial code below produces what I consider inconsistent results. Specifically there is a difference between results when the mod intrinsic function is presented with a plain argument (c) or with an argument which requires evaluation (a+b). The code works exactly as expected when compiled with gfortran on Linux, but gives confusing results with Silverfrost Fortran. I am using version 8.82.0
program main
double precision :: a,b,c,d
logical :: f1, f2
a = 1.0d0/3.0d0
b = 2.0d0/3.0d0
write(*,'(a,e15.10e1)') 'a = ',a
write(*,'(a,e15.10e1)') 'b = ',b
c = a + b
write(*,'(a,e15.10e1)') 'a + b = ',a+b
write(*,'(a,e15.10e1)') 'c = ',c
f1 = (c > 1.0d0)
f2 = (c < 1.0d0)
write(*,*) 'c is smaller than 1.0 = ',f2
write(*,*) 'c is larger than 1.0 = ',f1
f1 = (a+b > 1.0d0)
f2 = (a+b < 1.0d0)
write(*,*) 'a+b is smaller than 1.0 = ',f2
write(*,*) 'a+b is larger than 1.0 = ',f1
c = abs(mod(c,1.0d0))
d = abs(mod(a+b,1.0d0))
write(*,'(a,e15.10e1)') 'abs(mod(c,1.0)) = ',c
write(*,'(a,e15.10e1)') 'abs(mod(a+b,1.0)) = ',d
end 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2587 Location: Sydney

Posted: Mon Oct 07, 2024 12:27 pm Post subject: 


I changed your program to:
Code:  program main
use iso_fortran_env
implicit none
double precision :: a,b,c,d
logical :: f1, f2
write (*,*) 'Vern : ',compiler_version ()
a = 1.0d0/3.0d0
b = 2.0d0/3.0d0
write(*,'(a,e18.12)') 'a = ',a
write(*,'(a,e18.12)') 'b = ',b
c = a + b
write(*,'(a,e18.12)') 'a + b = ',a+b
write(*,'(a,e18.12)') 'c = ',c
f1 = (c > 1.0d0)
f2 = (c < 1.0d0)
write(*,*) 'c is smaller than 1.0 = ',f2
write(*,*) 'c is larger than 1.0 = ',f1
f1 = (a+b > 1.0d0)
f2 = (a+b < 1.0d0)
write(*,*) 'a+b is smaller than 1.0 = ',f2
write(*,*) 'a+b is larger than 1.0 = ',f1
c = abs(mod(c,1.0d0))
d = abs(mod(a+b,1.0d0))
write(*,'(a,e18.12)') 'abs(mod(c,1.0)) = ',c
write(*,'(a,e18.12)') 'abs(mod(a+b,1.0)) = ',d
end 
The results are:
Code:  Vern : FTN95 v9.04.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00
Vern : GCC version 14.1.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00

I can't see the confusion ?
perhaps you could use 1.0d0 in your mod functions ? 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 717 Location: Hamilton, Lanarkshire, Scotland.

Posted: Mon Oct 07, 2024 4:31 pm Post subject: 


Running John's program with FTN95 v9.03.0.
The 32 bit complier returns:
Code:  Vern : FTN95 v9.03.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = T
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.100000000000E+01 
The 64 bit complier returns:
Code:  Vern : FTN95 v9.03.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00 
Also
Code:  Vern : GCC version 12.3.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00 
Code:  Vern :
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel
(R) 64, Version 2021.4.0 Build 20210910_000000
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00 


Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2587 Location: Sydney

Posted: Tue Oct 08, 2024 3:52 am Post subject: 


OK,
A new program that demonstrates 32bit has different roundoff.
Code:  program main
use iso_fortran_env
implicit none
double precision :: a,b,c,d, one = 1
logical :: f1, f2
write (*,*) 'Version : ',compiler_version ()
write (*,*) 'Options : ',compiler_options ()
a = 1.0d0/3.0d0
b = 2.0d0/3.0d0
write(*,11) 'a = ',a
write(*,11) 'b = ',b
c = a + b
write(*,11) 'a + b = ',a+b
write(*,11) 'c = ',c
f1 = (c > one)
f2 = (c < one)
write(*,*) 'c is smaller than 1.0 = ',f2
write(*,*) 'c is larger than 1.0 = ',f1
f1 = (a+b > one)
f2 = (a+b < one)
write(*,*) 'a+b is smaller than 1.0 = ',f2
write(*,*) 'a+b is larger than 1.0 = ',f1
c = abs(mod(c,one))
d = abs(mod(a+b,one))
write(*,11) 'abs(mod(c,1.0)) = ',c
write(*,11) 'abs(mod(a+b,1.0)) = ',d
write(*,11) 'one(a+b) = ',one(a+b)
write(*,11) '(onea)b = ',(onea)b
write(*,11) 'oneab = ',oneab
11 format (a,e24.18)
end 
Results for revised program:
Code:  Version : FTN95 v9.04.0
Options : CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
a = 0.333333333333333315E+00
b = 0.666666666666666630E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = T
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.100000000000000000E+01
one(a+b) = 0.555111512312578271E16
(onea)b = 0.555111512312578271E16
oneab = 0.555111512312578271E16
Version : FTN95 v9.04.0
Options : 64;CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
a = 0.333333333333333303E+00
b = 0.666666666666666606E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000000000E+00
one(a+b) = 0.000000000000000000E+00
(onea)b = 0.111022302462515654E15
oneab = 0.111022302462515654E15
Version : GCC version 14.1.0
Options : cpp iprefix C:/Program Files (x86)/gcc_eq/gcc_14.1.0/bin/../lib/gcc/x86_64w64mingw32/14.1.0/ U_REENTRANT mtune=generic march=x8664 fnounderscoring fdollarok
a = 0.333333333333333315E+00
b = 0.666666666666666630E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
c is smaller than 1.0 = F
c is larger than 1.0 = F
a+b is smaller than 1.0 = F
a+b is larger than 1.0 = F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000000000E+00
one(a+b) = 0.000000000000000000E+00
(onea)b = 0.111022302462515654E15
oneab = 0.111022302462515654E15

The differing results demonstrate that these type of tests should account for roundoff uncertainty.
I am really not sure if 32bit uses 80bit registers for f2 = (a+b < one) ?
c (= a+b) would be rounded to a 64bit real. 

Back to top 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 8019 Location: Salford, UK

Posted: Tue Oct 08, 2024 8:59 am Post subject: 


For maximum precision you can use real*10 with Win32 rather than x64.
Code:  options(alt_kinds)
program main
real*10 :: a,b,c,d
logical :: f1, f2
a = 1.0_10/3.0_10
b = 2.0_10/3.0_10
write(*,'(a,e22.17e1)') 'a = ',a
write(*,'(a,e22.17e1)') 'b = ',b
c = a + b
write(*,'(a,e22.17e1)') 'a + b = ',a+b
write(*,'(a,e22.17e1)') 'c = ',c
f1 = (c > 1.0_10)
f2 = (c < 1.0_10)
write(*,*) 'c is smaller than 1.0 = ',f2
write(*,*) 'c is larger than 1.0 = ',f1
f1 = (a+b > 1.0_10)
f2 = (a+b < 1.0_10)
write(*,*) 'a+b is smaller than 1.0 = ',f2
write(*,*) 'a+b is larger than 1.0 = ',f1
c = abs(mod(c,1.0_10))
d = abs(mod(a+b,1.0_10))
write(*,'(a,e22.17e1)') 'abs(mod(c,1.0)) = ',c
write(*,'(a,e22.17e1)') 'abs(mod(a+b,1.0)) = ',d
end 
a and b are evaluated at compile time whilst c and d are evaluated at runtime.
The main point is that the values of a and b are inevitably subject to roundoff error. 

Back to top 


ljfarrugia
Joined: 06 Jun 2016 Posts: 37

Posted: Tue Oct 08, 2024 12:35 pm Post subject: 


Many thanks for all the comments. I kind of guessed that rounding errors were the cause of the discrepancy, but I didn't realise that 32bit executables were so subject to this issue.
I came across this issue when trying to understand why I was getting seriously differing results from a rather large program compiled under Windows & Linux. I normally have used 32 bit executables and in the past I havent noticed any great issues between the platforms
Any suggestions as the the best way forward? Just use 64 bit executables?
Louis 

Back to top 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 8019 Location: Salford, UK

Posted: Tue Oct 08, 2024 1:05 pm Post subject: 


64 bit executables are not generally the answer to issues relating to floating point precision.
For example 32 bit FTN95 provides for 80 bit floating point values as well as 32 and 64. Most Fortran compilers only provide 32 bit and 64 bit floating point values.
It is the programmer's responsibility to be aware/beware of calculations that are prone to roundoff error. 

Back to top 


Kenneth_Smith
Joined: 18 May 2012 Posts: 717 Location: Hamilton, Lanarkshire, Scotland.

Posted: Tue Oct 08, 2024 3:46 pm Post subject: 


This is a contrived example, but it illustrates Paul's point about 64 bit executables not being the "fix all" answer to floating point issues.
Code:  program p
use iso_fortran_env
implicit none
integer, parameter :: n = 5
real*8 :: a(n) = [1.d0,1.d16,1.d16,1.d0,1.d16]
real*8 :: acc
integer :: i(n)
integer :: j
write (*,*) 'Vern : ',compiler_version ()
write (*,*) 'Options : ',compiler_options ()
write(*,'(a)') 'Input array'
do j = 1, n, 1
write(*,'(G16.10)') a(j)
end do
write(*,*)
write(*,'(a,G16.10)') 'Intrinsic sum = ', sum(a)
call dsort@(i,abs(a),5) ! or qsortd from NSWC Library if not using FTN95
acc = 0.d0
do j = n, 1, 1
acc = acc + a(i(j))
end do
write(*,*)
write(*,'(a,G16.10)') 'Sorted sum = ',acc
end program

32 bit executable:
Code:  Vern : FTN95 v9.03.0
Options : ERROR_NUMBERS;LINK;NO_BANNER;UNLIMITED_ERRORS;VS7;
Input array
1.000000000
0.1000000000E15
0.1000000000E15
1.000000000
0.1000000000E15
Intrinsic sum = 0.2999268806E15
Sorted sum = 0.3000000000E15 
64 bit executable
Code:  Vern : FTN95 v9.03.0
Options : 64;ERROR_NUMBERS;LINK;NO_BANNER;UNLIMITED_ERRORS;VS7;
Input array
1.000000000
0.1000000000E15
0.1000000000E15
1.000000000
0.1000000000E15
Intrinsic sum = 0.1000000000E15
Sorted sum = 0.3000000000E15 


Back to top 


DanRRight
Joined: 10 Mar 2008 Posts: 2867 Location: South Pole, Antarctica

Posted: Wed Oct 09, 2024 7:36 am Post subject: 


I am curious why 32bit FTN95 still retain 80bit accuracy while 64bit one has only 64bit one? So both Intel and AMD processors have hardware for 80bit while 64bit FTN95 ignores it (as well as other Fortran compilers?)
Also how about SSE2 and AVX instructions with 64bit compilers. Do they still support 80bit accuracy? I've heard this opinion on the net. Also have heard that it's Windows which does not support double precision 80bit accuracy, Linux supports it (is this true?)
If anyone has some free time to investigate all that mess please try to make short confirming tests 

Back to top 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 8019 Location: Salford, UK

Posted: Wed Oct 09, 2024 9:15 am Post subject: 


64 bit applications can in theory access the the x87 FPU coprocessor which has 80 bit floating point precision. This coprocessor is separate from the the CPU that uses SSE etc. instructions and which basically does not extend to 80 bits.
The problem is that, at least under Windows, a massive time penalty is incurred when mixing x87 instructions with SSE etc. instructions. So compilers normally avoid this mix. 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2587 Location: Sydney

Posted: Wed Oct 09, 2024 1:54 pm Post subject: Re: 


DanRRight wrote:  Also how about SSE2 and AVX instructions with 64bit compilers. Do they still support 80bit accuracy? 
No ! this is the problem with AVX registers, they only support grouped 4/8 byte instructions.
As for 10byte reals; I do not know if they are hardware supported in Intel or Gfortran x64 compilers. I think the answer could be No or maybe!
Given how simd instructions are so preferred, it is too inefficient to target 10byte reals.
I think this thread is all about how 80bit registers are used in the 32bit FTN95 compiler. If the latest value is retained in an 80bit register for the next calculation, then the roundoff errors are much less. 

Back to top 




You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum

Powered by phpBB © 2001, 2005 phpBB Group
