!*******************************************************************************
!       Start of the header for a fortran source file for a subroutine
!       of the Free_EOS stellar interior equation of state code
!       Copyright (C) 2006 Alan W. Irwin
!
!       $Id: bfgs.f,v 1.12 2007/04/01 00:16:39 airwin Exp $
!
!       For the latest version of this source code, please contact
!       Alan W. Irwin
!       Department of Physics and Astronomy
!       University of Victoria,
!       Box 3055
!       Victoria, B.C., Canada
!       V8W 3P6
!       e-mail: irwin@beluga.phys.uvic.ca.
!
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!       End of the header for a fortran source file for a subroutine
!       of the Free_EOS stellar interior equation of state code
!*************************************************************

!      BFGS algorithm based on 
!      R. Fletcher, "Practical Methods of Optimization" 2nd ed.
!      John Wiley and Sons, 1987.
!      N.B. in following comments this reference is referred to as PMOO.


   module bfgs_irwin
   implicit none
   contains

      subroutine bfgs_constant_vector(n, constant, x)
!      set vector x(n) to constant
      implicit none
      integer*4 n,i
      real*8 constant, x(n)
      do i = 1,n
        x(i) = constant
      enddo
      end subroutine

      subroutine bfgs_set(status, name,
     &  n, x, xi, f, gradient, p)
!      status is a control string that takes the following values:
!        'function'.  Request function value to be supplied externally
!        'gradient'.  Request gradient to be supplied externally.
!           N.B. gradient evaluation is always done at xi of prior
!           function evaluation
!        'both'.  Request both function value and gradient to be supplied
!                 externally.
!        'continue'. Continue after call to this routine rather than
!                    returning immediately
!        'error....' some aspect of line search or bfgs update has not
!                    completed successfully.  The 'error' preamble is then
!                    followed with a specific error message.

!      initialize line search
      implicit none
      include 'bfgs.h'
      integer*4 n
      real*8 p(n), x(n), xi, f, gradient
      character*(*) status, name
      save
!      sanity checks:
      if(len(status).lt.8)
     &  stop 'bfgs_set: length of status string too small'
      if(len(name).lt.3)
     &  stop 'bfgs_set:length of name string too small'
      if(status(:4).ne.'both') then
        call dcopy(n, x, 1, xi, 1)
        status = 'both'
        name = 'set'
      else
!        initial preferred direction for the line search is in the
!        direction of the gradient, i.e, steepest ascent.  (Negative
!        sign of direction, i.e., initial steepest descent
!        figured out later.)
        call dcopy(n, gradient, 1, p, 1)
!        Set to true for Fletcher estimate of alpha1 and set to false for
!        GSL estimate of alpha1 (or Numerical Recipes) estimate of alpha1.
        fletcher_estimate = .true.
        if(fletcher_estimate) then
!         _use_large initial value of deltaf_previous to force an initial
!          step size of unity as per PMOO, 2.6.8.
          deltaf_previous = 1.d300
        else
!          set alpha_previous to 1.d0 to force an initial step size
!          of unity for the GSL (or Numerical Recipes) estimate.
          alpha_previous = 1.d0
        endif
        status = 'continue'
      endif
      end subroutine
      
      subroutine bfgs_take_step(n, x, p, step, lambda, x1, dx)
!      take step in the desired -lambda*p direction.
!      unchanged input:
!      x(n) is the current position vector
!      lambda*p(n) is the unit vector in the currently desired uphill
!      (not necessarily steepest ascent) direction.
!      output results:
!      dx(n) = -step*lambda*p
!      x1(n) = x + dx
      implicit none
      integer*4 n
      real*8 x(n), p(n), step, lambda, x1(n), dx(1)
      call bfgs_constant_vector(n, 0.d0, dx)
      call daxpy (n, -step*lambda, p, 1, dx, 1)
      call dcopy(n, x, 1, x1, 1)
      call daxpy (n, 1.d0, dx, 1, x1, 1)
      end subroutine

      real*8 function bfgs_cubic_minimum(
     &  a, fa, fprimea, b, fb, fprimeb, x_lo, x_hi)
!     _Use_cubic interpolation or extrapolation to calculate
!      minimum value within a specified range.

!      input variables:
!      a is first defined point
!      fa and fprimea are the function and derivative values at that point.
!      b is the second defined point
!      fb and fprimeb are the function and derivative values at that point.
!      x_lo through x_hi specify the range where the minimum value
!      is to be calculated.

!      The returned function value is the independent variable x where
!      the minimum of the interpolated/extrapolated function occurs
!      within [x_lo, x_hi].

      implicit none
      real*8 a, fa, fprimea, b, fb, fprimeb, x_lo, x_hi
      real*8 dx, q0, q1, q2, q3, z_lo, z_hi,
     &  q_lo, q_hi, z_stationary, q_stationary, minimum,
     &  radical, quadratic_factor
      integer*4 nsolution

      dx = b - a
