Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mod_root
21 : use const_def, only: dp, arg_not_provided
22 :
23 : implicit none
24 :
25 : contains
26 :
27 5 : real(dp) function do_safe_root_with_brackets(f, x1_in, x3_in, y1_in, y3_in, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
28 : use mod_brent, only: eval_brent_safe_zero
29 : interface
30 : ! f provides function values
31 : #include "num_root_fcn.dek"
32 : end interface
33 : real(dp), intent(in) :: x1_in, x3_in
34 : real(dp), intent(in) :: y1_in, y3_in ! f(x1) and f(x3)
35 : integer, intent(in) :: lipar, lrpar
36 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
37 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
38 : integer, intent(in) :: imax
39 : real(dp), intent(in) :: epsx, epsy
40 : integer, intent(out) :: ierr
41 :
42 : do_safe_root_with_brackets = eval_brent_safe_zero(x1_in, x3_in, 1d-14, &
43 2 : epsx, epsy, f, y1_in, y3_in, lrpar, rpar, lipar, ipar, ierr)
44 :
45 2 : end function do_safe_root_with_brackets
46 :
47 13630 : real(dp) function do_safe_root_with_guess(f, x_guess, dx, x1_in, x3_in, y1_in, y3_in, &
48 : newt_imax, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
49 : interface
50 : ! f provides function values and optional derivatives
51 : #include "num_root_fcn.dek"
52 : end interface
53 : real(dp), intent(in) :: x_guess ! initial guess for the root (required)
54 : real(dp), intent(in) :: dx, x1_in, x3_in, y1_in, y3_in
55 : ! dx is increment for searching for brackets (not used if x1 and x3 are given)
56 : ! x1 and x3 bracket solution (can be arg_not_provided)
57 : ! y1 and y3 are f(x1) and f(x3) (can be arg_not_provided)
58 : integer, intent(in) :: lipar, lrpar
59 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
60 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
61 : integer, intent(in) :: newt_imax, imax
62 : real(dp), intent(in) :: epsx, epsy
63 : ! stop searching when x is determined to within epsx
64 : ! or when abs(f(x)) is less than epsy
65 : integer, intent(out) :: ierr
66 :
67 : integer :: i
68 : logical :: have_x1, have_x3, have_y1, have_y3
69 13630 : real(dp) :: x, y, x1, x3, y1, y3, dydx, absy, absy_prev
70 :
71 : logical, parameter :: dbg = .false.
72 : include 'formats'
73 :
74 13630 : ierr = 0
75 13630 : do_safe_root_with_guess = 0
76 13630 : x1 = x1_in
77 13630 : x3 = x3_in
78 13630 : y1 = y1_in
79 13630 : y3 = y3_in
80 13630 : have_x1 = (x1 /= arg_not_provided)
81 13630 : have_x3 = (x3 /= arg_not_provided)
82 13630 : have_y1 = (y1 /= arg_not_provided)
83 13630 : have_y3 = (y3 /= arg_not_provided)
84 13630 : x = x_guess
85 13630 : absy_prev = 0
86 :
87 45550 : do i = 1, newt_imax
88 45548 : y = f(x, dydx, lrpar, rpar, lipar, ipar, ierr)
89 45548 : if (ierr /= 0) then
90 : if (dbg) then
91 : write (*, 3) 'x y dydx', ierr, i, x, y, dydx
92 : stop
93 : end if
94 0 : ierr = 0; exit ! try safe_root
95 : end if
96 : if (dbg) then
97 : write (*, 2) 'x y dydx', i, x, y, dydx
98 : end if
99 45548 : absy = abs(y)
100 : ! converged if abs(y) < epsy or abs(y/dydx) < epsx
101 : !write(*,2) 'y, epsy, y/dydx, epsx', i, y, epsy, y/dydx, epsx
102 45548 : if (absy < max(epsy, abs(dydx)*epsx)) then
103 : if (dbg) write (*, 1) 'converged', x
104 13627 : do_safe_root_with_guess = x
105 13627 : return
106 : end if
107 31921 : if (dydx == 0) exit ! try safe_root
108 31921 : if (i > 1 .and. absy > absy_prev) exit ! not converging
109 31920 : absy_prev = absy
110 31920 : if (y < 0) then
111 10118 : x1 = x; y1 = y; have_x1 = .true.
112 : else
113 21802 : x3 = x; y3 = y; have_x3 = .true.
114 : end if
115 31923 : x = x - y/dydx
116 : end do
117 :
118 3 : if (.not. (have_x1 .and. have_x3)) then
119 2 : call do_look_for_brackets(x_guess, dx, x1, x3, f, y1, y3, imax, lrpar, rpar, lipar, ipar, ierr)
120 2 : if (ierr /= 0) then
121 : if (dbg) then
122 : write (*, *) 'failed in do_look_for_brackets'
123 : stop 'do_safe_root_with_guess'
124 : end if
125 : return
126 : end if
127 : else
128 1 : if (.not. have_y1) then
129 1 : y1 = f(x1, dydx, lrpar, rpar, lipar, ipar, ierr)
130 1 : if (ierr /= 0) then
131 : if (dbg) then
132 : write (*, *) 'failed evaluating f(x1)'
133 : stop 'do_safe_root_with_guess'
134 : end if
135 : return
136 : end if
137 : end if
138 1 : if (.not. have_y3) then
139 1 : y3 = f(x3, dydx, lrpar, rpar, lipar, ipar, ierr)
140 1 : if (ierr /= 0) then
141 : if (dbg) then
142 : write (*, *) 'failed evaluating f(x3)'
143 : stop 'do_safe_root_with_guess'
144 : end if
145 : return
146 : end if
147 : end if
148 : end if
149 :
150 3 : do_safe_root_with_guess = do_safe_root_with_brackets(f, x1, x3, y1, y3, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
151 : if (ierr /= 0) then
152 : if (dbg) then
153 : write (*, *) 'failed in do_safe_root'
154 : stop 'do_safe_root_with_guess'
155 : end if
156 : return
157 : end if
158 :
159 : if (dbg) then
160 : write (*, 1) 'do_safe_root result:', do_safe_root_with_guess
161 : write (*, *)
162 : end if
163 :
164 : end function do_safe_root_with_guess
165 :
166 : ! safe_root requires bracketing values for the root.
167 : ! if you don't have them, you can use this routine to do a (not too dumb) search.
168 3 : subroutine do_look_for_brackets(x, dx, x1, x3, f, y1, y3, imax, lrpar, rpar, lipar, ipar, ierr)
169 : real(dp), intent(in) :: x, dx ! x is initial guess and dx is increment for searching
170 : real(dp), intent(out) :: x1, x3 ! bounds
171 : real(dp), intent(out) :: y1, y3 ! f(x1) and f(x3)
172 : interface
173 : ! f provides function values
174 : #include "num_root_fcn.dek"
175 : end interface
176 : integer, intent(in) :: imax
177 : integer, intent(in) :: lipar, lrpar
178 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
179 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
180 : integer, intent(out) :: ierr
181 :
182 3 : real(dp) :: jump, dfdx
183 : integer :: i
184 : logical :: move_x1, move_x3
185 :
186 : include 'formats'
187 :
188 3 : ierr = -1
189 3 : y1 = 0
190 3 : y3 = 0
191 :
192 3 : x1 = x; x3 = x
193 : ! after this, keep x1 < x3
194 :
195 : !write(*,2) 'do_look_for_brackets imax x dx', imax, x, dx
196 :
197 8 : do i = 1, imax
198 :
199 8 : jump = (2**i)*dx
200 :
201 8 : if (i == 1) then
202 3 : x1 = x1 - jump
203 3 : x3 = x3 + jump
204 : move_x1 = .true.
205 : move_x3 = .true.
206 5 : else if (y1 > 0) then ! both positive. move x for smaller
207 2 : if (y1 < y3) then
208 2 : x1 = x1 - jump
209 : move_x1 = .true.
210 : move_x3 = .false.
211 : else
212 0 : x3 = x3 + jump
213 : move_x3 = .true.
214 : move_x1 = .false.
215 : end if
216 : else ! both negative. move x for larger
217 3 : if (y1 > y3) then
218 0 : x1 = x1 - jump
219 : move_x1 = .true.
220 : move_x3 = .false.
221 : else
222 3 : x3 = x3 + jump
223 : move_x3 = .true.
224 : move_x1 = .false.
225 : end if
226 : end if
227 :
228 : if (move_x1) then
229 5 : y1 = f(x1, dfdx, lrpar, rpar, lipar, ipar, ierr)
230 5 : if (ierr /= 0) then
231 : !write(*,*) 'ierr from f(x1)'
232 3 : return
233 : end if
234 : end if
235 :
236 8 : if (move_x3) then
237 6 : y3 = f(x3, dfdx, lrpar, rpar, lipar, ipar, ierr)
238 6 : if (ierr /= 0) then
239 : !write(*,*) 'ierr from f(x3)'
240 : return
241 : end if
242 : end if
243 :
244 : !write(*,'(a,i4,4(3x,e18.8))') 'look_for_brackets', i, x1, y1, x3, y3
245 :
246 8 : if (y1*y3 <= 0) then
247 3 : ierr = 0
248 : !write(*,1) 'done do_look_for_brackets', y1, y3
249 3 : return
250 : end if
251 :
252 : end do
253 :
254 : !write(*,1) 'exit do_look_for_brackets'
255 :
256 : end subroutine do_look_for_brackets
257 :
258 0 : real(dp) function do_safe_root(f, x_guess, x1_in, x3_in, y1_in, y3_in, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
259 : use const_def, only: arg_not_provided
260 : use mod_brent, only: eval_brent_safe_zero
261 : interface
262 : include 'num_root_fcn.dek' ! f provides function values
263 : end interface
264 : real(dp), intent(in) :: x_guess, x1_in, x3_in
265 : real(dp), intent(in) :: y1_in, y3_in ! f(x1) and f(x3)
266 : integer, intent(in) :: lipar, lrpar
267 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
268 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
269 : integer, intent(in) :: imax
270 : real(dp), intent(in) :: epsx, epsy
271 : integer, intent(out) :: ierr
272 :
273 0 : do_safe_root = eval_brent_safe_zero(x1_in, x3_in, 1d-14, epsx, epsy, f, y1_in, y3_in, lrpar, rpar, lipar, ipar, ierr)
274 :
275 0 : end function do_safe_root
276 :
277 : end module mod_root
|
