linmin Function

private function linmin(mode, ax, bx, f, tol, a, b, d, e, p, q, r, u, v, w, x, m, fu, fv, fw, fx, tol1, tol2)

Linesearch without derivatives (used by slsqp if linesearch_mode=2). Returns the abscissa approximating the point where f attains a minimum.

purpose:

to find the argument linmin where the function f takes it's minimum on the interval ax, bx. It uses a combination of golden section and successive quadratic interpolation.

Reference

This function subprogram is a slightly modified version of the ALGOL 60 procedure localmin given in R.P. Brent: "Algorithms for minimization without derivatives", Prentice-Hall (1973).

History

  • Kraft, D., DFVLR - institut fuer dynamik der flugsysteme d-8031 oberpfaffenhofen
  • status: 31. august 1984
  • Jacob Williams, Jan 2016, Refactored into modern Fortran. Added saved variables as inouts to make the routine thread-safe.

Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: mode

controls reverse communication must be set to 0 initially, returns with intermediate values 1 and 2 which must not be changed by the user, ends with convergence with value 3.
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) :: f

function value at linmin which is to be brought in by reverse communication controlled by mode
real(kind=wp), intent(in) :: tol

desired length of interval of uncertainty of final result
real(kind=wp), intent(inout) :: a
real(kind=wp), intent(inout) :: b
real(kind=wp), intent(inout) :: d
real(kind=wp), intent(inout) :: e
real(kind=wp), intent(inout) :: p
real(kind=wp), intent(inout) :: q
real(kind=wp), intent(inout) :: r
real(kind=wp), intent(inout) :: u
real(kind=wp), intent(inout) :: v
real(kind=wp), intent(inout) :: w
real(kind=wp), intent(inout) :: x
real(kind=wp), intent(inout) :: m
real(kind=wp), intent(inout) :: fu
real(kind=wp), intent(inout) :: fv
real(kind=wp), intent(inout) :: fw
real(kind=wp), intent(inout) :: fx
real(kind=wp), intent(inout) :: tol1
real(kind=wp), intent(inout) :: tol2

Return Value real(kind=wp)


Called by

proc~~linmin~~CalledByGraph proc~linmin slsqp_core::linmin proc~slsqpb slsqp_core::slsqpb proc~slsqpb->proc~linmin proc~slsqp slsqp_core::slsqp proc~slsqp->proc~slsqpb proc~slsqp_wrapper slsqp_module::slsqp_solver%slsqp_wrapper proc~slsqp_wrapper->proc~slsqp

Source Code

    real(wp) function linmin(mode,ax,bx,f,tol,&
                                a,b,d,e,p,q,r,u,v,&
                                w,x,m,fu,fv,fw,fx,tol1,tol2)

    implicit none

    integer,intent(inout) :: mode  !! controls reverse communication
                                   !! must be set to 0 initially, returns with intermediate
                                   !! values 1 and 2 which must not be changed by the user,
                                   !! ends with convergence with value 3.
    real(wp) :: f                  !! function value at `linmin` which is to be brought in by
                                   !! reverse communication controlled by `mode`
    real(wp),intent(in) :: tol     !! desired length of interval of uncertainty of final result
    real(wp),intent(in) :: ax      !! left endpoint of initial interval
    real(wp),intent(in) :: bx      !! right endpoint of initial interval
    real(wp),intent(inout) :: a,b,d,e,p,q,r,u,v,w,x,m,fu,fv,fw,fx,tol1,tol2

    real(wp),parameter :: c = (3.0_wp-sqrt(5.0_wp))/2.0_wp  !! golden section ratio = `0.381966011`
    real(wp),parameter :: sqrteps = sqrt(eps)  !! square root of machine precision

    select case (mode)
    case(1)

        ! main loop starts here

        fx = f
        fv = fx
        fw = fv

    case(2)

        fu = f
        ! update a, b, v, w, and x
        if ( fu>fx ) then
            if ( u<x ) a = u
            if ( u>=x ) b = u
            if ( fu<=fw .or. abs(w-x)<=zero ) then
                v = w
                fv = fw
                w = u
                fw = fu
            else if ( fu<=fv .or. abs(v-x)<=zero .or. abs(v-w)<=zero ) then
                v = u
                fv = fu
            end if
        else
            if ( u>=x ) a = x
            if ( u<x ) b = x
            v = w
            fv = fw
            w = x
            fw = fx
            x = u
            fx = fu
        end if

    case default

        ! initialization
        a = ax
        b = bx
        e = zero
        v = a + c*(b-a)
        w = v
        x = w
        linmin = x
        mode = 1
        return

    end select

    m = 0.5_wp*(a+b)
    tol1 = sqrteps*abs(x) + tol
    tol2 = tol1 + tol1

    ! test convergence

    if ( abs(x-m)<=tol2-0.5_wp*(b-a) ) then
        ! end of main loop
        linmin = x
        mode = 3
    else
        r = zero
        q = r
        p = q
        if ( abs(e)>tol1 ) then
            ! fit parabola
            r = (x-w)*(fx-fv)
            q = (x-v)*(fx-fw)
            p = (x-v)*q - (x-w)*r
            q = q - r
            q = q + q
            if ( q>zero ) p = -p
            if ( q<zero ) q = -q
            r = e
            e = d
        end if

        ! is parabola acceptable
        if ( abs(p)>=0.5_wp*abs(q*r) .or. p<=q*(a-x) .or. p>=q*(b-x) ) then
            ! golden section step
            if ( x>=m ) e = a - x
            if ( x<m ) e = b - x
            d = c*e
        else
            ! parabolic interpolation step
            d = p/q
            ! f must not be evaluated too close to a or b
            if ( u-a<tol2 ) d = sign(tol1,m-x)
            if ( b-u<tol2 ) d = sign(tol1,m-x)
        end if

        ! f must not be evaluated too close to x
        if ( abs(d)<tol1 ) d = sign(tol1,d)
        u = x + d
        linmin = u
        mode = 2

    end if

    end function linmin