# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module bob_solver_mod
! 0 = W (F J m - d)
  use chain0_mod
  use solver_report_mod
  implicit none
  logical, parameter, private :: T = .true., F = .false.
  contains
  subroutine bob_solver( m,d, Fop, stepper, niter , Wop&
    &,Jop,m0,err,resd,mmov,rmov,verb,nfreq,rwt)
    optional :: Wop,Jop,m0,err,resd,mmov,rmov,verb,nfreq,rwt
    interface
!-------------------------- begin definitions -----------
      integer function Fop(adj,add,m,d)
        real::m(:),d(:)
        logical,intent(in)::adj,add
      end function
      integer function Wop(adj,add,m,d)
        real::m(:),d(:)
        logical,intent(in)::adj,add
      end function
      integer function Jop(adj,add,m,d)
        real::m(:),d(:)
        logical,intent(in)::adj,add
      end function
      integer function rwt(r,wt)
        real::r(:),wt(:)
      end function
      integer function stepper(forget,m,g,rr,gg)
        real, dimension(:) :: m,g,rr,gg
        logical :: forget
      end function
    end interface
    real, dimension(:), intent(in) :: d, m0
    integer, intent(in) :: niter
    logical, intent(in) :: verb
    real, dimension(:), intent(out) :: m,err, resd
    real, dimension(:,:), intent(out) :: rmov, mmov
    real, dimension(:) ,pointer:: g
    real, dimension(:), pointer :: rr, gg
    real, dimension(:), pointer :: tt
    real, dimension(:), pointer :: rd, gd, td
    real, dimension(:), pointer :: rm, gm, tm,wt
    integer :: iter, stat,nfreq
    logical :: forget,reweight
    allocate(g(size(m)))
    allocate(rr(size(d)),gg(size(d)))
    allocate(tt(size(d)+size(m)))
!this makes the implicit assumption that the Wop range is the same as its domain
    reweight=.false.
    if (present(nfreq)) then
      if (nfreq < niter) then
        if (.not. present(rwt)) then
          call seperr("if reweighting must supply rwt")
        end if
        reweight=.true.
        allocate(wt(size(d)))
        wt=1
      end if
    end if
    rd => rr(1:size(d))
    gd => gg(1:size(d))
    td => tt(1:size(d))
    tm => tt(1+size(d):)
    if (present( Wop)) then
      stat=Wop(F,F,-d,rd) ! begin initialization --------
    else
      rd = -d !Rd = -W d
    end if
    if (present( m0)) then
      m=m0 !m = m0
      if (present( Wop)) then
        call chain0(Wop,Fop,F,T,m,rd,td)
      else
        stat = Fop(F,T,m,rd ) !Rd+= WF m0
      end if
    else
      m=0
    end if
    forget = T
!-------------------------- begin iterations ------------
    do iter = 1,niter
      if (present(nfreq) .and. present(rwt)) then
        if (mod(iter,nfreq).eq.0 ) then
          stat=rwt(rd,wt)
        end if
      end if
      if (reweight) then
        rd=rd*wt
      end if
      if (present(Wop)) then
        call chain0(Wop,Fop,T,F,g,rd,td)
      else
        stat = Fop(T,F,g,rd ) !g = (WF)'Rd
      end if
      if (present(Jop)) then
        tm=g
        stat = Jop(F,F,tm, g )
      end if
!g = J g
      if (present(Wop)) then
        call chain0(Wop,Fop,F,F,g,gd,td)
      else
        stat = Fop(F,F,g,gd ) !Gd = (WF) g
      end if
      if (reweight) then
        gd=gd*wt
      end if
! write(0,*) "CHECU",sum(m),sum(g),sum(rr),sum(gg)
      stat = stepper(forget, m,g, rr,gg)
      !m+=dm; R+=dR
!write(0,*) "CHECK IT",rd(229),iter
! if(stat ==1) exit # got stuck descending
      if (present( mmov)) then
        mmov(:,iter) = m(:size(mmov,1)) ! report -----
      end if
      if (present( rmov)) then
        rmov(:,iter) = rd(:size(rmov,1))
      end if
      if (present( err )) then
        err( iter) = dot_product(rd,rd)
      end if
      if (present( verb)) then
        if (verb) then
          call solver_report(iter,m,g,rd)
        end if
      end if
      forget=F
    end do
    if (reweight) then
      deallocate(wt)
    end if
    if (present( resd)) then
      resd = rd
    end if
    deallocate(g,rr,gg,tt)
  end subroutine
end module
