line_search_interpolation Subroutine

   subroutine line_search_interpolation(ro, rl, ru, ri, fo, fl, fu, fi, &
                                        po, r, mode, mtyp, merr)

      real(wp) :: fo   !! value of the objective function for r=ro.
      real(wp) :: fl   !! value of the objective function for r=rl.
      real(wp) :: fu   !! value of the objective function for r=ru.
      real(wp) :: fi   !! value of the objective function for r=ri.
      real(wp) :: po   !! initial value of the directional derivative.
      real(wp) :: r    !! value of the stepsize parameter obtained.
      real(wp) :: rl   !! lower value of the stepsize parameter.
      real(wp) :: ru   !! upper value of the stepsize parameter.
      real(wp) :: ri   !! inner value of the stepsize parameter.
      real(wp) :: ro   !! initial value of the stepsize parameter.
      integer :: merr  !! error indicator. merr=0 for normal return.
      integer :: mode  !! mode of line search.
      integer :: mtyp  !! method selection
                       !!
                       !! * mtyp=1 - bisection.
                       !! * mtyp=2 - two point quadratic interpolation.
                       !! * mtyp=2 - three point quadratic interpolation.

      real(wp) :: ai, al, au, den, dis
      integer :: ntyp
      logical :: l1, l2

      real(wp), parameter :: zero = 0.0_wp
      real(wp), parameter :: half = 0.5_wp
      real(wp), parameter :: one = 1.0_wp
      real(wp), parameter :: two = 2.0_wp
      real(wp), parameter :: three = 3.0_wp
      real(wp), parameter :: c1l = 1.1_wp
      real(wp), parameter :: c1u = 1000.0_wp
      real(wp), parameter :: c2l = 1.0e-2_wp
      real(wp), parameter :: c2u = 0.9_wp
      real(wp), parameter :: c3l = 1.0e-1_wp

      merr = 0
      if (mode <= 0) return
      if (po >= zero) then
         merr = 2
         return
      elseif (ru <= rl) then
         merr = 3
         return
      end if
      l1 = rl <= ro
      l2 = ri <= rl
      main: do ntyp = mtyp, 1, -1
         if (ntyp == 1) then
            ! bisection
            if (mode == 1) then
               r = two*ru
               return
            elseif (ri - rl <= ru - ri) then
               r = half*(ri + ru)
               return
            else
               r = half*(rl + ri)
               return
            end if
         elseif (ntyp == mtyp .and. l1) then
            if (.not. l2) ai = (fi - fo)/(ri*po)
            au = (fu - fo)/(ru*po)
         end if
         if (l1 .and. (ntyp == 2 .or. l2)) then
            ! two point quadratic extrapolation or interpolation
            if (au >= one) cycle main
            r = half*ru/(one - au)
         elseif (.not. l1 .or. .not. l2 .and. ntyp == 3) then
            ! three point quadratic extrapolation or interpolation
            al = (fi - fl)/(ri - rl)
            au = (fu - fi)/(ru - ri)
            den = au - al
            if (den <= zero) cycle main
            r = ri - half*(au*(ri - rl) + al*(ru - ri))/den
         elseif (l1 .and. .not. l2 .and. ntyp == 4) then
            ! three point cubic extrapolation or interpolation
            dis = (ai - one)*(ru/ri)
            den = (au - one)*(ri/ru) - dis
            dis = au + ai - den - two*(one + dis)
            dis = den*den - three*dis
            if (dis < zero) cycle main
            den = den + sqrt(dis)
            if (den == zero) cycle main
            r = (ru - ri)/den
         else
            cycle main
         end if
         if (mode == 1 .and. r > ru) then
            ! extrapolation accepted
            r = max(r, c1l*ru)
            r = min(r, c1u*ru)
            return
         elseif (mode == 2 .and. r > rl .and. r < ru) then
            ! interpolation accepted
            if (ri == zero .and. ntyp /= 4) then
               r = max(r, rl + c2l*(ru - rl))
            else
               r = max(r, rl + c3l*(ru - rl))
            end if
            r = min(r, rl + c2u*(ru - rl))
            if (r /= ri) return
         end if
      end do main

   end subroutine line_search_interpolation