module row_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
  

  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 row_findmatch(data1,data2,targt,match1,match2,&
        & niter,mu,dtest,wght,nmem_in,nfreq_in)
    real, dimension(:), pointer     :: data1,data2,targt,wght
    integer                         :: niter,nmem_in,nfreq_in
	real                            :: mu
    type (filter)                   :: match1,match2
    logical                         :: init,dtest

    optional                        :: wght,nmem_in,nfreq_in

    integer                         :: stat,i,nmem,nfreq
    real                            :: sd,rnd

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

    if (present(nmem_in)) then
       call row_cgfindmatch(data1,data2,targt,match1,match2,niter,mu,dtest, &
            &nmem,nfreq)
    else
       if (present(wght)) then 
         call row_cgfindmatch(data1,data2,targt,match1,match2,niter,mu,dtest,&
           & weight=wght)
       else
         call row_cgfindmatch(data1,data2,targt,match1,match2,niter,mu,dtest)
       end if
    end if
  end subroutine row_findmatch

  !! Find match with conjugate gradients (slow)
  subroutine row_cgfindmatch(data1,data2,targt,match1,match2,niter,mu,dtest1,&
       & nmem,nfreq,weight)
    real, dimension(:), pointer         :: data1,data2,targt,weight
    integer                             :: nmem,nfreq
    type (filter)                       :: match1,match2
    integer                             :: niter
	real                                :: mu
	logical                             :: dtest1

    optional                            :: nmem,nfreq,weight

    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'

    ncoeffs1=size(match1%flt)
    ncoeffs2=size(match2%flt)
    allocate(coeffs(ncoeffs1+ncoeffs2))
    coeffs=0.

    call row_nchconest_init(data1,data2,match1,match2,ncoeffs1,ncoeffs2,mu)

if (dtest1) then
    write(0,*)'*****************************************'
    write(0,*)'Dot product test for row_nchconest_lop'
    call dot_est(row_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(row_nchconest_lop,cgstep,match%flt,targt,niter,wght=l1_priv, &
       call solver(row_nchconest_lop,cgstep,coeffs,targt,niter,x0=coeffs,&
            &     nmem=nmem,nfreq=nfreq)
    else

       if (present(weight)) then
            call solver(row_nchconest_lop,cgstep,coeffs,targt,niter, &
                 & x0=coeffs,wt=weight)
       else
            sverb=.true.
            call solver(row_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 row_cgfindmatch


end module row_findmatch_mod
