Line data Source code
1 : module test_brent
2 :
3 : use num_def
4 : use num_lib
5 : use math_lib
6 : use utils_lib, only: mesa_error
7 : use const_def, only: dp
8 :
9 : implicit none
10 :
11 : logical, parameter :: dbg = .false.
12 :
13 : contains
14 :
15 1 : subroutine do_test_brent
16 1 : write (*, *) 'test brent'
17 :
18 1 : call test_global_min_all
19 1 : call test_local_min_all
20 1 : call test_brent_zero
21 :
22 1 : end subroutine do_test_brent
23 :
24 1 : subroutine test_global_min_all
25 :
26 : !*****************************************************************************80
27 : !
28 : !! TEST_GLOMIN_ALL tests the Brent global minimizer on all test functions.
29 : !
30 : ! Licensing:
31 : !
32 : ! This code is distributed under the GNU LGPL license.
33 : !
34 : ! Modified:
35 : !
36 : ! 12 April 2008
37 : !
38 : ! Author:
39 : !
40 : ! John Burkardt
41 : !
42 : implicit none
43 :
44 : real(dp) :: a
45 : real(dp) :: b
46 : real(dp) :: c
47 : real(dp) :: e
48 : real(dp) :: m
49 : real(dp) :: machep
50 : real(dp) :: t
51 :
52 1 : write (*, '(a)') ' '
53 1 : write (*, '(a)') 'TEST_GLOMIN_ALL'
54 1 : write (*, '(a)') ' Test the Brent GLOMIN routine, which seeks'
55 1 : write (*, '(a)') ' a global minimizer of a function F(X)'
56 1 : write (*, '(a)') ' in an interval [A,B],'
57 1 : write (*, '(a)') ' given some upper bound M for F".'
58 :
59 1 : machep = epsilon(machep)
60 1 : e = sqrt(machep)
61 1 : t = sqrt(machep)
62 :
63 1 : a = 7.0D+00
64 1 : b = 9.0D+00
65 1 : c = (a + b)/2.0D+00
66 1 : m = 0.0D+00
67 :
68 : call test_glomin_one(a, b, c, m, machep, e, t, h_01, &
69 1 : 'h_01(x) = 2 - x')
70 :
71 : a = 7.0D+00
72 : b = 9.0D+00
73 : c = (a + b)/2.0D+00
74 1 : m = 100.0D+00
75 :
76 : call test_glomin_one(a, b, c, m, machep, e, t, h_01, &
77 1 : 'h_01(x) = 2 - x')
78 :
79 1 : a = -1.0D+00
80 1 : b = +2.0D+00
81 1 : c = (a + b)/2.0D+00
82 1 : m = 2.0D+00
83 :
84 : call test_glomin_one(a, b, c, m, machep, e, t, h_02, &
85 1 : 'h_02(x) = x * x')
86 :
87 : a = -1.0D+00
88 : b = +2.0D+00
89 : c = (a + b)/2.0D+00
90 1 : m = 2.1D+00
91 :
92 : call test_glomin_one(a, b, c, m, machep, e, t, h_02, &
93 1 : 'h_02(x) = x * x')
94 :
95 : a = -0.5D+00
96 : b = +2.0D+00
97 : c = (a + b)/2.0D+00
98 : m = 14.0D+00
99 :
100 : !call test_glomin_one ( a, b, c, m, machep, e, t, h_03, &
101 : ! 'h_03(x) = x^3 + x^2' )
102 :
103 : a = -0.5D+00
104 : b = +2.0D+00
105 : c = (a + b)/2.0D+00
106 : m = 28.0D+00
107 :
108 : !call test_glomin_one ( a, b, c, m, machep, e, t, h_03, &
109 : ! 'h_03(x) = x^3 + x^2' )
110 :
111 1 : a = -10.0D+00
112 1 : b = +10.0D+00
113 1 : c = (a + b)/2.0D+00
114 1 : m = 72.0D+00
115 :
116 : call test_glomin_one(a, b, c, m, machep, e, t, h_04, &
117 1 : 'h_04(x) = ( x + sin(x) ) * exp(-x*x)')
118 :
119 : a = -10.0D+00
120 : b = +10.0D+00
121 : c = (a + b)/2.0D+00
122 : m = 72.0D+00
123 :
124 : call test_glomin_one(a, b, c, m, machep, e, t, h_05, &
125 1 : 'h_05(x) = ( x - sin(x) ) * exp(-x*x)')
126 :
127 1 : end subroutine test_global_min_all
128 :
129 6 : subroutine test_glomin_one(a, b, c, m, machep, e, t, f, title)
130 : real(dp), intent(in) :: a, b, c, m, machep, e, t
131 : interface
132 : real(dp) function f(x)
133 : use const_def, only: dp
134 : implicit none
135 : real(dp), intent(in) :: x
136 : end function f
137 : end interface
138 :
139 6 : real(dp) :: fa
140 6 : real(dp) :: fb
141 6 : real(dp) :: fx
142 : character(len=*) :: title
143 : real(dp) :: x
144 : integer :: max_tries, ierr
145 : include 'formats'
146 :
147 6 : max_tries = 10000
148 : ierr = 0
149 6 : fx = brent_global_min(max_tries, a, b, c, m, machep, e, t, f, x, ierr)
150 6 : if (ierr /= 0) then
151 0 : write (*, *) 'failed in brent_global_min'
152 0 : call mesa_error(__FILE__, __LINE__)
153 : end if
154 6 : fa = f(a)
155 6 : fb = f(b)
156 :
157 6 : write (*, '(a)') ' '
158 6 : write (*, '(a)') trim(title)
159 6 : write (*, 1) 'a, b', a, b
160 6 : write (*, 1) 'fa, fb', fa, fb
161 :
162 6 : end subroutine test_glomin_one
163 :
164 31 : real(dp) function h_01(x)
165 : real(dp), intent(in) :: x
166 31 : h_01 = 2.0D+00 - x
167 31 : end function h_01
168 :
169 45 : real(dp) function h_02(x)
170 : real(dp), intent(in) :: x
171 45 : h_02 = x*x
172 45 : end function h_02
173 :
174 0 : real(dp) function h_03(x)
175 : real(dp), intent(in) :: x
176 0 : h_03 = x*x*(x + 1.0D+00)
177 0 : end function h_03
178 :
179 383 : real(dp) function h_04(x)
180 : real(dp), intent(in) :: x
181 383 : h_04 = (x + sin(x))*exp(-x*x)
182 383 : end function h_04
183 :
184 1135 : real(dp) function h_05(x)
185 : real(dp), intent(in) :: x
186 1135 : h_05 = (x - sin(x))*exp(-x*x)
187 1135 : end function h_05
188 :
189 1 : subroutine test_local_min_all
190 :
191 : !*****************************************************************************80
192 : !
193 : !! TEST_LOCAL_MIN_ALL tests Brent's local minimizer on all test functions.
194 : !
195 : ! Licensing:
196 : !
197 : ! This code is distributed under the GNU LGPL license.
198 : !
199 : ! Modified:
200 : !
201 : ! 12 April 2008
202 : !
203 : ! Author:
204 : !
205 : ! John Burkardt
206 : !
207 : implicit none
208 :
209 : real(dp) :: a
210 : real(dp) :: b
211 : real(dp) :: eps
212 : real(dp) :: t
213 :
214 1 : write (*, '(a)') ' '
215 1 : write (*, '(a)') 'TEST_LOCAL_MIN_ALL'
216 1 : write (*, '(a)') ' Test the Brent LOCAL_MIN routine, which seeks'
217 1 : write (*, '(a)') ' a local minimizer of a function F(X)'
218 1 : write (*, '(a)') ' in an interval [A,B].'
219 :
220 1 : eps = 10.0D+00*sqrt(epsilon(eps))
221 1 : t = eps
222 :
223 1 : a = 0.0D+00
224 1 : b = 3.141592653589793D+00
225 :
226 : call test_local_min_one(a, b, eps, t, g_01, &
227 1 : 'g_01(x) = ( x - 2 ) * ( x - 2 ) + 1')
228 :
229 : a = 0.0D+00
230 1 : b = 1.0D+00
231 :
232 : call test_local_min_one(a, b, eps, t, g_02, &
233 1 : 'g_02(x) = x * x + exp( - x )')
234 :
235 1 : a = -2.0D+00
236 1 : b = 2.0D+00
237 :
238 : call test_local_min_one(a, b, eps, t, g_03, &
239 1 : 'g_03(x) = x^4 + 2x^2 + x + 3')
240 :
241 1 : a = 0.0001D+00
242 1 : b = 1.0D+00
243 :
244 : call test_local_min_one(a, b, eps, t, g_04, &
245 1 : 'g_04(x) = exp( x ) + 1 / ( 100 x )')
246 :
247 1 : a = 0.0002D+00
248 1 : b = 2.0D+00
249 :
250 : call test_local_min_one(a, b, eps, t, g_05, &
251 1 : 'g_05(x) = exp( x ) - 2x + 1/(100x) - 1/(1000000x^2)')
252 :
253 1 : end subroutine test_local_min_all
254 :
255 5 : subroutine test_local_min_one(a, b, eps, t, f, title)
256 :
257 : !*****************************************************************************80
258 : !
259 : !! TEST_LOCAL_MIN_ONE tests Brent's local minimizer on one test function.
260 : !
261 : ! Licensing:
262 : !
263 : ! This code is distributed under the GNU LGPL license.
264 : !
265 : ! Modified:
266 : !
267 : ! 12 April 2008
268 : !
269 : ! Author:
270 : !
271 : ! John Burkardt
272 : !
273 : ! Parameters:
274 : !
275 : ! Input, real(dp) A, B, the endpoints of the interval.
276 : !
277 : ! Input, real(dp) EPS, a positive relative error tolerance.
278 : !
279 : ! Input, real(dp) T, a positive absolute error tolerance.
280 : !
281 : ! Input, external real(dp) F, the name of a user-supplied
282 : ! function, of the form "FUNCTION F ( X )", which evaluates the
283 : ! function whose local minimum is being sought.
284 : !
285 : ! Input, character ( LEN = * ) TITLE, a title for the problem.
286 : !
287 : implicit none
288 : real(dp), intent(in) :: a, b, eps, t
289 : interface
290 : real(dp) function f(x)
291 : use const_def, only: dp
292 : implicit none
293 : real(dp), intent(in) :: x
294 : end function f
295 : end interface
296 : character(len=*) :: title
297 :
298 5 : real(dp) :: fa
299 5 : real(dp) :: fb
300 5 : real(dp) :: fx
301 : real(dp) :: x
302 : integer :: max_tries, ierr
303 : include 'formats'
304 :
305 5 : max_tries = 10000
306 : ierr = 0
307 5 : fx = brent_local_min(max_tries, a, b, eps, t, f, x, ierr)
308 5 : if (ierr /= 0) then
309 0 : write (*, *) 'failed in brent_local_min'
310 0 : call mesa_error(__FILE__, __LINE__)
311 : end if
312 5 : fa = f(a)
313 5 : fb = f(b)
314 :
315 5 : write (*, '(a)') ' '
316 5 : write (*, '(a)') trim(title)
317 5 : write (*, 1) 'a, b', a, b
318 5 : write (*, 1) 'fa, fb', fa, fb
319 :
320 5 : end subroutine test_local_min_one
321 :
322 8 : real(dp) function g_01(x)
323 : real(dp), intent(in) :: x
324 8 : g_01 = (x - 2.0D+00)*(x - 2.0D+00) + 1.0D+00
325 8 : end function g_01
326 :
327 11 : real(dp) function g_02(x)
328 : real(dp), intent(in) :: x
329 11 : g_02 = x*x + exp(-x)
330 11 : end function g_02
331 :
332 13 : real(dp) function g_03(x)
333 : real(dp), intent(in) :: x
334 13 : g_03 = ((x*x + 2.0D+00)*x + 1.0D+00)*x + 3.0D+00
335 13 : end function g_03
336 :
337 15 : real(dp) function g_04(x)
338 : real(dp), intent(in) :: x
339 15 : g_04 = exp(x) + 0.01D+00/x
340 15 : end function g_04
341 :
342 12 : real(dp) function g_05(x)
343 : real(dp), intent(in) :: x
344 12 : g_05 = exp(x) - 2.0D+00*x + 0.01D+00/x - 0.000001D+00/x/x
345 12 : end function g_05
346 :
347 29 : real(dp) function f_01(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
348 : integer, intent(in) :: lrpar, lipar
349 : real(dp), intent(in) :: x
350 : real(dp), intent(out) :: dfdx
351 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
352 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
353 : integer, intent(out) :: ierr
354 29 : f_01 = sin(x) - 0.5D+00*x
355 29 : ierr = 0
356 29 : dfdx = 0
357 29 : end function f_01
358 :
359 10 : real(dp) function f_02(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
360 : integer, intent(in) :: lrpar, lipar
361 : real(dp), intent(in) :: x
362 : real(dp), intent(out) :: dfdx
363 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
364 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
365 : integer, intent(out) :: ierr
366 10 : f_02 = 2.0D+00*x - exp(-x)
367 10 : ierr = 0
368 10 : dfdx = 0
369 10 : end function f_02
370 :
371 15 : real(dp) function f_03(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
372 : integer, intent(in) :: lrpar, lipar
373 : real(dp), intent(in) :: x
374 : real(dp), intent(out) :: dfdx
375 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
376 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
377 : integer, intent(out) :: ierr
378 15 : f_03 = x*exp(-x)
379 15 : ierr = 0
380 15 : dfdx = 0
381 15 : end function f_03
382 :
383 18 : real(dp) function f_04(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
384 : integer, intent(in) :: lrpar, lipar
385 : real(dp), intent(in) :: x
386 : real(dp), intent(out) :: dfdx
387 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
388 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
389 : integer, intent(out) :: ierr
390 18 : f_04 = exp(x) - 1.0D+00/100.0D+00/x/x
391 18 : ierr = 0
392 18 : dfdx = 0
393 18 : end function f_04
394 :
395 18 : real(dp) function f_05(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
396 : integer, intent(in) :: lrpar, lipar
397 : real(dp), intent(in) :: x
398 : real(dp), intent(out) :: dfdx
399 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
400 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
401 : integer, intent(out) :: ierr
402 18 : f_05 = (x + 3.0D+00)*(x - 1.0D+00)*(x - 1.0D+00)
403 18 : ierr = 0
404 18 : dfdx = 0
405 18 : end function f_05
406 :
407 1 : subroutine test_brent_zero()
408 :
409 : !*****************************************************************************80
410 : !
411 : !! test_brent_zero tests Brent's zero finding routine on all test functions.
412 : !
413 : ! Licensing:
414 : !
415 : ! This code is distributed under the GNU LGPL license.
416 : !
417 : ! Modified:
418 : !
419 : ! 12 April 2008
420 : !
421 : ! Author:
422 : !
423 : ! John Burkardt
424 : !
425 : implicit none
426 :
427 : real(dp) :: a
428 : real(dp) :: b
429 : real(dp) :: machep
430 : real(dp) :: t
431 :
432 1 : machep = epsilon(machep)
433 1 : t = machep
434 :
435 1 : a = 1.0D+00
436 1 : b = 2.0D+00
437 :
438 : call test_zero_one(a, b, machep, t, f_01, &
439 1 : 'f_01(x) = sin ( x ) - x / 2')
440 :
441 1 : a = 0.0D+00
442 1 : b = 1.0D+00
443 :
444 : call test_zero_one(a, b, machep, t, f_02, &
445 1 : 'f_02(x) = 2 * x - exp( - x )')
446 :
447 1 : a = -1.0D+00
448 1 : b = 0.5D+00
449 :
450 : call test_zero_one(a, b, machep, t, f_03, &
451 1 : 'f_03(x) = x * exp( - x )')
452 :
453 1 : a = 0.0001D+00
454 1 : b = 20.0D+00
455 :
456 : call test_zero_one(a, b, machep, t, f_04, &
457 1 : 'f_04(x) = exp( x ) - 1 / ( 100 * x * x )')
458 :
459 1 : a = -5.0D+00
460 1 : b = 2.0D+00
461 :
462 : call test_zero_one(a, b, machep, t, f_05, &
463 1 : 'f_05(x) = (x+3) * (x-1) * (x-1)')
464 :
465 1 : end subroutine test_brent_zero
466 :
467 5 : subroutine test_zero_one(a, b, machep, t, f, title)
468 :
469 : !*****************************************************************************80
470 : !
471 : !! TEST_ZERO_ONE tests Brent's zero finding routine on one test function.
472 : !
473 : ! Licensing:
474 : !
475 : ! This code is distributed under the GNU LGPL license.
476 : !
477 : ! Modified:
478 : !
479 : ! 12 April 2008
480 : !
481 : ! Author:
482 : !
483 : ! John Burkardt
484 : !
485 : ! Parameters:
486 : !
487 : ! Input, real(dp) A, B, the two endpoints of the change of sign
488 : ! interval.
489 : !
490 : ! Input, real(dp) MACHEP, an estimate for the relative machine
491 : ! precision.
492 : !
493 : ! Input, real(dp) T, a positive error tolerance.
494 : !
495 : ! Input, external real(dp) F, the name of a user-supplied
496 : ! function, of the form "FUNCTION F ( X )", which evaluates the
497 : ! function whose zero is being sought.
498 : !
499 : ! Input, character ( LEN = * ) TITLE, a title for the problem.
500 : !
501 : implicit none
502 :
503 : interface
504 : include 'num_root_fcn.dek' ! f provides function values
505 : end interface
506 :
507 : real(dp) :: a
508 : real(dp) :: b
509 : real(dp) :: fa
510 : real(dp) :: fb
511 : real(dp) :: fz
512 : real(dp) :: machep
513 : real(dp) :: t
514 : character(len=*) :: title
515 : real(dp) :: z
516 : real(dp) :: dfdx
517 :
518 : integer, parameter :: lrpar = 0, lipar = 0
519 : integer :: ierr
520 : real(dp), target :: rpar_ary(lrpar)
521 : integer, target :: ipar_ary(lipar)
522 : real(dp), pointer :: rpar(:)
523 : integer, pointer :: ipar(:)
524 :
525 : include 'formats'
526 :
527 5 : rpar => rpar_ary
528 5 : ipar => ipar_ary
529 :
530 : ierr = 0
531 5 : fa = f(a, dfdx, lrpar, rpar, lipar, ipar, ierr)
532 5 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
533 5 : fb = f(b, dfdx, lrpar, rpar, lipar, ipar, ierr)
534 5 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
535 5 : z = brent_safe_zero(a, b, machep, t, 0d0, f, fa, fb, lrpar, rpar, lipar, ipar, ierr)
536 5 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
537 5 : fz = f(z, dfdx, lrpar, rpar, lipar, ipar, ierr)
538 5 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
539 :
540 5 : if (abs(fz) < 1d-14) return
541 :
542 0 : write (*, '(a)') ' '
543 0 : write (*, '(a)') trim(title)
544 0 : write (*, 1) 'a, z, b', a, z, b
545 0 : write (*, 1) 'fa, fz, fb', fa, fz, fb
546 0 : write (*, '(a)') ' '
547 :
548 : end subroutine test_zero_one
549 :
550 : end module test_brent
|
