module solver_sep_mod  
! 0 = F m - d
        use cgstep_module
        use chain0_mod 
        use ddot_mod
        use quantile_mod 
   
        implicit none
        logical, parameter, private     :: T = .true., F = .false.
        logical,            private     :: time
        real    , private               :: scale
        real    ,save                   :: tstart, tfinish, ttot, tinv, telapse, tleft, tnow, told
        integer ,save                   :: tarray(8), tarray1(8), tarray2(8), itniter, itleft
        real    ,save                   :: treshold,quantil
	character(len=256),private,save :: report_
        logical,            private     :: report_fil

        contains

!#################################    
  subroutine solver_tiny( m, d, Fop, stepper, niter, verb, resd, m0, timer)
        optional ::     verb, resd,m0, timer 
        interface  
                integer function Fop(adj,add,m,d); real,dimension(:) :: 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       ! data
        real, dimension(:),    intent(out) :: m       ! model
        real, dimension(:),    intent(in)  :: m0      ! initial model
        integer,               intent(in)  :: niter   ! number of iterations
        integer                            :: iter    ! iteration number    
        real, dimension(size( m))          :: g       ! gradient (dm)
        real, dimension(size( d)),target   :: rr      ! data residual ( vector)
        real, dimension(:),       pointer  :: rd      ! data residual (pointer)
        real, dimension(size( d)),target   :: gg      ! conj grad     ( vector)
        real, dimension(:),       pointer  :: gd      ! conj grad     (pointer)
        integer                            :: stat    ! status flag
        real, dimension(:),    intent(out) :: resd    ! residual 
        logical                            :: verb        
        logical                            :: timer

                if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif

                rd => rr(1:size( d))
                gd => gg(1:size( d))
                rd = -d                                         ! Rd =      - d
                m  = 0
                if (present( m0)) then 
                        m = m0                                  !  m =        m0
                        stat = Fop(.false.,.true.,m,rd)         ! Rd = Rd + F m0
                end if
                do iter = 1,niter  
                        stat = Fop( .true.,.false.,g,rd)        !  g =  F' Rd
                        stat = Fop(.false.,.false.,g,gd)        ! Gd =  F  g
                        stat = stepper(.false., m,g, rr,gg)     ! m = m+dm; R = R + dR
                        if (present( verb)) then; if (verb) then; call solver_report(iter,m,g,rd); end if; end if
                end do 
                if (present( resd)) then
                        resd = rd
                end if 
end subroutine 
!
!!#################################    
!  subroutine solver_tiny_freq( m, d, Fop, stepper, niter, verb, resd, m0, timer)
!        optional ::     verb, resd,m0, timer 
!        interface  
!                integer function Fop(adj,add,m,d); real,dimension(:) :: 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 
!        complex, dimension(:), intent(in)   :: d       ! data
!        real, dimension(:),    intent(out)  :: m       ! model
!        real, dimension(:),    intent(in)   :: m0      ! initial model
!        integer,               intent(in)   :: niter   ! number of iterations
!        integer                             :: iter    ! iteration number    
!        real, dimension(size( m))           :: g       ! gradient (dm)
!        complex, dimension(size( d)),target :: rr      ! data residual ( vector)
!        complex, dimension(:),       pointer:: rd      ! data residual (pointer)
!        complex, dimension(size( d)),target :: gg      ! conj grad     ( vector)
!        complex, dimension(:),       pointer:: gd   ! conj grad     (pointer)
!        integer                             :: stat    ! status flag
!        real, dimension(:),    intent(out)  :: resd    ! residual 
!        logical                             :: verb        
!        logical                             :: timer
!
!                if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif
!
!                rd => rr(1:size( d))
!                gd => gg(1:size( d))
!                rd = -d                                         ! Rd =      - d
!                m  = 0
!                if (present( m0)) then 
!                        m = m0                                  !  m =        m0
!                        stat = Fop(.false.,.true.,m,rd)         ! Rd = Rd + F m0
!                end if
!                do iter = 1,niter  
!                        stat = Fop( .true.,.false.,g,rd)        !  g =  F' Rd
!                        stat = Fop(.false.,.false.,g,gd)        ! Gd =  F  g
!                        stat = stepper(.false., m,g, rr,gg)     ! m = m+dm; R = R + dR
!                        if (present( verb)) then; if (verb) then; call solver_report(iter,m,g,rd); end if; end if
!                end do 
!                if (present( resd)) then
!                        resd = rd
!                end if 
!end subroutine 
!#################################    
! 0 = W (F J m - d)
! 0 =    A   m
subroutine solver_regs( m,d, Fop, Aop, stepper, nAop, niter,eps, Wop,Jop,m0,rm0,err,resd,resm,mmov,rmov,verb,timer)
        optional :: Wop,Jop,m0,rm0,err,resd,resm,mmov,rmov,verb, timer
        interface  
