        subroutine bfgsfm(fcn,nvar,np,p,h0,ftol,maxiter,
     *                  maxstep,fnew,iter,iprint,inform)

        implicit real*8(a-h,o-z)
        character*1 trans
        parameter(nvmax=100,eps=1.d-16,one=1.0,zero=0.0,trans='N')
        dimension p(np),pnew(nvmax),h0(np,np),g(nvmax),
     *            gnew(nvmax),h(nvmax,nvmax),xi(nvmax),work(nvmax),
     *            dx(nvmax),shy(nvmax),shys(nvmax,nvmax)
        external fcn
        sftol=sqrt(ftol)
        cftol=ftol**0.35d0
        ymin=-0.001
        inform=0
        jset=1

        call fcn(2,nvar,p,fp,g,0)
        do 10 i=1,nvar
            do 5 j=1,nvar
                h(i,j)=h0(i,j)
5           continue
10      continue

c---------------- set xi = -h*g
        call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,zero,xi,1)

        nfail=0
        jhct=0
        do 99 iter=1,maxiter
12          call bhhhstp(fcn,nvar,p,fp,xi,g,maxstep,s,fnew,iters,info)

c---------------- check for failed step; reset hessian & search direction
            if (fnew.gt.fp) then
                if (jset.eq.1) then
                    inform=1
                    fnew=fp
                    goto 999
                endif
                nfail=nfail+1
                jset=1
                do 15 i=1,nvar
                    do 14 j=1,nvar
                        h(i,j)=h0(i,j)
14                  continue
15              continue

                call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,zero,xi,1)
                goto 12
            else

c---------------- set pnew=p+s*xi and compute new gradient
            call dcopy(nvar,p,1,pnew,1)
            call daxpy(nvar,s,xi,1,pnew,1)
            call fcn(1,nvar,pnew,fjunk,gnew,0)

c------------------ print out results every iprint iteration
            itermod=mod(iter,iprint)
            if (itermod.eq.0) then
                imax=idamax(nvar,gnew,1)
                gnorm=ddot(nvar,gnew,1,gnew,1)
                write(*,21)iter,fnew,(fp-fnew)/abs(fp),gnorm,
     *                      nfail,jhct,jbr
21              format(1x,i8,2x,f12.8,2x,e13.6,2x,f12.8,2x,3i4)
                nfail=0
                jhct=0
		jbr=0
            endif

c----------------- Check if convergence criteria satisfied
            imax=idamax(nvar,gnew,1)
            gmax=abs(gnew(imax))
            if(gmax.lt.eps) then
                goto 999
            else
                if ( (fp-fnew).lt.ftol*(1.+fnew)) then
                    if(gmax.lt.cftol*(1.+abs(fnew))) then
                        imax=idamax(nvar,pnew,1)
                        pmax=abs(pnew(imax))
                        call dcopy(nvar,pnew,1,work,1)
                        call daxpy(nvar,-1.d0,p,1,work,1)
                        imax=idamax(nvar,work,1)
                        dpmax=abs(work(imax))
                        if (dpmax.lt.sftol*(1.+pmax)) then
                            goto 999
                        endif
                    endif
                endif
            endif

c--------------------- If not, prepare for the next iteration
c        set dx=pnew-p
            call dcopy(nvar,pnew,1,dx,1)
            call daxpy(nvar,-1.d0,p,1,dx,1)

c        set shy=dx-h*gnew+h*g
            call dcopy(nvar,dx,1,shy,1)
            call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,gnew,1,1.d0,shy,1)
            call dgemv(trans,nvar,nvar,1.d0,h,nvmax,g,1,1.d0,shy,1)

c        set shys=(shy*dx')+(shy*dx')'
            do 30 i=1,nvar
                do 25 j=1,nvar
                    shys(i,j)=shy(i)*dx(j) + shy(j)*dx(i)
25              continue
30          continue

c        set yshy=gnew'*shy-g'*shy
            y1=ddot(nvar,gnew,1,shy,1)
            y2=ddot(nvar,g,1,shy,1)
            yshy=y1-y2

c        set ys=gnew'*dx-g'*dx
            y1=ddot(nvar,gnew,1,dx,1)
            y2=ddot(nvar,g,1,dx,1)
            ys=y1-y2

c       update h if ys>0, otherwise reset h
            if (ys.gt.0.) then
                jset=0
                call dger(nvar,nvar,-yshy/ys,dx,1,dx,1,shys,nvmax)
                do 40 i=1,nvar
                    do 35 j=1,nvar
                        h(i,j)=h(i,j)+shys(i,j)/ys
35                  continue
40              continue
            else
                jset=1
                do 50 i=1,nvar
                    do 45 j=1,nvar
                        h(i,j)=h0(i,j)
45                  continue
50              continue
                jhct=jhct+1
            endif

c       update f,g,p, and search direction
            call dcopy(nvar,gnew,1,g,1)
            call dcopy(nvar,pnew,1,p,1)
            fp=fnew
            call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,0.d0,xi,1)

c------------------- If next step "points up", reset Hessian
            y1=ddot(nvar,xi,1,g,1)
            y2=ddot(nvar,xi,1,xi,1)
            if (y1.ge.ymin*y2) then
                jset=1
                do 60 i=1,nvar
                    do 55 j=1,nvar
                        h(i,j)=h0(i,j)
55                  continue
60              continue
                call dgemv(trans,nvar,nvar,-1.d0,h,nvmax,g,1,0.d0,xi,1)
                jhct=jhct+1
            endif

        endif

99      continue

999     return
        end

        subroutine bhhhstp(fcn,nvar,x0,f0,d,g,itermax,s,fs,
     *                      iters,info)
        implicit real*8(a-h,o-z)
        parameter(delta=0.25,pwr=0.618,one=1.0,nvmax=100)
        dimension x0(nvar),d(nvar),g(nvar),work(nvmax)

c--------------------- Initializations
        info=1
        iter=0
        iup=0
        idown=0
        factor=2.
        dg=ddot(nvar,d,1,g,1)
        downfact=dg*delta
        upfact=dg*(1.-delta)
        alam=1.0
        if (itermax.eq.0) itermax=25

        do 10 iter=1,itermax

c----------------- Set g= alam*d + x0, vof=fcn(g)
            call dcopy(nvar,x0,1,work,1)
            call daxpy(nvar,alam,d,1,work,1)
            call fcn(0,nvar,work,vof,g,1)

c------------------ vof above cone, hence reduce step size
            if ((vof-f0).gt.(downfact*alam)) then
                idown=1
                if (iup.eq.1) then
                    factor=factor**pwr
                    iup=0
                endif
                alam=alam/factor

c------------------ vof below cone, hence increase step size
            else if ((vof-f0).lt.(upfact*alam)) then
                iup=1
                if (idown.eq.1) then
                    factor=factor**pwr
                    idown=0
                endif
                alam=alam*factor

c------------------ vof within cone, hence stop iterations
            else
                info=0
                goto 20
            endif
10      continue

20      continue

c------------------- compute returned values
        iters=iter
        s=alam
        fs=vof
        return
        end

