# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module tomo_lin_cg_mod
  use bob_solver_mod
  use cgstep_mod
  use inv_basic_func_mod
  use reweight_mod
  use write_inv_params_mod
  use multiple_realizations_mod
  implicit none
  real, pointer,private,save :: dat(:)
  real, pointer,private,save :: mod(:)
  character(len=128), private, save :: inv_method
  integer, private, save :: niter_reg,niter_prec,nxtra
  real, private, save :: eps_reg,eps_prec
  logical, private,save :: verb
  logical,private, save :: write_xmovie,write_resid,write_rmovie
  logical, private,save :: styling_reg
  integer,private,save :: nreweight,nd1
  integer,private,save :: nm(3),nmod
  contains
  logical function init_lin_cg_inverse()
    real, pointer :: dat1(:),dat2(:),dat3(:),temp1(:),sc&
      &(:)
    type(sep3d) :: vel_s
    real :: om(3),dm(3)
    integer :: ierr
    init_lin_cg_inverse=.false.
    if (.not. tomo_grab_dat0(dat1)) then
      write(0,*) "      ERROR: grabbing data for tomograhy"
      return
    end if
    nd1=size(dat1)
    if (.not. add_realization_d_noise(dat1)) then
      write(0,*) "      ERROR: adding realization noise to data"
      return
    end if
    if (.not. vel_tomo_structure(vel_s)) then
      write(0,*) "      ERROR: grabbing vel_tomo_structure"
      return
    end if
    call set_n_o_d_3d(vel_s,nm,om,dm)
!INIT
    if (model_des_method(1:4).ne."NONE") then
      nmod=model_des_ndomain()
    else
      nmod=product(nm)
    end if
!i am making the assumption that my preconditioning operator operatres
! on the same space as my model
    if (.not. init_chain_extra(get_extra_nrange())) then
      write(0,*) "      ERROR: with init_chain_extra"
      return
    end if
    if (.not. get_model_styling_dat0(dat2)) then
      write(0,*) "      ERROR: getting model styling dat0"
      return
    end if
    if (.not. add_realization_m_noise(dat2)) then
      write(0,*) "      ERROR: adding realization noise to data"
      return
    end if
    if (extra_op_method(1:4).ne."NONE") then
      if (.not. get_extra_dat0(dat3)) then
        write(0,*) "      ERROR: getting extra dat0"
        return
      end if
      if (.not. add_realization_e_noise(dat3)) then
        write(0,*) "      ERROR: adding realization noise to extra"
        return
      end if
      allocate(dat(size(dat1)+size(dat2)+size(Dat3)+product(nm)))
      dat(1 :size(dat1)) =dat1
      dat(1+size(dat1) :size(dat1)+size(dat2))=dat2
      dat(1+size(dat1)+size(Dat2):size(dat)-product(nm))=dat3
      nxtra=size(dat3)
      deallocate(dat3)
    else
      allocate(dat(size(dat1)+size(dat2)+product(nm)))
      dat(1 :size(dat1)) =dat1
      dat(1+size(dat1) :size(dat1)+size(dat2))=dat2
      nxtra=0
    end if
    if (.not. setup_write_files(size(dat1),size(dat2),size(dat)-size&
      &(dat1)-size(dat2),niter_reg+niter_prec)) then
      write(0,*) "      ERROR: setting up cg_write files"
      return
    end if
    nmod=model_des_ndomain()
    if (.not. init_model_chains(nmod)) then
      write(0,*) "      ERROR: initializing model chains"
      return
    end if
    allocate(temp1(size(dat1)))
    temp1=1.
    allocate(sc(product(nm)))
    sc=0
    ierr=tomo_op(.true.,.true.,sc,temp1)
