module hconest
use adj_mod
use helix
  implicit none
    real, dimension (:), pointer , private :: x
    type( filter)           , private :: aa
    integer			 , private :: na, nFilters
		logical,save,private :: edgeson
contains

!--init
  subroutine hconest_init ( x_in, aa_in,edges)
    real, dimension (:), pointer	:: x_in
    type( filter), intent (in)	:: aa_in
		logical, optional :: edges

		if(present(edges)) then
      edgeson=edges
		else 
			edgeson=.false.
		end if
			
    x => x_in
    aa =  aa_in
    nFilters = product(aa%dataDims)/aa%filterJ; 
    na = size(aa%coeffs)/nFilters; 
    if( mod(product(aa%dataDims),aa%filterJ) .ne. 0) call erexit('hconest_lop: its my first day')
    if( mod(size(aa%coeffs),nFilters) .ne. 0) call erexit('hconest_lop: its my first day')
  end subroutine hconest_init

!--operator
  function hconest_lop (adj, add, a, y) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: a, y
    integer  ia, ix, iy


    call adjnull (adj,add,a,y)
    do ia = 1, na
      do iy = 1, size( y)
				if( edgeson) then
        	if (aa%edges( iy)) cycle  
        end if
        ix = iy - aa%lags(  na*((iy-1)/aa%filterJ) +ia)
	if ( ix < 1 ) then
		cycle
	end if
        if ( adj) then
          a( na*((iy-1)/aa%filterJ) + ia) = &
			a( na*((iy-1)/aa%filterJ) + ia) +  y( iy) * x( ix)
	!	a( ia) =a( ia) +  y( iy) * x( ix)
        else
          y( iy) =y( iy) +  a( na*((iy-1)/aa%filterJ) + ia) * x( ix)
	!	y( iy) =y( iy) +  a( ia) * x( ix)
        end if 
      end do 
    end do 
    stat = 0
  end function hconest_lop

!--close
  subroutine hconest_close ()
  end subroutine hconest_close
end module hconest