!-------------------------- begin definitions -----------
                integer function Fop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Aop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Wop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Jop(adj,add,m,d); real,dimension(:) :: 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,rm0
        integer,               intent(in)          :: niter, nAop
        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))                  :: g
        real, dimension(size( d) + nAop), target   :: rr, gg, tt
        real, dimension(:), pointer                :: rd, gd, td
        real, dimension(:), pointer                :: rm, gm, tm 
        integer                                    :: iter, stat
        logical                                    :: forget
        logical                                    :: timer

        if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif

        rd => rr(1:size(d)) ; gd => gg(1:size(d)) ; td => tt(1:size(d)); 
        rm => rr(1+size(d):); gm => gg(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( 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 
                stat = Aop(F,T,eps*m0,rm)                               !Rm += e A m0
        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);end if       !g  = (WF)'Rd 
                stat =  Aop(T,T,g,eps*rm)                                                                               !g +=  e A'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,F,F,g,gd,td); else;        stat=Fop(F,F,g,gd);end if       !Gd  = (WF) g
                stat =  Aop(F,F,eps*g,gm)                                                                               !Gm  =  e A g
                stat = stepper(forget, m,g, rr,gg)                                                                      !m+=dm; R+=dR
                if (stat .eq.1) exit                                                                                    ! got stuck descending
                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)) then; if (verb) then; call solver_report(iter,m,g,rd,rm); end if; end if
                forget=F 
        end do 
        if (present( resd)) resd = rd
        if (present( resm)) resm = rm(:size(resm))
end subroutine 
!################################# two regs
! 0 = W (F J m - d)
! 0 =    A   m
subroutine solver_regs2( m,d, Fop, Aop1, Aop2, stepper, nAop, niter,epsi,lmbd, Wop,Jop,m0,rm0,err,resd,resm,mmov,rmov,verb,timer)
        optional :: Wop,Jop,m0,rm0,err,resd,resm,mmov,rmov,verb, timer
        interface  