!      Find coefficients of Hermite polynomial on z which ranges from
!      0 to 1 within defined interval [a, b], which agrees with fa and
!      fprimea*dx at 0 and which agrees with fb and fprimeb*dx at 1.
!      z = (x - a)/dx or x = a + z*dx
!      f(z) = q0 + z*(q1 + z*(q2 + z*q3))
!      df(z) = q1 + 2.d0*q2*z + 3.d0*q3*z*z
      q0 = fa
      q1 = fprimea*dx
      q2 = 3.d0*(fb-fa) - (2.d0*fprimea + fprimeb)*dx
      q3 = (fprimea+fprimeb)*dx - 2.d0*(fb-fa)
!      N.B.
!      f(0) = q0 = fa
!      f(1) = q0 + q1 + q2 + q3 = fb
!      df(0) = q1 = fprimea*dx
!      df(1) = q1 + 2*q2 + 3*q3 = ((1-4+3)*fprimea + (-2+3)*fprimeb)*dx
!            = fprimeb*dx
!      z_lo and z_hi correspond to x_lo and x_hi.
      z_lo = (x_lo - a)/dx
      z_hi = (x_hi - a)/dx
      q_lo = q0 + z_lo*(q1 + z_lo*(q2 + z_lo*q3))
      q_hi = q0 + z_hi*(q1 + z_hi*(q2 + z_hi*q3))
      if(q_lo.lt.q_hi) then
        bfgs_cubic_minimum = x_lo
        minimum = q_lo
      else
        bfgs_cubic_minimum = x_hi
        minimum = q_hi
      endif
      if(q3.eq.0.d0) then
        nsolution = 1
        z_stationary = -q1*0.5d0/q2
      else
!        radical = b^2 - 4 a c
        radical = 4.d0*q2*q2 - 12.d0*q1*q3
        if(radical.lt.0.d0) then
!          no roots
          nsolution = 0
        elseif(radical.eq.0.d0) then
!          one root.
          nsolution = 1
          z_stationary = - q2/(3.d0*q3)
        else
!          two roots.
          nsolution = 2
          if(q2.gt.0) then
            quadratic_factor = 2.d0*q2 + sqrt(radical)
          else
            quadratic_factor = 2.d0*q2 - sqrt(radical)
          endif
          z_stationary = -quadratic_factor/(6.d0*q3)
        endif
      endif
      if(nsolution.ge.1) then
!        at least one possible solution calculated at z_stationary
        if((z_lo.le.z_stationary.and.z_stationary.le.z_hi).or.
     &      (z_hi.le.z_stationary.and.z_stationary.le.z_lo)) then
          q_stationary = q0 + z_stationary*(q1 + z_stationary*
     &      (q2 + z_stationary*q3))
          if(q_stationary.lt.minimum) then
            bfgs_cubic_minimum = a + z_stationary*dx
            minimum = q_stationary
          endif
        endif
        if(nsolution.eq.2) then
!          try alternate quadratic solution
          z_stationary = -2.d0*q1/quadratic_factor
          if((z_lo.le.z_stationary.and.z_stationary.le.z_hi).or.
     &        (z_hi.le.z_stationary.and.z_stationary.le.z_lo)) then
            q_stationary = q0 + z_stationary*(q1 + z_stationary*
     &        (q2 + z_stationary*q3))
            if(q_stationary.lt.minimum) then
              bfgs_cubic_minimum = a + z_stationary*dx
            endif
          endif
        endif
      endif
      end function

      real*8 function bfgs_quadratic_minimum(
     &  a, fa, fprimea, b, fb, x_lo, x_hi)
!     _Use_quadratic interpolation or extrapolation to calculate
!      minimum value within a specified range.

!      input variables:
!      a is first defined point
!      fa and fprimea are the function and derivative values at that point.
!      b is the second defined point
!      fb is the function value at that point.
!      x_lo through x_hi specify the range where the minimum value
!      is to be calculated.

!      The returned function value is the independent variable x where
!      the minimum of the interpolated/extrapolated function occurs
!      within [x_lo, x_hi].

      implicit none
      real*8 a, fa, fprimea, b, fb, x_lo, x_hi
      real*8 dx, q0, q1, q2, z_lo, z_hi,
     &  q_lo, q_hi, z_stationary, q_stationary, minimum
      integer*4 nsolution

      dx = b - a
!      Find coefficients of Hermite polynomial on z which ranges from
!      0 to 1 within defined interval [a, b], which agrees with fa and
!      fprimea*dx at 0 and which agrees with fb at 1.
!      z = (x - a)/dx or x = a + z*dx
!      f(z) = q0 + z*(q1 + z*q2)
!      df(z) = q1 + 2.d0*q2*z with stationary point z = -q1*0.5/q2
      q0 = fa
      q1 = fprimea*dx
      q2 = fb-fa-q1
