Line data Source code
1 : module mod_brent
2 : use const_def, only: dp
3 :
4 : implicit none
5 :
6 : contains
7 :
8 : ! this code was written by John Burkardt
9 : ! see his wonderful site for lots of numerical software.
10 : ! here's the link for the implementation that we use here.
11 : ! http://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.html
12 :
13 :
14 6 : real(dp) function eval_global_min(max_tries, a, b, c, m, machep, e, t, f, x, ierr)
15 : integer, intent(in) :: max_tries
16 : real(dp), intent(in) :: a, b, c, m, machep, e, t
17 : interface
18 : real(dp) function f(x)
19 : use const_def, only: dp
20 : implicit none
21 : real(dp), intent(in) :: x
22 : end function f
23 : end interface
24 : real(dp), intent(out) :: x
25 : integer, intent(out) :: ierr
26 :
27 : ! Parameters:
28 : !
29 : ! Input, real ( kind = 8 ) A, B, the endpoints of the interval.
30 : ! It must be the case that A < B.
31 : !
32 : ! Input, real ( kind = 8 ) C, an initial guess for the global
33 : ! minimizer. If no good guess is known, C = A or B is acceptable.
34 : !
35 : ! Input, real ( kind = 8 ) M, the bound on the second derivative.
36 : !
37 : ! Input, real ( kind = 8 ) MACHEP, an estimate for the relative machine
38 : ! precision.
39 : !
40 : ! Input, real ( kind = 8 ) E, a positive tolerance, a bound for the
41 : ! absolute error in the evaluation of F(X) for any X in [A,B].
42 : !
43 : ! Input, real ( kind = 8 ) T, a positive error tolerance.
44 : !
45 : ! Input, external real ( kind = 8 ) F, the name of a user-supplied
46 : ! function, of the form "FUNCTION F ( X )", which evaluates the
47 : ! function whose global minimum is being sought.
48 : !
49 : ! Output, real ( kind = 8 ) X, the estimated value of the abscissa
50 : ! for which F attains its global minimum value in [A,B].
51 : !
52 : ! Output, real ( kind = 8 ) GLOMIN, the value F(X).
53 : !
54 :
55 6 : real ( kind = 8 ) :: a0
56 6 : real ( kind = 8 ) :: a2
57 6 : real ( kind = 8 ) :: a3
58 6 : real ( kind = 8 ) :: d0
59 6 : real ( kind = 8 ) :: d1
60 6 : real ( kind = 8 ) :: d2
61 6 : real ( kind = 8 ) :: h
62 : integer ( kind = 4 ) :: k, iter
63 6 : real ( kind = 8 ) :: m2
64 6 : real ( kind = 8 ) :: p
65 6 : real ( kind = 8 ) :: q
66 6 : real ( kind = 8 ) :: qs
67 6 : real ( kind = 8 ) :: r
68 6 : real ( kind = 8 ) :: s
69 6 : real ( kind = 8 ) :: sc
70 6 : real ( kind = 8 ) :: y
71 6 : real ( kind = 8 ) :: y0
72 6 : real ( kind = 8 ) :: y1
73 6 : real ( kind = 8 ) :: y2
74 6 : real ( kind = 8 ) :: y3
75 6 : real ( kind = 8 ) :: yb
76 6 : real ( kind = 8 ) :: z0
77 6 : real ( kind = 8 ) :: z1
78 6 : real ( kind = 8 ) :: z2
79 :
80 6 : ierr = 0
81 6 : a0 = b
82 6 : x = a0
83 6 : a2 = a
84 6 : y0 = f ( b )
85 6 : yb = y0
86 6 : y2 = f ( a )
87 6 : y = y2
88 :
89 6 : if ( y0 < y ) then
90 : y = y0
91 : else
92 4 : x = a
93 : end if
94 :
95 6 : if ( m <= 0.0D+00 .or. b <= a ) then
96 : eval_global_min = y
97 : return
98 : end if
99 :
100 5 : m2 = 0.5D+00 * ( 1.0D+00 + 16.0D+00 * machep ) * m
101 :
102 5 : if ( c <= a .or. b <= c ) then
103 0 : sc = 0.5D+00 * ( a + b )
104 : else
105 5 : sc = c
106 : end if
107 :
108 5 : y1 = f ( sc )
109 5 : k = 3
110 5 : d0 = a2 - sc
111 5 : h = 9.0D+00 / 11.0D+00
112 :
113 5 : if ( y1 < y ) then
114 2 : x = sc
115 2 : y = y1
116 : end if
117 :
118 1084 : do iter = 1, max_tries+1
119 :
120 1084 : if (iter > max_tries) then
121 0 : ierr = -1
122 0 : exit
123 : end if
124 :
125 1084 : d1 = a2 - a0
126 1084 : d2 = sc - a0
127 1084 : z2 = b - a2
128 1084 : z0 = y2 - y1
129 1084 : z1 = y2 - y0
130 1084 : r = d1 * d1 * z0 - d0 * d0 * z1
131 1084 : p = r
132 1084 : qs = 2.0D+00 * ( d0 * z1 - d1 * z0 )
133 1084 : q = qs
134 :
135 1084 : if ( k < 1000000 .or. y2 <= y ) then
136 :
137 : do
138 :
139 1114 : if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) < &
140 : z2 * m2 * r * ( z2 * q - r ) ) then
141 470 : a3 = a2 + r / q
142 470 : y3 = f ( a3 )
143 :
144 470 : if ( y3 < y ) then
145 140 : x = a3
146 140 : y = y3
147 : end if
148 : end if
149 :
150 1114 : k = mod ( 1611 * k, 1048576 )
151 1114 : q = 1.0D+00
152 1114 : r = ( b - a ) * 0.00001D+00 * real ( k, kind = 8 )
153 :
154 1114 : if ( z2 <= r ) then
155 : exit
156 : end if
157 :
158 : end do
159 :
160 : else
161 :
162 43 : k = mod ( 1611 * k, 1048576 )
163 43 : q = 1.0D+00
164 43 : r = ( b - a ) * 0.00001D+00 * real ( k, kind = 8 )
165 :
166 47 : do while ( r < z2 )
167 :
168 4 : if ( q * ( r * ( yb - y2 ) + z2 * q * ( ( y2 - y ) + t ) ) < &
169 : z2 * m2 * r * ( z2 * q - r ) ) then
170 4 : a3 = a2 + r / q
171 4 : y3 = f ( a3 )
172 :
173 4 : if ( y3 < y ) then
174 0 : x = a3
175 0 : y = y3
176 : end if
177 : end if
178 :
179 4 : k = mod ( 1611 * k, 1048576 )
180 4 : q = 1.0D+00
181 4 : r = ( b - a ) * 0.00001D+00 * real ( k, kind = 8 )
182 :
183 : end do
184 :
185 : end if
186 :
187 1084 : r = m2 * d0 * d1 * d2
188 1084 : s = sqrt ( ( ( y2 - y ) + t ) / m2 )
189 1084 : h = 0.5D+00 * ( 1.0D+00 + h )
190 1084 : p = h * ( p + 2.0D+00 * r * s )
191 1084 : q = q + 0.5D+00 * qs
192 1084 : r = - 0.5D+00 * ( d0 + ( z0 + 2.01D+00 * e ) / ( d0 * m2 ) )
193 :
194 1084 : if ( r < s .or. d0 < 0.0D+00 ) then
195 1059 : r = a2 + s
196 : else
197 25 : r = a2 + r
198 : end if
199 :
200 1084 : if ( 0.0D+00 < p * q ) then
201 1053 : a3 = a2 + p / q
202 : else
203 31 : a3 = r
204 : end if
205 :
206 13 : do
207 :
208 1097 : a3 = max ( a3, r )
209 :
210 1097 : if ( b <= a3 ) then
211 6 : a3 = b
212 6 : y3 = yb
213 : else
214 1091 : y3 = f ( a3 )
215 : end if
216 :
217 1097 : if ( y3 < y ) then
218 27 : x = a3
219 27 : y = y3
220 : end if
221 :
222 1097 : d0 = a3 - a2
223 :
224 1097 : if ( a3 <= r ) then
225 : exit
226 : end if
227 :
228 41 : p = 2.0D+00 * ( y2 - y3 ) / ( m * d0 )
229 :
230 41 : if ( ( 1.0D+00 + 9.0D+00 * machep ) * d0 <= abs ( p ) ) then
231 : exit
232 : end if
233 :
234 40 : if ( 0.5D+00 * m2 * ( d0 * d0 + p * p ) <= &
235 : ( y2 - y ) + ( y3 - y ) + 2.0D+00 * t ) then
236 : exit
237 : end if
238 :
239 13 : a3 = 0.5D+00 * ( a2 + a3 )
240 41 : h = 0.9D+00 * h
241 :
242 : end do
243 :
244 1084 : if ( b <= a3 ) then
245 : exit
246 : end if
247 :
248 1079 : a0 = sc
249 1079 : sc = a2
250 1079 : a2 = a3
251 1079 : y0 = y1
252 1079 : y1 = y2
253 1084 : y2 = y3
254 :
255 : end do
256 :
257 : eval_global_min = y
258 :
259 : return
260 : end function eval_global_min
261 :
262 :
263 5 : real(dp) function eval_local_min(max_tries, a, b, eps, t, f, x, ierr)
264 : integer, intent(in) :: max_tries
265 : real(dp), intent(in) :: a, b, eps, t
266 : interface
267 : real(dp) function f(x)
268 : use const_def, only: dp
269 : implicit none
270 : real(dp), intent(in) :: x
271 : end function f
272 : end interface
273 : real(dp), intent(out) :: x
274 : integer, intent(out) :: ierr
275 :
276 : ! *****************************************************************************80
277 : !
278 : !! LOCAL_MIN seeks a local minimum of a function F(X) in an interval [A,B].
279 : !
280 : ! Discussion:
281 : !
282 : ! The method used is a combination of golden section search and
283 : ! successive parabolic interpolation. Convergence is never much slower
284 : ! than that for a Fibonacci search. If F has a continuous second
285 : ! derivative which is positive at the minimum (which is not at A or
286 : ! B), then convergence is superlinear, and usually of the order of
287 : ! about 1.324....
288 : !
289 : ! The values EPS and T define a tolerance TOL = EPS * abs ( X ) + T.
290 : ! F is never evaluated at two points closer than TOL.
291 : !
292 : ! If F is a unimodal function and the computed values of F are always
293 : ! unimodal when separated by at least SQEPS * abs ( X ) + (T/3), then
294 : ! LOCAL_MIN approximates the abscissa of the global minimum of F on the
295 : ! interval [A,B] with an error less than 3*SQEPS*abs(LOCAL_MIN)+T.
296 : !
297 : ! If F is not unimodal, then LOCAL_MIN may approximate a local, but
298 : ! perhaps non-global, minimum to the same accuracy.
299 : !
300 : ! Licensing:
301 : !
302 : ! This code is distributed under the GNU LGPL license.
303 : !
304 : ! Modified:
305 : !
306 : ! 13 April 2008
307 : !
308 : ! Author:
309 : !
310 : ! Original FORTRAN77 version by Richard Brent
311 : ! FORTRAN90 version by John Burkardt
312 : !
313 : ! Reference:
314 : !
315 : ! Richard Brent,
316 : ! Algorithms for Minimization Without Derivatives,
317 : ! Dover, 2002,
318 : ! ISBN: 0-486-41998-3,
319 : ! LC: QA402.5.B74.
320 : !
321 : ! Parameters:
322 : !
323 : ! Input, real ( kind = 8 ) A, B, the endpoints of the interval.
324 : !
325 : ! Input, real ( kind = 8 ) EPS, a positive relative error tolerance.
326 : ! EPS should be no smaller than twice the relative machine precision,
327 : ! and preferably not much less than the square root of the relative
328 : ! machine precision.
329 : !
330 : ! Input, real ( kind = 8 ) T, a positive absolute error tolerance.
331 : !
332 : ! Input, external real ( kind = 8 ) F, the name of a user-supplied
333 : ! function, of the form "FUNCTION F ( X )", which evaluates the
334 : ! function whose local minimum is being sought.
335 : !
336 : ! Output, real ( kind = 8 ) X, the estimated value of an abscissa
337 : ! for which F attains a local minimum value in [A,B].
338 : !
339 : ! Output, real ( kind = 8 ) LOCAL_MIN, the value F(X).
340 : !
341 :
342 5 : real ( kind = 8 ) :: c
343 5 : real ( kind = 8 ) :: d
344 5 : real ( kind = 8 ) :: e
345 5 : real ( kind = 8 ) :: fu
346 5 : real ( kind = 8 ) :: fv
347 5 : real ( kind = 8 ) :: fw
348 5 : real ( kind = 8 ) :: fx
349 5 : real ( kind = 8 ) :: m
350 5 : real ( kind = 8 ) :: p
351 5 : real ( kind = 8 ) :: q
352 5 : real ( kind = 8 ) :: r
353 5 : real ( kind = 8 ) :: sa
354 5 : real ( kind = 8 ) :: sb
355 5 : real ( kind = 8 ) :: t2
356 5 : real ( kind = 8 ) :: tol
357 5 : real ( kind = 8 ) :: u
358 5 : real ( kind = 8 ) :: v
359 5 : real ( kind = 8 ) :: w
360 : integer ( kind = 4 ) :: iter
361 :
362 5 : ierr = 0
363 :
364 : !
365 : ! C is the square of the inverse of the golden ratio.
366 : !
367 5 : c = 0.5D+00 * ( 3.0D+00 - sqrt ( 5.0D+00 ) )
368 :
369 5 : sa = a
370 5 : sb = b
371 5 : x = sa + c * ( sb - sa )
372 5 : fx = f ( x )
373 5 : w = x
374 5 : v = w
375 5 : e = 0.0D+00
376 5 : fw = fx
377 5 : fv = fw
378 5 : d = 0
379 :
380 49 : do iter = 1, max_tries+1
381 :
382 49 : m = 0.5D+00 * ( sa + sb )
383 49 : tol = eps * abs ( x ) + t
384 49 : t2 = 2.0D+00 * tol
385 : !
386 : ! Check the stopping criterion.
387 : !
388 49 : if ( abs ( x - m ) <= t2 - 0.5D+00 * ( sb - sa ) ) then
389 : exit
390 : end if
391 :
392 44 : if (iter > max_tries) then
393 0 : ierr = -1
394 0 : exit
395 : end if
396 : !
397 : ! Fit a parabola.
398 : !
399 44 : r = 0.0D+00
400 44 : q = r
401 44 : p = q
402 :
403 44 : if ( tol < abs ( e ) ) then
404 :
405 39 : r = ( x - w ) * ( fx - fv )
406 39 : q = ( x - v ) * ( fx - fw )
407 39 : p = ( x - v ) * q - ( x - w ) * r
408 39 : q = 2.0D+00 * ( q - r )
409 :
410 39 : if ( 0.0D+00 < q ) then
411 14 : p = - p
412 : end if
413 :
414 39 : q = abs ( q )
415 :
416 39 : r = e
417 39 : e = d
418 :
419 : end if
420 :
421 : if ( abs ( p ) < abs ( 0.5D+00 * q * r ) .and. &
422 44 : q * ( sa - x ) < p .and. &
423 : p < q * ( sb - x ) ) then
424 : !
425 : ! Take the parabolic interpolation step.
426 : !
427 32 : d = p / q
428 32 : u = x + d
429 : !
430 : ! F must not be evaluated too close to A or B.
431 : !
432 32 : if ( ( u - sa ) < t2 .or. ( sb - u ) < t2 ) then
433 :
434 5 : if ( x < m ) then
435 : d = tol
436 : else
437 3 : d = - tol
438 : end if
439 :
440 : end if
441 : !
442 : ! A golden-section step.
443 : !
444 : else
445 :
446 12 : if ( x < m ) then
447 6 : e = sb - x
448 : else
449 6 : e = sa - x
450 : end if
451 :
452 12 : d = c * e
453 :
454 : end if
455 : !
456 : ! F must not be evaluated too close to X.
457 : !
458 44 : if ( tol <= abs ( d ) ) then
459 41 : u = x + d
460 3 : else if ( 0.0D+00 < d ) then
461 2 : u = x + tol
462 : else
463 1 : u = x - tol
464 : end if
465 :
466 44 : fu = f ( u )
467 : !
468 : ! Update A, B, V, W, and X.
469 : !
470 49 : if ( fu <= fx ) then
471 :
472 26 : if ( u < x ) then
473 : sb = x
474 : else
475 11 : sa = x
476 : end if
477 :
478 26 : v = w
479 26 : fv = fw
480 26 : w = x
481 26 : fw = fx
482 26 : x = u
483 26 : fx = fu
484 :
485 : else
486 :
487 18 : if ( u < x ) then
488 : sa = u
489 : else
490 10 : sb = u
491 : end if
492 :
493 18 : if ( fu <= fw .or. w == x ) then
494 : v = w
495 : fv = fw
496 : w = u
497 : fw = fu
498 3 : else if ( fu <= fv .or. v == x .or. v == w ) then
499 3 : v = u
500 3 : fv = fu
501 : end if
502 :
503 : end if
504 :
505 : end do
506 :
507 5 : eval_local_min = fx
508 :
509 : return
510 : end function eval_local_min
511 :
512 :
513 10 : real(dp) function eval_brent_safe_zero ( a, b, machep, t, epsy, f, fa_in, fb_in, lrpar, rpar, lipar, ipar, ierr )
514 :
515 : ! *****************************************************************************80
516 : !
517 : !! seeks the root of a function F(X) in an interval [A,B].
518 : !
519 : ! Discussion:
520 : !
521 : ! The interval [A,B] must be a change of sign interval for F.
522 : ! That is, F(A) and F(B) must be of opposite signs. Then
523 : ! assuming that F is continuous implies the existence of at least
524 : ! one value C between A and B for which F(C) = 0.
525 : !
526 : ! The location of the zero is determined to within an accuracy
527 : ! of 6 * MACHEPS * abs ( C ) + 2 * T. or if abs(f(x)) < epsy.
528 : !
529 : ! Licensing:
530 : !
531 : ! This code is distributed under the GNU LGPL license.
532 : !
533 : ! Modified:
534 : !
535 : ! 12 April 2008
536 : !
537 : ! Author:
538 : !
539 : ! Original FORTRAN77 version by Richard Brent.
540 : ! FORTRAN90 version by John Burkardt.
541 : !
542 : ! Reference:
543 : !
544 : ! Richard Brent,
545 : ! Algorithms for Minimization Without Derivatives,
546 : ! Dover, 2002,
547 : ! ISBN: 0-486-41998-3,
548 : ! LC: QA402.5.B74.
549 : !
550 : ! Parameters:
551 : !
552 : ! Input, real ( kind = 8 ) A, B, the endpoints of the change of
553 : ! sign interval.
554 : !
555 : ! Input, real ( kind = 8 ) MACHEP, an estimate for the relative machine
556 : ! precision.
557 : !
558 : ! Input, real ( kind = 8 ) T, a positive error tolerance.
559 : !
560 : ! Input, external real ( kind = 8 ) F, the name of a user-supplied
561 : ! function, of the form "FUNCTION F ( X )", which evaluates the
562 : ! function whose zero is being sought.
563 : !
564 : ! Output, real ( kind = 8 ) ZERO, the estimated value of a zero of
565 : ! the function F.
566 : !
567 :
568 : interface
569 : #include "num_root_fcn.dek"
570 : end interface
571 : integer, intent(in) :: lipar, lrpar
572 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
573 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
574 : integer, intent(out) :: ierr
575 :
576 : real ( kind = 8 ) :: a
577 : real ( kind = 8 ) :: b
578 10 : real ( kind = 8 ) :: c
579 10 : real ( kind = 8 ) :: d
580 10 : real ( kind = 8 ) :: e
581 10 : real ( kind = 8 ) :: fa, fa_in
582 10 : real ( kind = 8 ) :: fb, fb_in
583 10 : real ( kind = 8 ) :: fc
584 10 : real ( kind = 8 ) :: m
585 : real ( kind = 8 ) :: machep
586 10 : real ( kind = 8 ) :: p
587 10 : real ( kind = 8 ) :: q
588 10 : real ( kind = 8 ) :: r
589 10 : real ( kind = 8 ) :: s
590 10 : real ( kind = 8 ) :: sa
591 : real ( kind = 8 ) :: sb
592 : real ( kind = 8 ) :: t, epsy
593 10 : real ( kind = 8 ) :: tol
594 10 : real ( kind = 8 ) :: dfdx
595 :
596 10 : ierr = 0
597 :
598 : !
599 : ! Make local copies of A and B.
600 : !
601 10 : sa = a
602 10 : sb = b
603 10 : fa = fa_in
604 : !fa = f ( sa, dfdx, lrpar, rpar, lipar, ipar, ierr )
605 : !if (ierr /= 0) return
606 10 : fb = fb_in
607 : !fb = f ( sb, dfdx, lrpar, rpar, lipar, ipar, ierr )
608 : !if (ierr /= 0) return
609 :
610 10 : c = sa
611 10 : fc = fa
612 10 : e = sb - sa
613 10 : d = e
614 :
615 : do
616 :
617 125 : if ( abs ( fc ) < abs ( fb ) ) then
618 :
619 35 : sa = sb
620 35 : sb = c
621 35 : c = sa
622 35 : fa = fb
623 35 : fb = fc
624 35 : fc = fa
625 :
626 : end if
627 :
628 125 : tol = 2.0D+00 * machep * abs ( sb ) + t
629 125 : m = 0.5D+00 * ( c - sb )
630 :
631 125 : if ( abs ( m ) <= tol .or. fb == 0.0D+00 ) then
632 : exit
633 : end if
634 :
635 118 : if ( abs ( e ) < tol .or. abs ( fa ) <= abs ( fb ) ) then
636 :
637 : e = m
638 : d = e
639 :
640 : else
641 :
642 115 : s = fb / fa
643 :
644 115 : if ( sa == c ) then
645 :
646 53 : p = 2.0D+00 * m * s
647 53 : q = 1.0D+00 - s
648 :
649 : else
650 :
651 62 : q = fa / fc
652 62 : r = fb / fc
653 62 : p = s * ( 2.0D+00 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0D+00 ) )
654 62 : q = ( q - 1.0D+00 ) * ( r - 1.0D+00 ) * ( s - 1.0D+00 )
655 :
656 : end if
657 :
658 115 : if ( 0.0D+00 < p ) then
659 48 : q = - q
660 : else
661 67 : p = - p
662 : end if
663 :
664 115 : s = e
665 115 : e = d
666 :
667 115 : if ( 2.0D+00 * p < 3.0D+00 * m * q - abs ( tol * q ) .and. &
668 : p < abs ( 0.5D+00 * s * q ) ) then
669 97 : d = p / q
670 : else
671 : e = m
672 : d = e
673 : end if
674 :
675 : end if
676 :
677 118 : sa = sb
678 118 : fa = fb
679 :
680 118 : if ( tol < abs ( d ) ) then
681 112 : sb = sb + d
682 6 : else if ( 0.0D+00 < m ) then
683 3 : sb = sb + tol
684 : else
685 3 : sb = sb - tol
686 : end if
687 :
688 118 : fb = f ( sb, dfdx, lrpar, rpar, lipar, ipar, ierr )
689 118 : if (ierr /= 0) return
690 118 : if (abs(fb) < epsy) exit
691 :
692 115 : if ( ( 0.0D+00 < fb .and. 0.0D+00 < fc ) .or. &
693 10 : ( fb <= 0.0D+00 .and. fc <= 0.0D+00 ) ) then
694 50 : c = sa
695 50 : fc = fa
696 50 : e = sb - sa
697 50 : d = e
698 : end if
699 :
700 : end do
701 :
702 10 : eval_brent_safe_zero = sb
703 :
704 10 : return
705 : end function eval_brent_safe_zero
706 :
707 : end module mod_brent
|
