module tomo_inverse_mod
use sep
use sep3d_struct_mod
use tomo_chain_mod
use bob_solver_mod
use cgstep_mod
use lbfgs_solver_mod
implicit none
logical,private, save :: prec,verb,dot,lbfgs
integer, private, save :: niter
real, private, save , pointer :: resid(:)
real(kind=dbl), pointer,save :: l_b(:),u_b(:)

contains

logical function init_solver(slow)
real :: slow(:)

init_solver=.false.

call from_param("niter",niter,50)
call from_param("prec",prec,.true.)
call from_param("dot",dot,.false.)
call from_param("verb",verb,.false.)
lbfgs=lbfgs_val()
if(lbfgs) then
 if(.not. lbfgs_init(slow)) then
    write(0,*) "trouble initializing lbfgs tomo params"
    return
 end if

end if


init_solver=.true.
end function

logical function lbfgs_init(slow) 
  integer, external :: auxin,sreed
  real, allocatable :: buf1(:),buf2(:)
  real :: slow(:)
  real    :: tempr
  integer :: n1,n
  
  lbfgs_init=.false.

  n=size(slow)
  allocate(u_b(n),l_b(n))
  allocate(buf2(n))
  if(-1/=auxin("vel_min")) then
    if(4*n /= sreed("vel_min",buf2,size(buf2)*4)) then
      write(0,*) "trouble reading vel_min"
      return
    end if
  else
    call from_param("vmin",tempr,1.5)
    buf2=tempr
  end if
  buf2=1./buf2
  u_b=buf2-slow
  if(-1/=auxin("vel_max")) then
    if(4*n /= sreed("vel_max",buf2,size(buf2)*4)) then
      write(0,*) "trouble reading vel_max"
      return
    end if
  else
    call from_param("vmax",tempr,7.)
    buf2=tempr
  end if
  buf2=1./buf2
  l_b=buf2-slow
  deallocate(buf2)
  lbfgs_init=.true.
end function


logical function run_solver(data,model)
real :: model(:),data(:)
real,allocatable :: dm(:),junk(:),junk2(:)
integer :: ierr

if(lbfgs) then
 run_solver=lbfgs_solver(data,model)
else 
 run_solver=cg_solver(data,model)
end if

end function


integer function tomo_calc(x,fnc,grad)
  real :: x(:),grad(:)
  real(kind=dbl) :: fnc
  real, allocatable :: rr(:),tm(:),t_d(:)
  integer :: ierr
                                                                                          
  allocate(rr(size(resid)),tm(size(grad)))
  rr=resid
  tm=0
!  ierr=prec_op(.true., .false., tm, x)
!  ierr=prec_op(.false., .false., tm, x)
  ierr=op(.false.,.true.,x  ,rr)
  fnc=0.5*(sum(dprod(rr,rr)))
  ierr=op(.true.,.false.,grad  ,rr)
!  ierr=prec_op(.true., .false., tm, grad)
!  ierr=prec_op(.false., .false., tm, grad)
                                                                                          
  deallocate(rr,tm)
  tomo_calc=0
end function




logical function lbfgs_solver(data,model)
real,allocatable :: dm(:),junk(:),junk2(:)
real,allocatable :: rr(:),r_m(:),t_d(:),s0(:)
  real :: model(:),data(:)
  real :: e
  real(kind=dbl) :: factr,pgtol
  integer,  allocatable :: nbd(:)
  integer :: iprint,nmetric,ierr
  integer, external :: auxin
                                                                                          
  lbfgs_solver=.false.
if(dot) then
  if(.not. chain_dot() ) call seperr("didn't pass dot test")
end if 

  iprint=0
  if(verb) iprint=1
  nmetric=5
  factr=1.0e-12
  pgtol=.0001

  allocate(resid(size(data)))
  resid=-data
  allocate(nbd(size(model)))

  nbd=2
  model=0.
  call solver_lbfgs(tomo_calc,model,niter,l_b,u_b,nbd,nmetric,iprint,&
    factr,pgtol)

lbfgs_solver=.true.
end function




logical function cg_solver(data,model)
real :: model(:),data(:)
real,allocatable :: dm(:),junk(:),junk2(:),junk3(:)
integer :: ierr
allocate(dm(size(model)))
allocate(junk(size(model)),junk3(size(model)))
!dot=.true.
if(dot) then
  if(.not. chain_dot() ) call seperr("didn't pass dot test")
end if 
    
allocate(junk2(size(data)))
call bob_solver( junk,data, op, cgstep, niter,resd=junk2,verb=verb)


ierr=m_p(.false.,.false.,junk,model)

ierr=op(.true.,.false.,junk3,junk2);
ierr=m_p(.false.,.false.,junk3,junk)
call srite("r_in_m.H",junk,size(junk)*4)



!call srite("aa.H",junk,size(junk)*4)
!call srite("aa.H",model,size(junk)*4)

call srite("res.H",data,size(data)*4)
call srite("res.H",-junk2,size(junk2)*4)
cg_solver=.true.
end function







end module
