module OW_invint
  use lint1
  use igrad1
  use cgstep_mod
  use random_mod
  use solver_reg_mod

  implicit none

  logical, save                          :: multreal = .false.
  integer, save                          :: seed = 2941
  real, save                             :: scale = 1.0
  real,dimension(:),private,save,pointer :: mask
contains

  ! --------------------------------------------------------------------------------
  subroutine init_random( mu, se, sc)
    logical                         :: mu
    integer                         :: se
    real                            :: sc

    seed = se;  scale = sc;  multreal = mu
  end subroutine init_random
  ! --------------------------------------------------------------------------------


  ! --------------------------------------------------------------------------------
  subroutine invint( niter, pair, nd, o1, d1, mm, eps)
    integer,               intent (in)   :: niter, nd
    real,                  intent (in)   :: o1, d1, eps
    real, dimension (:),   pointer       :: aa
    real, dimension (:,:), pointer       :: pair
    real, dimension (:),   intent (out)  :: mm
    optional    :: eps
    real, dimension (:),   pointer       :: ord, dd
    real                                 :: ep, rh
    integer                              :: nreg, nm, stat,i,m 
    real, dimension (:),   pointer       :: xnull, ones, fold

    ord => pair (1 : nd,1)
    dd  => pair (1 : nd,2)

    nm = size (mm)
    call lint1_init (o1, d1, ord)

    allocate( ones(nd), fold(nm),mask(nm)); ones = 1; fold = 0
!    stat = lint1_lop( .false., .false., ones, fold)

    mask=0
    do i=1,nd
       m=int((dd(i)-o1)/d1)+1.5
       mask(m)=1.
    end do

    ep = 1.; if (present (eps)) ep = eps

    allocate( xnull(nm) );  xnull = 0
    if( multreal) then
       call random_init(seed)
       call rand_real_array_1d( xnull)
       xnull=scale*(xnull-0.5)
    end if

!    mask=1
!    where(ord.eq.0.0)
!      mask = 0
!    end where

    call solver_reg(m=mm, d=ord, &
    Fop=lint1_lop, stepper=cgstep, niter=niter, &
    Aop=igrad1_lop, nAop=2, eps=ep , &
    Jop=identity_lop,rm0=xnull, verb=.false.)

    call cgstep_close ()
  end subroutine invint
! --------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function identity_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical                       :: adj,add
    call adjnull( adj, add, model, data)
    if (adj) then
      model = model + data*mask
    else
      data  = data + model*mask
    end if 
    identity_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
end module