!      z_lo and z_hi correspond to x_lo and x_hi.
      z_lo = (x_lo - a)/dx
      z_hi = (x_hi - a)/dx
      q_lo = q0 + z_lo*(q1 + z_lo*q2)
      q_hi = q0 + z_hi*(q1 + z_hi*q2)
      if(q_lo.lt.q_hi) then
        bfgs_quadratic_minimum = x_lo
        minimum = q_lo
      else
        bfgs_quadratic_minimum = x_hi
        minimum = q_hi
      endif
      if(q2.eq.0.d0) then
        nsolution = 0
      else
        nsolution = 1
        z_stationary = -q1*0.5d0/q2
      endif
      if(nsolution.eq.1) then
        if((z_lo.le.z_stationary.and.z_stationary.le.z_hi).or.
     &      (z_hi.le.z_stationary.and.z_stationary.le.z_lo)) then
          q_stationary = q0 + z_stationary*(q1 + z_stationary*q2)
          if(q_stationary.lt.minimum) then
            bfgs_quadratic_minimum = a + z_stationary*dx
          endif
        endif
      endif
      end function

      subroutine bfgs_linesearch(status, name,
     &  fbar, epsilon, n, p, x, xi, f, gradient, alpha_ratio, dx)
!      Do a linesearch using algorithm presented by PMOO, pp. 34-39

!      status is a control string that takes the following values:
!        'function'.  Request function value to be supplied externally
!        'gradient'.  Request gradient to be supplied externally.
!           N.B. gradient evaluation is always done at xi of prior
!           function evaluation
!        'both'.  Request both function value and gradient to be supplied
!                 externally.
!        'continue'. Continue after call to this routine rather than
!                    returning immediately
!        'error....' some aspect of line search or bfgs update has not
!                    completed successfully.  The 'error' preamble is then
!                    followed with a specific error message.

!      name is a control string set to 'linesearch'
!      when there is a 'function', 'gradient', or 'both' request or
!      when there is an error.

!      input quantities:
!      fbar is a user-specified minimum possible f used as per
!        PMOO, 2.6.1 to control the bracketing step size.  If actual f
!        values are <= fbar, then this routine stops with an error so
!        be realistic in how you specify fbar making it small enough
!        to avoid the error, but large enough to provide some control
!        over the maximum size of the bracketing step.
!      epsilon is the user-specified convergence criterion for f (not
!        relative f) discussed on p. 38 of PMOO.  It is used to terminate
!        the line-search routine to work around the case when round-off
!        errors are causing convergence problems close to the minimum.
!      p(n) is the unscaled line search direction.  The PMOO scaled s vector
!        is defined by
!        s = -(dir/pnorm)*p
!        where dir = 1 if p has an acute angle with the input gradient and
!        dir = -1 otherwise.  pnorm is the norm of p.  dir and pnorm are
!        calculated internally.

!      input and output quantities:
!      x(n) is the starting point of the line search on input and on output
!        is the ending point of the line search.
!      N.B. both f and gradient must be precalculated on input and calculated
!        on output.  That is:
!      f is f(x(n)) on both input and output.
!      gradient(n) is the gradient(x(n)) on both input and output.

!      output quantities:
!      alpha_ratio is the ratio of the calculated step size to the initial
!      estimate, alpha1 for the step size.
!      dx(n) is the vector of differences between the initial x and final x.
      implicit none
      include 'bfgs.h'
      integer*4 n
      real*8 fbar, epsilon,
     &  p(n), x(n), xi(n), f, gradient(n), alpha_ratio, dx(n)
!      dimensionless line search constants recommended by PMOO, page. 37
!      (superseded by PMOO, page 69 for tau2 = 0.05 rather than tau2 = 0.10).
      real*8 sigma, rho, tau1, tau2, tau3
!     _use_fairly accurate line search.
      parameter(sigma = 0.1d0)
!      rho must be less than or equal to sigma
      parameter(rho = 0.01d0)
!      jump size factor increase used in bracketing phase.
      parameter(tau1 = 9.d0)
!      0 < tau2 < tau3 <= 0.5 and tau2 <= sigma is advisable
      parameter(tau2 = 0.05d0)
      parameter(tau3 = 0.5d0)
      logical bracket
      real*8 alpha1, pnorm, dnrm2, pg, ddot, dir,
     &  f0, fprime0, mu,
     &  alphaim, fim, fprimeim, alphai, fi, fprimei,
     &  ai, fai, fprimeai, bi, fbi, fprimebi,
     &  dalpha, alphaip, bfgs_quadratic_minimum, bfgs_cubic_minimum
      logical acceptable, fprimebi_calc, ifcubic
!     _Use_cubic interpolation/extrapolation if fprimebi_calc is true
      parameter(ifcubic=.true.)
      real*8 aip, bip, fmin
      logical debug
      parameter(debug=.false.)
      character*(*) status, name
!      N.B. the *need logic replaces go to logic with typically three or four
!      if statements. However, I decided to_use_the more complicated *need
!      approach because the g77 compiler generated all sorts of warnings
!      about the go tos into the middle of control blocks leading
!      to uncertainty and doubt about whether the go to logic would actually
!      work correctly for most compilers/optimizers.

!      N.B. the various need variables and the not_need variables derived
!      from them are used to control re-entry into the routine with
!      the 'function', 'gradient', or 'both' data that are needed.
!      to understand the ordinary flow of this routine simply ignore
!      how the need and not_need variables change that flow.
!      To understand the re-entry flow of control, note that only
!      at most one of need1,...,need5 are true at any given time, and
!      for the case when one of those are true, you should be able to
!      swiftly follow the logic back to where the data are needed and
!      continue from there.
      logical need1, need2, need3, need4, need5,
     &  need12, need345,
     &  not_need12345, not_need345, not_need12, not_need2,
     &  not_need45, not_need3, not_need5
      data need1, need2, need3, need4, need5/5*.false./
      save
