- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hi all

According to the Intel documentation, the default I/O rounding mode is 'PROCESSOR_DEFINED'. For ifort, this corresponds to 'NEAREST' which conforms to the ISO/IEC/IEEE 60559:2011 standard specification for 'round to nearest, ties to even'.

The following program shows a different behavior:

```
program round
!
implicit none
!
integer, parameter :: sp = kind(0.0)
integer, parameter :: dp = kind(0.0d0)
integer, parameter :: qp = selected_real_kind(33)
!
character(len=:), allocatable :: inp_sp, inp_dp, inp_qp
real(sp) :: vio_dec_sp, vio_hex_sp, vas_dec_sp
real(dp) :: vio_dec_dp, vio_hex_dp, vas_dec_dp
real(qp) :: vio_dec_qp, vio_hex_qp, vas_dec_qp
!
write (*, '(/a/a)') '1. Single precision', &
'-------------------'
! The number 2**23 + 2**(-1) in decimal and hexadecimal:
inp_sp = '8388608.5 0x1.000001p+23'
read (inp_sp, *) vio_dec_sp, vio_hex_sp
vas_dec_sp = 8388608.5_sp
!
write (*, '(a)') 'Values from I/O:'
write (*, '(a,z8.8,2x,es16.9)') 'decimal : ', vio_dec_sp, vio_dec_sp
write (*, '(a,z8.8,2x,es16.9)') 'hexadecimal: ', vio_hex_sp, vio_hex_sp
!
write (*, '(/a)') 'Value from assignment:'
write (*, '(a,z8.8,2x,es16.9)') 'decimal : ', vas_dec_sp, vas_dec_sp
!
!
write (*, '(/a/a)') '2. Double precision', &
'-------------------'
! The number 2**52 + 2**(-1) in decimal and hexadecimal:
inp_dp = '4503599627370496.5 0x1.00000000000008p+52'
read (inp_dp, *) vio_dec_dp, vio_hex_dp
vas_dec_dp = 4503599627370496.5_dp
!
write (*, '(/a)') 'Values from I/O:'
write (*, '(a,z16.16,2x,es23.16)') 'decimal : ', vio_dec_dp, vio_dec_dp
write (*, '(a,z16.16,2x,es23.16)') 'hexadecimal: ', vio_hex_dp, vio_hex_dp
!
write (*, '(/a)') 'Value from assignment:'
write (*, '(a,z16.16,2x,es23.16)') 'decimal : ', vas_dec_dp, vas_dec_dp
!
!
write (*, '(/a/a)') '3. Quad precision', &
'-----------------'
! The number 2**112 + 2**(-1) in decimal and hexadecimal:
inp_qp = '5192296858534827628530496329220096.5 &
&0x1.00000000000000000000000000008p+112'
read (inp_qp, *) vio_dec_qp, vio_hex_qp
vas_dec_qp = 5192296858534827628530496329220096.5_qp
!
write (*, '(/a)') 'Values from I/O:'
write (*, '(a,z32.32,2x,es40.33)') 'decimal : ', vio_dec_qp, vio_dec_qp
write (*, '(a,z32.32,2x,es40.33)') 'hexadecimal: ', vio_hex_qp, vio_hex_qp
!
write (*, '(/a)') 'Value from assignment:'
write (*, '(a,z32.32,2x,es40.33)') 'decimal : ', vas_dec_qp, vas_dec_qp
!
end program round
```

Output from running the program on Windows 10 64bit with the ifort (2021.9.0) and the ifx (2023.1.0) compiler:

```
1. Single precision
-------------------
Values from I/O:
decimal : 4B000000 8.388608000E+06
hexadecimal: 4B000000 8.388608000E+06
Value from assignment:
decimal : 4B000001 8.388609000E+06
2. Double precision
-------------------
Values from I/O:
decimal : 4330000000000001 4.5035996273704970E+15
hexadecimal: 4330000000000000 4.5035996273704960E+15
Value from assignment:
decimal : 4330000000000001 4.5035996273704970E+15
3. Quad precision
-----------------
Values from I/O:
decimal : 406F0000000000000000000000000001 5.192296858534827628530496329220097E+33
hexadecimal: 406F0000000000000000000000000000 5.192296858534827628530496329220096E+33
Value from assignment:
decimal : 406F0000000000000000000000000001 5.192296858534827628530496329220097E+33
```

