# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!ds
!
!Usage:
!call ds (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Sat Apr 17 20:02:33 PDT 1999
!
!Purpose:
!
!
module ds
!based on HCL ds
  implicit none
  real,private,save :: minstep,maxstep,fncdecreasetol,slopedecreasetol
  real,private,save :: bracketincrease
  integer,private,save :: nmod,ndat,maxsample
  contains
  subroutine ds_init(min_step,max_step,max_sample,fnc_decrease_tol&
    &,slope_decrease_tol,bracket_increase)
    integer :: n_mod,n_dat
    integer, optional :: max_sample
    real, optional :: min_step,max_step,fnc_decrease_tol,slope_decrease_tol&
      &,bracket_increase
    if (present(min_step)) then
      minstep=min_step
    else
      minstep=1e-20
    end if
    if (present(max_step)) then
      maxstep=max_step
    else
      maxstep=1e20
    end if
    if (present(max_sample)) then
      maxsample=max_sample
    else
      maxsample=8
    end if
    if (present(fnc_decrease_tol)) then
      fncdecreasetol=fnc_decrease_tol
    else
      fncdecreasetol=.0001
    end if
    if (present(slope_decrease_tol)) then
      slopedecreasetol=slope_decrease_tol
    else
      slopedecreasetol=2.
    end if
    if (present(bracket_increase)) then
      bracketincrease=bracket_increase
    else
      bracketincrease=2.
    end if
  end subroutine
  subroutine ds_size(n_mod,n_dat)
    integer n_mod,n_dat
    nmod=n_mod
    ndat=n_dat
  end subroutine
  function ds_search_simple(oper,func,xcur,xcycle,pdir,grad,eval,reg&
    &,nreg)result(status)
    optional :: reg,nreg
    real eval
    interface
      integer function reg(adj,add,model,data) result(stat)
        integer :: stat
        logical :: adj ,add
        real :: model(:),data(:)
      end function
      integer function oper(adj,add,model,data) result(stat)
        integer :: stat
        logical :: adj ,add
        real :: model(:),data(:)
      end function
      real function func(oper,xval,grad,reg,nreg) result(val)
        optional :: reg,nreg,grad
        real :: val
        real :: grad(:),xval(:)
        interface
          integer function oper(adj,add,model,data) result(stat)
            integer :: stat
            logical :: adj ,add
            real :: model(:),data(:)
          end function
          integer function reg(adj,add,model,data) result(stat)
            integer :: stat
            logical :: adj ,add
            real :: model(:),data(:)
          end function
        end interface
      end function
    end interface
    real, intent(in) :: xcur(:),pdir(:)
    real, intent(out) :: xcycle(:),grad(:)
    real, intent(in) :: pdir
    real :: newtlen
    real :: tempr
    real :: initslope,minmu,maxmu,mu,mu0
    integer :: First,nreg,status
    logical :: exit
    real :: muprev,mudiff,muinc,mutemp,mulo
    real :: fhi,flo,fprev
    real :: newslope
    real :: t1,t2,t3,a,b,disc,f_x,f_nex,v1,v2,v3,w1,w2,w3
    integer :: isample
    newtlen=sqrt(sum(pdir**2))
!we don't believe in functionals not in the space
    if (newtlen>maxstep) then
      write(0,*) "newlen=",newtlen," max step=",maxstep
      call erexit("don't undestand newten length greater thanmax step"&
        &)
    end if
    if (present(nreg)) then
      f_x=func(oper,xcur,grad,reg,nreg)
    else
      f_x=func(oper,xcur,grad)
    end if
    initslope=dot_product(grad,pdir)
    if ( initslope>0.) then
! write(0,*) "function evaluation",f_x
! write(0,*) "x",xcur
! write(0,*) "grad pdir initslope",grad,pdir,initslope
      status=1
      return
    end if
    if (newtlen.eq.0) then
      status=3
      return
    end if
    minmu=minstep/newtlen
    maxmu=maxstep/newtlen
    mu=min(1.,maxmu)
    mu=max(mu,minmu)
    mu0=mu
    isample=0
    do while (isample< maxsample)
      isample= isample+1
      xcycle=xcur+mu*pdir
      if (present(nreg)) then
        f_nex=func(oper,xcycle,grad,reg,nreg)
      else
        f_nex=func(oper,xcycle,grad)
      end if
      write(0,*) "xnext 1",mu,f_nex
! write(0,*) "THE SECOND",xnext,grad
! write(0,*) "here is the test",f_nex,f_x+fncdecreasetol * mu * initslope,f_x,fncdecreasetol,mu,initslope
      if (f_nex <= f_x + fncdecreasetol * mu * initslope) then
!xnext satisfactory
! write(0,*) "xnext satisfactory"
        newslope=dot_product(pdir,grad)
        if (newslope<SlopeDecreaseTol*initslope ) then
! write(0,*) "insufficient decrease in slope",newslope
! write(0,*) "before conditional",mu,mu0,mu,maxmu
          if ( mu .eq. mu0 .and. mu < maxmu ) then
! This is the first step. Keep increasing the step
! length until there is simultaneously sufficient
! decrease in the function and sufficient increase
! in the gradient, or until the maximum step length
! is reached or until there is no longer sufficient
! decrease in the function.
! write(0,*) "here we go increase"
            do while ( (f_nex <= f_x+FncDecreaseTol*mu*initslope .and.&
              & newslope < SlopeDecreaseTol*initslope .and. mu <&
              & 0.99*maxmu .and. isample < MaxSample ) .or. f_nex>huge(f_nex))
              muprev = mu
              fprev = f_x
              if ( BracketIncrease*mu >= maxmu ) then
                mu = 0.99*maxmu
              else
                mu=mu*bracketincrease
              end if
!AAAAA
              xcycle=xcur+mu*pdir
              write(0,*) "xnext",mu
              isample= isample+1
              if (present(nreg)) then
                f_nex=func(oper,xcycle,grad,reg,nreg)
              else
                f_nex=func(oper,xcycle,grad)
              end if
! write(0,*) "xnext 1",mu,f_nex
! write(0,*) "reevaluate 1", f_nex
              if (f_nex <= f_x + FncDecreaseTol*mu*initslope) then
                newslope=dot_product(pdir,grad)
              end if
            end do
! write(0,*) "done with the increase"
          end if
! write(0,*) "now quadatic"
          if ( mu < mu0 .or. (mu > mu0 .and.f_nex > f_x+FncDecreaseTol&
            &*mu*initslope)) then
! There is an acceptable step length
! between mulo and the muhi; find it by using successive
! quadratic interpolation.
            mulo = min( mu,muprev )
            mudiff = abs( mu-muprev )
            if ( mu < mulo ) then
              flo = f_nex
              fhi = fprev
            else
              flo = fprev
              fhi = f_nex
            end if
! write(0,*) "acceptable mu stuff"
            do while ( newslope < SlopeDecreaseTol*initslope .and. mudiff&
              & >= minmu .and. isample < MaxSample )
! write(0,*) "in cubic while"
              if (2.0*(fhi-(flo+newslope*mudiff)).eq.0.) then
                muinc=.2*mudiff
              else
                muinc=(-newslope*mudiff*mudiff)/(2.0*(fhi-(flo+newslope&
                  &*mudiff)))
              end if
              if (muinc < 0.2*mudiff) then
                muinc = 0.2*mudiff
              end if
              mu = mulo + muinc
              xcycle=xcur+mu*pdir
! write(0,*) "xnext",mu
! write(0,*) "xnext 3",isample,mu,xnext
              if (present(nreg)) then
                f_nex=func(oper,xcycle,grad,reg,nreg)
              else
                f_nex=func(oper,xcycle,grad)
              end if
! write(0,*) "xnext 1",mu,f_nex
! write(0,*) "reevaluate 2", f_nex
              isample= isample+1
              if (f_nex > f_x+FncDecreaseTol*mu*initslope) then
                mudiff = muinc
                fhi = f_nex
              else
!AAAAA
                newslope=dot_product(pdir,grad)
                if (newslope < SlopeDecreaseTol*initslope) then
                  mulo = mu
                  mudiff = mudiff - muinc
                  flo = f_nex
                end if
              end if
            end do
          end if
          if ( newslope < SlopeDecreaseTol*initslope ) then
! write(0,*)" ",isample," << ", mu , " << mu << endl"
            write(0,*) "slope1"
            status=2
            return
! call erexit("Failed to satisfy SlopeDecreaseTol condition")
          else
            exit
          end if
        else
! write(0,*) "sufficient decrease"
          exit
        end if
      else
! write(0,*) "xnext not satisfactory"
!failed to decrease function enough
        if ( mu <= minmu ) then
          if ( f_nex < f_x ) then
! write(0,*)" << ",isample," << ", mu , " << mu << endl";
! call erexit("Failed to make sufficient progress (mu too small)")
            status=2
            write(0,*) "slope2"
            return
          else
! write(0,*)" << ",isample," << ", mu , " << mu << endl";
            status=2
            return
! call erexit("Failed to make sufficient progress (mu too small)")
          end if
        else
! write(0,*) "back it not satisfactory"
!Otherwise, perform cubic backtrack (quadratic on
! first step).
          if ( mu .eq. mu0 ) then
            if ( 2.0*(f_nex-f_x-mu*initslope).eq.0.) then
              mutemp = 0.5*mu
            else
              mutemp= -initslope/(2.0*(f_nex-f_x-mu*initslope))
            end if
          else
            exit=.false.
            t1 = f_nex-f_x-mu*initslope
            t2 = fprev-f_x-muprev*initslope
            if (mu-muprev.eq.0.) then
              exit=.true.
            else
              t3=1./(mu-muprev)
            end if
            if (.not. exit .and. mu.eq.0) then
              exit=.true.
            else
              v1=t1/mu**2
            end if
            if (.not. exit .and. muprev.eq.0.) then
              exit=.true.
            else
              v2=t2/muprev**2
            end if
            if (.not. exit) then
              a = t3*(v1-v2)
              b = t3*(mu*t2/(muprev*muprev)-muprev*t1/(mu*mu))
              disc = b*b-3.0*a*initslope
              if (.not. exit .and. b.eq.0) then
                exit=.true.
              else
                w1=-initslope/(2.*b)
              end if
              if (.not. exit .and. a.eq.0) then
                exit=.true.
              else
                w2=(sqrt(disc)-b)/(3.0*a)
              end if
            end if
            if (exit) then
              mutemp = 0.5*mu
            else if (a.eq.0.) then
              mutemp=w1
            else
              mutemp=w2
            end if
          end if
          mutemp = max( minmu,min(0.5*mu,mutemp ) )
        end if
        muprev = mu
        fprev = f_nex
        mu = max( 0.1*mu,mutemp )
      end if
    end do
    if (f_nex>f_x) then
! write(0,*) "this is what I am seeing",f_x,f_nex
! write(0,*) "slope4"
      status=2
      return
    end if
! call erexit("didn't reduce function,bad ")
    status=0
    eval=f_nex
! write(0,*) "sss",eval,mu
  end function
end module