!-------------------------- begin definitions -----------
                integer function Fop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Aop1(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Aop2(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Wop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Jop(adj,add,m,d); real,dimension(:) :: 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,rm0
        integer,               intent(in)          :: niter, nAop
        logical,               intent(in)          :: verb
        real, dimension(:),    intent(in)          :: epsi,lmbd
        real, dimension(:),    intent(out)         :: m,err, resd,resm
        real, dimension(:,:),  intent(out)         ::        rmov,mmov
        real, dimension(size( m))                  :: g
        real, dimension(size( d) + nAop), target   :: rr, gg, tt
        real, dimension(:), pointer                :: rd, gd, td
        real, dimension(:), pointer                :: rm, gm, tm 
        integer                                    :: iter, stat
        logical                                    :: forget
        logical                                    :: timer

        if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif

        rd => rr(1:size(d)) ; gd => gg(1:size(d)) ; td => tt(1:size(d)); 
        rm => rr(1+size(d):); gm => gg(1+size(d):); tm => tt(1+size(d):)
        if (present(Wop)) then; stat=Wop(F,F,-d,rd); else; rd = -d; endif !Rd  = -W  d
        rm = 0.
        if (present(rm0)) then; rm=rm0; endif                           !Rm  =    Rm0
        if (present( m0)) then; m=m0  ; 
                if (present(Wop)) then; call chain0(Wop,Fop,F,T,m,rd,td); else; stat=Fop(F,T,m,rd   ); endif !Rd +=  WF m0
                stat = Aop1(F,T,epsi*m0,rm)
                stat = Aop2(F,T,lmbd*m0,rm)                             !Rm2 += eps A1 m0 + lmb A2 m0
        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);end if       !g  = (WF)'Rd 
                stat =  Aop1(T,T,g,epsi*rm)                                                                             !   
                stat =  Aop2(T,T,g,lmbd*rm)                                                                             !g +=  (eps A1'Rm1 +  lmb A2'Rm2)    
                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);end if       !Gd = (WF) g
                stat =  Aop1(F,F,epsi*g,gm)                                                                             !       
                stat =  Aop2(F,T,lmbd*g,gm)                                                                             !Gm  = (eps A1 dm1 + lmb A2 dm2)
                stat = stepper(forget, m,g, rr,gg)                                                                      !m+=dm; R+=dR
                if (stat .eq.1) exit                                                                                    ! got stuck descending
                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)) then; if (verb) then; call solver_report(iter,m,g,rd,rm); end if; end if
                forget=F 
        end do 
        if (present( resd)) resd = rd
        if (present( resm)) resm = rm(:size(resm))
end subroutine 
!#################################    
! 0 = Wdiag Wop (Fop Jop m - d) and periodic restart
subroutine solver_irls( m,d, Fop, stepper, niter, Wop,Jop,Wdiag,Scale,m0,nmem,nfreq,err,resd,mmov,rmov,verb, timer)
        optional :: Wop,Jop,Wdiag,Scale,m0,nmem,nfreq,err,resd,mmov,rmov,verb,timer
        interface  
                integer function Fop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Aop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Wop(adj,add,m,d); real,dimension(:) :: m,d; logical,intent(in):: adj,add; end function 
                integer function Jop(adj,add,m,d); real,dimension(:) :: 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
                integer function Wdiag(res, w);    real::res(:),w(:); end function 
                integer function Scale(iter,niter,adj);integer::iter,niter; logical,intent(in)::adj; end function 		!scale preconditioner

        end interface 
        real, dimension(:),    intent(in)          :: d, m0
        integer,               intent(in)          :: niter, nmem, nfreq
        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,               intent(in)          :: verb
        logical                                    :: forget
        logical                           	   :: timer

        if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif
        rd => rr(1:size(d)); gd => gg(1:size(d)); td => tt(1:size(d)); tm => tt(1+size(d):)
        ! begin initialization --------
        if (present( Wop )) then; stat=Wop(F,F,-d,rd); else; rd = -d ;end if    !Rd = -W  d                     
        if (present(Wdiag)) allocate(wht(size(d))); 
        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);end if   !Rd+=  WF m0
        else
                m=0
        end if
        forget = F
