# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module wilson_bob
! Wilson's factorization
  use helicon
  use polydiv
  implicit none
  integer, private :: n
  real, dimension( :,:), allocatable, private :: auto, bb, cc, b, c
  contains
  subroutine wilson_bob_init( nmax,nthread)
    integer, intent (in) :: nmax,nthread
    n = nmax
    allocate ( auto( 2*n-1,nthread), bb( 2*n-1,nthread), cc( 2*n-1,nthread&
      &), b(n,nthread), c(n,nthread))
  end subroutine
  logical function wilson_bob_factor( niter, s0, ss, a0, aa, verb)
    integer, intent( in) :: niter ! Newton iterations
    real, intent( in) :: s0 ! autocorrelation zero lag
    type( filter), intent( in) :: ss
    ! autocorrelation, other lags
    real, intent( out) :: a0 ! factor, zero lag
    type( filter) :: aa ! factor, other lags
    logical, intent( in) :: verb
    optional :: verb
    real :: eps(niter),eold
    integer :: i, stat,Ithread
    ithread=1
    auto(:,ithread) = 0.
    auto( n,ithread) = s0
    b(1,:) =1. ! initialize
    auto( n+ss%lag,ithread) = ss%flt
    ! autocorrelation
    auto( n-ss%lag,ithread) = ss%flt
    ! symmetrize input auto.
    call helicon_init( aa) ! multiply polynoms
    call polydiv_init( 2*n-1, aa) ! divide polynoms
    eold=1.e10
    wilson_bob_factor=.true.
    do i = 1, niter
      stat= polydiv_lop(.false.,.false., auto(:,ithread), bb(:,ithread&
        &)) ! bb = S/A
      stat= polydiv_lop( .true.,.false., cc(:,ithread), bb(:,ithread&
        &)) ! cc = S/(AA')
      b( 2:n,ithread) = 0.5*( cc( n+1:2*n-1,ithread ) +cc( n-1:1 :&
        &-1,ithread)) / cc( n,ithread) ! b = plusside(1+cc)
      eps(i) = maxval( abs( b(2:n,ithread)))
! "L1 norm"
! if (verb) write (0,*) i, eps(i)
      stat= helicon_lop( .false., .false., b(:,ithread), c(:,ithread))
! c = A b
      aa%flt = c( 1+aa%lag,ithread)
      ! put on helix
      if ( eps(i) < epsilon( a0)) then
        exit ! convergence
      end if
      if (eold*10. < eps(i)) then
        wilson_bob_factor=.false.
        return
      end if
      if (i>4) then
        if (eps(i)*1.1 >eps(i-3)) then
          return
        end if
      end if
    end do
    a0 = sqrt( cc( n,ithread))
  end function
  subroutine wilson_bob_close ()
    deallocate( auto, bb, cc, b, c)
    call polydiv_close()
  end subroutine
end module