!      if(need1) then
!        go to 1
!      elseif(need2) then
!        go to 2
!      elseif(need3) then
!        go to 3
!      elseif(need4) then
!        go to 4
!      elseif(need5) then
!        go to 5
!      endif
      need12 = need1.or.need2
      need345 = need3.or.need4.or.need5
      not_need12345 = .not.(need12.or.need345)
      not_need345 = .not.need345
      not_need12 = .not.need12
      not_need2 = .not.need2
      not_need45 = .not.(need4.or.need5)
      not_need3 = .not.need3
      not_need5 = .not.need5
      if(not_need12345) then
        
!        sanity check
        if(n.gt.nmax_bfgs) stop 'bfgs_linesearch: n too large'
        if(len(status).lt.8)
     &    stop 'bfgs_linesearch: length of status string too small'
        if(len(name).lt.11)
     &    stop 'bfgs_linesearch: length of name string too small'
        
!        parameters of PMOO scaled s vector where s = -(dir/pnorm)*p
        pnorm = dnrm2(n, p, 1)
!        pg subsequently used to calculate fprime0
        pg = ddot(n, p, 1, gradient, 1)
        dir = sign(1.d0, pg)

!        early convergence test to avoid divide by zero.
!        Note, pnorm can be zero, if p is initialized to the gradient and
!        the gradient is zero.
        if(pnorm.eq.0.d0) then
          call bfgs_constant_vector(n, 0.d0, dx)
          alpha_ratio = 0.d0
          alpha_previous = 0.d0
          deltaf_previous = 0.d0
          status = 'continue'
          return
        endif
!        Initialize "0" variables:
        f0 = f
!        fprime = s dot gradient (PMOO, 1.2.6)
        fprime0 = -(dir/pnorm)*pg

!        early convergence test on fprime0 to avoid a divide by zero
        if(fprime0.ge.0.d0) then
          call bfgs_constant_vector(n, 0.d0, dx)
          alpha_ratio = 0.d0
          alpha_previous = 0.d0
          deltaf_previous = 0.d0
          status = 'continue'
          return
        endif
        if(f0.le.fbar) then
          write(0,*) 'bfgs_linesearch: ERROR f0 = f(alpha=0) '//
     &      '<= fbar, the minimum possible function value.'
          write(0,*) 'bfgs_linesearch: respecify fbar and try again.'
          stop
        else
!          maximum line search range from PMOO, 2.6.1
!          Note from above test this must always be positive
!          (unless underflow zeroes it)
          mu = (fbar - f0)/(rho*fprime0)
        endif
        
        if(fletcher_estimate) then
          if(2.d0*max(deltaf_previous,10.d0*epsilon).gt.
     &        -fprime0) then
!            -2 deltaf_previous/fprime0 will generally be close to
!            norm(x).  Thus, this branch generally taken when
!            norm(x)> 1.  alpha1 = 1 implies a relative change in x
!            of 1/norm(x).
!            According to PMOO discussion this gives rapid ultimate
!            convergence for Newton-like methods.
!            PMOO, 2.6.8
            alpha1 = 1.d0
          else
!            -2 deltaf_previous/fprime0 will generally be close to
!            norm(x).  Thus, this branch generally taken when
!            norm(x) < 1.  Thus, the alpha1 extimate below (if we ignore the
!            10.d0*epsilon roundoff error safeguard) implies
!            the relative change in x is generally close to unity for
!            this branch.
!            PMOO, 2.6.7.
            alpha1 = 2.d0*max(deltaf_previous,10.d0*epsilon)/
     &        (-fprime0)
          endif
        else
!          Estimate of initial step size based on best alpha from previous
!          iteration.  This is the GSL estimate.
!          In final super-linear convergence stage this
!          should always be an overestimate, but quadratic interpolation
!          should save the day to get an immediate reasonable minimum for
!          small enough alpha_previous
          alpha1 = alpha_previous
!          Numerical Recipes estimate
!          According to PMOO discussion this gives rapid ultimate convergence
!          for Newton-like methods.
          !alpha1 = 1.d0
        endif
        if(alpha1.le.0.d0) then
          write(0,*) 'alpha1 = ', alpha1
          stop 'bfgs_linesearch: bad initial estimate of step size'
        endif
!        Initialize "i-1" iteration variables (which have "im" suffix).
        alphaim = 0.d0
        fim = f0
        fprimeim = fprime0
!        note condition that fim <= f0+0.*rho*fprime0 is
!        automatically satisfied.
        fmin = fim
        if(debug) then
          write(0,*) "alpha, f(alpha) = ", alphaim, fim
          write(0,*) "alpha, f'(alpha) =", alphaim, fprimeim
        endif
        
!        Initialize "i" iteration variables (which have "i" suffix).
        alphai = alpha1
        
!        force at least one bracketing attempt.
        bracket = .false.
      endif  !not_need12345
      if(not_need345) then
        do while (need12.or..not.bracket)
          if(not_need12) then