!-------------------------- begin iterations ------------
        do iter = 1,niter  
	      	if(present(Scale))stat= Scale(iter,niter,T)	
                if (present(nmem)) forget = (iter>nmem)                         ! initialize forget
                if (present(Wdiag).and.forget) stat = Wdiag(rd,wht)             ! estimate weights
                if (present(Wdiag)) then                                                                                        !===>[dm = F'W'r]
                        rd = wht*rd; if (present(Wop)) then; call chain0(Wop,Fop,T,F,g,wht*rd,td); else; stat = Fop(T,F,g,wht*rd); endif ! g  = (WF)'Rd using diagonal weights
                else
                                     if (present(Wop)) then; call chain0(Wop,Fop,T,F,g,    rd,td); else; stat = Fop(T,F,g,    rd); endif ! g  = (WF)'Rd using only weighting operator
                end if
                if (forget.and.present(nfreq))forget=(mod(iter,nfreq).eq.0)     ! periodic restart
                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); endif !Gd = (WF) g        !===>[dr = WFdm]
                if (present(Wdiag)) gd=wht*gd
                stat = stepper(forget, m,g, rr,gg)                              !m+=dm; R+=dR                                   !===>[cgstep===>>>m=m+dm,r=r+dr]
                if (stat .eq.1) exit                                            !got stuck descending
                if (present(Wdiag)) then
                        rd = -d; 
                        stat = Fop(F,T,m,rd);
                        if (present( Wop)) then; td = rd; stat = Wop(F,F,td,rd); end if                                         !===>[r=Fm-d]
                end if
                if (present( mmov)) mmov(:,iter) =  m(:size(mmov,1))            
                if (present( rmov)) rmov(:,iter) = rd(:size(rmov,1))
                if (present( err )) err(  iter) = dot_product(rd,rd)
                if (present( verb)) then; if (verb) then; call solver_report(iter,m,g,rd); end if; end if			! report -----
        end do 
        if (present( resd))resd = rd
  end subroutine 
!#################################irls-antoine
  subroutine IRLS( oper, solv, x, dat, niter, nfreq, wght, err,Scale)
        optional ::        err,Scale
        interface; 
                integer function wght( res, w)                  ; real, dimension (:)   :: res, w  ;                                      end function 
                integer function oper( adj, add, x, dat)        ; logical, intent (in)  :: adj, add; real, dimension (:)  :: x, dat      ;end function 
                integer function solv( forget, x, g, rr, gg)    ; logical               :: forget  ;  real, dimension (:) :: x, g, rr, gg;end function
                integer function Scale(iter,niter,adj);integer::iter,niter; logical,intent(in)::adj; end function 		!scale preconditioner
        end interface 
        real, dimension (:), intent (in)  :: dat                ! data
        real, dimension (:), intent (out) :: x                  ! solution
        real, dimension (:), intent (out) :: err                ! error
        integer,             intent (in)  :: niter, nfreq       ! iterations
        real, dimension (size (x))        :: g                  ! gradient
        real, dimension (size (dat))      :: rr, gg, wt         ! res, CG, weight
        integer                           :: i, stat
        logical                           :: forget
        logical				  :: timer
        real                              :: f0

        rr = -dat

        x = 0.
        wt = 1.
        forget = T                                      ! initialize

        do i = 1, niter  
	      	if(present(Scale))stat= Scale(i,niter,T)
                if (forget) then                        !restart
                        rr = - dat
                        stat = oper( F, T, x, rr)       ! rr = Fx - d 
                        stat = wght( rr, wt)            ! compute weighting
                        rr = rr * wt
                end if    
                if(i.eq.1) f0 = sum(dprod(rr,rr))
                if (present(err)) then
                        err(i) = sum(dprod(rr,rr))/f0
                endif
                write(0,*) 'residual(',i,')=',sum(dprod(rr,rr))/f0
                stat = oper( T, F, g, wt*rr)            ! g  = F' W' rr
                stat = oper( F, F, g, gg)               ! gg = Fg
                gg = gg * wt      
                stat = solv( forget, x, g, rr, gg)      ! step in x and rr  
                forget = ( mod( i, nfreq) .eq. 0)       ! periodic restart                     
        end do
 end subroutine 
!#################################irls weights    
    subroutine IRLS_init (tr)
      real   :: tr
      treshold = tr
    end subroutine IRLS_init

    function weight_IRLS(res,w) result(stat)
       integer                    ::stat,i
       real, dimension(:)         ::res,w  
        do i=1, size(res)
                w(i) = 1./sqrt(sqrt(1.+res(i)*res(i)/(treshold*treshold)))
        end do
        stat = 0
    end function weight_IRLS

    subroutine l1_init (qnt)
      real   :: qnt
      quantil = qnt
    end subroutine l1_init

    integer function l1 (res, weight)
      real, dimension (:)  :: res, weight
      real                 :: rbar
      	 if(quantil==0)quantil=0.5
         rbar = quantile( int( quantil*size(res)), abs (res))         
         weight = 1. / sqrt( sqrt (1. + (res/rbar)**2))
         l1 = 0
  end function l1 
