Line data Source code
1 : ! Module : root_m
2 : ! Purpose : root finding algorithms
3 : !
4 : ! Copyright 2013-2025 Rich Townsend & The GYRE Team
5 : !
6 : ! This file is part of GYRE. GYRE is free software: you can
7 : ! redistribute it and/or modify it under the terms of the GNU General
8 : ! Public License as published by the Free Software Foundation, version 3.
9 : !
10 : ! GYRE is distributed in the hope that it will be useful, but WITHOUT
11 : ! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 : ! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
13 : ! License for more details.
14 : !
15 : ! You should have received a copy of the GNU General Public License
16 : ! along with this program. If not, see <http://www.gnu.org/licenses/>.
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 : module root_m
67 :
68 : ! Uses
69 :
70 : use forum_m, only: RD
71 :
72 : use ext_m
73 : use math_m
74 : use num_par_m
75 : use status_m
76 :
77 : use ISO_FORTRAN_ENV
78 :
79 : ! No implicit typing
80 :
81 : implicit none (type, external)
82 :
83 : ! Interfaces
84 :
85 :
86 : interface solve_root
87 : module procedure solve_root_r_
88 : end interface solve_root
89 :
90 : interface narrow_bracket
91 : module procedure narrow_bracket_r_
92 : end interface narrow_bracket
93 :
94 : interface expand_bracket
95 : module procedure expand_bracket_r_
96 : end interface expand_bracket
97 :
98 :
99 : interface solve_root
100 : module procedure solve_root_c_
101 : end interface solve_root
102 :
103 : interface narrow_bracket
104 : module procedure narrow_bracket_c_
105 : end interface narrow_bracket
106 :
107 : interface expand_bracket
108 : module procedure expand_bracket_c_
109 : end interface expand_bracket
110 :
111 :
112 : interface solve_root
113 : module procedure solve_root_xr_
114 : end interface solve_root
115 :
116 : interface narrow_bracket
117 : module procedure narrow_bracket_xr_
118 : end interface narrow_bracket
119 :
120 : interface expand_bracket
121 : module procedure expand_bracket_xr_
122 : end interface expand_bracket
123 :
124 :
125 : interface solve_root
126 : module procedure solve_root_xc_
127 : end interface solve_root
128 :
129 : interface narrow_bracket
130 : module procedure narrow_bracket_xc_
131 : end interface narrow_bracket
132 :
133 : interface expand_bracket
134 : module procedure expand_bracket_xc_
135 : end interface expand_bracket
136 :
137 :
138 : ! Access specifiers
139 :
140 : public :: solve_root
141 : public :: narrow_bracket
142 : public :: expand_bracket
143 :
144 : ! Default access
145 :
146 : private
147 :
148 : contains
149 :
150 :
151 72 : subroutine solve_root_r_(eval_func, x_a, x_b, x_tol, nm_p, x_root, stat, &
152 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
153 :
154 : interface
155 : subroutine eval_func(x, func, stat)
156 : use forum_m, only: RD
157 : use ext_m
158 : implicit none (type, external)
159 : real(RD), intent(in) :: x
160 : real(RD), intent(out) :: func
161 : integer, intent(out) :: stat
162 : end subroutine eval_func
163 : end interface
164 : real(RD), intent(in) :: x_a
165 : real(RD), intent(in) :: x_b
166 : real(RD), intent(in) :: x_tol
167 : class(num_par_t), intent(in) :: nm_p
168 : real(RD), intent(out) :: x_root
169 : integer, intent(out) :: stat
170 : integer, optional, intent(out) :: n_iter
171 : integer, optional, intent(in) :: n_iter_max
172 : logical, optional, intent(in) :: relative_tol
173 : real(RD), optional, intent(in) :: f_x_a
174 : real(RD), optional, intent(in) :: f_x_b
175 :
176 : real(RD) :: a
177 : real(RD) :: b
178 : real(RD) :: f_a
179 : real(RD) :: f_b
180 :
181 : ! Starting from the bracket [x_a,x_b], solve for a root of the
182 : ! function
183 :
184 36 : a = x_a
185 36 : b = x_b
186 :
187 36 : if (PRESENT(f_x_a)) then
188 0 : f_a = f_x_a
189 : else
190 36 : call eval_func(a, f_a, stat)
191 36 : if (stat /= STAT_OK) return
192 : endif
193 :
194 36 : if (PRESENT(f_x_b)) then
195 0 : f_b = f_x_b
196 : else
197 36 : call eval_func(b, f_b, stat)
198 36 : if (stat /= STAT_OK) return
199 : endif
200 :
201 : call narrow_bracket_r_(eval_func, a, b, x_tol, nm_p, stat, &
202 36 : n_iter, n_iter_max, relative_tol, f_a, f_b)
203 :
204 36 : x_root = b
205 :
206 : ! Finish
207 :
208 36 : return
209 :
210 : end subroutine solve_root_r_
211 :
212 :
213 136 : subroutine solve_root_xr_(eval_func, x_a, x_b, x_tol, nm_p, x_root, stat, &
214 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
215 :
216 : interface
217 : subroutine eval_func(x, func, stat)
218 : use forum_m, only: RD
219 : use ext_m
220 : implicit none (type, external)
221 : type(ext_rt), intent(in) :: x
222 : type(ext_rt), intent(out) :: func
223 : integer, intent(out) :: stat
224 : end subroutine eval_func
225 : end interface
226 : type(ext_rt), intent(in) :: x_a
227 : type(ext_rt), intent(in) :: x_b
228 : type(ext_rt), intent(in) :: x_tol
229 : class(num_par_t), intent(in) :: nm_p
230 : type(ext_rt), intent(out) :: x_root
231 : integer, intent(out) :: stat
232 : integer, optional, intent(out) :: n_iter
233 : integer, optional, intent(in) :: n_iter_max
234 : logical, optional, intent(in) :: relative_tol
235 : type(ext_rt), optional, intent(in) :: f_x_a
236 : type(ext_rt), optional, intent(in) :: f_x_b
237 :
238 : type(ext_rt) :: a
239 : type(ext_rt) :: b
240 : type(ext_rt) :: f_a
241 : type(ext_rt) :: f_b
242 :
243 : ! Starting from the bracket [x_a,x_b], solve for a root of the
244 : ! function
245 :
246 68 : a = x_a
247 68 : b = x_b
248 :
249 68 : if (PRESENT(f_x_a)) then
250 68 : f_a = f_x_a
251 : else
252 0 : call eval_func(a, f_a, stat)
253 0 : if (stat /= STAT_OK) return
254 : endif
255 :
256 68 : if (PRESENT(f_x_b)) then
257 68 : f_b = f_x_b
258 : else
259 0 : call eval_func(b, f_b, stat)
260 0 : if (stat /= STAT_OK) return
261 : endif
262 :
263 : call narrow_bracket_xr_(eval_func, a, b, x_tol, nm_p, stat, &
264 68 : n_iter, n_iter_max, relative_tol, f_a, f_b)
265 :
266 68 : x_root = b
267 :
268 : ! Finish
269 :
270 68 : return
271 :
272 : end subroutine solve_root_xr_
273 :
274 :
275 : !****
276 :
277 :
278 0 : subroutine solve_root_c_(eval_func, z_a, z_b, z_tol, nm_p, z_root, stat, &
279 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
280 :
281 : interface
282 : subroutine eval_func(x, func, stat)
283 : use forum_m, only: RD
284 : use ext_m
285 : implicit none (type, external)
286 : complex(RD), intent(in) :: x
287 : complex(RD), intent(out) :: func
288 : integer, intent(out) :: stat
289 : end subroutine eval_func
290 : end interface
291 : complex(RD), intent(in) :: z_a
292 : complex(RD), intent(in) :: z_b
293 : real(RD), intent(in) :: z_tol
294 : class(num_par_t), intent(in) :: nm_p
295 : complex(RD), intent(out) :: z_root
296 : integer, intent(out) :: stat
297 : integer, optional, intent(out) :: n_iter
298 : integer, optional, intent(in) :: n_iter_max
299 : logical, optional, intent(in) :: relative_tol
300 : complex(RD), optional, intent(in) :: f_z_a
301 : complex(RD), optional, intent(in) :: f_z_b
302 :
303 : complex(RD) :: a
304 : complex(RD) :: b
305 : complex(RD) :: f_a
306 : complex(RD) :: f_b
307 :
308 : ! Starting from the bracket [z_a,z_b], solve for a root of the
309 : ! function
310 :
311 0 : a = z_a
312 0 : b = z_b
313 :
314 0 : if (PRESENT(f_z_a)) then
315 0 : f_a = f_z_a
316 : else
317 0 : call eval_func(a, f_a, stat)
318 0 : if (stat /= STAT_OK) return
319 : endif
320 :
321 0 : if (PRESENT(f_z_b)) then
322 0 : f_b = f_z_b
323 : else
324 0 : call eval_func(b, f_b, stat)
325 0 : if (stat /= STAT_OK) return
326 : endif
327 :
328 : call narrow_bracket_c_(eval_func, a, b, z_tol, nm_p, stat, &
329 0 : n_iter, n_iter_max, relative_tol, f_a, f_b)
330 :
331 0 : z_root = b
332 :
333 : ! Finish
334 :
335 0 : return
336 :
337 : end subroutine solve_root_c_
338 :
339 :
340 0 : subroutine solve_root_xc_(eval_func, z_a, z_b, z_tol, nm_p, z_root, stat, &
341 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
342 :
343 : interface
344 : subroutine eval_func(x, func, stat)
345 : use forum_m, only: RD
346 : use ext_m
347 : implicit none (type, external)
348 : type(ext_ct), intent(in) :: x
349 : type(ext_ct), intent(out) :: func
350 : integer, intent(out) :: stat
351 : end subroutine eval_func
352 : end interface
353 : type(ext_ct), intent(in) :: z_a
354 : type(ext_ct), intent(in) :: z_b
355 : type(ext_rt), intent(in) :: z_tol
356 : class(num_par_t), intent(in) :: nm_p
357 : type(ext_ct), intent(out) :: z_root
358 : integer, intent(out) :: stat
359 : integer, optional, intent(out) :: n_iter
360 : integer, optional, intent(in) :: n_iter_max
361 : logical, optional, intent(in) :: relative_tol
362 : type(ext_ct), optional, intent(in) :: f_z_a
363 : type(ext_ct), optional, intent(in) :: f_z_b
364 :
365 : type(ext_ct) :: a
366 : type(ext_ct) :: b
367 : type(ext_ct) :: f_a
368 : type(ext_ct) :: f_b
369 :
370 : ! Starting from the bracket [z_a,z_b], solve for a root of the
371 : ! function
372 :
373 0 : a = z_a
374 0 : b = z_b
375 :
376 0 : if (PRESENT(f_z_a)) then
377 0 : f_a = f_z_a
378 : else
379 0 : call eval_func(a, f_a, stat)
380 0 : if (stat /= STAT_OK) return
381 : endif
382 :
383 0 : if (PRESENT(f_z_b)) then
384 0 : f_b = f_z_b
385 : else
386 0 : call eval_func(b, f_b, stat)
387 0 : if (stat /= STAT_OK) return
388 : endif
389 :
390 : call narrow_bracket_xc_(eval_func, a, b, z_tol, nm_p, stat, &
391 0 : n_iter, n_iter_max, relative_tol, f_a, f_b)
392 :
393 0 : z_root = b
394 :
395 : ! Finish
396 :
397 0 : return
398 :
399 : end subroutine solve_root_xc_
400 :
401 :
402 : !****
403 :
404 :
405 36 : subroutine narrow_bracket_r_(eval_func, x_a, x_b, x_tol, nm_p, stat, &
406 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
407 :
408 : interface
409 : subroutine eval_func (x, func, stat)
410 : use forum_m, only: RD
411 : use ext_m
412 : implicit none (type, external)
413 : real(RD), intent(in) :: x
414 : real(RD), intent(out) :: func
415 : integer, intent(out) :: stat
416 : end subroutine eval_func
417 : end interface
418 : real(RD), intent(inout) :: x_a
419 : real(RD), intent(inout) :: x_b
420 : real(RD), intent(in) :: x_tol
421 : class(num_par_t), intent(in) :: nm_p
422 : integer, intent(out) :: stat
423 : integer, optional, intent(out) :: n_iter
424 : integer, optional, intent(in) :: n_iter_max
425 : logical, optional, intent(in) :: relative_tol
426 : real(RD), optional, intent(inout) :: f_x_a
427 : real(RD), optional, intent(inout) :: f_x_b
428 :
429 : ! Narrow the bracket [x_a,x_b] on a root of the function
430 :
431 72 : select case (nm_p%r_root_solver)
432 : case ('BRENT')
433 36 : call narrow_brent_r_(eval_func, x_a, x_b, x_tol, stat, n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
434 : case default
435 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 233 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
436 0 : write(UNIT=ERROR_UNIT, FMT=*) 'invalid r_root_solver'
437 36 : error stop
438 : end select
439 :
440 : ! Finish
441 :
442 36 : return
443 :
444 : end subroutine narrow_bracket_r_
445 :
446 : !****
447 :
448 36 : subroutine narrow_brent_r_(eval_func, x_a, x_b, x_tol, stat, &
449 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
450 :
451 : interface
452 : subroutine eval_func(x, func, stat)
453 : use forum_m, only: RD
454 : use ext_m
455 : implicit none (type, external)
456 : real(RD), intent(in) :: x
457 : real(RD), intent(out) :: func
458 : integer, intent(out) :: stat
459 : end subroutine eval_func
460 : end interface
461 : real(RD), intent(inout) :: x_a
462 : real(RD), intent(inout) :: x_b
463 : real(RD), intent(in) :: x_tol
464 : integer, intent(out) :: stat
465 : integer, optional, intent(out) :: n_iter
466 : integer, optional, intent(in) :: n_iter_max
467 : logical, optional, intent(in) :: relative_tol
468 : real(RD), optional, intent(inout) :: f_x_a
469 : real(RD), optional, intent(inout) :: f_x_b
470 :
471 : logical :: relative_tol_
472 : real(RD) :: a
473 : real(RD) :: b
474 : real(RD) :: c
475 : real(RD) :: d
476 : real(RD) :: e
477 : real(RD) :: f_a
478 : real(RD) :: f_b
479 : real(RD) :: f_c
480 : real(RD) :: tol
481 : real(RD) :: m
482 : real(RD) :: p
483 : real(RD) :: q
484 : real(RD) :: r
485 : real(RD) :: s
486 : integer :: i_iter
487 :
488 36 : if (PRESENT(relative_tol)) then
489 0 : relative_tol_ = relative_tol
490 : else
491 : relative_tol_ = .FALSE.
492 : endif
493 :
494 : ! Narrow the bracket [x_a,x_b] on a root of the function
495 : ! using Brent's method [based on the ALGOL 60 routine 'zero'
496 : ! published in [Bre1973]
497 :
498 : ! Set up the initial state
499 :
500 36 : a = x_a
501 36 : b = x_b
502 :
503 36 : if (PRESENT(f_x_a)) then
504 36 : f_a = f_x_a
505 : else
506 0 : call eval_func(a, f_a, stat)
507 0 : if (stat /= STAT_OK) return
508 : endif
509 :
510 36 : if (PRESENT(f_x_b)) then
511 36 : f_b = f_x_b
512 : else
513 0 : call eval_func(b, f_b, stat)
514 0 : if (stat /= STAT_OK) return
515 : endif
516 :
517 36 : c = b
518 36 : f_c = f_b
519 :
520 : ! Confirm that the bracket contains a root
521 :
522 36 : if (.NOT. ((f_a >= 0._RD .AND. f_b <= 0._RD) .OR. (f_a <= 0._RD .AND. f_b >= 0._RD))) then
523 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 318 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
524 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion (f_a >= 0._RD .AND. f_b <= 0._RD) .OR. (f_a <= 0._RD .AND. f_b >= 0._RD) failed&
525 0 : & with message "'//'root is not bracketed'//'"'
526 0 : error stop
527 : end if
528 :
529 : ! Iterate until convergence to the desired tolerance, or the
530 : ! maximum number of iterations is exceeded
531 :
532 36 : i_iter = 0
533 :
534 36 : stat = STAT_OK
535 :
536 212 : iterate_loop : do
537 :
538 248 : i_iter = i_iter + 1
539 :
540 248 : if (PRESENT(n_iter_max)) then
541 0 : if (i_iter > n_iter_max) then
542 0 : stat = STAT_ITER_MAX
543 0 : exit iterate_loop
544 : endif
545 : endif
546 :
547 : ! Reorder c so that it has the opposite sign to b
548 :
549 248 : if (f_b > 0._RD .EQV. f_c > 0._RD) then
550 168 : c = a
551 168 : f_c = f_a
552 168 : d = b - a
553 168 : e = d
554 : endif
555 :
556 : ! Make sure that the function is smallest in magnitude
557 : ! at b
558 :
559 248 : if (abs(f_c) < abs(f_b)) then
560 56 : a = b
561 56 : b = c
562 56 : c = a
563 56 : f_a = f_b
564 56 : f_b = f_c
565 56 : f_c = f_a
566 : endif
567 :
568 248 : if (relative_tol_) then
569 0 : tol = (2._RD*EPSILON(0._RD) + x_tol)*abs(b)
570 : else
571 248 : tol = 2._RD*EPSILON(0._RD)*abs(b) + x_tol
572 : endif
573 :
574 248 : m = 0.5_RD*(c - b)
575 :
576 : ! Check for convergence
577 :
578 248 : if (abs(m) <= tol .OR. f_b == 0._RD) then
579 : exit iterate_loop
580 : endif
581 :
582 : ! See if bisection is forced
583 :
584 212 : if (abs(e) < tol .OR. abs(f_a) < abs(f_b)) then
585 :
586 : d = m
587 : e = d
588 :
589 : else
590 :
591 212 : s = f_b/f_a
592 :
593 212 : if (a == c) then
594 :
595 : ! Linear interpolation
596 :
597 132 : p = 2._RD*m*s
598 132 : q = 1._RD - s
599 :
600 : else
601 :
602 : ! Inverse quadratic interpolation
603 :
604 80 : q = f_a/f_c
605 80 : r = f_b/f_c
606 :
607 80 : p = s*(2._RD*m*q*(q - r) - (b - a)*(r - 1._RD))
608 80 : q = (q - 1._RD)*(r - 1._RD)*(s - 1._RD)
609 :
610 : endif
611 :
612 212 : if (p > 0._RD) then
613 132 : q = -q
614 : else
615 80 : p = -p
616 : endif
617 :
618 212 : s = e
619 212 : e = d
620 :
621 212 : if (2._RD*p < 3._RD*m*q - abs(tol*q) .AND. p < abs(0.5_RD*s*q)) then
622 208 : d = p/q
623 : else
624 : d = m
625 : e = d
626 : endif
627 :
628 : endif
629 :
630 : ! Store the old value of b in a
631 :
632 212 : a = b
633 212 : f_a = f_b
634 :
635 : ! Update b
636 :
637 212 : if (abs(d) > tol) then
638 176 : b = b + d
639 : else
640 36 : if(m > 0._RD) then
641 14 : b = b + tol
642 : else
643 22 : b = b - tol
644 : endif
645 : endif
646 :
647 212 : call eval_func(b, f_b, stat)
648 248 : if (stat /= STAT_OK) exit iterate_loop
649 :
650 : end do iterate_loop
651 :
652 : ! Store the results
653 :
654 36 : x_a = a
655 36 : x_b = b
656 :
657 36 : if (PRESENT(n_iter)) then
658 0 : n_iter = i_iter
659 : end if
660 :
661 36 : if (PRESENT(f_x_a)) f_x_a = f_a
662 36 : if (PRESENT(f_x_b)) f_x_b = f_b
663 :
664 : ! Finish
665 :
666 : return
667 :
668 : end subroutine narrow_brent_r_
669 :
670 : !****
671 :
672 0 : subroutine expand_bracket_r_(eval_func, x_a, x_b, stat, clamp_a, clamp_b, f_x_a, f_x_b)
673 :
674 : interface
675 : subroutine eval_func(x, func, stat)
676 : use forum_m, only: RD
677 : use ext_m
678 : implicit none (type, external)
679 : real(RD), intent(in) :: x
680 : real(RD), intent(out) :: func
681 : integer, intent(out) :: stat
682 : end subroutine eval_func
683 : end interface
684 : real(RD), intent(inout) :: x_a
685 : real(RD), intent(inout) :: x_b
686 : integer, intent(out) :: stat
687 : logical, optional, intent(in) :: clamp_a
688 : logical, optional, intent(in) :: clamp_b
689 : real(RD), optional, intent(out) :: f_x_a
690 : real(RD), optional, intent(out) :: f_x_b
691 :
692 : logical :: clamp_a_
693 : logical :: clamp_b_
694 : real(RD) :: f_a
695 : real(RD) :: f_b
696 : logical :: move_a
697 :
698 0 : if (PRESENT(clamp_a)) then
699 0 : clamp_a_ = clamp_a
700 : else
701 : clamp_a_ = .FALSE.
702 : endif
703 :
704 0 : if (PRESENT(clamp_b)) then
705 0 : clamp_b_ = clamp_b
706 : else
707 : clamp_b_ = .FALSE.
708 : endif
709 :
710 0 : if (.NOT. (.NOT. (clamp_a_ .AND. clamp_b_))) then
711 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 501 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
712 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion .NOT. (clamp_a_ .AND. clamp_b_) failed with message "'//'cannot clamp both&
713 0 : & points'//'"'
714 0 : error stop
715 : end if
716 :
717 0 : if (.NOT. (x_a /= x_b)) then
718 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 503 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
719 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion x_a /= x_b failed with message "'//'invalid initial bracket'//'"'
720 0 : error stop
721 : end if
722 :
723 : ! Expand the bracket [x_a,x_b] until it contains a root of the
724 : ! function
725 :
726 0 : call eval_func(x_a, f_a, stat)
727 0 : if (stat /= STAT_OK) return
728 :
729 0 : call eval_func(x_b, f_b, stat)
730 0 : if (stat /= STAT_OK) return
731 :
732 0 : stat = STAT_OK
733 :
734 : expand_loop : do
735 :
736 0 : if ((f_a > 0._RD .AND. f_b < 0._RD) .OR. &
737 : (f_a < 0._RD .AND. f_b > 0._RD)) exit expand_loop
738 :
739 0 : if (clamp_a_) then
740 : move_a = .FALSE.
741 0 : elseif (clamp_b_) then
742 : move_a = .TRUE.
743 : else
744 0 : move_a = abs(f_b) > abs(f_a)
745 : endif
746 :
747 0 : if (move_a) then
748 :
749 0 : x_a = x_a + GOLDEN_R*(x_a - x_b)
750 :
751 0 : call eval_func(x_a, f_a, stat)
752 0 : if (stat /= STAT_OK) exit expand_loop
753 :
754 : else
755 :
756 0 : x_b = x_b + GOLDEN_R*(x_b - x_a)
757 :
758 0 : call eval_func(x_b, f_b, stat)
759 0 : if (stat /= STAT_OK) exit expand_loop
760 :
761 : endif
762 :
763 : end do expand_loop
764 :
765 : ! Store the results
766 :
767 0 : if (PRESENT(f_x_a)) f_x_a = f_a
768 0 : if (PRESENT(f_x_b)) f_x_b = f_b
769 :
770 : ! Finish
771 :
772 : return
773 :
774 : end subroutine expand_bracket_r_
775 :
776 :
777 68 : subroutine narrow_bracket_xr_(eval_func, x_a, x_b, x_tol, nm_p, stat, &
778 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
779 :
780 : interface
781 : subroutine eval_func (x, func, stat)
782 : use forum_m, only: RD
783 : use ext_m
784 : implicit none (type, external)
785 : type(ext_rt), intent(in) :: x
786 : type(ext_rt), intent(out) :: func
787 : integer, intent(out) :: stat
788 : end subroutine eval_func
789 : end interface
790 : type(ext_rt), intent(inout) :: x_a
791 : type(ext_rt), intent(inout) :: x_b
792 : type(ext_rt), intent(in) :: x_tol
793 : class(num_par_t), intent(in) :: nm_p
794 : integer, intent(out) :: stat
795 : integer, optional, intent(out) :: n_iter
796 : integer, optional, intent(in) :: n_iter_max
797 : logical, optional, intent(in) :: relative_tol
798 : type(ext_rt), optional, intent(inout) :: f_x_a
799 : type(ext_rt), optional, intent(inout) :: f_x_b
800 :
801 : ! Narrow the bracket [x_a,x_b] on a root of the function
802 :
803 136 : select case (nm_p%r_root_solver)
804 : case ('BRENT')
805 68 : call narrow_brent_xr_(eval_func, x_a, x_b, x_tol, stat, n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
806 : case default
807 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 233 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
808 0 : write(UNIT=ERROR_UNIT, FMT=*) 'invalid r_root_solver'
809 68 : error stop
810 : end select
811 :
812 : ! Finish
813 :
814 68 : return
815 :
816 : end subroutine narrow_bracket_xr_
817 :
818 : !****
819 :
820 68 : subroutine narrow_brent_xr_(eval_func, x_a, x_b, x_tol, stat, &
821 : n_iter, n_iter_max, relative_tol, f_x_a, f_x_b)
822 :
823 : interface
824 : subroutine eval_func(x, func, stat)
825 : use forum_m, only: RD
826 : use ext_m
827 : implicit none (type, external)
828 : type(ext_rt), intent(in) :: x
829 : type(ext_rt), intent(out) :: func
830 : integer, intent(out) :: stat
831 : end subroutine eval_func
832 : end interface
833 : type(ext_rt), intent(inout) :: x_a
834 : type(ext_rt), intent(inout) :: x_b
835 : type(ext_rt), intent(in) :: x_tol
836 : integer, intent(out) :: stat
837 : integer, optional, intent(out) :: n_iter
838 : integer, optional, intent(in) :: n_iter_max
839 : logical, optional, intent(in) :: relative_tol
840 : type(ext_rt), optional, intent(inout) :: f_x_a
841 : type(ext_rt), optional, intent(inout) :: f_x_b
842 :
843 : logical :: relative_tol_
844 : type(ext_rt) :: a
845 : type(ext_rt) :: b
846 : type(ext_rt) :: c
847 : type(ext_rt) :: d
848 : type(ext_rt) :: e
849 : type(ext_rt) :: f_a
850 : type(ext_rt) :: f_b
851 : type(ext_rt) :: f_c
852 : type(ext_rt) :: tol
853 : type(ext_rt) :: m
854 : type(ext_rt) :: p
855 : type(ext_rt) :: q
856 : type(ext_rt) :: r
857 : type(ext_rt) :: s
858 : integer :: i_iter
859 :
860 68 : if (PRESENT(relative_tol)) then
861 0 : relative_tol_ = relative_tol
862 : else
863 : relative_tol_ = .FALSE.
864 : endif
865 :
866 : ! Narrow the bracket [x_a,x_b] on a root of the function
867 : ! using Brent's method [based on the ALGOL 60 routine 'zero'
868 : ! published in [Bre1973]
869 :
870 : ! Set up the initial state
871 :
872 68 : a = x_a
873 68 : b = x_b
874 :
875 68 : if (PRESENT(f_x_a)) then
876 68 : f_a = f_x_a
877 : else
878 0 : call eval_func(a, f_a, stat)
879 0 : if (stat /= STAT_OK) return
880 : endif
881 :
882 68 : if (PRESENT(f_x_b)) then
883 68 : f_b = f_x_b
884 : else
885 0 : call eval_func(b, f_b, stat)
886 0 : if (stat /= STAT_OK) return
887 : endif
888 :
889 68 : c = b
890 68 : f_c = f_b
891 :
892 : ! Confirm that the bracket contains a root
893 :
894 68 : if (.NOT. ((f_a >= 0._RD .AND. f_b <= 0._RD) .OR. (f_a <= 0._RD .AND. f_b >= 0._RD))) then
895 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 318 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
896 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion (f_a >= 0._RD .AND. f_b <= 0._RD) .OR. (f_a <= 0._RD .AND. f_b >= 0._RD) failed&
897 0 : & with message "'//'root is not bracketed'//'"'
898 0 : error stop
899 : end if
900 :
901 : ! Iterate until convergence to the desired tolerance, or the
902 : ! maximum number of iterations is exceeded
903 :
904 68 : i_iter = 0
905 :
906 68 : stat = STAT_OK
907 :
908 368 : iterate_loop : do
909 :
910 436 : i_iter = i_iter + 1
911 :
912 436 : if (PRESENT(n_iter_max)) then
913 436 : if (i_iter > n_iter_max) then
914 0 : stat = STAT_ITER_MAX
915 0 : exit iterate_loop
916 : endif
917 : endif
918 :
919 : ! Reorder c so that it has the opposite sign to b
920 :
921 436 : if (f_b > 0._RD .EQV. f_c > 0._RD) then
922 332 : c = a
923 332 : f_c = f_a
924 332 : d = b - a
925 332 : e = d
926 : endif
927 :
928 : ! Make sure that the function is smallest in magnitude
929 : ! at b
930 :
931 436 : if (abs(f_c) < abs(f_b)) then
932 104 : a = b
933 104 : b = c
934 104 : c = a
935 104 : f_a = f_b
936 104 : f_b = f_c
937 104 : f_c = f_a
938 : endif
939 :
940 436 : if (relative_tol_) then
941 0 : tol = (2._RD*EPSILON(0._RD) + x_tol)*abs(b)
942 : else
943 436 : tol = 2._RD*EPSILON(0._RD)*abs(b) + x_tol
944 : endif
945 :
946 436 : m = 0.5_RD*(c - b)
947 :
948 : ! Check for convergence
949 :
950 436 : if (abs(m) <= tol .OR. f_b == 0._RD) then
951 : exit iterate_loop
952 : endif
953 :
954 : ! See if bisection is forced
955 :
956 368 : if (abs(e) < tol .OR. abs(f_a) < abs(f_b)) then
957 :
958 6 : d = m
959 6 : e = d
960 :
961 : else
962 :
963 362 : s = f_b/f_a
964 :
965 362 : if (a == c) then
966 :
967 : ! Linear interpolation
968 :
969 268 : p = 2._RD*m*s
970 268 : q = 1._RD - s
971 :
972 : else
973 :
974 : ! Inverse quadratic interpolation
975 :
976 94 : q = f_a/f_c
977 94 : r = f_b/f_c
978 :
979 94 : p = s*(2._RD*m*q*(q - r) - (b - a)*(r - 1._RD))
980 94 : q = (q - 1._RD)*(r - 1._RD)*(s - 1._RD)
981 :
982 : endif
983 :
984 362 : if (p > 0._RD) then
985 190 : q = -q
986 : else
987 172 : p = -p
988 : endif
989 :
990 362 : s = e
991 362 : e = d
992 :
993 362 : if (2._RD*p < 3._RD*m*q - abs(tol*q) .AND. p < abs(0.5_RD*s*q)) then
994 362 : d = p/q
995 : else
996 0 : d = m
997 0 : e = d
998 : endif
999 :
1000 : endif
1001 :
1002 : ! Store the old value of b in a
1003 :
1004 368 : a = b
1005 368 : f_a = f_b
1006 :
1007 : ! Update b
1008 :
1009 368 : if (abs(d) > tol) then
1010 304 : b = b + d
1011 : else
1012 64 : if(m > 0._RD) then
1013 32 : b = b + tol
1014 : else
1015 32 : b = b - tol
1016 : endif
1017 : endif
1018 :
1019 368 : call eval_func(b, f_b, stat)
1020 368 : if (stat /= STAT_OK) exit iterate_loop
1021 :
1022 : end do iterate_loop
1023 :
1024 : ! Store the results
1025 :
1026 68 : x_a = a
1027 68 : x_b = b
1028 :
1029 68 : if (PRESENT(n_iter)) then
1030 68 : n_iter = i_iter
1031 : end if
1032 :
1033 68 : if (PRESENT(f_x_a)) f_x_a = f_a
1034 68 : if (PRESENT(f_x_b)) f_x_b = f_b
1035 :
1036 : ! Finish
1037 :
1038 : return
1039 :
1040 : end subroutine narrow_brent_xr_
1041 :
1042 : !****
1043 :
1044 0 : subroutine expand_bracket_xr_(eval_func, x_a, x_b, stat, clamp_a, clamp_b, f_x_a, f_x_b)
1045 :
1046 : interface
1047 : subroutine eval_func(x, func, stat)
1048 : use forum_m, only: RD
1049 : use ext_m
1050 : implicit none (type, external)
1051 : type(ext_rt), intent(in) :: x
1052 : type(ext_rt), intent(out) :: func
1053 : integer, intent(out) :: stat
1054 : end subroutine eval_func
1055 : end interface
1056 : type(ext_rt), intent(inout) :: x_a
1057 : type(ext_rt), intent(inout) :: x_b
1058 : integer, intent(out) :: stat
1059 : logical, optional, intent(in) :: clamp_a
1060 : logical, optional, intent(in) :: clamp_b
1061 : type(ext_rt), optional, intent(out) :: f_x_a
1062 : type(ext_rt), optional, intent(out) :: f_x_b
1063 :
1064 : logical :: clamp_a_
1065 : logical :: clamp_b_
1066 : type(ext_rt) :: f_a
1067 : type(ext_rt) :: f_b
1068 : logical :: move_a
1069 :
1070 0 : if (PRESENT(clamp_a)) then
1071 0 : clamp_a_ = clamp_a
1072 : else
1073 : clamp_a_ = .FALSE.
1074 : endif
1075 :
1076 0 : if (PRESENT(clamp_b)) then
1077 0 : clamp_b_ = clamp_b
1078 : else
1079 : clamp_b_ = .FALSE.
1080 : endif
1081 :
1082 0 : if (.NOT. (.NOT. (clamp_a_ .AND. clamp_b_))) then
1083 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 501 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1084 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion .NOT. (clamp_a_ .AND. clamp_b_) failed with message "'//'cannot clamp both&
1085 0 : & points'//'"'
1086 0 : error stop
1087 : end if
1088 :
1089 0 : if (.NOT. (x_a /= x_b)) then
1090 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 503 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1091 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion x_a /= x_b failed with message "'//'invalid initial bracket'//'"'
1092 0 : error stop
1093 : end if
1094 :
1095 : ! Expand the bracket [x_a,x_b] until it contains a root of the
1096 : ! function
1097 :
1098 0 : call eval_func(x_a, f_a, stat)
1099 0 : if (stat /= STAT_OK) return
1100 :
1101 0 : call eval_func(x_b, f_b, stat)
1102 0 : if (stat /= STAT_OK) return
1103 :
1104 0 : stat = STAT_OK
1105 :
1106 : expand_loop : do
1107 :
1108 0 : if ((f_a > 0._RD .AND. f_b < 0._RD) .OR. &
1109 : (f_a < 0._RD .AND. f_b > 0._RD)) exit expand_loop
1110 :
1111 0 : if (clamp_a_) then
1112 : move_a = .FALSE.
1113 0 : elseif (clamp_b_) then
1114 : move_a = .TRUE.
1115 : else
1116 0 : move_a = abs(f_b) > abs(f_a)
1117 : endif
1118 :
1119 0 : if (move_a) then
1120 :
1121 0 : x_a = x_a + GOLDEN_R*(x_a - x_b)
1122 :
1123 0 : call eval_func(x_a, f_a, stat)
1124 0 : if (stat /= STAT_OK) exit expand_loop
1125 :
1126 : else
1127 :
1128 0 : x_b = x_b + GOLDEN_R*(x_b - x_a)
1129 :
1130 0 : call eval_func(x_b, f_b, stat)
1131 0 : if (stat /= STAT_OK) exit expand_loop
1132 :
1133 : endif
1134 :
1135 : end do expand_loop
1136 :
1137 : ! Store the results
1138 :
1139 0 : if (PRESENT(f_x_a)) f_x_a = f_a
1140 0 : if (PRESENT(f_x_b)) f_x_b = f_b
1141 :
1142 : ! Finish
1143 :
1144 : return
1145 :
1146 : end subroutine expand_bracket_xr_
1147 :
1148 :
1149 : !****
1150 :
1151 :
1152 0 : subroutine narrow_bracket_c_(eval_func, z_a, z_b, z_tol, nm_p, stat, &
1153 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1154 :
1155 : interface
1156 : subroutine eval_func (z, func, stat)
1157 : use forum_m, only: RD
1158 : use ext_m
1159 : implicit none (type, external)
1160 : complex(RD), intent(in) :: z
1161 : complex(RD), intent(out) :: func
1162 : integer, intent(out) :: stat
1163 : end subroutine eval_func
1164 : end interface
1165 : complex(RD), intent(inout) :: z_a
1166 : complex(RD), intent(inout) :: z_b
1167 : real(RD), intent(in) :: z_tol
1168 : class(num_par_t), intent(in) :: nm_p
1169 : integer, intent(out) :: stat
1170 : integer, optional, intent(out) :: n_iter
1171 : integer, optional, intent(in) :: n_iter_max
1172 : logical, optional, intent(in) :: relative_tol
1173 : complex(RD), optional, intent(inout) :: f_z_a
1174 : complex(RD), optional, intent(inout) :: f_z_b
1175 :
1176 : ! Narrow the bracket [z_a,z_b] toward a root of the function
1177 :
1178 0 : select case (nm_p%c_root_solver)
1179 : case ('SECANT')
1180 : call narrow_secant_c_(eval_func, z_a, z_b, z_tol, stat, &
1181 0 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1182 : case ('RIDDERS')
1183 : call narrow_ridders_c_(eval_func, z_a, z_b, z_tol, stat, &
1184 0 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1185 : case default
1186 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 598 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1187 0 : write(UNIT=ERROR_UNIT, FMT=*) 'invalid c_root_solver'
1188 0 : error stop
1189 : end select
1190 :
1191 : ! Finish
1192 :
1193 0 : return
1194 :
1195 : end subroutine narrow_bracket_c_
1196 :
1197 : !****
1198 :
1199 0 : subroutine narrow_secant_c_(eval_func, z_a, z_b, z_tol, stat, &
1200 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1201 :
1202 : interface
1203 : subroutine eval_func(z, func, stat)
1204 : use forum_m, only: RD
1205 : use ext_m
1206 : implicit none (type, external)
1207 : complex(RD), intent(in) :: z
1208 : complex(RD), intent(out) :: func
1209 : integer, intent(out) :: stat
1210 : end subroutine eval_func
1211 : end interface
1212 : complex(RD), intent(inout) :: z_a
1213 : complex(RD), intent(inout) :: z_b
1214 : real(RD), intent(in) :: z_tol
1215 : integer, intent(out) :: stat
1216 : integer, optional, intent(out) :: n_iter
1217 : integer, optional, intent(in) :: n_iter_max
1218 : logical, optional, intent(in) :: relative_tol
1219 : complex(RD), optional, intent(inout) :: f_z_a
1220 : complex(RD), optional, intent(inout) :: f_z_b
1221 :
1222 : logical :: relative_tol_
1223 : complex(RD) :: a
1224 : complex(RD) :: b
1225 : complex(RD) :: c
1226 : complex(RD) :: f_a
1227 : complex(RD) :: f_b
1228 : complex(RD) :: f_c
1229 : integer :: i_iter
1230 : complex(RD) :: f_dz
1231 : complex(RD) :: rho
1232 : real(RD) :: tol
1233 :
1234 0 : if (PRESENT(relative_tol)) then
1235 0 : relative_tol_ = relative_tol
1236 : else
1237 : relative_tol_ = .FALSE.
1238 : endif
1239 :
1240 : ! Narrow the bracket [z_a,z_b] toward a root of the function !
1241 : ! using the secant method
1242 :
1243 : ! Set up the initial state
1244 :
1245 0 : a = z_a
1246 0 : b = z_b
1247 :
1248 0 : if (PRESENT(f_z_a)) then
1249 0 : f_a = f_z_a
1250 : else
1251 0 : call eval_func(a, f_a, stat)
1252 0 : if (stat /= STAT_OK) return
1253 : endif
1254 :
1255 0 : if (PRESENT(f_z_b)) then
1256 0 : f_b = f_z_b
1257 : else
1258 0 : call eval_func(b, f_b, stat)
1259 0 : if (stat /= STAT_OK) return
1260 : endif
1261 :
1262 0 : if (abs(f_a) < abs(f_b)) then
1263 :
1264 0 : c = a
1265 0 : a = b
1266 0 : b = c
1267 :
1268 0 : f_c = f_a
1269 0 : f_a = f_b
1270 0 : f_b = f_c
1271 :
1272 : endif
1273 :
1274 : ! Iterate until convergence to the desired tolerance, or the
1275 : ! maximum number of iterations is exceeded
1276 :
1277 0 : i_iter = 0
1278 :
1279 0 : stat = STAT_OK
1280 :
1281 0 : iterate_loop : do
1282 :
1283 0 : if (f_b == 0._RD) exit iterate_loop
1284 :
1285 0 : i_iter = i_iter + 1
1286 :
1287 0 : if (PRESENT(n_iter_max)) then
1288 0 : if (i_iter > n_iter_max) then
1289 0 : stat = STAT_ITER_MAX
1290 0 : exit iterate_loop
1291 : endif
1292 : endif
1293 :
1294 : ! Calculate the correction
1295 :
1296 0 : f_dz = f_b*(b - a)
1297 :
1298 0 : rho = f_b - f_a
1299 :
1300 : ! Check for a singular correction
1301 :
1302 0 : if (abs(b*rho) < 8._RD*EPSILON(0._RD)*abs(f_dz)) then
1303 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 713 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1304 0 : write(UNIT=ERROR_UNIT, FMT=*) 'singular correction in secant'
1305 0 : error stop
1306 : endif
1307 :
1308 : ! Update the root
1309 :
1310 0 : a = b
1311 0 : f_a = f_b
1312 :
1313 0 : b = b - f_dz/rho
1314 0 : call eval_func(b, f_b, stat)
1315 0 : if (stat /= STAT_OK) exit iterate_loop
1316 :
1317 : ! Check for convergence
1318 :
1319 0 : if (relative_tol_) then
1320 0 : tol = (4._RD*EPSILON(0._RD) + z_tol)*abs(b)
1321 : else
1322 0 : tol = 4._RD*EPSILON(0._RD)*abs(b) + z_tol
1323 : endif
1324 :
1325 0 : if (abs(b - a) <= tol) exit iterate_loop
1326 :
1327 : end do iterate_loop
1328 :
1329 : ! Store the results
1330 :
1331 0 : z_a = a
1332 0 : z_b = b
1333 :
1334 0 : if (PRESENT(n_iter)) n_iter = i_iter
1335 :
1336 0 : if (PRESENT(f_z_a)) f_z_a = f_a
1337 0 : if (PRESENT(f_z_b)) f_z_b = f_b
1338 :
1339 : ! Finish
1340 :
1341 : end subroutine narrow_secant_c_
1342 :
1343 : !****
1344 :
1345 0 : subroutine narrow_ridders_c_(eval_func, z_a, z_b, z_tol, stat, &
1346 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1347 :
1348 : interface
1349 : subroutine eval_func (z, func, stat)
1350 : use forum_m, only: RD
1351 : use ext_m
1352 : implicit none (type, external)
1353 : complex(RD), intent(in) :: z
1354 : complex(RD), intent(out) :: func
1355 : integer, intent(out) :: stat
1356 : end subroutine eval_func
1357 : end interface
1358 : complex(RD), intent(inout) :: z_a
1359 : complex(RD), intent(inout) :: z_b
1360 : real(RD), intent(in) :: z_tol
1361 : integer, intent(out) :: stat
1362 : integer, optional, intent(out) :: n_iter
1363 : integer, optional, intent(in) :: n_iter_max
1364 : logical, optional, intent(in) :: relative_tol
1365 : complex(RD), optional, intent(inout) :: f_z_a
1366 : complex(RD), optional, intent(inout) :: f_z_b
1367 :
1368 : logical :: relative_tol_
1369 : complex(RD) :: a
1370 : complex(RD) :: b
1371 : complex(RD) :: c
1372 : complex(RD) :: f_a
1373 : complex(RD) :: f_b
1374 : complex(RD) :: f_c
1375 : integer :: i_iter
1376 : complex(RD) :: exp_Q_p
1377 : complex(RD) :: exp_Q_m
1378 : complex(RD) :: exp_Q
1379 : complex(RD) :: f_dz
1380 : complex(RD) :: rho
1381 : real(RD) :: tol
1382 :
1383 0 : if (PRESENT(relative_tol)) then
1384 0 : relative_tol_ = relative_tol
1385 : else
1386 : relative_tol_ = .FALSE.
1387 : endif
1388 :
1389 0 : if (.NOT. (z_a /= z_b)) then
1390 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 797 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1391 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion z_a /= z_b failed with message "'//'invalid initial bracket'//'"'
1392 0 : error stop
1393 : end if
1394 :
1395 : ! Narrow the bracket [z_a,z_b] toward a root of the function using
1396 : ! a complex Ridders' method (with secant updates, rather than
1397 : ! regula falsi)
1398 :
1399 : ! Set up the initial state
1400 :
1401 0 : a = z_a
1402 0 : b = z_b
1403 :
1404 0 : if (PRESENT(f_z_a)) then
1405 0 : f_a = f_z_a
1406 : else
1407 0 : call eval_func(a, f_a, stat)
1408 0 : if (stat /= STAT_OK) return
1409 : endif
1410 :
1411 0 : if (PRESENT(f_z_b)) then
1412 0 : f_b = f_z_b
1413 : else
1414 0 : call eval_func(b, f_b, stat)
1415 0 : if (stat /= STAT_OK) return
1416 : endif
1417 :
1418 0 : if (abs(f_a) < abs(f_b)) then
1419 :
1420 : c = a
1421 0 : a = b
1422 0 : b = c
1423 :
1424 0 : f_c = f_a
1425 0 : f_a = f_b
1426 0 : f_b = f_c
1427 :
1428 : endif
1429 :
1430 : ! Iterate until convergence to the desired tolerance, or the
1431 : ! maximum number of iterations is exceeded
1432 :
1433 0 : i_iter = 0
1434 :
1435 0 : stat = STAT_OK
1436 :
1437 0 : iterate_loop : do
1438 :
1439 0 : if (f_b == 0._RD) exit iterate_loop
1440 :
1441 0 : i_iter = i_iter + 1
1442 :
1443 0 : if (PRESENT(n_iter_max)) then
1444 0 : if (i_iter > n_iter_max) then
1445 0 : stat = STAT_ITER_MAX
1446 0 : exit iterate_loop
1447 : end if
1448 : endif
1449 :
1450 : ! Calculate the mid-point values
1451 :
1452 0 : c = 0.5_RD*(a + b)
1453 :
1454 0 : call eval_func(c, f_c, stat)
1455 0 : if (stat /= STAT_OK) exit iterate_loop
1456 :
1457 : ! Solve for the re-scaling exponential
1458 :
1459 0 : exp_Q_p = (f_c + sqrt(f_c*f_c - f_a*f_b))/f_b
1460 0 : exp_Q_m = (f_c - sqrt(f_c*f_c - f_a*f_b))/f_b
1461 :
1462 0 : if (abs(exp_Q_p-1._RD) < abs(exp_Q_m-1._RD)) then
1463 : exp_Q = exp_Q_p
1464 : else
1465 0 : exp_Q = exp_Q_m
1466 : endif
1467 :
1468 : ! Apply the secant method to the re-scaled problem
1469 :
1470 0 : f_dz = f_b*(exp_Q*exp_Q)*(b - a)
1471 :
1472 0 : rho = f_b*(exp_Q*exp_Q) - f_a
1473 :
1474 : ! Check for a singular correction
1475 :
1476 0 : if (abs(b*rho) < 8._RD*EPSILON(0._RD)*abs(f_dz)) then
1477 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 881 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1478 0 : write(UNIT=ERROR_UNIT, FMT=*) 'singular correction in secant'
1479 0 : error stop
1480 : endif
1481 :
1482 : ! Update the root
1483 :
1484 0 : a = b
1485 0 : f_a = f_b
1486 :
1487 0 : b = b - f_dz/rho
1488 0 : call eval_func(b, f_b, stat)
1489 0 : if (stat /= STAT_OK) exit iterate_loop
1490 :
1491 : ! Check for convergence
1492 :
1493 0 : if (relative_tol_) then
1494 0 : tol = (4._RD*EPSILON(0._RD) + z_tol)*abs(b)
1495 : else
1496 0 : tol = 4._RD*EPSILON(0._RD)*abs(b) + z_tol
1497 : endif
1498 :
1499 0 : if (abs(b - a) <= tol) exit iterate_loop
1500 :
1501 : end do iterate_loop
1502 :
1503 : ! Store the results
1504 :
1505 0 : z_a = a
1506 0 : z_b = b
1507 :
1508 0 : if (PRESENT(n_iter)) n_iter = i_iter
1509 :
1510 0 : if (PRESENT(f_z_a)) f_z_a = f_a
1511 0 : if (PRESENT(f_z_b)) f_z_b = f_b
1512 :
1513 : ! Finish
1514 :
1515 : return
1516 :
1517 : end subroutine narrow_ridders_c_
1518 :
1519 : !****
1520 :
1521 0 : subroutine expand_bracket_c_(eval_func, z_a, z_b, f_z_tol, stat, &
1522 : clamp_a, clamp_b, relative_tol, f_z_a, f_z_b)
1523 :
1524 : interface
1525 : subroutine eval_func(z, func, stat)
1526 : use forum_m, only: RD
1527 : use ext_m
1528 : implicit none (type, external)
1529 : complex(RD), intent(in) :: z
1530 : complex(RD), intent(out) :: func
1531 : integer, intent(out) :: stat
1532 : end subroutine eval_func
1533 : end interface
1534 : complex(RD), intent(inout) :: z_a
1535 : complex(RD), intent(inout) :: z_b
1536 : real(RD), intent(in) :: f_z_tol
1537 : integer, intent(out) :: stat
1538 : complex(RD), optional, intent(out) :: f_z_a
1539 : complex(RD), optional, intent(out) :: f_z_b
1540 : logical, optional, intent(in) :: clamp_a
1541 : logical, optional, intent(in) :: clamp_b
1542 : logical, optional, intent(in) :: relative_tol
1543 :
1544 : logical :: relative_tol_
1545 : logical :: clamp_a_
1546 : logical :: clamp_b_
1547 : complex(RD) :: f_a
1548 : complex(RD) :: f_b
1549 : real(RD) :: tol
1550 : logical :: move_a
1551 :
1552 0 : if (PRESENT(clamp_a)) then
1553 0 : clamp_a_ = clamp_a
1554 : else
1555 : clamp_a_ = .FALSE.
1556 : endif
1557 :
1558 0 : if (PRESENT(clamp_b)) then
1559 0 : clamp_b_ = clamp_b
1560 : else
1561 : clamp_b_ = .FALSE.
1562 : endif
1563 :
1564 0 : if (.NOT. (.NOT. (clamp_a_ .AND. clamp_b_))) then
1565 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 966 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1566 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion .NOT. (clamp_a_ .AND. clamp_b_) failed with message "'//'cannot clamp both&
1567 0 : & points'//'"'
1568 0 : error stop
1569 : end if
1570 :
1571 0 : if (PRESENT(relative_tol)) then
1572 0 : relative_tol_ = relative_tol
1573 : else
1574 : relative_tol_ = .FALSE.
1575 : endif
1576 :
1577 0 : if (.NOT. (z_a /= z_b)) then
1578 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 974 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1579 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion z_a /= z_b failed with message "'//'invalid initial bracket'//'"'
1580 0 : error stop
1581 : end if
1582 :
1583 : ! Expand the bracket [z_a,z_b] until the difference between f(z_a)
1584 : ! and f(z_b) exceeds the tolerance
1585 :
1586 0 : call eval_func(z_a, f_a, stat)
1587 0 : if (stat /= STAT_OK) return
1588 :
1589 0 : call eval_func(z_b, f_b, stat)
1590 0 : if (stat /= STAT_OK) return
1591 :
1592 0 : stat = STAT_OK
1593 :
1594 : expand_loop : do
1595 :
1596 0 : if (relative_tol_) then
1597 0 : tol = (4._RD*EPSILON(0._RD) + f_z_tol)*max(abs(f_a), abs(f_b))
1598 : else
1599 0 : tol = 4._RD*EPSILON(0._RD)*max(abs(f_a), abs(f_b)) + f_z_tol
1600 : endif
1601 :
1602 0 : if (abs(f_a - f_b) > tol) exit expand_loop
1603 :
1604 0 : if (clamp_a_) then
1605 : move_a = .FALSE.
1606 0 : elseif (clamp_b_) then
1607 : move_a = .TRUE.
1608 : else
1609 0 : move_a = abs(f_b) > abs(f_a)
1610 : endif
1611 :
1612 0 : if (move_a) then
1613 0 : z_a = z_a + GOLDEN_R*(z_a - z_b)
1614 0 : call eval_func(z_a, f_a, stat)
1615 0 : if (stat /= STAT_OK) exit expand_loop
1616 : else
1617 0 : z_b = z_b + GOLDEN_R*(z_b - z_a)
1618 0 : call eval_func(z_b, f_b, stat)
1619 0 : if (stat /= STAT_OK) exit expand_loop
1620 : endif
1621 :
1622 : end do expand_loop
1623 :
1624 : ! Store f_a and f_b
1625 :
1626 0 : if (PRESENT(f_z_a)) f_z_a = f_a
1627 0 : if (PRESENT(f_z_b)) f_z_b = f_b
1628 :
1629 : ! Finish
1630 :
1631 : return
1632 :
1633 : end subroutine expand_bracket_c_
1634 :
1635 :
1636 0 : subroutine narrow_bracket_xc_(eval_func, z_a, z_b, z_tol, nm_p, stat, &
1637 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1638 :
1639 : interface
1640 : subroutine eval_func (z, func, stat)
1641 : use forum_m, only: RD
1642 : use ext_m
1643 : implicit none (type, external)
1644 : type(ext_ct), intent(in) :: z
1645 : type(ext_ct), intent(out) :: func
1646 : integer, intent(out) :: stat
1647 : end subroutine eval_func
1648 : end interface
1649 : type(ext_ct), intent(inout) :: z_a
1650 : type(ext_ct), intent(inout) :: z_b
1651 : type(ext_rt), intent(in) :: z_tol
1652 : class(num_par_t), intent(in) :: nm_p
1653 : integer, intent(out) :: stat
1654 : integer, optional, intent(out) :: n_iter
1655 : integer, optional, intent(in) :: n_iter_max
1656 : logical, optional, intent(in) :: relative_tol
1657 : type(ext_ct), optional, intent(inout) :: f_z_a
1658 : type(ext_ct), optional, intent(inout) :: f_z_b
1659 :
1660 : ! Narrow the bracket [z_a,z_b] toward a root of the function
1661 :
1662 0 : select case (nm_p%c_root_solver)
1663 : case ('SECANT')
1664 : call narrow_secant_xc_(eval_func, z_a, z_b, z_tol, stat, &
1665 0 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1666 : case ('RIDDERS')
1667 : call narrow_ridders_xc_(eval_func, z_a, z_b, z_tol, stat, &
1668 0 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1669 : case default
1670 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 598 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1671 0 : write(UNIT=ERROR_UNIT, FMT=*) 'invalid c_root_solver'
1672 0 : error stop
1673 : end select
1674 :
1675 : ! Finish
1676 :
1677 0 : return
1678 :
1679 : end subroutine narrow_bracket_xc_
1680 :
1681 : !****
1682 :
1683 0 : subroutine narrow_secant_xc_(eval_func, z_a, z_b, z_tol, stat, &
1684 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1685 :
1686 : interface
1687 : subroutine eval_func(z, func, stat)
1688 : use forum_m, only: RD
1689 : use ext_m
1690 : implicit none (type, external)
1691 : type(ext_ct), intent(in) :: z
1692 : type(ext_ct), intent(out) :: func
1693 : integer, intent(out) :: stat
1694 : end subroutine eval_func
1695 : end interface
1696 : type(ext_ct), intent(inout) :: z_a
1697 : type(ext_ct), intent(inout) :: z_b
1698 : type(ext_rt), intent(in) :: z_tol
1699 : integer, intent(out) :: stat
1700 : integer, optional, intent(out) :: n_iter
1701 : integer, optional, intent(in) :: n_iter_max
1702 : logical, optional, intent(in) :: relative_tol
1703 : type(ext_ct), optional, intent(inout) :: f_z_a
1704 : type(ext_ct), optional, intent(inout) :: f_z_b
1705 :
1706 : logical :: relative_tol_
1707 : type(ext_ct) :: a
1708 : type(ext_ct) :: b
1709 : type(ext_ct) :: c
1710 : type(ext_ct) :: f_a
1711 : type(ext_ct) :: f_b
1712 : type(ext_ct) :: f_c
1713 : integer :: i_iter
1714 : type(ext_ct) :: f_dz
1715 : type(ext_ct) :: rho
1716 : type(ext_rt) :: tol
1717 :
1718 0 : if (PRESENT(relative_tol)) then
1719 0 : relative_tol_ = relative_tol
1720 : else
1721 : relative_tol_ = .FALSE.
1722 : endif
1723 :
1724 : ! Narrow the bracket [z_a,z_b] toward a root of the function !
1725 : ! using the secant method
1726 :
1727 : ! Set up the initial state
1728 :
1729 0 : a = z_a
1730 0 : b = z_b
1731 :
1732 0 : if (PRESENT(f_z_a)) then
1733 0 : f_a = f_z_a
1734 : else
1735 0 : call eval_func(a, f_a, stat)
1736 0 : if (stat /= STAT_OK) return
1737 : endif
1738 :
1739 0 : if (PRESENT(f_z_b)) then
1740 0 : f_b = f_z_b
1741 : else
1742 0 : call eval_func(b, f_b, stat)
1743 0 : if (stat /= STAT_OK) return
1744 : endif
1745 :
1746 0 : if (abs(f_a) < abs(f_b)) then
1747 :
1748 0 : c = a
1749 0 : a = b
1750 0 : b = c
1751 :
1752 0 : f_c = f_a
1753 0 : f_a = f_b
1754 0 : f_b = f_c
1755 :
1756 : endif
1757 :
1758 : ! Iterate until convergence to the desired tolerance, or the
1759 : ! maximum number of iterations is exceeded
1760 :
1761 0 : i_iter = 0
1762 :
1763 0 : stat = STAT_OK
1764 :
1765 0 : iterate_loop : do
1766 :
1767 0 : if (f_b == 0._RD) exit iterate_loop
1768 :
1769 0 : i_iter = i_iter + 1
1770 :
1771 0 : if (PRESENT(n_iter_max)) then
1772 0 : if (i_iter > n_iter_max) then
1773 0 : stat = STAT_ITER_MAX
1774 0 : exit iterate_loop
1775 : endif
1776 : endif
1777 :
1778 : ! Calculate the correction
1779 :
1780 0 : f_dz = f_b*(b - a)
1781 :
1782 0 : rho = f_b - f_a
1783 :
1784 : ! Check for a singular correction
1785 :
1786 0 : if (abs(b*rho) < 8._RD*EPSILON(0._RD)*abs(f_dz)) then
1787 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 713 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1788 0 : write(UNIT=ERROR_UNIT, FMT=*) 'singular correction in secant'
1789 0 : error stop
1790 : endif
1791 :
1792 : ! Update the root
1793 :
1794 0 : a = b
1795 0 : f_a = f_b
1796 :
1797 0 : b = b - f_dz/rho
1798 0 : call eval_func(b, f_b, stat)
1799 0 : if (stat /= STAT_OK) exit iterate_loop
1800 :
1801 : ! Check for convergence
1802 :
1803 0 : if (relative_tol_) then
1804 0 : tol = (4._RD*EPSILON(0._RD) + z_tol)*abs(b)
1805 : else
1806 0 : tol = 4._RD*EPSILON(0._RD)*abs(b) + z_tol
1807 : endif
1808 :
1809 0 : if (abs(b - a) <= tol) exit iterate_loop
1810 :
1811 : end do iterate_loop
1812 :
1813 : ! Store the results
1814 :
1815 0 : z_a = a
1816 0 : z_b = b
1817 :
1818 0 : if (PRESENT(n_iter)) n_iter = i_iter
1819 :
1820 0 : if (PRESENT(f_z_a)) f_z_a = f_a
1821 0 : if (PRESENT(f_z_b)) f_z_b = f_b
1822 :
1823 : ! Finish
1824 :
1825 : end subroutine narrow_secant_xc_
1826 :
1827 : !****
1828 :
1829 0 : subroutine narrow_ridders_xc_(eval_func, z_a, z_b, z_tol, stat, &
1830 : n_iter, n_iter_max, relative_tol, f_z_a, f_z_b)
1831 :
1832 : interface
1833 : subroutine eval_func (z, func, stat)
1834 : use forum_m, only: RD
1835 : use ext_m
1836 : implicit none (type, external)
1837 : type(ext_ct), intent(in) :: z
1838 : type(ext_ct), intent(out) :: func
1839 : integer, intent(out) :: stat
1840 : end subroutine eval_func
1841 : end interface
1842 : type(ext_ct), intent(inout) :: z_a
1843 : type(ext_ct), intent(inout) :: z_b
1844 : type(ext_rt), intent(in) :: z_tol
1845 : integer, intent(out) :: stat
1846 : integer, optional, intent(out) :: n_iter
1847 : integer, optional, intent(in) :: n_iter_max
1848 : logical, optional, intent(in) :: relative_tol
1849 : type(ext_ct), optional, intent(inout) :: f_z_a
1850 : type(ext_ct), optional, intent(inout) :: f_z_b
1851 :
1852 : logical :: relative_tol_
1853 : type(ext_ct) :: a
1854 : type(ext_ct) :: b
1855 : type(ext_ct) :: c
1856 : type(ext_ct) :: f_a
1857 : type(ext_ct) :: f_b
1858 : type(ext_ct) :: f_c
1859 : integer :: i_iter
1860 : type(ext_ct) :: exp_Q_p
1861 : type(ext_ct) :: exp_Q_m
1862 : type(ext_ct) :: exp_Q
1863 : type(ext_ct) :: f_dz
1864 : type(ext_ct) :: rho
1865 : type(ext_rt) :: tol
1866 :
1867 0 : if (PRESENT(relative_tol)) then
1868 0 : relative_tol_ = relative_tol
1869 : else
1870 : relative_tol_ = .FALSE.
1871 : endif
1872 :
1873 0 : if (.NOT. (z_a /= z_b)) then
1874 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 797 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1875 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion z_a /= z_b failed with message "'//'invalid initial bracket'//'"'
1876 0 : error stop
1877 : end if
1878 :
1879 : ! Narrow the bracket [z_a,z_b] toward a root of the function using
1880 : ! a complex Ridders' method (with secant updates, rather than
1881 : ! regula falsi)
1882 :
1883 : ! Set up the initial state
1884 :
1885 0 : a = z_a
1886 0 : b = z_b
1887 :
1888 0 : if (PRESENT(f_z_a)) then
1889 0 : f_a = f_z_a
1890 : else
1891 0 : call eval_func(a, f_a, stat)
1892 0 : if (stat /= STAT_OK) return
1893 : endif
1894 :
1895 0 : if (PRESENT(f_z_b)) then
1896 0 : f_b = f_z_b
1897 : else
1898 0 : call eval_func(b, f_b, stat)
1899 0 : if (stat /= STAT_OK) return
1900 : endif
1901 :
1902 0 : if (abs(f_a) < abs(f_b)) then
1903 :
1904 0 : c = a
1905 0 : a = b
1906 0 : b = c
1907 :
1908 0 : f_c = f_a
1909 0 : f_a = f_b
1910 0 : f_b = f_c
1911 :
1912 : endif
1913 :
1914 : ! Iterate until convergence to the desired tolerance, or the
1915 : ! maximum number of iterations is exceeded
1916 :
1917 0 : i_iter = 0
1918 :
1919 0 : stat = STAT_OK
1920 :
1921 0 : iterate_loop : do
1922 :
1923 0 : if (f_b == 0._RD) exit iterate_loop
1924 :
1925 0 : i_iter = i_iter + 1
1926 :
1927 0 : if (PRESENT(n_iter_max)) then
1928 0 : if (i_iter > n_iter_max) then
1929 0 : stat = STAT_ITER_MAX
1930 0 : exit iterate_loop
1931 : end if
1932 : endif
1933 :
1934 : ! Calculate the mid-point values
1935 :
1936 0 : c = 0.5_RD*(a + b)
1937 :
1938 0 : call eval_func(c, f_c, stat)
1939 0 : if (stat /= STAT_OK) exit iterate_loop
1940 :
1941 : ! Solve for the re-scaling exponential
1942 :
1943 0 : exp_Q_p = (f_c + sqrt(f_c*f_c - f_a*f_b))/f_b
1944 0 : exp_Q_m = (f_c - sqrt(f_c*f_c - f_a*f_b))/f_b
1945 :
1946 0 : if (abs(exp_Q_p-1._RD) < abs(exp_Q_m-1._RD)) then
1947 0 : exp_Q = exp_Q_p
1948 : else
1949 0 : exp_Q = exp_Q_m
1950 : endif
1951 :
1952 : ! Apply the secant method to the re-scaled problem
1953 :
1954 0 : f_dz = f_b*(exp_Q*exp_Q)*(b - a)
1955 :
1956 0 : rho = f_b*(exp_Q*exp_Q) - f_a
1957 :
1958 : ! Check for a singular correction
1959 :
1960 0 : if (abs(b*rho) < 8._RD*EPSILON(0._RD)*abs(f_dz)) then
1961 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 881 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
1962 0 : write(UNIT=ERROR_UNIT, FMT=*) 'singular correction in secant'
1963 0 : error stop
1964 : endif
1965 :
1966 : ! Update the root
1967 :
1968 0 : a = b
1969 0 : f_a = f_b
1970 :
1971 0 : b = b - f_dz/rho
1972 0 : call eval_func(b, f_b, stat)
1973 0 : if (stat /= STAT_OK) exit iterate_loop
1974 :
1975 : ! Check for convergence
1976 :
1977 0 : if (relative_tol_) then
1978 0 : tol = (4._RD*EPSILON(0._RD) + z_tol)*abs(b)
1979 : else
1980 0 : tol = 4._RD*EPSILON(0._RD)*abs(b) + z_tol
1981 : endif
1982 :
1983 0 : if (abs(b - a) <= tol) exit iterate_loop
1984 :
1985 : end do iterate_loop
1986 :
1987 : ! Store the results
1988 :
1989 0 : z_a = a
1990 0 : z_b = b
1991 :
1992 0 : if (PRESENT(n_iter)) n_iter = i_iter
1993 :
1994 0 : if (PRESENT(f_z_a)) f_z_a = f_a
1995 0 : if (PRESENT(f_z_b)) f_z_b = f_b
1996 :
1997 : ! Finish
1998 :
1999 : return
2000 :
2001 : end subroutine narrow_ridders_xc_
2002 :
2003 : !****
2004 :
2005 0 : subroutine expand_bracket_xc_(eval_func, z_a, z_b, f_z_tol, stat, &
2006 : clamp_a, clamp_b, relative_tol, f_z_a, f_z_b)
2007 :
2008 : interface
2009 : subroutine eval_func(z, func, stat)
2010 : use forum_m, only: RD
2011 : use ext_m
2012 : implicit none (type, external)
2013 : type(ext_ct), intent(in) :: z
2014 : type(ext_ct), intent(out) :: func
2015 : integer, intent(out) :: stat
2016 : end subroutine eval_func
2017 : end interface
2018 : type(ext_ct), intent(inout) :: z_a
2019 : type(ext_ct), intent(inout) :: z_b
2020 : type(ext_rt), intent(in) :: f_z_tol
2021 : integer, intent(out) :: stat
2022 : type(ext_ct), optional, intent(out) :: f_z_a
2023 : type(ext_ct), optional, intent(out) :: f_z_b
2024 : logical, optional, intent(in) :: clamp_a
2025 : logical, optional, intent(in) :: clamp_b
2026 : logical, optional, intent(in) :: relative_tol
2027 :
2028 : logical :: relative_tol_
2029 : logical :: clamp_a_
2030 : logical :: clamp_b_
2031 : type(ext_ct) :: f_a
2032 : type(ext_ct) :: f_b
2033 : type(ext_rt) :: tol
2034 : logical :: move_a
2035 :
2036 0 : if (PRESENT(clamp_a)) then
2037 0 : clamp_a_ = clamp_a
2038 : else
2039 : clamp_a_ = .FALSE.
2040 : endif
2041 :
2042 0 : if (PRESENT(clamp_b)) then
2043 0 : clamp_b_ = clamp_b
2044 : else
2045 : clamp_b_ = .FALSE.
2046 : endif
2047 :
2048 0 : if (.NOT. (.NOT. (clamp_a_ .AND. clamp_b_))) then
2049 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 966 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
2050 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion .NOT. (clamp_a_ .AND. clamp_b_) failed with message "'//'cannot clamp both&
2051 0 : & points'//'"'
2052 0 : error stop
2053 : end if
2054 :
2055 0 : if (PRESENT(relative_tol)) then
2056 0 : relative_tol_ = relative_tol
2057 : else
2058 : relative_tol_ = .FALSE.
2059 : endif
2060 :
2061 0 : if (.NOT. (z_a /= z_b)) then
2062 0 : write(UNIT=ERROR_UNIT, FMT=*) 'ABORT at line 974 of /home/runner/work/mesa/mesa/build/gyre/src/src/math/root_m.fypp'
2063 0 : write(UNIT=ERROR_UNIT, FMT=*) 'assertion z_a /= z_b failed with message "'//'invalid initial bracket'//'"'
2064 0 : error stop
2065 : end if
2066 :
2067 : ! Expand the bracket [z_a,z_b] until the difference between f(z_a)
2068 : ! and f(z_b) exceeds the tolerance
2069 :
2070 0 : call eval_func(z_a, f_a, stat)
2071 0 : if (stat /= STAT_OK) return
2072 :
2073 0 : call eval_func(z_b, f_b, stat)
2074 0 : if (stat /= STAT_OK) return
2075 :
2076 0 : stat = STAT_OK
2077 :
2078 : expand_loop : do
2079 :
2080 0 : if (relative_tol_) then
2081 0 : tol = (4._RD*EPSILON(0._RD) + f_z_tol)*max(abs(f_a), abs(f_b))
2082 : else
2083 0 : tol = 4._RD*EPSILON(0._RD)*max(abs(f_a), abs(f_b)) + f_z_tol
2084 : endif
2085 :
2086 0 : if (abs(f_a - f_b) > tol) exit expand_loop
2087 :
2088 0 : if (clamp_a_) then
2089 : move_a = .FALSE.
2090 0 : elseif (clamp_b_) then
2091 : move_a = .TRUE.
2092 : else
2093 0 : move_a = abs(f_b) > abs(f_a)
2094 : endif
2095 :
2096 0 : if (move_a) then
2097 0 : z_a = z_a + GOLDEN_R*(z_a - z_b)
2098 0 : call eval_func(z_a, f_a, stat)
2099 0 : if (stat /= STAT_OK) exit expand_loop
2100 : else
2101 0 : z_b = z_b + GOLDEN_R*(z_b - z_a)
2102 0 : call eval_func(z_b, f_b, stat)
2103 0 : if (stat /= STAT_OK) exit expand_loop
2104 : endif
2105 :
2106 : end do expand_loop
2107 :
2108 : ! Store f_a and f_b
2109 :
2110 0 : if (PRESENT(f_z_a)) f_z_a = f_a
2111 0 : if (PRESENT(f_z_b)) f_z_b = f_b
2112 :
2113 : ! Finish
2114 :
2115 : return
2116 :
2117 : end subroutine expand_bracket_xc_
2118 :
2119 :
2120 : end module root_m
|
