# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
module solver_irls_mod
! 0 = Wdiag Wop (Fop Jop m - d) and periodic restart
  use chain0_mod
  use solver_report_mod
  implicit none
  logical, parameter, private :: T = .true., F = .false.
  contains
  subroutine solver_irls( m,d, Fop, stepper, niter , Wop&
    &,Jop,Wdiag,m0,nmem,nfreq,err,resd,mmov,rmov,verb)
    optional :: Wop,Jop,Wdiag,m0,nmem,nfreq,err,resd,mmov,rmov,verb
    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 Wdiag(res, w)
        real::res(:),w(:)
      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, nmem, nfreq
    logical, intent(in) :: verb
    real, dimension(:), intent(out) :: m,err, resd
    real, dimension(:,:), intent(out) :: rmov, mmov
    real, dimension(size(m)) :: g
    real, dimension(size(d)), target :: rr, gg
    real, dimension(size(d)+size(m)), target :: tt
    real, dimension(:), pointer :: rd, gd, td
    real, dimension(:), pointer :: rm, gm, tm, wht
    integer :: iter, stat
    logical :: forget
    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(Wdiag)) then
      allocate(wht(size(d)))
    end if
    wht=1.
    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 = F
!-------------------------- begin iterations ------------
    do iter = 1,niter
      if (present(nmem)) then
        forget = (iter>nmem) ! initialize forget
      end if
      if (present(Wdiag).and.forget) then
        stat = Wdiag(rd,wht) ! estimate weights
      end if
      if (present(Wdiag)) then
        rd = wht*rd
        if (present(Wop)) then
          call chain0(Wop,Fop,T,F,g,wht*rd,td)
        else
          stat = Fop(T,F,g,wht*rd ) ! g = (WF)'Rd
        end if
      else
        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
      end if
      if (forget.and.present(nfreq)) then
        forget=(mod(iter,nfreq).eq.0) ! periodic restart
      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 (present(Wdiag)) then
        gd=wht*gd
      end if
      stat = stepper(forget, m,g, rr,gg)
      !m+=dm; R+=dR
      if (stat .eq.1) then
        exit ! got stuck descending
      end if
      if (present(Wdiag)) then
        rd = -d
        stat = Fop(F,T,m,rd)
        if (present( Wop)) then
          td = rd
          stat = Wop(F,F,td,rd)
        end if
      end if
      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
    end do
    if (present( resd)) then
      resd = rd
    end if
  end subroutine
end module