module shaping_solver_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 shaping_solver( m,d, Fop,sop, niter,np, eps,Wop&
    &,m0,resd,verb)
    optional :: Wop,m0,resd,verb
    interface  
!-------------------------- begin definitions -----------
      integer function sop(adj,add,m,d)
        real::m(:),d(:)
        logical,intent(in)::adj,add
      end function 
      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 
    end interface 
    real, dimension(:),    intent(in)          :: d, m0
    integer,               intent(in)          :: niter
    logical,               intent(in)          :: verb
    real, dimension(:),    intent(out)         :: m,resd
    real, dimension(:)              ,pointer:: r
    integer                                    :: dp,np
    double precision                          :: gn,gnp,alpha,beta,g0,dg,r0,b0
    real, dimension(:), pointer                ::  gp,sp,sx,sr,tt,p,gx,gr
    integer                                    :: iter, stat
    logical                                    :: forget
    real :: eps

!this makes the implicit assumption that the Wop range is the same as its domain
    allocate(r(size(d)))

    if (present( Wop)) then
      stat=Wop(F,F,-d,r) ! begin initialization --------
    else
      r = -d                                           !Rd = -W  d
    end if 
!-------------------------- begin iterations ------------

    
    dg=0;g0=0;b0=0;gnp=0
    allocate(gp(np),p(np),gx(size(m)),gr(size(d)))
    allocate(sp(np),sx(size(m)),sr(size(d)))
    if(present(wop)) allocate(tt(size(d)))
    p=0.
    m=0

call srite("resd.H",r,size(r)*4)
    do iter = 1,niter  
!write(0,*) "XXxx",size(r)
      gp=eps*eps*p
      gx=-eps*eps*m
      if (present(Wop)) then
        call chain0(Wop,Fop,T,T,gx,r,tt)
      else
        stat = Fop(T,T,gx,r)   !g  = (WF)'Rd
      end if 
      stat=sop(.true.,.true.,gp,gx)
      stat=sop(.false.,.false.,gp,gx)
      if (present(Wop)) then
        call chain0(Wop,Fop,F,F,gx,gr,tt)
      else
        stat = Fop(F,F,gx,gr)   !g  = (WF)'Rd
      end if 
      gn=norm_d(gp)

      if(iter==1) then
        g0=gn 
        b0=abs(gn+eps*eps*(norm_d(gr)-norm_d(gx)))
        sp=gp
        sx=gx
        sr=gr
      else
        alpha=gn/gnp
        dg=gn/g0
        sp=gp+alpha*sp
        sx=gx+alpha*sx
        sr=gr+alpha*sr
      end if
      beta=norm_d(sr)+eps*eps*(norm_d(sp)-norm_d(sx))
      alpha=-gn/beta
      write(0,*) "iteration",iter,norm_d(r),dg,alpha,beta
      p=p+alpha*sp
      m=m+alpha*sx
      r=r+alpha*sr
      gnp=gn
    end do
!call srite("resd.H",r,size(r)*4)
    if (present( resd)) then
      resd = r
    end if
    deallocate(r,gp,p,gx,sp,sx,sr,gr)
    if(present(wop)) deallocate(tt)
  end subroutine 
  double precision function norm_d(vec)
    real :: vec(:)
    double precision v
    integer :: i
    norm_d=0
    do i=1,size(vec)
     v=vec(i)
      norm_d=norm_d +v*v
    end do
  end function
end module 