!#################################   
  subroutine solver_prec( oper, solv, prec, nprec, x, dat, niter, eps, p0, timer)
    optional :: p0,timer
    interface  
      integer function oper( adj, add, x, dat);logical, intent (in)     :: adj, add; real, dimension (:)        :: x, dat       ; end function 
      integer function solv( forget, x, g, rr, gg); logical             :: forget; real, dimension (:)          :: x, g, rr, gg ; end function 
      integer function prec( adj, add, x, dat);    logical, intent (in) :: adj, add; real, dimension (:)        :: x, dat       ; end function
    end interface 
    real, dimension (:), intent (in)  :: dat, p0           ! data, initial
    real, dimension (:), intent (out) :: x                 ! solution
    real,                intent (in)  :: eps               ! scaling
    integer,             intent (in)  :: niter, nprec      ! iterations, size
    real, dimension (nprec)           :: g, p              ! gradient, precond
    real, dimension (size (dat) + nprec), target :: rr, gg ! residual, conj grad
    real, dimension (:), pointer      :: rd, rm, gd, gm
    integer                           :: i, stat1, stat2, stat
    logical                           :: timer
        if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif

    rm => rr (1 : nprec) 
    rd => rr (1 + nprec :)          ! model and data resids
    gm => gg (1 : nprec) 
    gd => gg (1 + nprec :)          ! model and data grads
    rd = -dat                       ! initialize r_d
    if ( present( p0)) then
      stat1 = prec( F, F, p0, x)       ! x = S p
      stat2 = oper( F, T, x, rd)       ! r_d = L x - dat
      p = p0 
      rm = p0*eps
! start with p0
    else
      p = 0. 
      rm = 0.
    end if 
! start with zero
    do i = 1, niter  
      stat2 = oper( T, F, x, rd)    
      stat1 = prec( T, F, g, x)        ! g = S' L' r_d
      g = g + eps*rm                   ! g = S' L' r_d + eps I r_m
      stat1 = prec( F, F, g, x)     
      stat2 = oper( F, F, x, gd)       ! g_d = L S g
      gm = eps*g                       ! g_m = eps I g
      stat = solv (F, p, g, rr, gg)    ! step in p and rr 
    end do 
    stat1 = prec( F, F, p, x)           ! x = S p
  end subroutine 
!#################################   
! 0 = W (F S J p - d)
  subroutine solver_prc( m,d, Fop, Sop, stepper, nSop, niter,eps   ,  &
    &           Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer)
    optional :: Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer
    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
    logical                                    :: timer

    if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif
    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):)
! begin initialization --------
    if (present( Wop)) then; stat=Wop(F,F,-d,rd) ; else; rd = -d; end if !Rd  = -W    d
    rm = 0.
    if (present(rm0)) then; rm=rm0; end if      !Rm  =      Rm0
    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   );end if    !Rd +=  WFS  p0
      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   ); end if    !g  = (WFS)'Rd
      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   ); end if    !Gd = (WFS) g
      gm   =  eps*g                                                                                                             !Gm =   e I g
      stat = stepper(forget, p,g, rr,gg)                                                                                        !m+=dm; R+=dR
      if (stat .eq.1) then; exit; end if                                                                                        ! got stuck descending
      stat = Sop(F,F,p,m)                                                                                                       !m  = S p
      if (present( mmov)) then; mmov(:,iter) =  m(:size(mmov,1)); end if                                                        ! report -----
      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))resd = rd
    if (present( resm))resm = rm(:size(resm))
  end subroutine 