!            All code in this loop follows pseudo-code in PMOO, 2.6.2.
            
!            Compute new trial point at alphai corresponding to
!            xi = x + alphai * s
            call bfgs_take_step(n, x, p, alphai, dir/pnorm, xi, dx)
!            EVALUATE f(alphai) = function at xi
            status = 'function'
            name = 'linesearch'
            need1 = .true.
            return
          endif  !not_need12
          if(need1) then
! 1          continue
            need1 = .false.
            need12 = need1.or.need2
            not_need12345 = .not.(need12.or.need345)
            not_need12 = .not.need12
            fi = f
            if(fi.le.f0+alphai*rho*fprime0) fmin = min(fi, fmin)
            if(debug) then
              write(0,*) "alpha, f(alpha) = ", alphai, fi
            endif
            if(fi.le.fbar) then
              write(0,*) 'bfgs_linesearch: ERROR fi = f(alphai) '//
     &          '<= fbar, the minimum possible function value.'
              write(0,*) 'bfgs_linesearch: respecify fbar and '//
     &          'try again.'
              stop
            endif
!            N.B. in my copy of PMO, there is a misprint in the 2.6.2
!            pseudo-code where the rho factor in the next uncommented line is
!            ignored.  But the rho factor must be there if the conditions
!            given by 2.6.3 are to be satisfied at the end of the bracket
!            loop.
          endif  !need1
          if(not_need2.and.
     &        (fi.gt.f0+alphai*rho*fprime0.or.fi.ge.fim)) then
            ai = alphaim
            fai = fim
            fprimeai = fprimeim
            bi = alphai
            fbi = fi
!            At this stage in bracketing iteration gradient for xi
!            is _not_ known.
            fprimebi_calc = .false.
            bracket = .true.
          else
            if(not_need2) then
!              EVALUATE gradient of function at same xi as previous.
              status = 'gradient'
              name = 'linesearch'
              need2 = .true.
              return
            endif  !not_need2
            if(need2) then
! 2            continue
              need2 = .false.
              need12 = need1.or.need2
              not_need12345 = .not.(need12.or.need345)
              not_need12 = .not.need12
              not_need2 = .not.need2
!              from PMOO, 1.2.6
!              df(alphai) = s dot gradient
              fprimei = (-dir/pnorm)*ddot(n, p, 1, gradient, 1)
              if(debug) then
                write(0,*) "alpha, f'(alpha) =", alphai, fprimei
              endif
!              good converged solution to line search if xi satisfies
!              two-sided test (PMOO, equation 2.5.6)
              if(abs(fprimei).le.-sigma*fprime0) then
!                acceptable point found.
!                output x = xi
                call dcopy(n, xi, 1, x, 1)
!                output f = f(xi)
                f = fi
!                gradient at xi already stored ready for output.
                alpha_ratio = alphai/alpha1
                alpha_previous = alphai
                deltaf_previous = f0-fi
                status = 'continue'
                return
              endif
              if(fprimei.ge.0.d0) then
                ai = alphai
                fai = fi
                fprimeai = fprimei
                bi = alphaim
!                fprimeim always known.
                fprimebi_calc = .true.
                fprimebi = fprimeim
                fbi = fim
                bracket = .true.
              else
                dalpha = alphai - alphaim
                if(mu.le.alphai + dalpha) then
                  alphaip = mu
                else
!                  extrapolation to find possible bracket in the range
!                  [alphai + dalpha, min(mu, alphai + tau*dalpha)]
                  if(ifcubic) then
                    alphaip = bfgs_cubic_minimum(
     &                alphaim, fim, fprimeim, alphai, fi, fprimei,
     &                alphai + dalpha, min(mu, alphai + tau1*dalpha))
                  else
                    alphaip = bfgs_quadratic_minimum(
     &                alphaim, fim, fprimeim, alphai, fi,
     &                alphai + dalpha, min(mu, alphai + tau1*dalpha))
                  endif
                endif
!                prepare for the next bracket attempt
                alphaim = alphai
                fim = fi
                fprimeim = fprimei
                alphai = alphaip
              endif
            else
              stop 'bfgs_linesearch: bad logic1'
            endif  !need2
          endif
        enddo
!        found bracket where the following conditions occur:
!        (i) ai is the current best trial point (least f) that
!          satisfies PMOO, 2.5.1, i.e., fai = f(ai) <= f(0) + ai*rho*f'(0)
!        (ii) fprimeai = f'(ai) has been evaluated and satisfies
!          (bi-ai)*f'(ai) <0 but |f'(ai)| > -sigma*f'(0)
!        (iii) bi satisfies either fbi = f(bi) > f(0) + bi*rho*f'(0) or
!          fbi = f(bi) >= f(ai) or both conditions are true.
!        PMOO, Lemma 2.6.1:
!        if sigma >= rho then such [ai, bi] brackets contain an interval of
!        acceptable alpha points such that
!        (i) f(alpha) <= f(0) + alpha*rho*f'(0)  PMO, 2.5.1.
!        (ii) |f'(alpha)'| <= -sigma*f'(0) PMO, 2.5.6 (two sided condition).

