program main
  use base_lapc2d_mod
  use base_lapc3d_mod
  use posmva_struct_mod
  use posmva_kernel_mod
  use cgstep_mod
  use solver_smp_mod
  use solver_reg_mod
  use solver_irls_mod
  use irls
  use laplac2
  type(posmva_par)                :: par
  integer                         :: stat, nAop
  logical                         :: adj, add, bkg, dottest
  real, dimension(:), allocatable :: dS, dI
  real, dimension(:), allocatable :: dat1, dat2, mod1, mod2

  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("adj"    , adj    , .false.)
  call from_param("dottest", dottest, .false.)

  stat = posmva_params_init(par)

  allocate(dS(par%ax%n*par%ay%n*par%az%n))
  allocate(dI(par%ax%n*par%ay%n*par%az%n))

  stat = posmva_arrays_init()
  stat = posmva_bgslow_init()
  stat = posmva_wfdexp_init()

  write(0,*)"-----------------------------------------------"
  write(0,*)"      computing background wavefield..."
  write(0,*)"-----------------------------------------------"

  stat = posmva_bgwfld_oper()


  if (dottest) then
     write(0,*)"-----------------------------------------------"
     write(0,*)"      dot-product test...                      "
     write(0,*)"-----------------------------------------------"

     allocate(mod1(par%ax%n*par%ay%n*par%az%n), mod2(par%ax%n*par%ay%n*par%az%n))
     allocate(dat1(par%ax%n*par%ay%n*par%az%n), dat2(par%ax%n*par%ay%n*par%az%n))
     call random_number(dat1)
     call random_number(mod2)

     stat =  posmva_tomgph_oper(.true. , .false., mod1, dat1)
     stat =  posmva_tomgph_oper(.false., .false., mod2, dat2)

     write(0,*)"dot(1)=", sum(mod1*mod2)
     write(0,*)"dot(2)=", sum(dat1*dat2)

  else
 
     if (adj) then
        call houaxis(par%dS_, par%arx, 1)
        call houaxis(par%dS_, par%ary, 2)
        call houaxis(par%dS_, par%az , 3)
        !call sep_read (dI(par%jx1:par%jx2,par%jy1:par%jy2,:), par%dI_)
        call sep_read( dI, par%dI_)

        !stat = posmva_tomgph_oper(adj, .false., dS, dI)

        stat = base_lapc3d_init(par%ax%n,par%ay%n,par%az%n)
        nAop = par%ax%n*par%ay%n*par%az%n
        call laplac2_init(par%ax%n,par%az%n)

        !call solver_smp(m=dS, d=dI, Fop=posmva_tomgph_oper, stepper=cgstep, niter=50, verb=.true.)

        !call solver_reg(m=dS, d=dI, Fop=posmva_tomgph_oper, Aop=laplac2_lop, nAop=nAop, stepper=cgstep, niter=30, eps=20000., verb=.true.)

        call solver_irls(m=dS, d=dI, Fop=posmva_tomgph_oper, stepper=cgstep, Wdiag=cauchy, &
                         niter=30, verb=.true.)


        !call sep_write(dS(par%jx1:par%jx2,par%jy1:par%jy2,:), par%dS_)
        call sep_write(dS, par%dS_)
     else
        call houaxis(par%dI_, par%arx, 1)
        call houaxis(par%dI_, par%ary, 2)
        call houaxis(par%dI_, par%az , 3)
        call sep_read (dS, par%dS_)
        stat = posmva_tomgph_oper(adj, .false., dS, dI)
        call sep_write(dI, par%dI_)
     end if

  end if

  write(0,*)"------------------------------------------------"
  write(0,*)"       FINISHED!!!                              "
  write(0,*)"------------------------------------------------"

end program
