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

!--init
  subroutine hconest_init ( x_in, aa_in)
    real, dimension (:), pointer	:: x_in
    type( filter), intent (in)	:: aa_in
    x => x_in
    aa =  aa_in
    nFilters = product(aa%dataDims)/aa%filterJ; write(0,*)nFilters
    na = size(aa%coeffs)/nFilters; write(0,*)na
    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)
    write(0,*)size(y)*na
    do ia = 1, na
!     write(0,*)size(aa%lags),size(a)
!     write(0,*)   na*((iy-1)/aa%filterJ) +ia
!     write(0,*) aa%lags(  na*((iy-1)/aa%filterJ) +ia)
!     do iy = 1  + aa%lags(  na*((iy-1)/aa%filterJ) +ia), size( y)
!      do iy = 1  + aa%lags( ia), size( y)
      do iy = 1, size( y)
        if ( .not. aa%edges( iy)) then
          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
	write(0,*) adj,"minval",minval(x),minval(x)
	write(0,*) adj,"maxval",maxval(x),maxval(y)
  end function hconest_lop

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