module solver_smp_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 solver_smp( m,d, Fop, stepper, niter   ,             Wop&
    &,Jop,m0,err,resd,mmov,rmov,verb)
    optional :: Wop,Jop,m0,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 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(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
    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( 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(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 
      stat = stepper(forget, m,g, rr,gg)                  
      !m+=dm; R+=dR
      if (stat .eq.1) then
        exit ! got stuck descending
      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
      forget=F
    end do 
    if (present( resd)) then
      resd = rd
    end if
  end subroutine 
end module 