!sc=1./(abs(sc)+.01)
    if (.not. init_punish_change(sc)) then
      write(0,*) "      ERROR: trouble initializing punish change"
      return
    end if
    deallocate(temp1,sc)
    dat(size(dat)-product(nm)+1:)=return_change_sum()
    deallocate(dat1,dat2)
    call add_mem(size(dat))
    allocate(mod(nmod))
    mod=0.
    init_lin_cg_inverse=.true.
  end function
  logical function lin_cg_inverse()
    real, pointer :: temp(:)
    real, allocatable :: junk(:),junk2(:)
    integer :: ierr
    lin_cg_inverse=.false.
    if (niter_prec >0) then
      if (.not. do_prec()) then
        write(0,*) "      ERROR: setting to preconditioned problem"
        return
      end if
      if (.not. iterate_it(niter_prec)) then
        write(0,*) "      ERROR: iterating preconditioning"
        return
      end if
    end if
    if (niter_reg>0) then
      if (.not. do_reg()) then
        write(0,*) "      ERROR: setting to preconditioned problem"
        return
      end if
      if (.not. iterate_it(niter_reg)) then
        write(0,*) "      ERROR: iterating preconditioning"
        return
      end if
    end if
    allocate(temp(product(nm)))
    ierr=model_transform_op(.false.,.false.,mod,temp)
! call srite("mine.H",temp,size(temp)*4)
    ierr=extra_op(.false.,.false.,temp,mod)
! call srite("mtest2.H",mod,size(mod)*4)
    if (.not. add_change(temp,.true.)) then
      write(0,*) "      ERROR: setting change"
      return
    end if
    deallocate(temp)
    lin_cg_inverse=.true.
  end function
  logical function iterate_it(niter)
    integer :: niter
    integer :: npt1,npt2
    real, pointer :: resid(:),rmovie(:,:),xmovie(:,:)
    iterate_it=.false.
    if (write_resid) then
      allocate(resid(size(dat)))
    else
      nullify(resid)
    end if
    if (write_rmovie) then
      allocate(rmovie(size(dat),niter))
    else
      nullify(rmovie)
    end if
    if (write_xmovie) then
      allocate(xmovie(size(mod),niter))
    else
      nullify(xmovie)
    end if
    write(0,*) "SIZE OF DATA",size(dat)
    call srite("dat_test.H",dat,size(dat)*4)
    if (write_resid) then
      if (write_rmovie) then
        if (write_xmovie) then
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,mmov=xmovie,resd=resid,rmov=rmovie,verb=verb)
        else
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,resd=resid,rmov=rmovie,verb=verb)
        end if
      else
        if (write_xmovie) then
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,mmov=xmovie,resd=resid,verb=verb)
        else
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,resd=resid,verb=verb)
        end if
      end if
    else
      if (write_rmovie) then
        if (write_xmovie) then
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,mmov=xmovie,rmov=rmovie,verb=verb)
        else
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,rmov=rmovie,verb=verb)
        end if
      else
        if (write_xmovie) then
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op ,mmov=xmovie,verb=verb)
        else
          call bob_solver(mod,dat,total_op,cgstep,niter,nfreq=nreweight&
            &,rwt=reweight_op,verb=verb)
        end if
      end if
    end if
    if (.not. write_out_solve_info(resid,rmovie,xmovie)) then
      write(0,*) "      ERROR: writing out relevant info"
      return
    end if
    write(0,*) "RESIDE SIZE",size(resid)
    call srite("dat_test.H",resid,size(resid)*4)
    if (associated(rmovie)) then
      deallocate(rmovie)
    end if
    if (associated(xmovie)) then
      deallocate(xmovie)
    end if
    if (associated(resid)) then
      deallocate(resid)
    end if
    iterate_it=.true.
  end function
  logical function check_lin_cg_pars()
    character(len=50) :: temp
    check_lin_cg_pars=.false.
    call from_param("verb",verb,.false.)
    call from_param("niter_reg",niter_reg,0)
    call from_param("niter_prec",niter_prec,0)
    if (niter_reg.eq.0 .and. niter_prec.eq.0) then
      write(0,*) "    ERROR-Having 0 iterations of preconditioning and&
        & regularization isn't logical"
      return
    end if
    if (.not. check_reweight_params(niter_prec+niter_reg,nreweight))&
      & then
      write(0,*) "      ERROR: grabbing weighting params"
      return
    end if
    if (.not. check_write_inv_pars(niter_prec+niter_reg,write_resid&
      &,write_rmovie,write_xmovie)) then
      write(0,*) "      ERROR: writing inverse parameters"
      return
    end if
    check_lin_cg_pars=.true.
  end function
end module
