module solver_H_mod ! 0 = G m - d
  implicit none
  contains
  subroutine solver_tiny( m,d, Fop, niter)
    interface
      integer function Fop(adj,add,m,d)
        real,dimension(:) :: m,d
        logical,intent(in):: adj,add
      end function
    end interface
    real, dimension(:),    intent(in)  :: d         ! data
    real, dimension(:),    intent(out) :: m         ! model
    integer,               intent(in)  :: niter     ! number of iterations
    integer                            :: iter      ! iteration number
    real, dimension(size( m))          :: p         ! p
    real, dimension(size( m))          :: ppp         ! ppp
    real, dimension(size( m))          :: pppp         ! ppp
    real, dimension(size( m))          :: qqqq         ! ppp
    real, dimension(size( m))          :: c         ! c
    real, dimension(size( m))          :: r         ! residual
    real, dimension(size( m))          :: q         ! q
    integer                            :: stat ,tt,tt1,tt2     ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    c   =-d                                       ! c   =- d
    m   =0.                                       ! x0=0
    r   =c                                        ! r=c for x0=0
    write(0,*)"          iter       dat res     p*q       "
    do iter = 1,niter     
      p =-r+beta*p                             
      stat = Fop(.false.,.false.,p,q)          ! q =  G  p data space 
      pp=dp_ale(p,q)                           ! model times data space
      rr=dp_ale(r,r)                           ! data space residual
      if(rr < .01) then
        write(0,*) "residual small"
        return 
      end if
      alpha=rr/pp
      m=m+alpha*p                              ! model update

!      stat = Fop(.false.,.false.,m,ppp)        ! r = r - G m
!      r=c+ppp
      r=r+alpha*q                              ! data space residual
      rrr=dp_ale(r,r)
      beta=rrr/rr
      if (beta > 1.) then
        write(0,*) "residual negative"
        return
      endif
      pppp=alpha*p
      qqqq=alpha*q
      tt=dp_ale(pppp,pppp)
      tt1=dp_ale(p,p)
      tt2=dp_ale(m,m)

      write(0,*)"",iter,int(rr),int(pp) 
   end do
  end subroutine

  subroutine solver_precond( m,d, Fop,Pop, niter)
    interface
      integer function Fop(adj,add,m,d)
        real,dimension(:) :: m,d
        logical,intent(in):: adj,add
      end function
      integer function Pop(adj,add,m,d)
        real,dimension(:) :: m,d
        logical,intent(in):: adj,add
      end function
    end interface
    real, dimension(:),    intent(in)  :: d         ! data
    real, dimension(:),    intent(out) :: m         ! model
    integer,               intent(in)  :: niter     ! number of iterations
    integer                            :: iter      ! iteration number
    real, dimension(size( m))          :: p         ! p
    real, dimension(size( m))          :: ppp       ! ppp
    real, dimension(size( m))          :: c         ! c
    real, dimension(size( m))          :: r         ! residual
    real, dimension(size( m))          :: wr        ! residual
    real, dimension(size( m))          :: q         ! q
    integer                            :: stat      ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real			       :: eps

    eps=10000.
    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    c   =-d                                       ! c   =- d
    m   =0.                                       ! x0=0
    r   =c                                        ! r=c for x0=0
    wr  =0.
    ppp=0.
    write(0,*)"          iter       dat res     p*q"
    do iter = 1,niter     
      stat = Pop(.true.,.false.,wr,eps*r)
      p =-wr+beta*p      
      r=c+ppp
      rr=dp_ale(r,wr)
      stat = Fop(.false.,.false.,p,q)          ! q =  G  p 
      pp=dp_ale(p,q)
      alpha=rr/pp
      m=m+alpha*p
      stat = Fop(.false.,.false.,m,ppp)          ! r = r - G m
      r=c+ppp
!      r=r+alpha*q
      stat = Pop(.true.,.false.,wr,eps*r)
      rrr=dp_ale(r,wr)
      beta=rrr/rr
      if (beta > 1.) then
        write(0,*) "residual negative"
        return
      endif
     write(0,*)"",iter,int(rr),int(pp)
   end do
  end subroutine


  real function dp_ale (a,b)
       real, dimension(:) :: a,b 
       double precision  :: aa,bb,dp
       integer           :: i

      dp=0.
      do i=1,size(a)
        aa=a(i);bb=b(i)
        dp =dp+aa*bb
      end do
      dp_ale=dp
  end function dp_ale

end module