As you can see, the hexadecimal values are all correctly rounded. From the decimal values, only the single precision value is correctly rounded. The other precisions do not honor the 'ties to even' rule.

Another point is rounding on assignment. None of the assigned decimal values is correctly rounded according to the 'round to nearest, ties to even' rule.

I have not found any information about the rounding rules for assignment in the Fortran standard or in the Intel documentation. But I think the rule used by most compilers is 'round to nearest, ties to even'.

Questions:

Is this behavior intended, or is it eventually a bug?

Are there any compiler options to change the behavior?

Thanks for any insight.

--

Thomas

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I have not studied this carefully enough, so help me out.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:

8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

```
gfortran -O0 -g wnew.f90
rwgreen@orcsle153:~/quad$ ./a.out
1. Single precision
-------------------
At line 18 of file wnew.f90
Fortran runtime error: Bad real number in item 2 of list input
Error termination. Backtrace:
#0 0x7fc2f4c23a12 in ???
#1 0x7fc2f4c24509 in ???
#2 0x7fc2f4c2505f in ???
#3 0x7fc2f4e62c73 in ???
#4 0x7fc2f4e6403d in ???
#5 0x7fc2f4e64c8a in ???
#6 0x4013f2 in round
at /home/quad/wnew.f90:18
#7 0x402181 in main
at /home/quad/wnew.f90:59
```

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

```
read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp
```

to explicit format edit descriptors

```
read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp
```

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct.

I expect that the output from the decimal floating point input values is the same as the output from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

```
read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp
```

to explicit format edit descriptors

```
read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp
```

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct. I expect the output from the decimal floating point input values to be the same as the output from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

```
read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp
```

to explicit format edit descriptors

```
read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp
```

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct. I expect the same output from the decimal floating point input values as from from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

--

Thomas

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Please excuse the multiple replies. The first two replies where marked as spam and didn't showed up in the forum first. I observed that my Pi-Hole blocked some URLs referenced from the Intel forum. I put them on the whitelist and replied a third time and this time it was successful. Unfortunately the first two replies showed now up too.

Never less, I made a different example to show what I think is a problem:

```
program round_c
!
use, intrinsic :: iso_c_binding
!
implicit none
!
integer, parameter :: dp = kind(0.0d0)
!
interface
function strtod(str, endptr) bind(c)
use, intrinsic :: iso_c_binding
implicit none
real(c_double) :: strtod
character(kind=c_char) :: str(*)
type(c_ptr), optional :: endptr
end function strtod
end interface
!
character(len=:), allocatable :: inp
real(dp) :: v1, v2
!
! The number 2**52 + 2**(-1) in decimal:
inp = '4503599627370496.5'
!
read (inp, '(f18.0)') v1
!
v2 = strtod(inp//c_null_char)
!
write (*, '(/a/a)') 'Results', &
'-------'
write (*, '(a,z16.16,2x,es23.16)') 'internal read : ', v1, v1
write (*, '(a,z16.16,2x,es23.16)') 'C strtod() : ', v2, v2
write (*, '(a,a)') 'Are they equal? : ', trim(merge('Yes', 'No ', v1 == v2))
!
end program round_c
```

Output from ifort (2021.9.0) and ifx (2023.1.0):

```
Results
-------
internal read : 4330000000000001 4.5035996273704970E+15
C strtod() : 4330000000000000 4.5035996273704960E+15
Are they equal? : No
```

Output from gfortran 12.2

```
Results
-------
internal read : 4330000000000000 4.5035996273704960E+15
C strtod() : 4330000000000000 4.5035996273704960E+15
Are they equal? : Yes
```

I would expect that the Intel compilers produce the same result as the gfortran compiler.

Can someone else reproduce these results?

Thanks

--

Thomas

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