module solver_irls_mod {                        # 0 = Wdiag Wop (Fop Jop m - d) and periodic restart
  use chain0_mod + solver_report_mod
  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}
       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 Wdiag(res, w)     {real::res(:),w(:)}
       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, 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)) stat=Wop(F,F,-d,rd) # begin initialization --------
    else rd = -d                                           #Rd = -W  d
    if(present(Wdiag))allocate(wht(size(d)))
    wht=1.
    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 = F; #-------------------------- begin iterations ------------
    do iter = 1,niter {
       if(present(nmem)) forget = (iter>nmem)                     # initialize forget
       if(present(Wdiag).and.forget) stat = Wdiag(rd,wht)         # estimate weights
       if(present(Wdiag)) {
          rd = wht*rd	
          if(present(Wop)) call chain0(Wop,Fop,T,F,g,wht*rd,td) 
          else             stat =          Fop(T,F,g,wht*rd   )   # g  = (WF)'Rd
       } else {
          if(present(Wop)) call chain0(Wop,Fop,T,F,g,    rd,td) 
          else             stat =          Fop(T,F,g,    rd   )   # g  = (WF)'Rd
       }
       if(forget.and.present(nfreq)) forget=(mod(iter,nfreq)==0)  # periodic restart	
       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 
       if(present(Wdiag)) gd=wht*gd
       stat = stepper(forget, m,g, rr,gg)                  #m+=dm; R+=dR
       if(stat ==1) exit # got stuck descending
       if(present(Wdiag)) { 
            rd = -d
            stat = Fop(F,T,m,rd)
            if(present( Wop)){ 
            td = rd
            stat = Wop(F,F,td,rd)} 
       }
       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)}
    }
    if(present( resd)) resd = rd
  }
}