automatic_reblocking_analysis Subroutine

    subroutine automatic_reblocking_analysis(this, blocklength, corrlength, filename, tPrint, iValue)
        ! Performs automatic reblocking (Flyvbjerg and Petersen)
        ! plus some crude selection scheme
        ! General routine, does not require global data
        ! iValue indicates the value which is currently being blocked:
        ! 1 = denominantor
        ! 2 = real numerator
        ! 3 = shift
        ! 4 = imaginary numerator
        integer, intent(in) :: iValue
        real(dp), intent(in) :: this(:)
        character(len=*), intent(in) :: filename
        logical, intent(in) :: tPrint
        integer :: length
        integer, intent(in) :: blocklength ! this is the minimum step size (e.g. 2)
        integer, intent(out) :: corrlength ! this is the correlation length in iterations (e.g. 2**N)
        integer :: i
        real(dp) :: mean, error, eie
        real(dp), allocatable :: that(:) ! stores this array
        integer :: blocking_events
        real(dp), allocatable :: mean_array(:), error_array(:), eie_array(:)
        integer :: which_element, iunit
        real(dp) :: final_error
        character(len=*), parameter :: t_r = "automatic_reblocking_analysis"

        length = size(this, 1)
        allocate(that(length))
        that = this
        blocking_events = int(log(real(size(that), dp)) / log(2.0_dp) - 1.0_dp)
        ! Yuck! This justs finds the
        ! expected number of reblocking analyses automatically
        allocate(mean_array(blocking_events + 1))
        allocate(error_array(blocking_events + 1))
        allocate(eie_array(blocking_events + 1))

        iunit = 0
        if (tPrint) then
            iunit = get_free_unit()
            open(iunit, file=filename)
            write(iunit, "(A)") "# No.Blocks      Mean       Error        Error_in_Error"
        end if

        call analyze_data(that, mean, error, eie)
        if (tPrint) write(iunit, *) size(that), mean, error, eie
        mean_array(1) = mean
        error_array(1) = error
        eie_array(1) = eie
        do i = 1, blocking_events
            call reblock_data(that, blocklength)
            call analyze_data(that, mean, error, eie)
            mean_array(i + 1) = mean
            error_array(i + 1) = error
            eie_array(i + 1) = eie
            if (tPrint) write(iunit, *) size(that), mean, error, eie
        end do
        call check_reblock_monotonic_inc(error_array, tPrint, iValue)
        call find_max_error(error_array, final_error, which_element)
        corrlength = blocklength**(which_element - 1)
        if (tPrint) then
            if (iValue == 1) then
                write(stdout, "(A,I7)") "Number of blocks assumed for calculation of error in projected energy denominator: ", &
                    length / corrlength
            else if (iValue == 2) then
                write(stdout, "(A,I7)") "Number of blocks assumed for calculation of error in projected energy numerator: ", &
                    length / corrlength
            else if (iValue == 3) then
                write(stdout, "(A,I7)") "Number of blocks assumed for calculation of error in shift: ", &
                    length / corrlength
            else if (iValue == 4) then
                write(stdout, "(A,I7)") "Number of blocks assumed for calculation of error in Im projected energy numerator: ", &
                    length / corrlength
            else
                call stop_all(t_r, "Error in iValue")
            end if
        end if
        if (errordebug > 0) then
            write(stdout, *) "Mean", mean_array(1)
            write(stdout, *) "Final error", final_error, "number of blocks", length / blocklength**(which_element - 1)
            write(stdout, *) "Correlation length", corrlength
        end if
        if (tPrint) close(iunit)

        deallocate(mean_array, error_array, eie_array)

    end subroutine automatic_reblocking_analysis