module solver_prc_mod{                         # 0 = W (F S J p - d)  
  use chain0_mod + solver_report_mod           # 0 =      I   p
  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}
       integer function Sop(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, 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)) stat=Wop(F,F,-d,rd) # begin initialization --------
    else rd = -d                                                #Rd  = -W    d
    rm = 0.; if(present(rm0)) rm=rm0			        #Rm  =      Rm0
    if(present( p0)){ p=p0			                # p  =       p0
       if(present( Wop)) call chain0(Wop,Fop,Sop,F,T,p,rd,tm,td)
       else              call chain0(    Fop,Sop,F,T,p,rd,tm   )#Rd +=  WFS  p0
       rm = rm + eps*p                                          #Rm +=   e I p0
    } else p=0
    forget = T; #-------------------------- begin iterations ------------
    do iter = 1,niter {
       if(present(Wop)) call chain0(Wop,Fop,Sop,T,F,g,rd,tm,td)
       else             call chain0(    Fop,Sop,T,F,g,rd,tm   ) #g  = (WFS)'Rd
       g = g + eps*rm                                           #g +=   e I'Rm
       if(present(Jop)){ tm=g;         stat=Jop(F,F,tm,g      )}#g  =     J g
       if(present(Wop)) call chain0(Wop,Fop,Sop,F,F,g,gd,tm,td)
       else             call chain0(    Fop,Sop,F,F,g,gd,tm   ) #Gd = (WFS) g
       gm   =  eps*g                                            #Gm =   e I g
       stat = stepper(forget, p,g, rr,gg)                       #m+=dm; R+=dR
       if(stat ==1) exit # got stuck descending
       stat = Sop(F,F,p,m)                                      #m  = S p
       if(present( mmov)) mmov(:,iter) =  m(:size(mmov,1)) # report -----
       if(present( rmov)) rmov(:,iter) = rr(:size(rmov,1))
       if(present( err ))  err(  iter) = dot_product(rd,rd)
       if(present( verb)){ if(verb) call solver_report(iter,m,g,rd,rm)}
       forget=F
    }
    if(present( resd)) resd = rd
    if(present( resm)) resm = rm(:size(resm))
  }
}