module solver_H_mod ! 0 = G m - d
  implicit none
  contains

  subroutine solver_tiny( m,d, Fop, niter,beta_max,Wop,m_mov,r_mov,r_nor)
  optional ::m_mov,r_mov,r_nor,Wop
    interface
      integer function Fop(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
    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
    real, dimension(size( m))          :: t
    integer                            :: stat ,tt,tt1,tt2     ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    integer, dimension (:)             :: r_nor
    real,    dimension (:,:)           :: r_mov
    real                               :: beta_max
    write(0,*)"beta_max",beta_max

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0

    if (present(Wop)) then
       stat = Wop(.false.,.false.,-d,c)
    else
       c   =-d                                       ! c   =- d
    end if

    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
      t = q 
      if (present(Wop)) stat = Wop(.false.,.false.,t,q)   ! q = W q
      pp=dp_ale(p,q)                           ! 
      rr=dp_ale(r,r)                           ! 
      alpha=rr/pp
      m=m+alpha*p                              ! model update
      r=r+alpha*q                              ! data space residual
      rrr=dp_ale(r,r)
      beta=rrr/rr

      write(0,*)"",iter,rr,pp,beta

      if (beta > beta_max) then
        beta=beta_max
      endif 

      if(present(m_mov)) m_mov (:,iter) = m(1:size(m_mov,1))
      if(present(r_mov)) r_mov (:,iter) = r(1:size(r_mov,1))
      if(present(r_nor)) r_nor (iter)   = int(rrr)
   end do

  end subroutine

  subroutine solver_precond1( m,d, Fop,Pop,niter,beta_max,Wop,mo,m_mov,r_mov,r_nor)
  optional ::Wop,mo,m_mov,r_mov,r_nor
    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
      integer function Wop(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))          :: mo        ! initial model
    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
    real, dimension(size( m))          :: t
    integer                            :: stat      ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov
    integer, dimension (:)             :: r_nor
    real                               :: beta_max
    write(0,*)"beta_max",beta_max
    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0

    if (present(Wop)) then
       stat = Wop(.false.,.false.,-d,c)
    else
       c   =-d                                       ! c   =- d
    end if

    m   =0.                                       ! x0=0
    r   =c                                        ! r=c for x0=0
    wr  =0.
    if (present( mo)) then 
      m = mo
      stat = Fop(.false.,.true.,m,r   )
    endif

    ppp=0.
    write(0,*)"          iter       dat res     p*q"
    do iter = 1,niter     
      stat = Pop(.true.,.false.,wr,r)          ! wr(k)= W^(-1)r(k)
      p =-wr+beta*p                            !  p(k)=-wr(k)+beta(k-1)*p(k-1)
      rr=dp_ale(r,wr)
      stat = Fop(.false.,.false.,p,q)          ! q =  G  p
      t = q
      if (present(Wop)) stat = Wop(.false.,.false.,t,q)
      pp=dp_ale(p,q)
      alpha=rr/pp
      m=m+alpha*p
      r=r+alpha*q
      stat = Pop(.true.,.false.,wr,r)
      rrr=dp_ale(r,wr)
      beta=rrr/rr

      write(0,*)"",iter,rr,pp,beta

      if (beta > beta_max) then
        beta=beta_max
      endif

      if(present(m_mov)) m_mov (:,iter) = m(1:size(m_mov,1))
      if(present(r_mov)) r_mov (:,iter) = r(1:size(r_mov,1))
      if(present(r_nor)) r_nor (iter)   = int(rrr)
   end do

  end subroutine

  real function dp_ale (a,b)
       real, dimension(:) :: a,b 
       double precision  :: dp,aa,bb
       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