!#################################   
! 0 = W (F S J p - d) using eps vector
  subroutine solver_prc2( m,d, Fop, Sop, stepper, nSop, niter,eps   ,  &
    &           Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer)
    optional :: Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer
    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, dimension(:),    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
    logical                                    :: timer

    if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif
    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):)
! begin initialization --------
    if (present( Wop)) then; stat=Wop(F,F,-d,rd) ; else; rd = -d; end if !Rd  = -W    d
    rm = 0.
    if (present(rm0)) then; rm=rm0; end if      !Rm  =      Rm0
    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   );end if    !Rd +=  WFS  p0
      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   ); end if    !g  = (WFS)'Rd
      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   ); end if    !Gd = (WFS) g
      gm   =  eps*g                                                                                                             !Gm =   e I g
      stat = stepper(forget, p,g, rr,gg)                                                                                        !m+=dm; R+=dR
      if (stat .eq.1) then; exit; end if                                                                                        ! got stuck descending
      stat = Sop(F,F,p,m)                                                                                                       !m  = S p
      if (present( mmov)) then; mmov(:,iter) =  m(:size(mmov,1)); end if                                                        ! report -----
      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))resd = rd
    if (present( resm))resm = rm(:size(resm))
  end subroutine 
!#################################   
! 0 = W (F S J p - d) using eps vector
  subroutine solver_prc_scale( m,d, Fop, Sop, stepper, nSop, niter,eps   ,  &
    &           Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer,Scale,report_in)      
    optional :: Wop,Jop,p0,rm0,err,resd,resm,mmov,rmov,verb,timer,Scale,report_in
    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
      integer function Scale(iter,niter,adj);integer::iter,niter; logical,intent(in)::adj; end function 		!scale preconditioner
    end interface 
    real, dimension(:),    intent(in)          :: d, p0,rm0
    integer,               intent(in)          :: niter, nSop 
    logical,               intent(in)          :: verb
    real, dimension(:),    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
    logical                                    :: timer
    character(len=256)                         :: report_in

    if (present(timer))then; time=timer ; itniter=niter; else; time=F; endif;
    if (present(report_in))then; report_=report_in ; report_fil=T; else; report_fil=F; endif;

    rr=0.; gg=0; tt=0.;
    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):)
! begin initialization --------
    if (present( Wop)) then; stat=Wop(F,F,-d,rd) ; else; rd = -d; end if !Rd  = -W    d
    rm = 0.
    if (present(rm0)) then; rm=rm0; end if      !Rm  =      Rm0
    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   );end if    !Rd +=  WFS  p0
      rm = rm + eps*p                           !Rm +=   e I p0
    else
      p=0
    end if 
    forget = T
!-------------------------- begin iterations ------------
    do iter = 1,niter  
      if(present(Scale))stat= Scale(iter,niter,T)	
      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   ); end if    !g  = (WFS)'Rd
      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   ); end if    !Gd = (WFS) g
      gm   =  eps*g                                                                                                             !Gm =   e I g
      stat = stepper(forget, p,g, rr,gg)                                                                                        !m+=dm; R+=dR
      if (stat .eq.1) then; exit; end if                                                                                        ! got stuck descending
      stat = Sop(F,F,p,m)                                                                                                       !m  = S p
      if (present( mmov)) then; mmov(:,iter) =  m(:size(mmov,1)); end if                                                        ! report -----
      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))resd = rd
    if (present( resm))resm = rm(:size(resm))
  end subroutine 