!       _Use_sectioning (while preserving above properties) to find
!        an acceptable alpha
        acceptable = .false.
      endif  !not_need345
      do while(need345.or..not.acceptable)
        if(not_need345) then
          if(debug) then
            write(0,*) 'ai, fai, fprimeai = ', ai, fai, fprimeai
            write(0,*) 'bi, fbi = ', bi, fbi
            write(0,*) 'fai.eq.fmin', fai.eq.fmin
            write(0,*) 'fai.le.f0+ai*rho*fprime0',
     &        fai.le.f0+ai*rho*fprime0
            write(0,*) '(bi-ai)*fprimeai.lt.0.d0.and.'//
     &        'abs(fprimeai).gt.-sigma*fprime0 =',
     &        (bi-ai)*fprimeai.lt.0.d0.and.
     &        abs(fprimeai).gt.-sigma*fprime0
            write(0,*) 'fbi.gt.f0+bi*rho*fprime0.or.fbi.ge.fai = ',
     &        fbi.gt.f0+bi*rho*fprime0.or.fbi.ge.fai
          endif
          if(.not.(
     &      (fai.eq.fmin).and.
     &      (fai.le.f0+ai*rho*fprime0).and.
     &      ((bi-ai)*fprimeai.lt.0.d0).and.
     &      (abs(fprimeai).gt.-sigma*fprime0).and.
     &      (fbi.gt.f0+bi*rho*fprime0.or.fbi.ge.fai)))
     &      stop 'bfgs_linesearch: internal bracketing logic error'
          dalpha = bi - ai
!          interpolation to find acceptable point in the
!          range [ai + tau2*dalpha, bi - tau3*dalpha].
          if(ifcubic.and.fprimebi_calc) then
            alphai = bfgs_cubic_minimum(
     &        ai, fai, fprimeai, bi, fbi, fprimebi,
     &        ai + tau2*dalpha, bi - tau3*dalpha)
          else
            alphai = bfgs_quadratic_minimum(
     &        ai, fai, fprimeai, bi, fbi,
     &        ai + tau2*dalpha, bi - tau3*dalpha)
          endif
        endif  !not_need345
        if(not_need45) then
!          terminate using epsilon condition on f (not relative f)
!          to avoid round-off error problems (see discussion on
!          p. 38 of PMOO).  Because assertion (bi-ai)*fprimeai.lt.0.d0
!          above is always correct, the left-hand side of the following
!          condition must always be positive.
          if(need3.or.(ai-alphai)*fprimeai.le.epsilon) then
            if(not_need3) then
