program main
  use base_struct_mod
  use base_weight_mod
  use base_mskopr_mod
  use hes_convol_mod2
  use hes_reglar_mod
  use solver_smp_mod
  use solver_reg_mod
  use solver_prc_mod
  use solver_irls_mod
  use cgstep_mod
  use cdstep_mod
  use lsqr_mod
  use laplac2
  use irls
  use sep
  use solver_H_mod
  use base_solver_mod
  implicit none

  type(hes_grid)                    :: hgrid
  real, dimension(:)  , allocatable :: modl
  real, dimension(:)  , allocatable :: data
  real, dimension(:)  , allocatable :: mask
  real, dimension(:)  , allocatable :: wght
  real, dimension(:)  , allocatable :: epsm
  real, dimension(:,:), allocatable :: mmov
  real, dimension(:,:), allocatable :: rmov
  real, dimension(:)  , allocatable :: fmov
  logical                           :: adj, search, vareps
  character(len=256)                :: modl_, data_, hess_, mask_, wght_
  character(len=256)                :: mmov_, rmov_, fmov_ 
  integer                           :: stat, ndim, iter, niter, solver, reglar
  real                              :: imgsign, beta_max, eps 
  double precision                  :: t1, t2

  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()

  modl_ = "modl"
  data_ = "data"
  hess_ = "hess"
  mask_ = 'mask'
  wght_ = 'wght'
  mmov_ = 'mmov'
  rmov_ = 'rmov'
  fmov_ = 'fmov'
 
  call from_param("adj"     , adj     , .true.)
  call from_param("imgsign" , imgsign ,    +1.)
  call from_param("niter"   , niter   ,    50 )
  call from_param("search"  , search  , .true.)
  call from_param("vareps"  , vareps  , .true.)
  call from_param("solver"  , solver  ,      1)
  call from_param("reglar"  , reglar  ,      1)
  call from_param("beta_max", beta_max,  0.001) 
  call from_param("eps"     ,      eps,  0.001)
  
  call auxinout(modl_)
  call auxinout(data_)
  call auxin   (hess_)
  call auxin   (mask_)
  call auxin   (wght_)
  call auxout  (mmov_)
  call auxout  (rmov_)
  call auxout  (fmov_)
 
  call hinaxis(hess_, hgrid%ax_modl, 1) 
  call hinaxis(hess_, hgrid%ay_modl, 2)
  call hinaxis(hess_, hgrid%az_modl, 3)
  call hinaxis(hess_, hgrid%ax_offd, 4)
  call hinaxis(hess_, hgrid%ay_offd, 5)
  call hinaxis(hess_, hgrid%az_offd, 6)

  ndim  = hgrid%ax_modl%n*hgrid%ay_modl%n*hgrid%az_modl%n

  allocate(modl(ndim))
  allocate(data(ndim))
  allocate(mask(ndim))
  allocate(wght(ndim))
  allocate(epsm(ndim))
  allocate(mmov(ndim,niter))
  allocate(rmov(ndim,niter))
  allocate(fmov(niter))

  modl = 0.; data = 0.; mask = 0.; wght = 0.
  mmov = 0.; rmov = 0.; fmov = 0.; epsm = 0.

  if (adj) then
     call houaxis(modl_, hgrid%ax_modl, 1)
     call houaxis(modl_, hgrid%ay_modl, 2)
     call houaxis(modl_, hgrid%az_modl, 3)

     call houaxis(mmov_, hgrid%ax_modl, 1)
     call houaxis(mmov_, hgrid%ay_modl, 2)
     call houaxis(mmov_, hgrid%az_modl, 3)
     call auxputch("n4", "i", niter, mmov_)

     call houaxis(rmov_, hgrid%ax_modl, 1)
     call houaxis(rmov_, hgrid%ay_modl, 2)
     call houaxis(rmov_, hgrid%az_modl, 3)
     call auxputch("n4", "i", niter, rmov_)

     call auxputch("n1", "i", niter, fmov_)

     call sep_read(data, data_) 
     call sep_read(wght, wght_)

     data = data * imgsign / maxval(abs(data))

     wght = 1. / wght

     if (vareps) then
        epsm = abs(wght) / maxval(abs(wght)) * eps
     else
        epsm = eps
     end if

  else
     call houaxis(data_, hgrid%ax_modl, 1)
     call houaxis(data_, hgrid%ay_modl, 2)
     call houaxis(data_, hgrid%az_modl, 3)
     call sep_read(modl, modl_)
  end if

  stat = hes_convol_init(hgrid, hess_)

  if (adj) then


     call sep_read(mask, mask_)

     stat = base_weight_init(wght)
     stat = base_mskopr_init(mask)
     stat = hes_reglar_init(hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n,reglar)

     select case (solver)


     !call solver_irls(m=modl, d=data, Fop=hes_convol_oper, Wop=base_mskopr_oper, stepper=cgstep, Wdiag=cauchy, &
     !                 niter=niter, mmov=mmov, rmov=rmov, err=fmov, verb=.true.)

     !call solver_prc(m=modl, d=data, Fop=hes_convol_oper, Sop=base_weight_oper, nSop=ndim, Wop=base_mskopr_oper, &
     !                stepper=cgstep, niter=niter, mmov=mmov, rmov=rmov, eps=0.1, verb=.true.)

     case (1)
        call base_solver_cglpd_smp(m=modl, d=data, Fop=hes_convol_oper, Wop=base_mskopr_oper, &
                                   niter=niter, search=search, beta_max=beta_max, mmov=mmov, rmov=rmov, fmov=fmov)

     case (2)
        call base_solver_cglpd_reg(m=modl, d=data, Fop=hes_convol_oper, Wop=base_mskopr_oper, Aop=hes_reglar_oper, &
                                   niter=niter, eps=epsm, search=search, beta_max=beta_max, mmov=mmov, rmov=rmov, fmov=fmov)

     case (3)
        call base_solver_cglpd_prc(m=modl, d=data, Fop=hes_convol_oper, Sop=base_weight_oper, Wop=base_mskopr_oper, &
                                   niter=niter, search=search, beta_max=beta_max, mmov=mmov, rmov=rmov, fmov=fmov)

     case (4)
        call solver_reg(m=modl, d=data, Fop=hes_convol_oper, Aop=hes_reglar_oper, Wop=base_mskopr_oper, stepper=cgstep, &
                        nAop=ndim, eps=eps, niter=niter, mmov=mmov, rmov=rmov, verb=.true.)

     !call base_solver_anneal(func=base_objfnc, Fop=hes_convol_oper, xmin=modl, dobs=data, Tred=0.618, KTmax=3, KTacc_max=500, Smax=2, &
     !                        Tmin=0.01, Omin=0.01, wt=1., Wop=base_mskopr_oper)
     end select
     stat = base_weight_clos()
     stat = base_mskopr_clos()

  else
     stat = hes_convol_oper(adj,.false.,modl,data)
  end if

  stat = hes_convol_clos()

  if (adj) then
     call sep_write(modl, modl_)
     call sep_write(mmov, mmov_)
     call sep_write(rmov, rmov_)
     call sep_write(fmov, fmov_)
  else
     call sep_write(data, data_)
  end if

  deallocate(modl, data)

end program
