module iwilson { # Wilson's factorization - inverse square root use helicon use polydiv integer, private :: n real, dimension( :), allocatable, private :: auto, bb, cc, b, c contains subroutine iwilson_init( nmax) { integer, intent (in) :: nmax; n = nmax allocate ( auto( 2*n-1), bb( 2*n-1), cc( 2*n-1), b(n), c(n)) } subroutine iwilson_factor( niter, s0, ss, a0, aa) { 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 real :: eps integer :: i, stat real, dimension (size (aa%flt)) :: rand # auto = 0.; auto( n) = s0; b( 1) =1. # initialize # auto( n+ss%lag) = ss%flt # autocorrelation # auto( n-ss%lag) = ss%flt # symmetrize input auto. call polydiv_init(2*n-1, aa) auto = 0.; auto( n) = 1.; b( 1) = 1. stat= polydiv_lop(.false.,.false., auto, bb) # bb = S*A stat= polydiv_lop( .true.,.false., cc, bb) # cc = S*AA' auto = cc call random_number (rand) aa%flt = aa%flt*(1.+0.25*rand) call helicon_init( aa) # multiply polynoms do i = 1, niter { stat= helicon_lop(.false.,.false., auto, bb) # bb = S*A stat= helicon_lop( .true.,.false., cc, bb) # cc = S*AA' b( 2:n) = -0.5*( cc( n+1:2*n-1 ) + cc( n-1:1 :-1)) / cc( n) # b = plusside(3-cc) eps = maxval( abs( b(2:n))); write (0,*) i, eps # "L1 norm" if( eps < epsilon( a0)) break # convergence stat= helicon_lop( .false., .false., b, c) # c = A b aa%flt = c( 1+aa%lag) # put on helix } a0 = sqrt( cc( n)) } subroutine iwilson_close () { deallocate( auto, bb, cc, b, c); } }