module solver_smp_mod {                        # 0 = W (F J m - d)
  use chain0_mod + solver_report_mod
  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}
       integer function Wop(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function Jop(adj,add,m,d){real::m(:),d(:);logical,intent(in)::adj,add}
       integer function stepper(forget,m,g,rr,gg) {
         real, dimension(:)  ::        m,g,rr,gg    
         logical             :: forget            }
    }
    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)) stat=Wop(F,F,-d,rd) # begin initialization --------
    else rd = -d                                           #Rd = -W  d
    if(present( m0)){ m=m0                                 #m  =     m0
	if(present( Wop)) call chain0(Wop,Fop,F,T,m,rd,td)
	else             stat =           Fop(F,T,m,rd   ) #Rd+=  WF m0
    } else m=0
    forget = T; #-------------------------- begin iterations ------------
    do iter = 1,niter {
       if(present(Wop)) call chain0(Wop,Fop,T,F,g,rd,td) 
       else             stat =          Fop(T,F,g,rd   )   #g  = (WF)'Rd
       if(present(Jop)){ tm=g;  stat  = Jop(F,F,tm, g  )}  #g  =   J  g
       if(present(Wop)) call chain0(Wop,Fop,F,F,g,gd,td) 
       else             stat =          Fop(F,F,g,gd   )   #Gd = (WF) g 
       stat = stepper(forget, m,g, rr,gg)                  #m+=dm; R+=dR
       if(stat ==1) exit # got stuck descending
       if(present( mmov)) mmov(:,iter) =  m(:size(mmov,1)) # report -----
       if(present( rmov)) rmov(:,iter) = rd(:size(rmov,1))
       if(present( err ))  err(  iter) = dot_product(rd,rd)
       if(present( verb)){ if(verb) call solver_report(iter,m,g,rd)}
       forget=F
    }
    if(present( resd)) resd = rd
  }
}