# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
module solver_prc_mod
! 0 = W (F S J p - d)
  use chain0_mod
  use solver_report_mod ! 0 = I p
  implicit none
  logical, parameter, private :: T = .true., F = .false.
  contains
  subroutine solver_prc( m,d, Fop, Sop, stepper, nSop, niter,eps , &
    & Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb)
    optional :: Wop,Jop,p0,rm0,err,resd,resm,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 Sop(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 stepper(forget,m,g,rr,gg)
        real, dimension(:) :: m,g,rr,gg
        logical :: forget
      end function
    end interface
    real, dimension(:), intent(in) :: d, p0,rm0
    integer, intent(in) :: niter, nSop
    logical, intent(in) :: verb
    real, intent(in) :: eps
    real, dimension(:), intent(out) :: m,err, resd,resm
    real, dimension(:,:), intent(out) :: rmov,mmov
    real, dimension(size( m)) :: p , g
    real, dimension(size( d) + nSop), target :: rr, gg, tt
    real, dimension(:), pointer :: rd, gd, td
    real, dimension(:), pointer :: rm, gm, tm
    integer :: iter, stat
    logical :: forget
    rd => rr(1:size(d))
    rm => rr(1+size(d):)
    gd => gg(1:size(d))
    gm => 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
    rm = 0.
    if (present(rm0)) then
      rm=rm0 !Rm = Rm0
    end if
    if (present( p0)) then
      p=p0 ! p = p0
      if (present( Wop)) then
        call chain0(Wop,Fop,Sop,F,T,p,rd,tm,td)
      else
        call chain0( Fop,Sop,F,T,p,rd,tm )!Rd += WFS p0
      end if
      rm = rm + eps*p
      !Rm += e I p0
    else
      p=0
    end if
    forget = T
!-------------------------- begin iterations ------------
    do iter = 1,niter
      if (present(Wop)) then
        call chain0(Wop,Fop,Sop,T,F,g,rd,tm,td)
      else
        call chain0( Fop,Sop,T,F,g,rd,tm ) !g = (WFS)'Rd
      end if
      g = g + eps*rm
      !g += e I'Rm
      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,Sop,F,F,g,gd,tm,td)
      else
        call chain0( Fop,Sop,F,F,g,gd,tm ) !Gd = (WFS) g
      end if
      gm = eps*g
      !Gm = e I g
      stat = stepper(forget, p,g, rr,gg)
      !m+=dm; R+=dR
      if (stat .eq.1) then
        exit ! got stuck descending
      end if
      stat = Sop(F,F,p,m)
      !m = S p
      if (present( mmov)) then
        mmov(:,iter) = m(:size(mmov,1)) ! report -----
      end if
      if (present( rmov)) then
        rmov(:,iter) = rr(: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,rm)
        end if
      end if
      forget=F
    end do
    if (present( resd)) then
      resd = rd
    end if
    if (present( resm)) then
      resm = rm(:size(resm))
    end if
  end subroutine
end module