!              recalculate everything for alphai = ai since that is less
!              expensive than carrying copies along for the entire line
!              search (and for _every_ line search until function convergence
!              _might_ finally trigger the above if statement.
              alphai = ai
!              Compute trial point at ai corresponding to
!              xi = x + alphai * s
              call bfgs_take_step(n, x, p, alphai, dir/pnorm, xi, dx)
!              EVALUATE f(ai) = function at xi.  This should be identical to
!              fai, but we recalculate here to set up the gradient calculation
!              to satisfy the condition (quite useful for calling routines)
!              that the gradient calculation always occurs at the same xi as
!              the preceeding function evaluation. Also EVALUATE gradient of
!              function at xi.
              status = 'both'
              name = 'linesearch'
              need3 = .true.
              return
            endif  !not_need3
            if(need3) then
! 3            continue
              need3 = .false.
              need345 = need3.or.need4.or.need5
              not_need12345 = .not.(need12.or.need345)
              not_need345 = .not.need345
              not_need3 = .not.need3
              fi = f
!              from PMOO, 1.2.6
!              df(alphai) = s dot gradient
              fprimei = (-dir/pnorm)*ddot(n, p, 1, gradient, 1)
!              output x = xi
              call dcopy(n, xi, 1, x, 1)
!              output f = f(xi)
              f = fi
!              gradient at x = xi already stored ready for output.
              alpha_ratio = alphai/alpha1
              alpha_previous = alphai
              deltaf_previous = f0-fi
              status = 'continue'
              return
            else
              stop 'bfgs_linesearch: bad logic2'
            endif  !need3
          endif
!          Compute new trial point at alphai corresponding to
!          xi = x + alphai * s
          call bfgs_take_step(n, x, p, alphai, dir/pnorm, xi, dx)
!          EVALUATE f(alphai) = function at xi
          status = 'function'
          name = 'linesearch'
          need4 = .true.
          return
        endif  !not_need45
        if(not_need5) then
! 4        continue
          need4 = .false.
          need345 = need3.or.need4.or.need5
          not_need12345 = .not.(need12.or.need345)
          not_need345 = .not.need345
          not_need45 = .not.(need4.or.need5)
          fi = f
          if(fi.le.f0+alphai*rho*fprime0) fmin = min(fi, fmin)
          if(debug) then
            write(0,*) "alpha, f(alpha) =", alphai, fi
          endif
        endif  !not_need5
        if(not_need5.and.(
     &      (fi.gt.f0 + rho*alphai*fprime0).or.(fi.ge.fai))) then
          aip = ai
!          fai and fprimeai unchanged.
          bip = alphai
          fbi = fi
!          dont know gradient at xi at this point in sectioning
!          iteration.
          fprimebi_calc = .false.
        else
          if(not_need5) then
!            EVALUATE gradient of function at same xi as previous.
            status = 'gradient'
            name = 'linesearch'
            need5 = .true.
            return
          endif
          if(need5) then
! 5          continue
            need5 = .false.
            need345 = need3.or.need4.or.need5
            not_need12345 = .not.(need12.or.need345)
            not_need345 = .not.need345
            not_need45 = .not.(need4.or.need5)
            not_need5 = .not.need5
!            from PMOO, 1.2.6
!            df(alphai) = s dot gradient
            fprimei = (-dir/pnorm)*ddot(n, p, 1, gradient, 1)
            if(debug) then
              write(0,*) "alpha, f'(alpha) =", alphai, fprimei
            endif
            if(abs(fprimei).le.-sigma*fprime0) then
!              acceptable point found.
!              output x = xi
              call dcopy(n, xi, 1, x, 1)
!              output f = f(xi)
              f = fi
!              gradient at xi already stored ready for output.
              alpha_ratio = alphai/alpha1
              alpha_previous = alphai
              deltaf_previous = f0-fi
              status = 'continue'
              return
            endif
            aip = alphai
            if(dalpha*fprimei.ge.0.d0) then
              bip = ai
              fbi = fai
              fprimebi_calc = .true.
              fprimebi = fprimeai
            else
              bip = bi
!              fbi, fprimebi_calc, and fprimebi (if defined) unchanged
            endif
            fai = fi
            fprimeai = fprimei
          else
            stop 'bfgs_linesearch: bad logic 3'
          endif  !need5
        endif
        ai = aip
        bi = bip
      enddo
      end subroutine

      subroutine bfgs_iterate(status, name,
     &  fbar, epsilon, n, x, xi, f, gradient, alpha_ratio, dx)
!      do one iteration of the BFGS minimization algorithm consisting
!      of a line search plus BFGS update.

!      status is a control string that takes the following values:
!        'start'.  Set externally to initialize bfgs iteration.
!        'function'.  Request function value to be supplied externally
!        'gradient'.  Request gradient to be supplied externally.
!           N.B. gradient evaluation is always done at xi of prior
!           function evaluation
!        'both'.  Request both function value and gradient to be supplied
!                 externally.
!        'complete'. line search + bfgs update has completed successfully.
!        'error....' some aspect of line search or bfgs update has not
!                    completed successfully.  The 'error' preamble is then
!                    followed with a specific error message.
                
!        Note, all values of status other than 'start' are set internally by
!        the BFGS routines and passed to the calling routine to control what
!        information is supplied on the next call to bfgs_iterate.

!      name is a control string set to the name of the appropriate bfgs
!      routine (with the bfgs_ prefix removed!) when there is a 'function',
!      'gradient', or 'both' request or when there is an error.

!      input quantities:
!      fbar is a user-specified minimum possible f used as per
!        PMOO, 2.6.1 to control the bracketing step size.  If actual f
!        values are <= fbar, then this routine stops with an error so
!        be realistic in how you specify fbar making it small enough
!        to avoid the error, but large enough to provide some control
!        over the maximum size of the bracketing step.
!      epsilon is the user-specified convergence criterion for f (not
!        relative f) discussed on p. 38 of PMOO.  It is used to terminate
!        the line-search routine to work around the case when round-off
!        errors are causing convergence problems close to the minimum.

!      input and output quantities:
!      x(n) is the starting point of the line search on input and on output
!        is the ending point of the line search.
!      xi(n) is trial values of x for external evaluations.
!        for status(:8).eq.'complete' result, x and xi are identical.
!      N.B. both f and gradient must be precalculated on input and calculated
!        on output.  That is:
!      f is f(xi(n)) on both input and output.
!      gradient(n) is the gradient(xi(n)) on both input and output.

!      output quantity:
!      alpha_ratio is the ratio of the calculated step size to the initial
!      estimate, alpha1 for the step size.
!      dx(n) is the vector of differences between the initial x and final x.

!      internal quantity:
!      p(nmax_bfgs) is the unscaled line search direction.
!      The PMOO scaled s vector is defined by
!        s = -(dir/pnorm)*p
!        where dir = 1 if p has an acute angle with the input gradient and
!        dir = -1 otherwise.  pnorm is the norm of p.  dir and pnorm are
!        calculated internally.

      implicit none
      include 'bfgs.h'
      integer*4 n
      real*8 fbar, epsilon,
     &  x(n), xi(n), f, gradient(n), alpha_ratio, dx(n)
      real*8
     &  x0(nmax_bfgs), dx0(nmax_bfgs),
     &  g0(nmax_bfgs), dg0(nmax_bfgs), p(nmax_bfgs),
     &  ddot, dnrm2, dxg, dgg, dxdg, dgnorm, b, a
      character*(*) status, name
      logical initialized, recall
      data initialized/.false./
      save
      
!      sanity check
      if(n.gt.nmax_bfgs) stop 'bfgs_iterate: n too large'
      if(len(status).lt.8)
     &  stop 'bfgs_iterate: length of status string too small'
      if(len(name).lt.10)
     &  stop 'bfgs_iterate:length of name string too small'
      if(status(:5).eq.'error')
     &  stop 'bfgs_iterate: attempt to re-enter with error condition'

      recall =
     &  status(:4).eq.'both'.or.
     &  status(:8).eq.'function'.or.
     &  status(:8).eq.'gradient'

      if(status(:5).eq.'start'.and..not.recall) then
        call bfgs_set(status, name, n, x, xi, f, gradient, p)
        if(status(:8).ne.'continue') then
          return
        else
          stop 'bfgs_iterate: bad logic1'
        endif
      elseif(recall.and.name(:3).eq.'set') then
        call bfgs_set(status, name, n, x, xi, f, gradient, p)
        if(status(:8).eq.'continue') then
          initialized = .true.
          status = "complete"
          return
        else
          stop 'bfgs_iterate: bad logic2'
        endif
      endif

!      further sanity check:
      if(.not.initialized)
     &  stop 'bfgs_iterate: not properly initialized'

      if(.not.recall) then
!        save old values.
        call dcopy(n, x, 1, x0, 1)
        call dcopy(n, gradient, 1, g0, 1)
        call bfgs_linesearch(status, name,
     &    fbar, epsilon, n, p, x, xi, f, gradient, alpha_ratio, dx)
        if(status(:8).ne.'continue') then
          return
        else
          stop 'bfgs_iterate: bad logic3'
        endif
      elseif(name(:10).eq.'linesearch') then
        call bfgs_linesearch(status, name,
     &    fbar, epsilon, n, p, x, xi, f, gradient, alpha_ratio, dx)
        if(status(:8).ne.'continue') return
      else
        stop 'bfgs_iterate: bad logic4'
      endif
      
!      This is the BFGS update:
!      dp = g1 - A dx - B dg
!      A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg
!      B = dx.g/dx.dg

!      dx0 = x - x0
      call dcopy(n, x, 1, dx0, 1)
      call daxpy(n, -1.d0, x0, 1, dx0, 1)

!        dg0 = gradient - g0
      call dcopy(n, gradient, 1, dg0, 1)
      call daxpy(n, -1.d0, g0, 1, dg0, 1)

      dxg = ddot(n, dx0, 1, gradient, 1)
      dgg = ddot(n, dg0, 1, gradient, 1)
      dxdg = ddot(n, dx0, 1, dg0, 1)

      dgnorm = dnrm2(n, dg0, 1)

!      Guard against zero step with either dx0 or dg0 or both equal to the
!      zero vector.
      if(dxdg.eq.0.d0) return

!      B = dx.g/dx.dg
      b = dxg/dxdg

!      A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg
      a = -(1.d0 + dgnorm*dgnorm/dxdg)*b + dgg/dxdg

!      dp = g1 - A dx - B dg
      call dcopy(n, gradient, 1, p, 1)
      call daxpy(n, -a, dx0, 1, p, 1)
      call daxpy(n, -b, dg0, 1, p, 1)
      status = 'complete'
      end subroutine

!==============================================================================

     SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
!     .. Scalar Arguments ..
      DOUBLE PRECISION DA
      INTEGER INCX,INCY,N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION DX(*),DY(*)
!     ..
!
!  Purpose
!  =======
!
!     constant times a vector plus a vector.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!
!     .. Local Scalars ..
      INTEGER I,IX,IY,M,MP1
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MOD
!     ..
      IF (N.LE.0) RETURN
      IF (DA.EQ.0.0d0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
          DY(IY) = DY(IY) + DA*DX(IX)
          IX = IX + INCX
          IY = IY + INCY
   10 CONTINUE
      RETURN
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 M = MOD(N,4)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          DY(I) = DY(I) + DA*DX(I)
   30 CONTINUE
      IF (N.LT.4) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
          DY(I) = DY(I) + DA*DX(I)
          DY(I+1) = DY(I+1) + DA*DX(I+1)
          DY(I+2) = DY(I+2) + DA*DX(I+2)
          DY(I+3) = DY(I+3) + DA*DX(I+3)
   50 CONTINUE
      RETURN
      END SUBROUTINE
!=====================================================

      SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
!     .. Scalar Arguments ..
      INTEGER INCX,INCY,N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION DX(*),DY(*)
!     ..
!
!  Purpose
!  =======
!
!     copies a vector, x, to a vector, y.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!
!     .. Local Scalars ..
      INTEGER I,IX,IY,M,MP1
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MOD
!     ..
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
          DY(IY) = DX(IX)
          IX = IX + INCX
          IY = IY + INCY
   10 CONTINUE
      RETURN
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 M = MOD(N,7)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
          DY(I) = DX(I)
   30 CONTINUE
      IF (N.LT.7) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
          DY(I) = DX(I)
          DY(I+1) = DX(I+1)
          DY(I+2) = DX(I+2)
          DY(I+3) = DX(I+3)
          DY(I+4) = DX(I+4)
          DY(I+5) = DX(I+5)
          DY(I+6) = DX(I+6)
   50 CONTINUE
      RETURN
      END SUBROUTINE
!==================================================

   end module
