1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
! { dg-do run }
! { dg-additional-options "-mfp-rounding-mode=d" { target alpha*-*-* } }
use, intrinsic :: ieee_features, only : ieee_rounding
use, intrinsic :: ieee_arithmetic
implicit none
interface check_equal
procedure check_equal_float, check_equal_double
end interface
interface check_not_equal
procedure check_not_equal_float, check_not_equal_double
end interface
interface divide
procedure divide_float, divide_double
end interface
real :: sx1, sx2, sx3
double precision :: dx1, dx2, dx3
type(ieee_round_type) :: mode
! We should support at least C float and C double types
if (ieee_support_rounding(ieee_nearest)) then
if (.not. ieee_support_rounding(ieee_nearest, 0.)) STOP 1
if (.not. ieee_support_rounding(ieee_nearest, 0.d0)) STOP 2
end if
! The initial rounding mode should probably be NEAREST
! (at least on the platforms we currently support)
if (ieee_support_rounding(ieee_nearest, 0.)) then
call ieee_get_rounding_mode (mode)
if (mode /= ieee_nearest) STOP 3
end if
if (ieee_support_rounding(ieee_up, sx1) .and. &
ieee_support_rounding(ieee_down, sx1) .and. &
ieee_support_rounding(ieee_nearest, sx1) .and. &
ieee_support_rounding(ieee_to_zero, sx1)) then
sx1 = 1
sx2 = 3
sx1 = divide(sx1, sx2, ieee_up)
sx3 = 1
sx2 = 3
sx3 = divide(sx3, sx2, ieee_down)
call check_not_equal(sx1, sx3)
call check_equal(sx3, nearest(sx1, -1.))
call check_equal(sx1, nearest(sx3, 1.))
call check_equal(1./3., divide(1., 3., ieee_nearest))
call check_equal(-1./3., divide(-1., 3., ieee_nearest))
call check_equal(divide(3., 7., ieee_to_zero), &
divide(3., 7., ieee_down))
call check_equal(divide(-3., 7., ieee_to_zero), &
divide(-3., 7., ieee_up))
end if
if (ieee_support_rounding(ieee_up, dx1) .and. &
ieee_support_rounding(ieee_down, dx1) .and. &
ieee_support_rounding(ieee_nearest, dx1) .and. &
ieee_support_rounding(ieee_to_zero, dx1)) then
dx1 = 1
dx2 = 3
dx1 = divide(dx1, dx2, ieee_up)
dx3 = 1
dx2 = 3
dx3 = divide(dx3, dx2, ieee_down)
call check_not_equal(dx1, dx3)
call check_equal(dx3, nearest(dx1, -1.d0))
call check_equal(dx1, nearest(dx3, 1.d0))
call check_equal(1.d0/3.d0, divide(1.d0, 3.d0, ieee_nearest))
call check_equal(-1.d0/3.d0, divide(-1.d0, 3.d0, ieee_nearest))
call check_equal(divide(3.d0, 7.d0, ieee_to_zero), &
divide(3.d0, 7.d0, ieee_down))
call check_equal(divide(-3.d0, 7.d0, ieee_to_zero), &
divide(-3.d0, 7.d0, ieee_up))
end if
contains
real function divide_float (x, y, rounding) result(res)
use, intrinsic :: ieee_arithmetic
real, intent(in) :: x, y
type(ieee_round_type), intent(in) :: rounding
type(ieee_round_type) :: old
call ieee_get_rounding_mode (old)
call ieee_set_rounding_mode (rounding)
res = x / y
call ieee_set_rounding_mode (old)
end function
double precision function divide_double (x, y, rounding) result(res)
use, intrinsic :: ieee_arithmetic
double precision, intent(in) :: x, y
type(ieee_round_type), intent(in) :: rounding
type(ieee_round_type) :: old
call ieee_get_rounding_mode (old)
call ieee_set_rounding_mode (rounding)
res = x / y
call ieee_set_rounding_mode (old)
end function
subroutine check_equal_float (x, y)
real, intent(in) :: x, y
if (x /= y) then
print *, x, y
STOP 4
end if
end subroutine
subroutine check_equal_double (x, y)
double precision, intent(in) :: x, y
if (x /= y) then
print *, x, y
STOP 5
end if
end subroutine
subroutine check_not_equal_float (x, y)
real, intent(in) :: x, y
if (x == y) then
print *, x, y
STOP 6
end if
end subroutine
subroutine check_not_equal_double (x, y)
double precision, intent(in) :: x, y
if (x == y) then
print *, x, y
STOP 7
end if
end subroutine
end
|
