module solver_dp ! 0 = G m - d

  use DataModule, only : dp

  implicit none
  contains

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

    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.
    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 
      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,int(rr),int(pp)
     write(0,*)" alpha    beta  ",alpha, beta

      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do
  end subroutine

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

end module
