module array_findmatch_mod
  
  use helix
  use helicon
  use polydiv
  use row_nchconest
  use nchelicon
  use findmatch_mod
  
  use nhelix
  use ncnhconest
  use ncnhelicon
  
  use cgstep_mod
  use solver2_mod
  use verbosity_mod
  use array_nchconest_mod
  

  implicit none
  
  integer, private                     :: nf,ndata
  type (filter), private               :: acorr
  real, private, dimension(:), pointer :: xcorr,space,matrix

contains

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Find a single matched-filter     !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! Minimizes residual:
  ! wght ( data * match - targt ) = 0

  subroutine array_findmatch(data1,data2,targt,match1,match2,&
        & niter,mu,eps1,eps2,eps3,dtest,wght,nmem_in,nfreq_in)

    !! subroutine arguments
    real, dimension(:), pointer     :: data1,data2,targt,wght
    integer                         :: niter,nmem_in,nfreq_in
    logical                         :: dtest
	real                            :: mu
	real                            :: eps1,eps2,eps3
    type (filter)                   :: match1,match2
    integer                         :: nmem,nfreq
    optional                        :: wght,nmem_in,nfreq_in

    !! internal variables
    logical                         :: i,init,stat
    real                            :: sd,rnd

    if (present(nfreq_in)) then 
       nmem=nmem_in
       nfreq=nfreq_in
    end if

    !! call conjugate gradient with array operator 
    if (present(nmem_in)) then
       call array_cgfindmatch(data1,data2,targt,match1,match2,niter,mu,&
            & eps1,eps2,eps3,dtest,nmem,nfreq)
    else
        if (present(wght)) then 
            call array_cgfindmatch(data1,data2,targt,match1,match2,niter,mu, &
                & eps1,eps2,eps3,dtest,weight=wght)
        else
            call array_cgfindmatch(data1,data2,targt,match1,match2,niter,mu, &
                & eps1,eps2,eps3,dtest)
        end if
    end if

  end subroutine array_findmatch



  !! Find match with conjugate gradients (slow)
  subroutine array_cgfindmatch(data1,data2,targt,match1,match2,niter, &
        & mu,eps1,eps2,eps3,dtest1,nmem,nfreq,weight)

    !! subroutine arguments
    real, dimension(:), pointer         :: data1,data2,targt,weight
    integer                             :: nmem,nfreq
    type (filter)                       :: match1,match2
    integer                             :: niter
	real                                :: mu
    real                                :: eps1,eps2,eps3
	logical                             :: dtest1
    optional                            :: nmem,nfreq,weight

    !! internal variables
    integer                   :: stat
    integer                   :: ncoeffs1
    integer                   :: ncoeffs2
	real,dimension(:),pointer :: coeffs
	logical                   :: sverb

    if (verb>1) write(0,*) ' Using CG to solve system- there may be better way'

    !! allocate space for filter coefficients
    ncoeffs1=size(match1%flt)
    ncoeffs2=size(match2%flt)
    allocate(coeffs(ncoeffs1+ncoeffs2))
    coeffs=0.

    !! initialize helical convolution operator
    call array_nchconest_init(data1,data2,match1,match2,ncoeffs1,ncoeffs2, &
         mu,eps1,eps2,eps3)

    !! if requested, compute the dot product of the helical conv operator
if (dtest1) then
    write(0,*)'*****************************************'
    write(0,*)'Dot product test for array_nchconest_lop'
    call dot_est(array_nchconest_lop, size(coeffs), size(targt))
    write(0,*)'*****************************************'
end if

    if (present(nfreq)) then
       write(0,*) 'Using L1 norm with nonlinear CG'
       write(0,*) 'nmem',nmem
       write(0,*) 'nfreq',nfreq
       call solver(array_nchconest_lop,cgstep,coeffs,targt,niter,x0=coeffs,&
            &     nmem=nmem,nfreq=nfreq)
    else

       if (present(weight)) then
            call solver(array_nchconest_lop,cgstep,coeffs,targt,niter,x0=coeffs,wt=weight)
       else
            sverb=.true.
            call solver(array_nchconest_lop,cgstep,coeffs,targt,niter,verb=sverb)
       end if
    end if

    !! reassign coeffcients to the corresponding filters
    match1%flt(:)=coeffs(1:ncoeffs1)
    match2%flt(:)=coeffs(ncoeffs1+1:ncoeffs1+ncoeffs2)
    deallocate(coeffs)

    call cgstep_close()

  end subroutine array_cgfindmatch


end module array_findmatch_mod
