fmin Function

private function fmin(me, ax, bx, tol) result(xmin)

An approximation x to the point where f attains a minimum on the interval (ax,bx) is determined.

the method used is a combination of golden section search and successive parabolic interpolation. convergence is never much slower than that for a fibonacci search. if f has a continuous second derivative which is positive at the minimum (which is not at ax or bx), then convergence is superlinear, and usually of the order of about 1.324.

the function f is never evaluated at two points closer together than epsabs(fmin) + (tol/3), where eps is approximately the square root of the relative machine precision. if f is a unimodal function and the computed values of f are always unimodal when separated by at least epsabs(x) + (tol/3), then fmin approximates the abcissa of the global minimum of f on the interval ax,bx with an error less than 3epsabs(fmin) + tol. if f is not unimodal, then fmin may approximate a local, but perhaps non-global, minimum to the same accuracy.

this function subprogram is a slightly modified version of the algol 60 procedure localmin given in richard brent, algorithms for minimization without derivatives, prentice - hall, inc. (1973).

See also

[1] http://www.netlib.org/fmm/fmin.f

Type Bound

brent_class

Arguments

Type IntentOptional Attributes Name
class(brent_class), intent(inout) :: me
real(kind=wp), intent(in) :: ax

left endpoint of initial interval

real(kind=wp), intent(in) :: bx

right endpoint of initial interval

real(kind=wp), intent(in) :: tol

desired length of the interval of uncertainty of the final result (>=0)

Return Value real(kind=wp)

abcissa approximating the point where f attains a minimum


Called by

proc~~fmin~~CalledByGraph proc~fmin brent_module::brent_class%fmin proc~brent_test brent_module::brent_test proc~brent_test->proc~fmin

Source Code

    function fmin(me,ax,bx,tol) result(xmin)

    implicit none

    class(brent_class),intent(inout) :: me
    real(wp),intent(in) :: ax   !! left endpoint of initial interval
    real(wp),intent(in) :: bx   !! right endpoint of initial interval
    real(wp),intent(in) :: tol  !! desired length of the interval of
                                !! uncertainty of the final result (>=0)
    real(wp)            :: xmin !! abcissa approximating the point where
                                !! f attains a minimum

    real(wp) :: a,b,d,e,xm,p,q,r,tol1,tol2,u,v,w
    real(wp) :: fu,fv,fw,fx,x
    real(wp) :: abs,sqrt,sign

    real(wp),parameter :: c = (three-sqrt(five))/two  !! squared inverse of golden ratio
    real(wp),parameter :: half = 0.5_wp
    real(wp),parameter :: eps = sqrt(epsilon(one))

    ! initialization

    a = ax
    b = bx
    v = a + c*(b - a)
    w = v
    x = v
    e = zero
    fx = me%f(x)
    fv = fx
    fw = fx

    do    !  main loop starts here

        xm = half*(a + b)
        tol1 = eps*abs(x) + tol/three
        tol2 = two*tol1

        !  check stopping criterion

        if (abs(x - xm) <= (tol2 - half*(b - a))) exit

        ! is golden-section necessary

        if (abs(e) <= tol1) then

            !  a golden-section step

            if (x >= xm) then
                e = a - x
            else
                e = b - x
            end if
            d = c*e

        else

            !  fit parabola

            r = (x - w)*(fx - fv)
            q = (x - v)*(fx - fw)
            p = (x - v)*q - (x - w)*r
            q = two*(q - r)
            if (q > zero) p = -p
            q =  abs(q)
            r = e
            e = d

            !  is parabola acceptable

            if ((abs(p) >= abs(half*q*r)) .or. (p <= q*(a - x)) .or. (p >= q*(b - x))) then

                !  a golden-section step

                if (x >= xm) then
                    e = a - x
                else
                    e = b - x
                end if
                d = c*e

            else

                !  a parabolic interpolation step

                d = p/q
                u = x + d

                !  f must not be evaluated too close to ax or bx

                if (((u - a) < tol2) .or. ((b - u) < tol2)) d = sign(tol1, xm - x)

            end if

        end if

        !  f must not be evaluated too close to x

        if (abs(d) >= tol1) then
            u = x + d
        else
            u = x + sign(tol1, d)
        end if
        fu = me%f(u)

        !  update a, b, v, w, and x

        if (fu <= fx) then
            if (u >= x) a = x
            if (u < x) b = x
            v = w
            fv = fw
            w = x
            fw = fx
            x = u
            fx = fu
        else
            if (u < x) a = u
            if (u >= x) b = u
            if (fu <= fw .or. w == x) then
                v = w
                fv = fw
                w = u
                fw = fu
            else if (fu <= fv .or. v == x .or. v == w ) then
                v = u
                fv = fu
            end if
        end if

    end do    !  end of main loop

    xmin = x

    end function fmin