Line data Source code
1 :
2 : module accurate_sum
3 :
4 : use const_def, only: qp
5 :
6 : implicit none
7 :
8 : private
9 : public :: accurate_real, neumaier_sum, operator(+), operator(-), assignment(=), &
10 : operator(<), operator(>), operator(*), operator(/)
11 :
12 : ! Type for easily using Neumaier's summation algorithm in place of normal addition.
13 : type accurate_real
14 : real(qp) :: sum
15 : real(qp) :: compensator
16 :
17 : contains
18 :
19 : procedure :: value
20 :
21 : end type accurate_real
22 :
23 : interface operator(+)
24 : procedure add_accurate_real
25 : procedure add_accurate_real_num
26 : procedure add_accurate_real_num_rev
27 : end interface operator(+)
28 :
29 : interface operator(-)
30 : procedure sub_accurate_real
31 : procedure sub_accurate_real_num
32 : procedure sub_accurate_real_num_rev
33 : end interface operator(-)
34 :
35 : interface operator(*)
36 : procedure mult_acc_real
37 : procedure mult_real_acc
38 : procedure mult_acc_acc
39 : end interface operator(*)
40 :
41 : interface operator(/)
42 : procedure div_acc_real
43 : procedure div_real_acc
44 : procedure div_acc_acc
45 : end interface operator(/)
46 :
47 : interface assignment(=)
48 : procedure set_accurate_real_num
49 : procedure set_num_accurate_real
50 : procedure set_acc_acc
51 : end interface assignment(=)
52 :
53 : interface operator(<)
54 : procedure acc_less_than_num
55 : procedure num_less_than_acc
56 : end interface operator(<)
57 :
58 : interface operator(>)
59 : procedure acc_greater_than_num
60 : procedure num_greater_than_acc
61 : end interface operator(>)
62 :
63 : contains
64 :
65 : ! Helper method for evaluating an accurate_real.
66 0 : function value(this) result(res)
67 : ! Inputs
68 : class(accurate_real), intent(in) :: this
69 : real(qp) :: res
70 :
71 0 : res = this%sum + this%compensator
72 0 : end function value
73 :
74 : ! Performs one step of Neumaier's summation algorithm
75 : ! (see doi:10.1002/zamm.19740540106)
76 : ! NOTE: Do not use --ffast-math or the like with this routine.
77 : ! The compiler will optimize away the trick which preserves
78 : ! accuracy.
79 0 : subroutine neumaier_sum(sum, compensator, summand)
80 : ! Inputs
81 : real(qp) :: sum, compensator, summand
82 :
83 : ! Intermediates
84 0 : real(qp) :: provisional
85 :
86 0 : provisional = sum + summand
87 0 : if (abs(sum) >= abs(summand)) then
88 0 : compensator = compensator + (sum - provisional) + summand
89 : else
90 0 : compensator = compensator + (summand - provisional) + sum
91 : end if
92 0 : sum = provisional
93 :
94 0 : end subroutine neumaier_sum
95 :
96 : ! The remaining are helper methods which overload +,-,*,/,=,<,>
97 : ! to work with accurate_real and real(qp) numbers interchangeably.
98 :
99 : ! *************
100 0 : type(accurate_real) function mult_acc_acc(op1, op2) result(ret)
101 : type(accurate_real), intent(in) :: op1
102 : type(accurate_real), intent(in) :: op2
103 :
104 0 : ret%sum = op1%sum*op2%sum
105 0 : ret%compensator = op1%compensator*op2%sum
106 :
107 0 : call neumaier_sum(ret%sum, ret%compensator, op1%sum*op2%compensator)
108 0 : call neumaier_sum(ret%sum, ret%compensator, op1%compensator*op2%compensator)
109 :
110 0 : end function mult_acc_acc
111 :
112 0 : type(accurate_real) function mult_acc_real(op1, op2) result(ret)
113 : type(accurate_real), intent(in) :: op1
114 : real(qp), intent(in) :: op2
115 :
116 0 : ret%sum = op1%sum*op2
117 0 : ret%compensator = op1%compensator*op2
118 :
119 0 : end function mult_acc_real
120 :
121 0 : type(accurate_real) function mult_real_acc(op1, op2) result(ret)
122 : real(qp), intent(in) :: op1
123 : type(accurate_real), intent(in) :: op2
124 :
125 0 : ret%sum = op2%sum*op1
126 0 : ret%compensator = op2%compensator*op1
127 :
128 0 : end function mult_real_acc
129 :
130 : ! ///////////
131 0 : type(accurate_real) function div_acc_acc(op1, op2) result(ret)
132 : type(accurate_real), intent(in) :: op1
133 : type(accurate_real), intent(in) :: op2
134 :
135 0 : ret%sum = op1%value()/op2%value()
136 0 : ret%compensator = 0
137 :
138 0 : end function div_acc_acc
139 :
140 0 : type(accurate_real) function div_acc_real(op1, op2) result(ret)
141 : type(accurate_real), intent(in) :: op1
142 : real(qp), intent(in) :: op2
143 :
144 0 : ret%sum = op1%sum/op2
145 0 : ret%compensator = op1%compensator/op2
146 :
147 0 : end function div_acc_real
148 :
149 0 : type(accurate_real) function div_real_acc(op1, op2) result(ret)
150 : real(qp), intent(in) :: op1
151 : type(accurate_real), intent(in) :: op2
152 :
153 0 : ret%sum = op1/op2%value()
154 0 : ret%compensator = 0
155 :
156 0 : end function div_real_acc
157 :
158 : ! ============
159 0 : subroutine set_accurate_real_num(this, new)
160 : type(accurate_real), intent(out) :: this
161 : real(qp), intent(in) :: new
162 0 : this%sum = new
163 0 : this%compensator = 0
164 0 : end subroutine set_accurate_real_num
165 :
166 0 : subroutine set_num_accurate_real(this, new)
167 : real(qp), intent(out) :: this
168 : type(accurate_real), intent(in) :: new
169 0 : this = new%value()
170 0 : end subroutine set_num_accurate_real
171 :
172 0 : subroutine set_acc_acc(this, new)
173 : type(accurate_real), intent(out) :: this
174 : type(accurate_real), intent(in) :: new
175 0 : this%sum = new%sum
176 0 : this%compensator = new%compensator
177 0 : end subroutine set_acc_acc
178 :
179 : ! <<<<<<<<<<<
180 0 : logical function acc_less_than_num(acc, num) result(ret)
181 : type(accurate_real), intent(in) :: acc
182 : real(qp), intent(in) :: num
183 0 : if (acc%value() < num) then
184 : ret = .true.
185 : else
186 0 : ret = .false.
187 : end if
188 0 : end function acc_less_than_num
189 :
190 0 : logical function num_less_than_acc(num, acc) result(ret)
191 : type(accurate_real), intent(in) :: acc
192 : real(qp), intent(in) :: num
193 0 : if (num < acc%value()) then
194 : ret = .true.
195 : else
196 0 : ret = .false.
197 : end if
198 0 : end function num_less_than_acc
199 :
200 : ! >>>>>>>>>>>
201 0 : logical function acc_greater_than_num(acc, num) result(ret)
202 : type(accurate_real), intent(in) :: acc
203 : real(qp), intent(in) :: num
204 0 : if (acc%value() > num) then
205 : ret = .true.
206 : else
207 0 : ret = .false.
208 : end if
209 0 : end function acc_greater_than_num
210 :
211 0 : logical function num_greater_than_acc(num, acc) result(ret)
212 : type(accurate_real), intent(in) :: acc
213 : real(qp), intent(in) :: num
214 0 : if (num > acc%value()) then
215 : ret = .true.
216 : else
217 0 : ret = .false.
218 : end if
219 0 : end function num_greater_than_acc
220 :
221 : ! +++++++++++
222 0 : type(accurate_real) function add_accurate_real(op1, op2) result(ret)
223 : type(accurate_real), intent(in) :: op1
224 : type(accurate_real), intent(in) :: op2
225 :
226 0 : ret%sum = op1%sum
227 0 : ret%compensator = op1%compensator
228 0 : call neumaier_sum(ret%sum, ret%compensator, op2%sum)
229 0 : call neumaier_sum(ret%sum, ret%compensator, op2%compensator)
230 :
231 0 : end function add_accurate_real
232 :
233 0 : type(accurate_real) function add_accurate_real_num(op1, op2) result(ret)
234 : type(accurate_real), intent(in) :: op1
235 : real(qp), intent(in) :: op2
236 :
237 0 : ret%sum = op1%sum
238 0 : ret%compensator = op1%compensator
239 0 : call neumaier_sum(ret%sum, ret%compensator, op2)
240 :
241 0 : end function add_accurate_real_num
242 :
243 0 : type(accurate_real) function add_accurate_real_num_rev(op1, op2) result(ret)
244 : real(qp), intent(in) :: op1
245 : type(accurate_real), intent(in) :: op2
246 :
247 0 : ret%sum = op2%sum
248 0 : ret%compensator = op2%compensator
249 0 : call neumaier_sum(ret%sum, ret%compensator, op1)
250 :
251 0 : end function add_accurate_real_num_rev
252 :
253 : ! -----------
254 0 : type(accurate_real) function sub_accurate_real(op1, op2) result(ret)
255 : type(accurate_real), intent(in) :: op1
256 : type(accurate_real), intent(in) :: op2
257 :
258 0 : ret%sum = op1%sum
259 0 : ret%compensator = op1%compensator
260 0 : call neumaier_sum(ret%sum, ret%compensator, -op2%sum)
261 0 : call neumaier_sum(ret%sum, ret%compensator, -op2%compensator)
262 :
263 0 : end function sub_accurate_real
264 :
265 0 : type(accurate_real) function sub_accurate_real_num(op1, op2) result(ret)
266 : type(accurate_real), intent(in) :: op1
267 : real(qp), intent(in) :: op2
268 :
269 0 : ret%sum = op1%sum
270 0 : ret%compensator = op1%compensator
271 0 : call neumaier_sum(ret%sum, ret%compensator, -op2)
272 :
273 0 : end function sub_accurate_real_num
274 :
275 0 : type(accurate_real) function sub_accurate_real_num_rev(op1, op2) result(ret)
276 : real(qp), intent(in) :: op1
277 : type(accurate_real), intent(in) :: op2
278 :
279 0 : ret%sum = op1
280 0 : ret%compensator = -op2%compensator
281 0 : call neumaier_sum(ret%sum, ret%compensator, -op2%sum)
282 :
283 0 : end function sub_accurate_real_num_rev
284 :
285 0 : end module accurate_sum
|
