program main
  use sep
  use solver_lbfgsb_mod
  use wemva_struct_mod
  use wemva_kerpix_incore_mod
  use wemva_offdso_incore_mod
  use wemva_objdso_incore_mod
  use triangle2

  implicit none

  type(img_grid)                              :: igrid
  type(vel_grid)                              :: vgrid
  type(wemva_ftag)                            :: ftag_
  type(axis)                                  :: axk, ayk, azk, aiter
  real            , dimension(:), allocatable :: bslw, dimg, dslw, dspl, temp
  integer                                     :: rect1, rect2, rect3

  real                                        :: bvel_aver
  double precision                            :: ax, bx, cx
  double precision                            :: fa, fb, fc
  double precision                            :: tol, xmin, fmin
  integer                                     :: niter, iter, ITMAX, stat
  character(len=6)                            :: lnsrch

  character(len=256)                          :: mmov_

! L-BFGS solver
  character*60                                :: task, csave
  logical                                     :: lsave(4), lbfgsb
  integer                                     :: n, m, iprint, isave(44)
  integer         , dimension(:), allocatable :: nbd, iwa
  double precision                            :: f, factr, pgtol, dsave(29)
  double precision, dimension(:), allocatable :: x, l, u, g, wa

  mmov_ = "mmov"

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()
  
  call from_param("lnsrch", lnsrch, "brent")
  call from_param("ITMAX" , ITMAX , 10     )

  call from_param("niter", niter, 5) 
  call from_param("rect1", rect1, 5)
  call from_param("rect2", rect2, 5)
  call from_param("rect3", rect3, 5)

  call from_param("oxk", axk%o, 0.)
  call from_param("dxk", axk%d, 1.)
  call from_param("nxk", axk%n, 1 )

  call from_param("oyk", ayk%o, 0.)
  call from_param("dyk", ayk%d, 1.)
  call from_param("nyk", ayk%n, 1 )

  call from_param("ozk", azk%o, 0.)
  call from_param("dzk", azk%d, 1.)
  call from_param("nzk", azk%n, 1 )


  stat = wemva_kerpix_incore_init(vgrid, igrid, ftag_)
  stat = wemva_offdso_incore_init(igrid%ax%n ,igrid%ay%n ,igrid%az%n ,igrid%ahx%n,igrid%ahy%n,igrid%ahz%n, &
                                  igrid%ahx%o,igrid%ahy%o,igrid%ahz%o,igrid%ahx%d,igrid%ahy%d,igrid%ahz%d)
  stat = base_spln3d_init(vgrid%ax,vgrid%ay,vgrid%az,axk,ayk,azk)
  call triangle2_init(rect1,rect3,vgrid%ax%n,vgrid%az%n)
 
  allocate(dimg(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))
  allocate(temp(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))
  allocate(bslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(dslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(dspl(axk%n*ayk%n*azk%n))

  n = vgrid%ax%n*vgrid%ay%n*vgrid%az%n
  m = 5

  allocate(x(n), l(n), u(n), g(n), nbd(n))
  allocate(wa(2*m*n+4*n+12*m*m+12*m), iwa(3*n))

  dimg = 0.
  temp = 0.
  bslw = 0. 
  dslw = 0.
  dspl = 0.

  call sseek(ftag_%bslw_, 0, 0)
  call sep_read(bslw, ftag_%bslw_)

  call auxinout(ftag_%dimg_)
  call wemva_histo_imag(igrid, ftag_%dimg_)
  call wemva_fio_imag_init(igrid, ftag_%dimg_)

  aiter%o = 1
  aiter%d = 1
  aiter%n = niter

  call auxinout(mmov_)

  call houaxis(mmov_, vgrid%ax, 1)
  call houaxis(mmov_, vgrid%ay, 2)
  call houaxis(mmov_, vgrid%az, 3)
  call houaxis(mmov_, aiter   , 4)

  iprint = 1

  factr  = 1.0d+7
  pgtol  = 1.0d-5

  task   = 'START'
  lbfgsb = .true.
  iter   = 0

  l      = dble(1./5500.)
  nbd    = 1

  x      = dble(bslw)
  f      = dble(0.)
  g      = dble(0.)

  do while (lbfgsb .and. iter<=niter)
     lbfgsb = .false.
     call setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, &
                  task, iprint, csave, lsave, isave, dsave)

      

     if (task(1:2) .eq. 'FG') then

         bslw = real(x)
         call sseek(ftag_%bslw_, 0, 0)
         call sep_write(bslw, ftag_%bslw_)

         stat = wemva_objdso_incore_oper(.false.,.false.,bslw,dimg)

         stat = wemva_offdso_incore_oper(.true. ,.false.,temp,dimg)

         call sseek(ftag_%dimg_, 0, 0)
         call sep_write(temp, ftag_%dimg_)

         stat = wemva_kerpix_incore_oper(.true. ,.false.,dslw,temp)

         stat = base_spln3d_oper(.true. , .false., dspl, dslw)
         stat = base_spln3d_oper(.false., .false., dspl, dslw)

         dslw = dslw/maxval(abs(dslw))
         
         call sseek(ftag_%dslw_, 0, 0)
         call sep_write(dslw, ftag_%dslw_)

         f = dble(sum(dimg*dimg))

         g = dslw

         lbfgsb = .true.

      endif

      if (task(1:5) .eq. 'NEW_X')  then
         iter = iter + 1
         lbfgsb = .true.
         bslw = real(x)
         call sseek(ftag_%bslw_, 0, 0)
         call sep_write(bslw, ftag_%bslw_)
         call sep_write(bslw, mmov_)
      end if
  end do

  
end program
