program main
  use sep
  use wemva_struct_mod
  use wemva_kerpix_mod
  use wemva_offdso_mod
  use solver_lbfgsb_mod
  use triangle2

  implicit none

  type(img_grid)                              :: igrid
  type(vel_grid)                              :: vgrid
  type(wemva_ftag)                            :: ftag_
  type(axis)                                  :: axk, ayk, azk
  double precision, dimension(:), allocatable :: x, l, u, g, wa
  real            , dimension(:), allocatable :: bslw, dimg, dslw, dspl
  integer         , dimension(:), allocatable :: nbd, iwa
  double precision                            :: f, factr, pgtol
  integer                                     :: rect1, rect2, rect3
  integer                                     :: nmax, mmax, iprint
  character(len=60)                           :: task, csave
  logical                                     :: lsave(4)
  integer                                     :: isave(44)
  double precision                            :: dsave(29)
  character(len=256)                          :: temp_

  integer :: stat, niter, iter, ix, iy, iz, ihx, ihy, ihz, n, m, nm, nh
  real    :: hcut, hord
  logical :: lbfgs_step

  temp_ = "temp"

  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("niter", niter, 5) 
  call from_param("rect1", rect1, 10)
  call from_param("rect3", rect3, 10)

  call from_param("hcut" , hcut , 0.)
  call from_param("hord" , hord , 1.)

  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_init(vgrid, igrid, ftag_)
  stat = wemva_offdso_init(igrid%ax%n ,igrid%ay%n ,igrid%az%n ,igrid%ahx%n,igrid%ahy%n,&
                           igrid%ahx%o,igrid%ahy%o,igrid%ahx%d,igrid%ahy%d, hcut, hord)
  stat = base_spln3d_init(vgrid%ax,vgrid%ay,vgrid%az,axk,ayk,azk)
  call triangle2_init(rect1,rect3,vgrid%ax%n,vgrid%az%n)
 
  n = vgrid%ax%n*vgrid%ay%n*vgrid%az%n
  m = 7

  iprint = 1
  factr  = 1.0d+1
  pgtol  = 1.0d-6
  
  allocate(dimg(igrid%ahx%n*igrid%ahy%n*igrid%ax%n*igrid%ay%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))

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

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

  x   = dble(bslw) 
  l   = dble(1./3000.)
  u   = dble(1./1000.)
  nbd = 2

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

  call auxinout(temp_)
  call wemva_histo_imag(igrid, temp_)
  call wemva_fio_imag_init(igrid, temp_)

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

  do while (lbfgs_step .and. iter<=niter)

     lbfgs_step = .false.

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

     call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint, &
                 csave,lsave,isave,dsave)

     if (task(1:2)=='FG') then

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

        !!.. compute the current objfnc and gradient

        stat = wemva_bimpix_oper(         .false., ftag_%bslw_, ftag_%bimg_)          !!.. bI
        stat = wemva_offdso_oper(.false., .false., ftag_%bimg_, ftag_%dimg_)          !!.. D bI
        stat = wemva_offdso_oper(.true. , .false., temp_      , ftag_%dimg_)          !!.. D'D bI
        stat = wemva_kerpix_oper(.true. , .false., ftag_%dslw_, temp_      )          !!.. F'D'D bI

        call sseek(ftag_%dslw_, 0, 0)
        call sep_read(dslw, ftag_%dslw_)

        stat = base_spln3d_oper(.true. , .false., dspl, dslw)
        stat = base_spln3d_oper(.false., .false., dspl, dslw)
 
!        stat = triangle2_lop(.true., .false., bslw, dslw) 
!        dslw = bslw

        g =  dslw / maxval(abs(dslw)) 

        call sseek(ftag_%dslw_, 0, 0)
        call sep_write(dslw, ftag_%dslw_)

        f = dble(0.)

        call sseek(ftag_%dimg_, 0, 0)
        do iz=1, igrid%az%n
           call sep_read(dimg, ftag_%dimg_)
           f = f + dot_product(dimg,dimg)
        end do

        write(0,*)"f =", f

        lbfgs_step = .true.

     end if

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

  end do

  
end program
