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