module solver_H_mod ! 0 = G m - d
  implicit none
  contains
  subroutine solver_tiny( m,d, Fop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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
    real,    dimension (:,:)           :: m_mov
    real,    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
    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)                           ! 
      rr=dp_ale(r,r)                           ! 
!      if (rr > rrr) then
!        alpha=.01*rr/pp
!      else
        alpha=rr/pp
!      endif
!      alpha=min(.01,max(-.01,alpha))
      m=m+alpha*p                              ! model update
      r=r+alpha*q                              ! data space residual
      rrr=dp_ale(r,r)
      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

  subroutine solver_tiny_w( m,d, Fop, niter,r1,m_mov,r_mov)
  optional ::m_mov,r_mov
    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))          :: r1         ! residual
    real, dimension(size( m))          :: q         ! q
    integer                            :: stat ,tt,tt1,tt2     ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    c   =-d                                       ! c   =- d
    m   =0.                                       ! x0=0
    r   =c*r1                                        ! 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)                           ! 
      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,int(rr),int(pp) 
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do

  end subroutine

  subroutine solver_tiny_square( m,d, Fop,mo, niter,m_mov,r_mov)
  optional ::m_mov,r_mov,mo
    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))          :: mo         ! p
    real, dimension(size( m))          :: pt         ! 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))          :: rt         ! residual
    real, dimension(size( m))          :: q         ! q
    integer                            :: stat ,tt,tt1,tt2     ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    stat= Fop(.true.,.false.,c,-d)                ! c   =- G'd
    m   =0.                                       ! x0=0
    r   =c                                        ! r=c for x0=0
    pt=0
    if (present( mo)) then 
      m = mo
      stat = Fop(.false.,.false.,m,rt   )
      stat = Fop(.true. ,.true. ,r,rt   )
    endif
    write(0,*)"          iter       dat res     p*q       "
    do iter = 1,niter     
      p =-r+beta*p                             
      stat = Fop(.false.,.false.,p,pt)         ! pt =  G  p  
      stat = Fop(.true. ,.false.,q,pt)         ! q  =  G' pt = G'G p  
      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,int(rr),int(pp) 
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do

  end subroutine


  subroutine solver_tiny_recomp( m,d, Fop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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))          :: 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,pppp,qqqq
    real,    dimension (:,:)           :: m_mov
    real,    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
    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 
      qqqq=dp_ale(q,q)
      pppp=dp_ale(p,p)
      pp=dp_ale(p,q)                           ! 
      rr=dp_ale(r,r)                           ! 
      alpha=rr/pp
      m=m+alpha*p                              ! model update
      stat = Fop(.false.,.false.,m,ppp)        ! r = c - G m
      r=c+ppp                             ! data space residual
      rrr=dp_ale(r,r)
      beta=rrr/rr
      write(0,*)"",iter,int(rr),int(pp)
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do

  end subroutine

  subroutine solver_conj_dir( m,d, Fop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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
    real,    dimension (:,:)           :: m_mov
    real,    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
    write(0,*)"          iter       dat res     p*q       "
    do iter = 1,niter     
      p =-r                             
      stat = Fop(.false.,.false.,p,q)          ! q =  G  p data space 
      pp=dp_ale(p,q)                           ! pTGp 
      rr=dp_ale(r,r)                           ! rTr 
      alpha=rr/pp
      m=m+alpha*p                              ! model update
      r=r+alpha*q                              ! data space residual

      write(0,*)"",iter,int(rr),int(pp), alpha
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do

!      stat = Fop(.false.,.false.,m,ppp)        ! r = r - G m
!      r=c+ppp

  end subroutine


  subroutine solver_precond( m,d, Fop,Pop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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,    dimension (:,:)           :: m_mov
    real,    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.
    ppp=0.
    write(0,*)"          iter       dat res     p*q"
    do iter = 1,niter     
      stat = Pop(.true.,.false.,wr,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
!      r=r+alpha*q
      stat = Fop(.false.,.false.,m,ppp)          ! r = r - G m
      r=c+ppp
      stat = Pop(.true.,.false.,wr,r)
      rrr=dp_ale(r,wr)
      beta=rrr/rr

     write(0,*)"",iter,int(rr),int(pp)
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do
  end subroutine

   subroutine solver_precond_conj_dir( m,d, Fop,Pop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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,    dimension (:,:)           :: m_mov
    real,    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.
    ppp=0.
    write(0,*)"          iter       dat res     p*q"
    do iter = 1,niter     
      stat = Pop(.true.,.false.,wr,r)
      p =-wr      
      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,r)
!      rrr=dp_ale(r,wr)
!      beta=rrr/rr
!      if (beta > 1.) then
!        write(0,*) "residual negative"
!        return
!      endif
     write(0,*)"",iter,int(1000*rr),int(1000*pp)
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do
  end subroutine

  subroutine solver_precond1( m,d, Fop,Pop,mo, niter,m_mov,r_mov,beta_max)
  optional ::m_mov,r_mov,mo
    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))          :: 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
    integer                            :: stat      ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov
    real                               :: beta_max
    write(0,*)"beta_max",beta_max
    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
      if (beta > beta_max) then
        beta=0.01*beta
      endif

     write(0,*)"",iter,int(rrr),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

  subroutine solver_precond1_square( m,d, Fop,Pop,mo,niter,m_mov,r_mov)
  optional ::m_mov,r_mov,mo
    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))          :: mo        ! initial model
    real, dimension(size( m))          :: c         ! c
    real, dimension(size( m))          :: r         ! residual
    real, dimension(size( m))          :: rt         ! residual
    real, dimension(size( m))          :: wr        ! residual
    real, dimension(size( m))          :: q         ! q
    real, dimension(size( m))          :: pt         ! q
    integer                            :: stat      ! status flag
    double precision                   :: beta,alpha,rr,rrr,pp
    real,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    stat= Fop(.true.,.false.,c,-d)                ! c   =- G'd
    m   =0.                                       ! x0=0
    r   =c                                        ! r=c for x0=0
    wr  =0.
    rt  =0.
    ppp=0.
    if (present( mo)) then 
      m = mo
      stat = Fop(.false.,.false.,m,rt   )
      stat = Fop(.true. ,.true. ,r,rt   )
    endif
    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,pt)        ! pt =  G  p  
      stat = Fop(.true.,.false.,q,pt)          ! q =  G' pt =G'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)
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do
  end subroutine



  subroutine solver_precond2( m,d, Fop,Pop,mo, niter,r1,m_mov,r_mov,beta_max)
  optional ::m_mov,r_mov,mo
    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))          :: 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))          :: r1         ! 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,    dimension (:,:)           :: m_mov
    real,    dimension (:,:)           :: r_mov
    real                               :: beta_max

    beta=0.                                       ! beta=0
    p   =0.                                       ! p   =0
    c   =-d                                       ! c   =- d
    m   =0.                                       ! x0=0
    r   =c*r1                                        ! r=c for x0=0
    if (present( mo)) then 
      m = mo
      stat = Fop(.false.,.true.,m,r   )
    endif
    wr  =0.
    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) preconditioning
      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
      if (beta > beta_max) then
        beta=0.01*beta
      endif

     write(0,*)"",iter,int(rrr),int(pp)
      m_mov (:,iter) = m(1:size(m_mov,1))
      r_mov (:,iter) = r(1:size(r_mov,1))
   end do
  end subroutine

  subroutine solver_precond_conj_dir1( m,d, Fop,Pop, niter,m_mov,r_mov)
  optional ::m_mov,r_mov
    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,    dimension (:,:)           :: m_mov
    real,    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.
    ppp=0.
    write(0,*)"          iter       dat res     p*q"
    do iter = 1,niter     
      
      stat = Pop(.true.,.false.,wr,r)          ! wr(k)= Wr(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

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