module array_findnsmatch_mod
  
  use nhelix
  use row_ncnhconest
  use helicon
  use ncnhelicon
  
  use cgstep_mod
  use solver2_mod
  use verbosity_mod
  use dotest
  use array_ncnhconest_mod
  
  implicit none
  
contains

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Non-stationary matched filtering !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !! Finds non-stationary matched-filter
  subroutine array_findnsmatch(data1,data2,targt,match1,match2,&
        & reg,mode,niter,eps,mu,eps1,eps2,eps3,dtest)
    real                                  :: eps
    integer                               :: niter    
    logical                               :: dtest
    real, dimension(:), pointer           :: data1,data2,targt
    type (nfilter)                        :: match1,match2
    type ( filter)                        :: reg
    integer                               :: mode
    real                                  :: mu
    real                                  :: eps1,eps2,eps3

    real, dimension(:), pointer :: coeffs
    integer :: stat,nhmax1,nhmax2,i
    logical :: sverb
    integer    :: ncoeffs1,ncoeffs2,nl

    sverb=(verb>2)
    nhmax1=size(match1%hlx(1)%lag)  ! assumes all filters are same length
    nhmax2=size(match2%hlx(1)%lag)  ! assumes all filters are same length

    write(0,*) '###epsilon=',eps,'niter=',niter,'mu=',mu
    write(0,*) '###eps1=',eps,'eps2=',eps2,'eps3=',eps3

    ncoeffs1=nhmax1*size(match1%hlx)
    ncoeffs2=nhmax2*size(match2%hlx)

    allocate (coeffs(ncoeffs1+ncoeffs2))
    coeffs=0.

    ! Regularized solver
    if (verb>2) write(0,*) '  using regularized conjugate gradients'

    call helicon_init(reg)   ! initialize regularizer

    call array_ncnhconest_init(data1,data2,match1,match2,nhmax1,nhmax2,&
        & ncoeffs1,ncoeffs2,mu,eps1,eps2,eps3) ! initialize operators

if (dtest) then
    write(0,*)'*****************************************'
	write(0,*)'size(data1)=',size(data1)
	write(0,*)'size(data2)=',size(data2)
	write(0,*)'size(target)=',size(targt)
    write(0,*)'Dot product test for array_ncnhconest_lop'
    call dot_est(array_ncnhconest_lop, size(coeffs), size(targt))
    write(0,*)'*****************************************'
end if

    ! solve problem
    call solver_reg(array_ncnhconest_lop,cgstep,helicon_lop,size(coeffs), &
            & coeffs,targt,niter,eps,x0=coeffs,verb=sverb)

    !! Copy coefficients into NS filter type
    do i=1,size(match1%hlx)
       match1%hlx(i)%flt(1:nhmax1)=coeffs(1+(i-1)*nhmax1:i*nhmax1)
    end do
    nl=size(match1%hlx)*nhmax1 !! last coeff from match1
    do i=1,size(match2%hlx)
       match2%hlx(i)%flt(1:nhmax2)=coeffs(nl+1+(i-1)*nhmax2:i*nhmax2)
    end do

    deallocate(coeffs)
    call cgstep_close()

  end subroutine array_findnsmatch

end module array_findnsmatch_mod
