SUBROUTINE PrintShiftBlocking(Iter)
use util_mod, only: get_free_unit
INTEGER :: i, NoBlocks, Iter, NoBlockSizes, NoContrib, iunit
real(dp) :: MeanShift, MeanShiftSqrd, StandardDev, Error, ErrorinError
!First find out how many blocks would have been formed with the number of iterations actually performed.
NoContrib = 0
NoContrib = ((Iter - StartShiftBlockIter) / StepsSft) + 1
NoBlockSizes = 0
NoBlockSizes = FLOOR((LOG10(REAL(NoContrib - 1, dp))) / (LOG10(2.0_dp)))
iunit = get_free_unit()
open(iunit, file='SHIFTBLOCKINGANALYSIS', status='unknown')
write(iunit, '(I4,A,I4)') NoBlockSizes, ' blocks were formed with sizes from 1 to ', (2**(NoBlockSizes))
write(iunit, '(3A16,5A20)') '1.Block No.', '2.Block Size ', '3.No. of Blocks', '4.Mean Shift', &
'5.Mean Shift^2', '6.SD', '7.Error', '8.ErrorinError'
do i = 0, NoBlockSizes - 1
! First need to find out how many blocks of this particular size contributed to the final sum in BlockSum.
! NoContrib is the total number of contributions to the blocking throughout the simulation.
NoBlocks = 0
NoBlocks = FLOOR(REAL(NoContrib / (2**i), dp))
! This finds the lowest integer multiple of 2**i (the block size).
MeanShift = ShiftBlockSum(i) / REAL(NoBlocks, dp)
MeanShiftSqrd = ShiftBlockSqrdSum(i) / REAL(NoBlocks, dp)
StandardDev = SQRT(MeanShiftSqrd - (MeanShift**2))
IF (abs(StandardDev) < 1.0e-12_dp) THEN
Error = 0.0_dp
ErrorinError = 0.0_dp
ELSE
Error = StandardDev / SQRT(REAL(NoBlocks - 1, dp))
ErrorinError = Error / SQRT(2.0_dp * (REAL(NoBlocks - 1, dp)))
end if
!This is from the blocking paper, and indicates the error in the blocking error, due to the limited number of blocks available.
write(iunit, '(3I16,5F20.10)') i, (2**i), NoBlocks, MeanShift, MeanShiftSqrd, StandardDev, Error, ErrorinError
end do
close(iunit)
END SUBROUTINE PrintShiftBlocking