module  lbfgs_solver_mod
use cp_lbfgs_db
implicit none


contains

!subroutine solver_lbfgs(oper,f_g_calc, x, dat, niter, lb,ub,nmetric,verb,factr,pgtol)
!  optional  :: lb,ub,nmetric,verb,factr,pgtol
!  interface
!    function f_g_calc (x, dat,  fnc, grad ) result (stat)
!         integer              :: stat
!      real, dimension (:)  :: x, dat,fnc,grad
!    end function f_g_calc
!  end interface
!  interface
!    function oper (adj, add, x, dat) result (stat)
!         integer              :: stat
!      logical, intent (in) :: adj, add
!      real, dimension (:)  :: x, dat
!    end function oper
!  end interface
!  real      :: x(:),dat(:),lb(:),ub(:)
!  integer   :: niter,nmetric,nmet,iprint
!  real(kind=dbl),pointer :: model(:),l_b(:),u_b(:),grad(:),wa(:)
!  real          ,pointer :: g(:),rr(:)
!  real(kind=dbl)         :: f,factr,ftr,ptol,pgtol,dsave(29),fct
!  integer, pointer       :: nbd(:),iwa(:)
!  logical                :: verb,v
!  character(len=60)      :: task, csave
!  logical                :: lsave(4),cont,good
!  integer                :: isave(44),iter,fnc_calc
!
!
!  allocate(model(size(x)))
!
!  !number of metrics
!  nmet=5
!  if(present(nmetric))nmet=nmetric
!  
!  !lower bounds
!  allocate(l_b(size(model)),u_b(size(model)),nbd(size(model)))
!  allocate(grad(size(model)),g(size(model)))
!  allocate(rr(size(dat)))
!  nbd=0
!  if(present(lb))then
!    l_b=lb
!    if(present(ub))then
!      u_b=ub
!      nbd=2
!    else
!      nbd=1
!    end if
!  else if(present(ub)) then
!    nbd=3
!    u_b=ub
!  end if
!  ftr=1e1
!  if(present(factr)) ftr=factr 
!  ptol=1e-6
!  if(present(pgtol)) ptol=pgtol 
!  
!  allocate(wa(2*nmet*size(model)+4*size(model)+12*nmet*nmet+12*nmet))
!  allocate(iwa(3*size(model)))
!
!  iprint=-1
!  v=.false.
!  if(present(verb)) then
!    if(verb) iprint=0
!    v=verb
!  end if
!
!
!
!
!  !calculate initial residual
!  rr=-data   
!  ierr=oper(.false.,.false.,x,rr) 
!
!
!  fct=0.5*sum(dprod(rr,rr)) !initial functional
!
!  iter=1
!  task="START"
!  good=.true.;fnc_calc=0
!  do while(iter < niter .and. good) 
!    good=.false.
!
!    if(task(1:2) == "FG" .or. task(1:4)=="START") then
!      fnc_calc=0
!      ierr=oper(.true. ,.false.,g,rr) !calculate gradient
!
!      grad=g  !double precission the gradient
!      call setulb(n123,mmem,model,l_b,u_b,nbd,fct,grad,ftr,ptol,wa,iwa,task,&
!       iprint,csave,lsave,isave,dsave)
!      x=model !single precision the model
!
!      !calculate the uppdated residual
!      rr=-data   
!      ierr=oper(.false.,.false.,x,rr) 
!
!      fct=0.5*sum(dprod(rr,rr)) !update functional
!
!      good=.true.
!      fnc_calc=fnc_calc+1
!    end if
!
!    if (task(1:5)=="NEW_X") then
!      good=.false.
!      fnc_calc=0
!      if(v) write(0,*) "Iter ",iter,fnc_calc,"=fnc_calc"
!      iter=iter+1
!    end if
!  end do
!  deallocate(model,l_b,u_b,nbd,rr,g,grad,wa,iwa)
!end subroutine

subroutine solver_lbfgs(f_g_calc, x, niter, l_b,u_b,nbd,nmetric,iprint,factr,pgtol)
  interface
    function f_g_calc (x,   fnc, grad ) result (stat)
       integer, parameter :: dbl=SELECTED_REAL_KIND(15,307)
         integer              :: stat
      real, dimension (:)  :: x, grad
      real(kind=dbl) :: fnc
    end function f_g_calc
  end interface
  real      :: x(:)
  integer   :: niter,nmetric,nmet,iprint,ierr,n123
  real(kind=dbl),pointer :: model(:),grad(:),wa(:)
  real(kind=dbl) :: l_b(:),u_b(:)
  real          ,pointer :: g(:),rr(:)
  real(kind=dbl)         :: f,factr,pgtol,dsave(29),fct
  integer, pointer       :: iwa(:)
  integer  :: nbd(:)
  character(len=60)      :: task, csave
  logical                :: lsave(4),cont,good
  integer                :: isave(44),iter,fnc_calc


  allocate(model(size(x)))

  !number of metrics
  
  !lower bounds
  allocate(grad(size(model)),g(size(model)))
  
  allocate(wa(2*nmetric*size(model)+4*size(model)+12*nmetric*nmetric+12*nmetric))
  allocate(iwa(3*size(model)))
 


  !calculate initial residual
  ierr=f_g_calc(x,fct,g)
  model=x



  iter=1
  task="START"
  good=.true.;fnc_calc=0
  do while(iter <= niter .and. good) 
    good=.false.

!  write(0,*) 'iter=',iter,'fct',fct
    grad=g  !double precission the gradient

    call setulb(size(model),nmetric,model,l_b,u_b,nbd,fct,grad,factr,pgtol,wa,iwa,task,&
     iprint,csave,lsave,isave,dsave)
    x=model !single precision the model

    if(task(1:2) == "FG" ) then
      ierr=f_g_calc(x,fct,g)
      good=.true.
      fnc_calc=fnc_calc+1
    else if (task(1:5)=="NEW_X") then
      good=.true.
      if(iprint >=0) write(0,*) "Iter ",iter,fnc_calc,"=fnc_calc",fct
      fnc_calc=0
      iter=iter+1
    else 
      write(0,*) "EXIT PROBLEM ",task
    end if
  end do
  deallocate(model,g,grad,wa,iwa)
end subroutine
 
end module
