Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementation of procedures in [pm_mathRoot](@ref pm_mathRoot).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #define CHECK_LF(LF) \
28 : CHECK_ASSERTION(__LINE__, abs(LF - lf) < abstol, SK_"@setRoot(): The condition `abs(getFunc(lb), lf) < abstol` must hold. lb, getFunc(lb), lf, abstol = "//getStr([lb, LF, lf, abstol]))
29 : #define CHECK_UF(UF) \
30 : CHECK_ASSERTION(__LINE__, abs(UF - uf) < abstol, SK_"@setRoot(): The condition `abs(getFunc(ub), uf) < abstol` must hold. ub, getFunc(ub), uf, abstol = "//getStr([ub, UF, uf, abstol]))
31 : #define CHECK_BRACKET \
32 : CHECK_ASSERTION(__LINE__, sign(lf, 1._RKC) /= sign(uf, 1._RKC), \
33 : SK_"@setRoot(): The condition `sign(lf, 1.) /= sign(uf, 1.)` must hold. lf, uf = "//getStr([lf, uf])) ! fpp
34 : #define CHECK_LB_LT_UB \
35 : CHECK_ASSERTION(__LINE__, lb < ub, SK_"@setRoot(): The condition `lb < ub` must hold. lb, ub = "//getStr([lb, ub])) ! fpp
36 : #define CHECK_ZERO_LT_ABSTOL \
37 : CHECK_ASSERTION(__LINE__, 0._RKC < abstol, SK_"@setRoot(): The condition `0. < abstol` must hold. abstol = "//getStr(abstol)) ! fpp
38 : #define CHECK_ABSTOL_LT_LB_UB_DIFF \
39 : CHECK_ASSERTION(__LINE__, abstol < ub - lb, \
40 : SK_"@setRoot(): The condition `abstol < ub - lb` must hold. abstol, ub, lb = "//getStr([abstol, ub, lb])) ! fpp
41 : real(RKC), parameter :: EPS10 = 10 * epsilon(0._RKC)
42 : #if Fixed_ENABLED
43 : #define CHECK_ZERO_LT_NITER
44 : integer(IK), parameter :: NITER = ceiling(50 * precision(0._RKC) / log10(2._RKC))
45 : #elif Niter_ENABLED
46 : #define CHECK_ZERO_LT_NITER \
47 : CHECK_ASSERTION(__LINE__, 0_IK < niter, SK_"@setRoot(): The condition `0 < niter` must hold. niter = "//getStr(niter)) ! fpp
48 : #elif !getRoot_ENABLED
49 : #error "Unrecognized interface."
50 : #endif
51 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 : #if getRoot_ENABLED && Def_ENABLED
53 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54 :
55 96006 : root = getRoot(brent, getFunc, lb, ub, abstol, neval, niter)
56 :
57 : !%%%%%%%%%%%%%%
58 : #elif getRoot_ENABLED
59 : !%%%%%%%%%%%%%%
60 :
61 : real(RKC) :: abstol_def, lf, uf
62 : integer(IK) :: neval_def
63 137633 : if (present(abstol)) then
64 68801 : abstol_def = abstol
65 : else
66 68832 : abstol_def = epsilon(0._RKC)**.8 * (abs(lb) + abs(ub))
67 : end if
68 : #if Newton_ENABLED || Halley_ENABLED || Schroder_ENABLED
69 19209 : lf = getFunc(lb, 0_IK)
70 19209 : uf = getFunc(ub, 0_IK)
71 19209 : if (present(init)) then
72 9600 : root = init
73 : else
74 9609 : root = 0.5_RKC * (lb + ub)
75 : end if
76 : #elif False_ENABLED || Bisection_ENABLED || Secant_ENABLED || Brent_ENABLED || Ridders_ENABLED || TOMS748_ENABLED
77 118424 : lf = getFunc(lb)
78 118424 : uf = getFunc(ub)
79 : #else
80 : #error "Unrecognized interface."
81 : #endif
82 137633 : if (present(niter)) then
83 100800 : call setRoot(method, getFunc, root, lb, ub, lf, uf, abstol_def, neval_def, niter)
84 : else
85 36833 : call setRoot(method, getFunc, root, lb, ub, lf, uf, abstol_def, neval_def)
86 : end if
87 137633 : if (present(neval)) then
88 100830 : neval = neval_def + 2_IK
89 36803 : elseif (neval_def < 0_IK) then
90 : error stop MODULE_NAME//SK_"@getRoot(): Failed to converge in search for the root of the specified function." ! LCOV_EXCL_LINE
91 : end if
92 :
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 : #elif setRoot_ENABLED && False_ENABLED
95 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 :
97 : real(RKC) :: diff, interval, func
98 4006 : if (abs(lf) <= EPS10) then
99 906 : neval = 0_IK
100 906 : root = lb
101 3100 : elseif (abs(uf) <= EPS10) then
102 982 : neval = 0_IK
103 982 : root = ub
104 2118 : elseif (abs(ub - lb) < abstol) then
105 0 : root = .5_RKC * (ub + lb)
106 0 : neval = 0_IK
107 : else
108 6354 : CHECK_BRACKET
109 6354 : CHECK_LB_LT_UB
110 2118 : CHECK_ZERO_LT_ABSTOL
111 10590 : CHECK_LF(getFunc(lb))
112 10590 : CHECK_UF(getFunc(ub))
113 8472 : CHECK_ABSTOL_LT_LB_UB_DIFF
114 2118 : if (0._RKC <= lf) then
115 1076 : diff = lb
116 1076 : lb = ub
117 1076 : ub = diff
118 : diff = lf
119 : lf = uf
120 : uf = diff
121 : end if
122 2118 : interval = ub - lb
123 56536 : do neval = 1_IK, niter
124 56536 : root = lb + interval * lf / (lf - uf)
125 56536 : func = getFunc(root)
126 56536 : if (func < 0._RKC) then
127 56468 : diff = lb - root
128 56468 : lb = root
129 : lf = func
130 : else
131 68 : diff = ub - root
132 68 : ub = root
133 : uf = func
134 : end if
135 56536 : interval = ub - lb
136 56536 : if (abs(diff) < abstol .or. func == 0._RKC) return
137 : end do
138 : ! Error occurred.
139 0 : root = .5_RKC * (lb + ub) ! for the sake of defining `root` on output (important for tests).
140 0 : neval = -neval
141 : end if
142 :
143 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144 : #elif setRoot_ENABLED && Bisection_ENABLED
145 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
146 :
147 : real(RKC) :: delta, center
148 4006 : if (abs(lf) <= EPS10) then
149 950 : neval = 0_IK
150 950 : root = lb
151 3056 : elseif (abs(uf) <= EPS10) then
152 1030 : neval = 0_IK
153 1030 : root = ub
154 2026 : elseif (abs(ub - lb) < abstol) then
155 0 : root = .5_RKC * (ub + lb)
156 0 : neval = 0_IK
157 : else
158 6078 : CHECK_BRACKET
159 6078 : CHECK_LB_LT_UB
160 1010 : CHECK_ZERO_LT_NITER
161 2026 : CHECK_ZERO_LT_ABSTOL
162 10130 : CHECK_LF(getFunc(lb))
163 10130 : CHECK_UF(getFunc(ub))
164 8104 : CHECK_ABSTOL_LT_LB_UB_DIFF
165 2026 : if (lf < 0._RKC) then
166 970 : root = lb
167 970 : delta = ub - lb
168 : else
169 1056 : root = ub
170 1056 : delta = lb - ub
171 : end if
172 78204 : do neval = 3_IK, NITER + 2_IK
173 78204 : delta = delta * 0.5_RKC
174 78204 : center = root + delta
175 78204 : uf = getFunc(center)
176 78204 : if (uf <= 0._RKC) root = center
177 78204 : if (abs(delta) < abstol .or. uf == 0._RKC) return
178 : end do
179 : ! Error occurred.
180 0 : root = .5_RKC * (lb + ub) ! for the sake of defining `root` on output (important for tests).
181 0 : neval = -neval
182 : end if
183 :
184 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185 : #elif setRoot_ENABLED && Secant_ENABLED
186 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
187 :
188 : real(RKC) :: x, delta
189 4006 : if (abs(lf) <= EPS10) then
190 1056 : neval = 0_IK
191 1056 : root = lb
192 2950 : elseif (abs(uf) <= EPS10) then
193 968 : neval = 0_IK
194 968 : root = ub
195 1982 : elseif (abs(ub - lb) < abstol) then
196 0 : root = .5_RKC * (ub + lb)
197 0 : neval = 0_IK
198 : else
199 5946 : CHECK_BRACKET
200 5946 : CHECK_LB_LT_UB
201 988 : CHECK_ZERO_LT_NITER
202 1982 : CHECK_ZERO_LT_ABSTOL
203 9910 : CHECK_LF(getFunc(lb))
204 9910 : CHECK_UF(getFunc(ub))
205 7928 : CHECK_ABSTOL_LT_LB_UB_DIFF
206 1982 : if (abs(lf) < abs(uf)) then
207 836 : root = lb
208 836 : x = ub
209 : delta = lf
210 : lf = uf
211 : uf = delta
212 : else
213 1146 : x = lb
214 1146 : root = ub
215 : end if
216 17206 : do neval = 1_IK, niter
217 17206 : delta = (x - root) * uf / (uf - lf)
218 : x = root
219 : lf = uf
220 17206 : root = root + delta
221 17206 : uf = getFunc(root)
222 17206 : if (abs(delta) < abstol .or. uf == 0._RKC) return
223 : end do
224 : ! Error occurred.
225 0 : root = .5_RKC * (lb + ub) ! for the sake of defining `root` on output (important for tests).
226 0 : neval = -neval
227 : end if
228 :
229 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
230 : #elif setRoot_ENABLED && Brent_ENABLED
231 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232 :
233 : real(RKC) :: halfAbsTol
234 : real(RKC) :: middle, interval, funMid, tol1, halfint, rnew, p, q, r, ratio
235 : real(RKC) , parameter :: EPS = epsilon(lb)
236 103212 : if (abs(lf) <= EPS10) then
237 26170 : neval = 0_IK
238 26170 : root = lb
239 77042 : elseif (abs(uf) <= EPS10) then
240 25518 : neval = 0_IK
241 25518 : root = ub
242 51524 : elseif (abs(ub - lb) < abstol) then
243 0 : root = .5_RKC * (ub + lb)
244 0 : neval = 0_IK
245 : else
246 154572 : CHECK_BRACKET
247 154572 : CHECK_LB_LT_UB
248 41740 : CHECK_ZERO_LT_NITER
249 51524 : CHECK_ZERO_LT_ABSTOL
250 257620 : CHECK_LF(getFunc(lb))
251 257620 : CHECK_UF(getFunc(ub))
252 206096 : CHECK_ABSTOL_LT_LB_UB_DIFF
253 51524 : halfAbsTol = 0.5_RKC * abstol
254 51524 : root = ub
255 : middle = ub
256 : funMid = uf
257 267808 : do neval = 0, niter
258 : ! adjust the search interval, if needed.
259 235840 : if ((funMid < 0._RKC .and. uf < 0._RKC) .or. (0._RKC < uf .and. 0._RKC < funMid)) then ! same sign
260 156781 : middle = lb
261 : funMid = lf
262 156781 : interval = root - lb
263 : rnew = interval
264 : end if
265 235840 : if (abs(funMid) < abs(uf)) then
266 39534 : lb = root
267 39534 : root = middle
268 : middle = lb
269 : lf = uf
270 : uf = funMid
271 : funMid = lf
272 : end if
273 235840 : halfint = .5_RKC * (middle - root)
274 235840 : tol1 = 2._RKC * EPS * abs(root) + halfAbsTol
275 235840 : if (abs(halfint) < tol1 .or. uf == 0._RKC) return
276 : ! See if a bisection is needed.
277 216284 : if ((abs(rnew) >= tol1) .and. (abs(lf) > abs(uf))) then
278 216284 : ratio = uf / lf ! 50
279 216284 : if (lb /= middle) then ! Try inverse quadratic interpolation.
280 68720 : q = lf / funMid
281 68720 : r = uf / funMid
282 68720 : p = ratio * (2._RKC * halfint * q * (q - r) - (root - lb) * (r - 1._RKC))
283 68720 : q = (q - 1._RKC) * (r - 1._RKC) * (ratio - 1._RKC)
284 : else ! Try linear interpolation.
285 147564 : p = 2._RKC * halfint * ratio
286 147564 : q = 1._RKC - ratio
287 : end if
288 216284 : if (0._RKC < p) then ! Check inbounds.
289 109119 : q = -q
290 : else
291 107165 : p = -p
292 : end if
293 216284 : if (((2._RKC * p) >= (3._RKC * halfint * q - abs(tol1 * q))) .or. (p >= abs(0.5_RKC * rnew * q))) then ! Interpolation failed, use bisection.
294 : interval = halfint
295 : rnew = interval
296 : else ! Accept interpolation.
297 : rnew = interval
298 199032 : interval = p / q
299 : end if
300 : else ! Bounds decreasing too slowly, use bisection.
301 : interval = halfint
302 : rnew = interval
303 : end if
304 : ! Move last best guess to a.
305 216284 : lb = root ! 110
306 : lf = uf
307 : ! Evaluate new trial root.
308 216284 : if (tol1 < abs(interval)) then
309 207104 : root = root + interval
310 9180 : elseif (0._RKC < halfint) then
311 4346 : root = root + tol1
312 : else
313 4834 : root = root - tol1
314 : end if
315 248252 : uf = getFunc(root)
316 : end do
317 : ! Error occurred.
318 31968 : root = .5_RKC * (lb + root) ! for the sake of defining `root` on output (important for tests).
319 31968 : neval = -neval
320 : end if
321 :
322 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
323 : #elif setRoot_ENABLED && Ridders_ENABLED
324 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
325 :
326 : real(RKC) :: halfint, rold, fold, func, multiplicity
327 4006 : neval = 0_IK
328 4006 : if (abs(lf) <= EPS10) then
329 1064 : root = lb
330 2942 : elseif (abs(uf) <= EPS10) then
331 1022 : root = ub
332 1920 : elseif (abs(ub - lb) < abstol) then
333 0 : root = .5_RKC * (ub + lb)
334 : else
335 5760 : CHECK_BRACKET
336 5760 : CHECK_LB_LT_UB
337 957 : CHECK_ZERO_LT_NITER
338 1920 : CHECK_ZERO_LT_ABSTOL
339 9600 : CHECK_LF(getFunc(lb))
340 9600 : CHECK_UF(getFunc(ub))
341 7680 : CHECK_ABSTOL_LT_LB_UB_DIFF
342 : !rnew = huge(rold)**.66
343 : do
344 12906 : halfint = 0.5_RKC * (ub - lb)
345 12906 : root = lb + halfint
346 12906 : func = getFunc(root)
347 12906 : neval = neval + 1_IK
348 12906 : if (niter < neval) exit
349 12906 : multiplicity = sqrt(func**2 - lf * uf)
350 : if (multiplicity == 0._RKC) return ! LCOV_EXCL_LINE ! This should rarely occur.
351 12906 : rold = root
352 : fold = func
353 12906 : root = rold + halfint * sign(fold, lf - uf) / multiplicity
354 12906 : if (abs(root - rold) < abstol) return ! This should never occur in the first iteration.
355 12582 : func = getFunc(root)
356 12582 : neval = neval + 1_IK
357 12582 : if (niter < neval) exit
358 : if (func == 0._RKC) return ! LCOV_EXCL_LINE
359 11112 : if (sign(fold, func) /= fold) then
360 7768 : lb = rold
361 : lf = fold
362 : uf = func
363 7768 : ub = root
364 3344 : elseif (sign(lf, func) /= lf) then
365 1644 : ub = root
366 : uf = func
367 1700 : elseif (sign(uf, func) /= uf) then
368 1700 : lb = root
369 : lf = func
370 : else
371 0 : error stop MODULE_NAME//SK_"@setRootRidders(): Internal library error detected. This condition should never occur."
372 : end if
373 11112 : if (abs(ub - lb) < abstol) exit
374 : end do
375 126 : if (niter < neval) then
376 0 : neval = -neval ! Error occurred.
377 : else
378 126 : root = .5_RKC * (lb + ub)
379 : end if
380 : end if
381 :
382 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
383 : #elif setRoot_ENABLED && TOMS748_ENABLED
384 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
385 :
386 : real(RKC) , parameter :: EPS5 = 5 * epsilon(0._RKC), SMALL = tiny(0._RKC) / epsilon(0._RKC) !* 32
387 : real(RKC) :: rmid, rupp, fupp, rold, fold, rnew, fnew, interval
388 4006 : neval = 0_IK
389 4006 : if (abs(lf) <= EPS10) then
390 988 : root = lb
391 3018 : elseif (abs(uf) <= EPS10) then
392 1004 : root = ub
393 2014 : elseif (abs(ub - lb) < abstol) then
394 0 : root = .5_RKC * (ub + lb)
395 : else
396 6042 : CHECK_BRACKET
397 6042 : CHECK_LB_LT_UB
398 1004 : CHECK_ZERO_LT_NITER
399 2014 : CHECK_ZERO_LT_ABSTOL
400 10070 : CHECK_LF(getFunc(lb))
401 10070 : CHECK_UF(getFunc(ub))
402 8056 : CHECK_ABSTOL_LT_LB_UB_DIFF
403 : !fnew = 1.e5_RKC
404 : !rnew = 1.e5_RKC
405 : !fold = 1.e5_RKC
406 : blockConvergence: block
407 : ! On the first step we take lb secant step.
408 2014 : rmid = lb - lf * (ub - lb) / (uf - lf)
409 2014 : if (rmid <= lb + abs(lb) * EPS5 .or. ub - abs(ub) * EPS5 <= rmid) rmid = .5_RKC * (lb + ub)
410 2014 : call setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
411 2014 : neval = neval + 1_IK
412 : ! Take a quadratic interpolation on the second step.
413 2014 : call setPolIntQuad(rmid, lb, ub, rold, lf, uf, fold, nstep = 2_IK)
414 2014 : rnew = rold
415 2014 : fnew = fold
416 2014 : neval = neval + 1_IK
417 2014 : call setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
418 7958 : loopConvergence: do
419 3522 : interval = ub - lb
420 3522 : if (niter < neval .or. lf == 0._RKC .or. abs(ub - lb) < abstol) exit blockConvergence
421 : ! Starting with the third step taken use either quadratic or cubic interpolation.
422 : ! Cubic interpolation requires that all four function values lf, uf, fold, and fnew are distinct.
423 : ! Should that not be the case, take a quadratic step instead.
424 : if ((abs(lf - uf) < SMALL) .or. & ! LCOV_EXCL_LINE
425 : (abs(lf - fold) < SMALL) .or. & ! LCOV_EXCL_LINE
426 : (abs(lf - fnew) < SMALL) .or. & ! LCOV_EXCL_LINE
427 : (abs(uf - fold) < SMALL) .or. & ! LCOV_EXCL_LINE
428 : (abs(uf - fnew) < SMALL) .or. & ! LCOV_EXCL_LINE
429 : (abs(fold - fnew) < SMALL)) then
430 0 : call setPolIntQuad(rmid, lb, ub, rold, lf, uf, fold, nstep = 2_IK)
431 : else
432 3520 : call setPolIntCubic(rmid, lb, ub, rold, rnew, lf, uf, fold, fnew)
433 : end if
434 : ! Re-bracket and check for termination.
435 3520 : rnew = rold
436 3520 : fnew = fold
437 3520 : neval = neval + 1_IK
438 3520 : call setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
439 3520 : if (niter < neval .or. lf == 0._RKC .or. abs(ub - lb) < abstol) exit blockConvergence
440 : ! Next interpolated step.
441 : if ((abs(lf - uf) < SMALL) .or. & ! LCOV_EXCL_LINE
442 : (abs(lf - fold) < SMALL) .or. & ! LCOV_EXCL_LINE
443 : (abs(lf - fnew) < SMALL) .or. & ! LCOV_EXCL_LINE
444 : (abs(uf - fold) < SMALL) .or. & ! LCOV_EXCL_LINE
445 : (abs(uf - fnew) < SMALL) .or. & ! LCOV_EXCL_LINE
446 : (abs(fold - fnew) < SMALL)) then
447 0 : call setPolIntQuad(rmid, lb, ub, rold, lf, uf, fold, nstep = 3_IK)
448 : else
449 2636 : call setPolIntCubic(rmid, lb, ub, rold, rnew, lf, uf, fold, fnew)
450 : end if
451 : ! Bracket, check termination condition, and update rnew.
452 2636 : neval = neval + 1_IK
453 2636 : call setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
454 2636 : if (niter < neval .or. lf == 0._RKC .or. abs(ub - lb) < abstol) exit blockConvergence
455 : ! Take a double-length secant step.
456 1802 : if (abs(lf) < abs(uf)) then
457 : rupp = lb
458 : fupp = lf
459 : else
460 : rupp = ub
461 : fupp = uf
462 : end if
463 1802 : rmid = rupp - 2 * fupp * (ub - lb) / (uf - lf)
464 1802 : if (abs(rmid - rupp) > .5_RKC * (ub - lb)) rmid = lb + .5_RKC * (ub - lb)
465 : ! Bracket and check termination condition.
466 1802 : rnew = rold
467 1802 : fnew = fold
468 1802 : neval = neval + 1_IK
469 1802 : call setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
470 1802 : if (niter < neval .or. lf == 0._RKC .or. abs(ub - lb) < abstol) exit blockConvergence
471 1508 : if (2 * (ub - lb) < interval) cycle
472 : ! Take an additional bisection step due to slow convergence.
473 0 : rnew = rold
474 0 : fnew = fold
475 0 : neval = neval + 1_IK
476 0 : call setBracket(getFunc, lb, ub, lb + .5_RKC * (ub - lb), lf, uf, rold, fold)
477 2014 : if (niter < neval .or. lf == 0._RKC .or. abs(ub - lb) < abstol) exit blockConvergence
478 : end do loopConvergence
479 : end block blockConvergence
480 2014 : if (niter < neval) neval = -neval
481 2014 : if (abs(lf) < EPS10) ub = lb
482 2014 : if (abs(uf) < EPS10) lb = ub
483 2014 : root = .5_RKC * (ub + lb)
484 : end if
485 :
486 : contains
487 :
488 11986 : subroutine setBracket(getFunc, lb, ub, rmid, lf, uf, rold, fold)
489 : ! Given a point `rmid` inside the existing enclosing interval `[lb, ub]`:
490 : ! 1. set lb = rmid if `getFunc(rmid) == 0`, otherwise,
491 : ! 2. find the new enclosing interval, either `[lb, rmid]` or `[rmid, ub]`,
492 : ! and set `rold` and `fold` to the point that has just been removed from the interval.
493 : ! In other words `rold` is the third best guess to the root.
494 : real(RKC) , parameter :: EPS2 = 2 * epsilon(0._RKC)
495 : real(RKC), intent(inout) :: lb, ub, lf, uf
496 : real(RKC), intent(out) :: rold, fold
497 : procedure(real(RKC)) :: getFunc
498 : real(RKC), value :: rmid
499 : real(RKC) :: fmid
500 : ! Adjust the location of `rmid` accordingly if the `[lb, ub]` is small, or `rmid` is too close to one interval end.
501 11986 : if (ub - lb < 2 * EPS2 * lb) then
502 0 : rmid = lb + .5_RKC * (ub - lb)
503 11986 : elseif (rmid <= lb + abs(lb) * EPS2) then
504 96 : rmid = lb + abs(lb) * EPS2
505 11890 : elseif (ub - abs(ub) * EPS2 <= rmid) then
506 114 : rmid = ub - abs(ub) * EPS2
507 : end if
508 11986 : fmid = getFunc(rmid)
509 11986 : if (fmid == 0._RKC) then
510 1528 : rold = 0._RKC
511 1528 : fold = 0._RKC
512 1528 : lf = 0._RKC
513 1528 : lb = rmid
514 1528 : return
515 : end if
516 : ! Update the interval.
517 10458 : if ((lf < 0._RKC .and. 0._RKC < fmid) .or. (fmid < 0._RKC .and. 0._RKC < lf)) then
518 5246 : rold = ub
519 5246 : fold = uf
520 5246 : ub = rmid
521 5246 : uf = fmid
522 : else
523 5212 : rold = lb
524 5212 : fold = lf
525 5212 : lb = rmid
526 5212 : lf = fmid
527 : end if
528 : end subroutine
529 :
530 : ! Perform quadratic interpolation to determine the next point,
531 : ! by taking `niter - neval` Newton steps to find the location of the quadratic polynomial.
532 : ! `rold`, the third best approximation to the root, after `lb` and `ub`, must lie outside of the interval `[lb, ub]`.
533 : ! Fall back to a secant step should the result be out of range.
534 : ! Start by obtaining the coefficients of the quadratic polynomial.
535 : ! Compute division without undue over/under-flow.
536 2616 : pure subroutine setPolIntQuad(rmid, lb, ub, rold, lf, uf, fold, nstep)
537 : real(RKC), parameter :: HUGE_RKC = huge(0._RKC)
538 : real(RKC), intent(in) :: lb, ub, rold, lf, uf, fold
539 : integer(IK), intent(in) :: nstep
540 : real(RKC), intent(out) :: rmid
541 : integer(IK) :: istep
542 : real(RKC) :: ratio1, ratio2, numer, denom, divres
543 : #define SET_DIV(DIVRES,NUMER,DENOM,OFLOW) \
544 : if (abs(DENOM) < 1._RKC) then; if (abs((DENOM) * HUGE_RKC) <= abs(NUMER)) then; \
545 : DIVRES = OFLOW; else; DIVRES = (NUMER) / (DENOM); end if; else; DIVRES = (NUMER) / (DENOM); end if;
546 2616 : SET_DIV(ratio1,uf - lf,ub - lb,HUGE_RKC)
547 2616 : SET_DIV(ratio2,fold - uf,rold - ub,HUGE_RKC)
548 2616 : SET_DIV(ratio2,ratio2 - ratio1,rold - lb,0._RKC)
549 2616 : if (ratio2 == 0._RKC) then ! Failed to determine coefficients. Try a secant step.
550 : call setPolIntSecant(rmid, lb, ub, lf, uf)
551 2 : return
552 : end if
553 : ! Determine the starting point of the Newton steps.
554 2614 : if ((ratio2 < 0._RKC .and. lf < 0._RKC) .or. (0._RKC < ratio2 .and. 0._RKC < lf)) then
555 1364 : rmid = lb
556 : else
557 1250 : rmid = ub
558 : end if
559 : ! Take the Newton steps.
560 8444 : do istep = 1, nstep
561 5830 : numer = lf + (ratio1 + ratio2 * (rmid - ub)) * (rmid - lb)
562 5830 : denom = ratio1 + ratio2 * (2 * rmid - lb - ub)
563 5830 : SET_DIV(divres,numer,denom,1 + rmid - lb)
564 8444 : rmid = rmid - divres
565 : end do
566 2614 : if (rmid <= lb .or. ub <= rmid) call setPolIntSecant(rmid, lb, ub, lf, uf) ! Failed. Try a secant step.
567 : #undef SET_DIV
568 : end subroutine
569 :
570 : ! Perform standard secant interpolation of `[lb, ub]` given function evaluations `(lf, uf)`.
571 : ! Perform a bisection if secant interpolation is very close to either `lb` or `ub`.
572 : ! This interpolation is needed when at least one other form of interpolation has already failed,
573 : ! implying that the function is unlikely to be smooth with a root near `lb` or `ub`.
574 : pure subroutine setPolIntSecant(rmid, lb, ub, lf, uf)
575 : real(RKC), intent(in) :: lb, ub, lf, uf
576 : real(RKC), intent(out) :: rmid
577 2 : rmid = lb - lf * (ub - lb) / (uf - lf)
578 2 : if (rmid <= lb + abs(lb) * EPS5 .or. ub - abs(ub) * EPS5 <= rmid) rmid = .5_RKC * (lb + ub)
579 : end subroutine
580 :
581 : ! Use inverse cubic interpolation of getFunc(x) at points `[lb, ub, rold, rnew]` to obtain an approximate root of getFunc(x).
582 : ! `rold` and `rnew` lie out of `[lb, ub]` and are the third and forth best approximations to the root that we have found so far.
583 : ! Fall back to quadratic interpolation in case of an erroneous result out of `[lb, ub]`.
584 6156 : pure subroutine setPolIntCubic(rmid, lb, ub, rold, rnew, lf, uf, fold, fnew)
585 : real(RKC), intent(in) :: lb, ub, rold, rnew, lf, uf, fold, fnew
586 : real(RKC), intent(out) :: rmid
587 : real(RKC) :: q11, q21, q31, d21, d31, q22, q32, d32, q33
588 6156 : q11 = (rold - rnew) * fold / (fnew - fold)
589 6156 : q21 = (ub - rold) * uf / (fold - uf)
590 6156 : q31 = (lb - ub) * lf / (uf - lf)
591 6156 : d21 = (ub - rold) * fold / (fold - uf)
592 6156 : d31 = (lb - ub) * uf / (uf - lf)
593 :
594 6156 : q22 = (d21 - q11) * uf / (fnew - uf)
595 6156 : q32 = (d31 - q21) * lf / (fold - lf)
596 6156 : d32 = (d31 - q21) * fold / (fold - lf)
597 6156 : q33 = (d32 - q22) * lf / (fnew - lf)
598 6156 : rmid = q31 + q32 + q33 + lb
599 6156 : if (rmid <= lb .or. ub <= rmid) call setPolIntQuad(rmid, lb, ub, rold, lf, uf, fold, nstep = 3_IK) ! Out of bounds step. Us quadratic interpolation.
600 6156 : end subroutine
601 :
602 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
603 : #elif setRoot_ENABLED && Newton_ENABLED
604 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
605 :
606 : integer(IK) :: iter
607 : real(RKC) :: func, grad, interval, func2grad, temp
608 7206 : neval = 0_IK
609 7206 : if (abs(lf) <= EPS10) then
610 2000 : root = lb
611 5206 : elseif (abs(uf) <= EPS10) then
612 1722 : root = ub
613 3484 : elseif (abs(ub - lb) < abstol) then
614 0 : root = .5_RKC * (ub + lb)
615 : else
616 10452 : CHECK_BRACKET
617 10452 : CHECK_LB_LT_UB
618 1739 : CHECK_ZERO_LT_NITER
619 3484 : CHECK_ZERO_LT_ABSTOL
620 17420 : CHECK_LF(getFunc(lb, 0_IK))
621 17420 : CHECK_UF(getFunc(ub, 0_IK))
622 13936 : CHECK_ABSTOL_LT_LB_UB_DIFF
623 3484 : if (root <= lb .or. ub <= root) root = .5_RKC * (lb + ub)
624 3484 : if (0._RKC <= lf) then
625 1464 : root = lb
626 1464 : lb = ub
627 1464 : ub = root
628 : end if
629 3484 : interval = abs(ub - lb)
630 : func2grad = interval
631 3484 : neval = neval + 1_IK
632 3484 : func = getFunc(root, 0_IK)
633 19016 : do iter = 1, niter
634 19016 : neval = neval + 1_IK
635 19016 : grad = getFunc(root, 1_IK)
636 19016 : if (0._RKC < ((root - ub) * grad - func) * ((root - lb) * grad - func) .or. abs(interval * grad) < abs(2._RKC * func)) then
637 : interval = func2grad
638 172 : func2grad = .5_RKC * (ub - lb)
639 172 : root = lb + func2grad
640 172 : if (lb == root) return
641 : else
642 : interval = func2grad
643 18844 : func2grad = func / grad
644 : temp = root
645 18844 : root = root - func2grad
646 18844 : if (temp == root) return
647 : end if
648 16680 : if (abs(func2grad) < abstol) return
649 15532 : neval = neval + 1_IK
650 15532 : func = getFunc(root, 0_IK)
651 15532 : if (func < 0._RKC) then
652 245 : lb = root
653 : else
654 15287 : ub = root
655 : end if
656 : end do
657 : ! Error occurred.
658 0 : root = .5_RKC * (lb + ub) ! for the sake of defining `root` on output (important for tests).
659 0 : neval = -neval
660 : end if
661 :
662 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
663 : #elif setRoot_ENABLED && (Halley_ENABLED || Schroder_ENABLED)
664 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
665 :
666 : logical(LK) :: out_of_bounds_sentry, varset
667 : real(RKC) :: func, grad, hess, rold, fold, fmax, fmin, delta, delta1, delta2, temp, diff, ratio
668 14412 : neval = 0_IK
669 14412 : if (abs(lf) <= EPS10) then
670 3496 : root = lb
671 10916 : elseif (abs(uf) <= EPS10) then
672 3738 : root = ub
673 7178 : elseif (abs(ub - lb) < abstol) then
674 0 : root = .5_RKC * (ub + lb)
675 : else
676 21534 : CHECK_BRACKET
677 21534 : CHECK_LB_LT_UB
678 3583 : CHECK_ZERO_LT_NITER
679 7178 : CHECK_ZERO_LT_ABSTOL
680 35890 : CHECK_LF(getFunc(lb, 0_IK))
681 35890 : CHECK_UF(getFunc(ub, 0_IK))
682 28712 : CHECK_ABSTOL_LT_LB_UB_DIFF
683 : !print *, "lb, root, ub", lb, root, ub
684 7178 : if (root <= lb .or. ub <= root) root = .5_RKC * (lb + ub)
685 : out_of_bounds_sentry = .false._LK
686 7178 : delta = huge(delta)
687 : fmax = 0._RKC
688 : fmin = 0._RKC
689 7178 : rold = root
690 7178 : func = 0._RKC
691 : delta1 = delta
692 : delta2 = delta
693 7178 : neval = 0_IK
694 : do
695 26752 : fold = func
696 : delta2 = delta1
697 26752 : delta1 = delta
698 26752 : neval = neval + 3_IK
699 26752 : if (niter < neval) exit
700 26752 : func = getFunc(root, 0_IK)
701 26752 : grad = getFunc(root, 1_IK)
702 26752 : hess = getFunc(root, 2_IK)
703 29006 : if (func == 0._RKC) return
704 21828 : if (grad == 0._RKC) then ! Handle zero derivative and hessian.
705 0 : if (fold == 0._RKC) then
706 : ! First iteration: pretend that we had a previous one at either `lb` or `ub`.
707 0 : if (root == lb) then
708 0 : rold = ub
709 : else
710 0 : rold = lb
711 : end if
712 0 : neval = neval + 1_IK
713 0 : if (niter < neval) exit
714 0 : fold = getFunc(rold, 0_IK)
715 0 : delta = rold - root
716 : end if
717 0 : if ((fold < 0._RKC .and. 0._RKC < func) .or. (func < 0._RKC .and. 0._RKC < fold)) then
718 : ! Crossed over. Move in opposite direction to last step.
719 0 : if (delta < 0._RKC) then
720 0 : delta = .5_RKC * (root - lb)
721 : else
722 0 : delta = .5_RKC * (root - ub)
723 : end if
724 : else ! Move in same direction as last step.
725 0 : if (delta < 0._RKC) then
726 0 : delta = .5_RKC * (root - ub)
727 : else
728 0 : delta = .5_RKC * (root - lb)
729 : end if
730 : end if
731 21828 : elseif (hess /= 0._RKC) then
732 : #if Halley_ENABLED || Schroder_ENABLED
733 21828 : ratio = 2 * func
734 21828 : temp = 2 * grad - func * (hess / grad)
735 21828 : varset = abs(temp) < 1._RKC
736 21828 : if (varset) then
737 0 : varset = abs(temp) * huge(temp) <= abs(ratio)
738 0 : if (varset) delta = func / grad ! Possible overflow, use Newton step.
739 : end if
740 21828 : if (.not. varset) delta = ratio / temp ! Use Halley step.
741 : #elif Schroder_ENABLED
742 : ! The Schroder stepper appears unstable for certain class of test problems where the Halley method does.
743 : ! For now, Schroder is switched to Halley, but the origins of this instability must be explored.
744 : ratio = func / grad
745 : varset = root /= 0._RKC
746 : if (varset) then
747 : varset = abs(ratio / root) < 0.1_RKC
748 : if (varset) then
749 : delta = ratio + (hess / (2 * grad)) * ratio**2
750 : ! Fall back to Newton iteration if hess contribution is larger than grad.
751 : varset = 0._RKC <= delta * ratio
752 : end if
753 : end if
754 : ! Fall back to Newton iteration.
755 : if (.not. varset) delta = ratio
756 : #else
757 : #error "Unrecognized interface."
758 : #endif
759 21828 : if (delta * grad / func < 0._RKC) then
760 : ! The Newton and Halley steps disagree about which way we should move,
761 : ! likely due to cancelation error in the calculation of the Halley step,
762 : ! or else the derivatives are so small that their values are basically trash.
763 : ! Move in the direction indicated by a Newton step, but by no more than twice the
764 : ! current rold value, otherwise the search can jump out of bounds.
765 0 : delta = func / grad
766 0 : if (2 * abs(rold) < abs(delta)) delta = sign(2._RKC, delta) * abs(rold)
767 : end if
768 : else
769 0 : delta = func / grad
770 : end if
771 21828 : temp = abs(delta / delta2)
772 21828 : if (0.8_RKC < temp .and. temp < 2._RKC) then
773 : ! The last two steps did not converge.
774 92 : if (0._RKC < delta) then
775 42 : delta = .5_RKC * (root - lb)
776 : else
777 50 : delta = .5_RKC * (root - ub)
778 : end if
779 92 : if (root /= 0._RKC .and. root < abs(delta)) delta = sign(.9_RKC, delta) * abs(root) ! protect against huge jumps.
780 : ! reset delta2 so that this branch will *not* be taken on the next iteration.
781 92 : delta1 = delta * 3._RKC
782 : delta2 = delta1
783 : end if
784 21828 : rold = root
785 21828 : root = root - delta
786 : ! Check for out of bounds step.
787 21828 : if (root < lb) then
788 22 : varset = abs(lb) < 1._RKC .and. 1._RKC < abs(root)
789 : if (varset) then
790 0 : varset = huge(root) / abs(root) < abs(lb)
791 0 : if (varset) diff = 1000._RKC
792 : end if
793 22 : if (.not. varset) then
794 22 : varset = abs(lb) < 1._RKC
795 22 : if (varset) then
796 0 : varset = huge(lb) * abs(lb) < abs(root)
797 0 : if (varset) then
798 0 : if (lb < 0._RKC .neqv. root < 0._RKC) then
799 : diff = -huge(diff)
800 : else
801 : diff = +huge(diff)
802 : end if
803 : end if
804 : end if
805 22 : if (.not. varset) diff = root / lb
806 : end if
807 22 : if (abs(diff) < 1._RKC) diff = 1 / diff
808 22 : if (0._RKC < diff .and. diff < 3._RKC .and. .not. out_of_bounds_sentry) then
809 : ! Only a small out of bounds step, lets assume that the root is probably approximately at `lb`.
810 : out_of_bounds_sentry = .true. ! only take this branch once!
811 22 : delta = 0.99_RKC * (rold - lb)
812 22 : root = rold - delta
813 : else
814 0 : root = .5_RKC * (lb + ub)
815 0 : if (abs(lb - ub) < 2 * spacing(root)) return
816 0 : call setBracketTowardMin(delta, getFunc, rold, func, lb, ub, neval, niter)
817 0 : if (niter < neval) exit
818 0 : root = rold - delta
819 0 : rold = lb
820 0 : cycle
821 : end if
822 21806 : elseif (ub < root) then
823 4 : varset = abs(ub) < 1._RKC .and. 1._RKC < abs(root)
824 : if (varset) then
825 0 : varset = huge(root) / abs(root) < abs(ub)
826 0 : if (varset) diff = 1000._RKC
827 : end if
828 4 : if (.not. varset) diff = root / ub
829 4 : if (abs(diff) < 1._RKC) diff = 1._RKC / diff
830 4 : if (0._RKC < diff .and. diff < 3._RKC .and. .not. out_of_bounds_sentry) then
831 : ! Only a small out of bounds step, assume that the root is approximately at `lb`.
832 : out_of_bounds_sentry = .true. ! Take this branch only once.
833 0 : delta = 0.99_RKC * (rold - ub)
834 0 : root = rold - delta
835 : else
836 4 : root = .5_RKC * (lb + ub)
837 4 : if (abs(lb - ub) < 2 * spacing(root)) return
838 4 : call setBracketTowardMax(delta, getFunc, rold, func, lb, ub, neval, niter)
839 4 : if (niter < neval) exit
840 4 : root = rold - delta
841 4 : rold = lb
842 4 : cycle
843 : end if
844 : end if
845 : ! update brackets:
846 21824 : if (0._RKC < delta) then
847 10755 : ub = rold
848 : fmax = func
849 : else
850 11069 : lb = rold
851 : fmin = func
852 : end if
853 : ! Sanity check that we bracket the rold:
854 21824 : if (0._RKC < fmax * fmin) exit ! no root, possibly a minimum around root.
855 21824 : if (abs(delta) < abs(root * abstol)) return
856 : end do
857 : ! Error occurred.
858 0 : root = .5_RKC * (lb + ub) ! for the sake of defining `root` on output (important for tests).
859 0 : neval = -neval
860 : end if
861 :
862 : contains
863 :
864 4 : subroutine setBracketTowardMax(delta, getFunc, rold, func, lb, ub, neval, niter)
865 : ! Move rold towards ub until we bracket the root, updating `lb` and `ub` on the fly.
866 : real(RKC), intent(inout) :: rold, lb, ub
867 : integer(IK), intent(inout) :: neval
868 : integer(IK), intent(in) :: niter
869 : real(RKC), intent(out) :: delta
870 : procedure(real(RKC)) :: getFunc
871 : real(RKC), intent(in) :: func
872 : real(RKC) :: guess0, funcurrent
873 : integer(IK) :: multiplier
874 : multiplier = 2_IK
875 4 : funcurrent = func
876 4 : guess0 = rold
877 4 : if (abs(lb) < abs(ub)) then
878 0 : loopSearch1: do
879 0 : if (funcurrent < 0._RKC .neqv. func < 0._RKC) exit loopSearch1
880 0 : lb = rold
881 0 : rold = rold * multiplier
882 0 : if (ub < rold) then
883 0 : rold = ub
884 0 : funcurrent = -funcurrent ! There must be a change of sign.
885 0 : exit loopSearch1
886 : end if
887 0 : neval = neval + 1_IK
888 0 : if (niter < neval) return
889 0 : funcurrent = getFunc(rold, 0_IK)
890 0 : multiplier = multiplier * 2_IK
891 : end do loopSearch1
892 : else ! If `lb` and `ub` are negative, divide to head towards `ub`.
893 3 : loopSearch2: do
894 7 : if (funcurrent < 0._RKC .neqv. func < 0._RKC) exit loopSearch2
895 4 : lb = rold
896 4 : rold = rold / multiplier
897 4 : if (ub < rold) then
898 1 : rold = ub
899 1 : funcurrent = -funcurrent ! There must be a change of sign.
900 1 : exit loopSearch2
901 : end if
902 3 : neval = neval + 1_IK
903 3 : if (niter < neval) return
904 3 : funcurrent = getFunc(rold, 0_IK)
905 6 : multiplier = multiplier * 2_IK
906 : end do loopSearch2
907 : end if
908 4 : if (neval < niter) then
909 4 : ub = rold
910 4 : if (16_IK < multiplier) then
911 0 : call setBracketTowardMin(delta, getFunc, rold, funcurrent, lb, ub, neval, niter)
912 0 : delta = delta + (guess0 - rold)
913 0 : return
914 : end if
915 : end if
916 4 : delta = guess0 - .5_RKC * (lb + ub)
917 : end subroutine
918 :
919 0 : subroutine setBracketTowardMin(delta, getFunc, rold, func, lb, ub, neval, niter)
920 : ! Move rold towards lb until we bracket the root, updating `lb` and `ub` on the fly.
921 : real(RKC), intent(inout) :: rold, lb, ub
922 : integer(IK), intent(inout) :: neval
923 : integer(IK), intent(in) :: niter
924 : real(RKC), intent(out) :: delta
925 : procedure(real(RKC)) :: getFunc
926 : real(RKC), intent(in) :: func
927 : real(RKC) :: guess0, funcurrent
928 : integer(IK) :: multiplier
929 : multiplier = 2_IK
930 0 : funcurrent = func
931 0 : guess0 = rold
932 0 : if (abs(lb) < abs(ub)) then
933 0 : loopSearch1: do
934 0 : if (funcurrent < 0._RKC .neqv. func < 0._RKC) exit loopSearch1
935 0 : ub = rold
936 0 : rold = rold / multiplier
937 0 : if (rold < lb) then
938 0 : rold = lb
939 0 : funcurrent = -funcurrent ! There must be a change of sign.
940 0 : exit loopSearch1
941 : end if
942 0 : neval = neval + 1_IK
943 0 : if (niter < neval) return
944 0 : funcurrent = getFunc(rold, 0_IK)
945 0 : multiplier = multiplier * 2_IK
946 : end do loopSearch1
947 : else ! If `lb` and `ub` are negative, multiply to head towards `lb`.
948 0 : loopSearch2: do
949 0 : if (funcurrent < 0._RKC .neqv. func < 0._RKC) exit loopSearch2
950 0 : ub = rold
951 0 : rold = rold * multiplier
952 0 : if (rold < lb) then
953 0 : rold = lb
954 0 : funcurrent = -funcurrent ! There must be a change of sign.
955 0 : exit loopSearch2
956 : end if
957 0 : neval = neval + 1_IK
958 0 : if (niter < neval) return
959 0 : funcurrent = getFunc(rold, 0_IK)
960 0 : multiplier = multiplier * 2_IK
961 : end do loopSearch2
962 : end if
963 0 : if (neval < niter) then
964 0 : lb = rold
965 0 : if (16_IK < multiplier) then
966 0 : call setBracketTowardMax(delta, getFunc, rold, funcurrent, lb, ub, neval, niter)
967 0 : delta = delta + (guess0 - rold)
968 0 : return
969 : end if
970 : end if
971 0 : delta = guess0 - .5_RKC * (lb + ub)
972 : end subroutine
973 :
974 : #else
975 : !%%%%%%%%%%%%%%%%%%%%%%%%
976 : #error "Unrecognized interface."
977 : !%%%%%%%%%%%%%%%%%%%%%%%%
978 : #endif
979 : #undef CHECK_ABSTOL_LT_LB_UB_DIFF
980 : #undef CHECK_ZERO_LT_ABSTOL
981 : #undef CHECK_ZERO_LT_NITER
982 : #undef CHECK_LB_LT_UB
983 : #undef CHECK_BRACKET
984 : #undef CHECK_LF
985 : #undef CHECK_UF
|