!#################################   
!  subroutine solver_prec_sergey (oper, solv, prec, nprec, x, dat, niter, eps &
!  ,             x0,nmem,nfreq,err,res,xp,nloper,xmov,rmov,wght,verb,known,wt,mwt) 
!    optional :: x0,nmem,nfreq,err,res,xp,nloper,xmov,rmov,wght,verb,known,wt,mwt
!    interface
!       function wght (res, w) result (stat); integer              :: stat; real, dimension (:)  :: res, w                                       ; end function wght
!       function solv (forget, x, g, rr, gg) result (stat); integer:: stat; logical              :: forget;  real, dimension  (:) :: x, g, rr, gg; end function solv
!       function oper (adj, add, x, dat) result (stat); integer    :: stat; logical, intent (in) :: adj, add; real, dimension (:)  :: x, dat     ; end function oper
!       function prec (adj, add, x, dat) result (stat); integer    :: stat; logical, intent (in) :: adj, add; real, dimension (:)  :: x, dat     ; end function prec
!       function nloper (adj, add, x, dat) result (stat); integer  :: stat; logical, intent (in) :: adj, add; real, dimension (:)  :: x, dat     ; end function nloper
!    end interface
!    real,    dimension (:),   intent (in)           :: dat, x0, wt, mwt
!    logical, dimension (:),   intent (in)           :: known
!    real,    dimension (:),   intent (out)          :: x, err, res, xp
!    integer,                  intent (in)           :: niter, nmem, nfreq, nprec
!    real,    dimension (:,:), intent (out)          :: xmov, rmov   
!    real,                     intent (in)           :: eps
!    logical,                  intent (in)           :: verb
!
!    real,    dimension (size (dat) + nprec), target :: p, g
!    real,    dimension (size (dat))                 :: rr, gg
!    real,    dimension (:), pointer                 :: pm, pd, gm, gd, wht
!    integer                                         :: iter, stat, nx, nd
!    integer                                         :: dprr, dppd, dppm, dpgm
!    logical                                         :: forget, v
!
!!    nd = size (dat)
!!    nx = size (x)
!!
!!    pm => p (1 : nprec) ; pd => p (1 + nprec:) ; pd = 0.
!!    gm => g (1 : nprec) ; gd => g (1 + nprec:)
!!
!!    if (present (wt) .or. present (wght)) then; allocate (wht (nd)); if (present (wt)) then; wht = wt; else; wht = 1.; end if; end if
!!
!!    rr = -dat
!!    if (present (x0)) then
!!       pm = x0 
!!       if (present (nloper)) then
!!          if (present (mwt)) then; call chain (nloper, prec, F, T, pm*mwt, rr, x); else; call chain (nloper, prec, F, T, pm, rr, x); end if
!!       else 
!!          if (present (mwt)) then; call chain (  oper, prec, F, T, pm*mwt, rr, x); else; call chain (  oper, prec, F, T, pm, rr, x); end if
!!       end if
!!    else
!!       pm = 0.
!!    end if
!!
!!    v = F ; if (present (verb)) v = verb
!!    forget = F
!!    do iter = 1, niter
!!       if (present (nmem )) forget = (iter > nmem) ! restart
!!       if (present (wght) .and. forget) stat = wght (rr, wht)
!!       if (present (wght) .or. present (wt)) then; rr = eps*pd + wht*rr ; call chain (oper, prec, T, F, gm, rr*wht, x) ; else; call chain (oper, prec, T, F, gm, rr,    x); end if
!!       if (present (mwt)) gm = gm*mwt
!!       gd = eps*rr
!!       if (present (known)) where (known) gm = 0.
!!       if (present (mwt)) then
!!          call chain (oper, prec, F, F, gm*mwt, gg,   x)
!!       else
!!          call chain (oper, prec, F, F, gm, gg,   x)
!!       end if
!!       if (present (wght) .or. present (wt)) then; gg = eps*gd + wht*gg; else; gg = eps*gd + gg; end if
!!       if (forget .and. present (nfreq)) forget = (mod (iter, nfreq) == 0) ! periodic restart
!!       if (v) then
!!          dprr = int (dot_product (rr, rr)/nd)
!!          dppd = int (dot_product (pd, pd)/nd)
!!          dppm = int (dot_product (pm, pm)/nprec)
!!          dpgm = int (dot_product (gm, gm)/nprec)
!!          write (0,*) "iteration ", iter, " res", dprr, " prec dat", dppd, " prec mod", dppm, " grad", dpgm
!!       end if
!!       stat = solv (forget, p, g, rr, gg)
!!       if (present (nloper)) then
!!          rr = eps*pd -dat
!!          if (present (mwt)) then; call chain (nloper, prec, F, T, pm*mwt, rr, x); else; call chain (nloper, prec, F, T, pm, rr, x); end if
!!       else if (present (wght) .or. present (wt)) then
!!          rr =        -dat
!!          if (present (mwt)) then; call chain (  oper, prec, F, T, pm*mwt, rr, x); else; call chain (  oper, prec, F, T, pm, rr, x); end if
!!       else if (present (xmov) .or. iter == niter) then
!!          if (present (mwt)) then 
!!             stat = prec (F, F, pm*mwt, x)
!!          else
!!             stat = prec (F, F, pm, x)
!!          end if
!!       end if
!!       if (present (xmov)) xmov (:,iter) =  x (1 : size (xmov, 1))
!!       if (present (rmov)) rmov (:,iter) =  pd * eps
!!       if (present (err)) err (iter) = dot_product (rr,rr) 
!!    end do
!!
!!    if (present (xp)) then
!!       if (present (mwt)) then
!!          xp = pm*mwt
!!       else
!!          xp = pm
!!       end if
!!    end if
!!    if (present (res)) res = rr
!!    if (present (wght) .and. .not. present (wt)) deallocate (wht)
!  end subroutine solver_prec_sergey
!#################################    
subroutine solver_report(iter,x,g,rd,rm)
        integer,            intent(in) :: iter 
        real, dimension(:), intent(in) :: x,g
        real, dimension(:), pointer    :: rd,rm
        optional                       :: rm
        double precision               :: dd,dm,dx,dg
        dd=dot_product(rd,rd)
        dx=dot_product( x, x) 
        dg=dot_product( g, g)
        dm=0.
        if (present(rm)) then; dm=dot_product(rm,rm); end if
        if (iter.eq.1) then;  scale=max(dd,dm)/10000.; write(0,*) "DOT PRODUCT VALUES SCALED BY ",scale; end if


        if(time)then
                call date_and_time(values=tarray)
                        if (iter==1) then
                                tstart  = (tarray (5) * 3600 + tarray (6) * 60 + tarray (7) + 0.001 * tarray (8))/60; 
                                told=tstart; 
                                if(report_fil)open(unit=3600, file=report_, status='unknown')
                        endif
                        tnow    = (tarray (5) * 3600 + tarray (6) * 60 + tarray (7) + 0.001 * tarray (8))/60
                        tinv    = (tnow-told); told=told+tinv
                        telapse = (tnow - tstart)+tinv
                        itleft  = (itniter-iter+1)
                        tleft   =  itleft*(telapse/iter)
                                      write (0,'(5(A11,1x,I12),A35,3F9.3,A5)')     "iteration:", iter , &
                                                        "| res dat:", int(dd/scale), &
                                                        "| res mod:", int(dm/scale), &
                                                        "| mod:", int(dx/scale), &
                                                        "| grad:", int(dg/scale), & 
                                                        "| Time: tintv, telapse, tleft: ", tinv, telapse, tleft,trim('mins')
                        if(report_fil)write (3600,'(5(A11,1x,I12),A35,3F9.3,A5)')     "iteration:", iter , &
                                                        "| res dat:", int(dd/scale), &
                                                        "| res mod:", int(dm/scale), &
                                                        "| mod:", int(dx/scale), &
                                                        "| grad:", int(dg/scale), & 
                                                        "| Time: tintv, telapse, tleft: ", tinv, telapse, tleft,trim('mins')
        else
                                      write (0,'(5(A11,1x,I12),A35,3F9.3,A5)') "iteration:", iter , &
                                                        "| res dat:", int(dd/scale), &
                                                        "| res mod:", int(dm/scale), &
                                                        "| mod:", int(dx/scale), &
                                                        "| grad:", int(dg/scale) 
        endif

        
end subroutine  solver_report


end module 
