module row_findnsmatch_mod
  
  use nhelix
!  use ncnhconest1
!  use ncnhconest2
    use row_ncnhconest
  use helicon
  use ncnhelicon
  use findmatch_mod
  
  use cgstep_mod
  use solver2_mod
  use verbosity_mod
    use dotest
  
  implicit none
  
contains

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

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

    real, dimension(:), pointer :: coeffs
    integer :: stat,nhmax1,nhmax2,i,iouter
    logical :: sverb
    integer	:: ncoeffs1,ncoeffs2,nl
    real, dimension(:), pointer :: tmp1,tmp2    
    real :: eps1,mu1

    if (nouter>1) then
	    allocate (tmp1(size(data1)))
	    allocate (tmp2(size(data2)))
	end if

    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

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

    allocate (coeffs(ncoeffs1+ncoeffs2))

    ! Regularized solver
    if (verb>2) write(0,*) '  using regularized conjugate gradients'
    call helicon_init(reg)                       ! initialize regularizer

    coeffs=0.
    do iouter=1,nouter
		tmp1=0.
		tmp2=0.
	    eps1=eps
		!mu=1.0/mu

write(0,*)'iouter=',iouter,'mu=',mu
        call row_ncnhconest_init(data1,data2,match1,match2,nhmax1,nhmax2,&
    	     & ncoeffs1,ncoeffs2,mu) ! initialize operators

if (dtest .and. iouter==1) then
    write(0,*)'*****************************************'
    write(0,*)'Dot product test for row_ncnhconest_lop'
    call dot_est(row_ncnhconest_lop, size(coeffs), size(targt))
    write(0,*)'*****************************************'
end if

        ! solve problem
        !if (iouter==1) eps1=0.
        call solver_reg(row_ncnhconest_lop,cgstep,helicon_lop,size(coeffs), &
            & coeffs,targt,niter,eps1,x0=coeffs,verb=sverb)
write(0,*) '***PASSEI***'
write(0,*)'iouter=',iouter
        !! Copy coefficients into NS filter type
        do i=1,size(match1%hlx)
!			write(0,*)'i1=',1+(i-1)*nhmax1,'i2=',i*nhmax1
            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)
!			write(0,*)'i1=',nl+1+(i-1)*nhmax2,'i2=',nl+i*nhmax2
            match2%hlx(i)%flt(1:nhmax2)=coeffs(nl+1+(i-1)*nhmax2:nl+i*nhmax2)
        end do
 
		if (iouter<nouter) then
            !! recompute estimate of primaries and multiples
            call applynsmatch(data1,tmp1,match1,.false.)
            call applynsmatch(data2,tmp2,match2,.false.)
		    data2=tmp2
		    data1=mu*tmp1
		    !data1=targt-tmp2
		    !data2=targt-mu*tmp1
			coeffs=0
        end if
write(0,*)'Despues de recomputar primarios y multiples'
	
    end do

    deallocate(coeffs)
    if (nouter>1) deallocate(tmp1,tmp2)
    call cgstep_close()

  end subroutine row_findnsmatch

end module row_findnsmatch_mod
