util_mod_numerical.F90 Source File

module util_mod_numerical_int32

    ! This module deals with templated things that must be numbers.
    !
    ! Particularly, +, -, <, >, <=, >= must be defined

    use constants
    implicit none
    interface binary_search_first_ge
        module procedure binary_search_first_ge_int32
    end interface

    interface stats_out
        module procedure stats_out_int32
    end interface

contains
    pure function binary_search_first_ge_int32 (arr, val, ind_change, offset) result(pos)

        ! Find the first element in an array which is >= val, in an array
        ! which has been sorted. Note that we only use the < comparison.
        !
        ! If there is no such element, the function returns -1.
        !
        ! If ind_change is specified, then the values from this index onwards
        ! are modified by the value in offset. This allows us to essentially
        ! remove an element from an underlying list that this is the cumulative
        ! sum of.

        integer(kind=int32), intent(in) :: arr(:)
        integer(kind=int32), intent(in) :: val
        integer, intent(in), optional :: ind_change
        integer(kind=int32), intent(in), optional :: offset
        integer :: pos

        integer :: hi, lo, ind_change_
        integer(kind=int32) :: tmp

        ! The search range
        lo = lbound(arr, 1)
        hi = ubound(arr, 1)

        ! The change/offset values
        if (present(ind_change)) then
            ind_change_ = ind_change
        else
            ind_change_ = hi + 1
        end if

        ! Test if such an element exists
        tmp = arr(hi)
        if (hi >= ind_change_) tmp = tmp + offset
        if (tmp < val) then
            pos = -1
            return
        endif

        do while (hi /= lo)
            pos = int(real(hi + lo,sp) / 2_sp)

            !if (arr(pos) >= val) then
            tmp = arr(pos)
            if (pos >= ind_change_) tmp = tmp + offset
            if (.not.(tmp < val)) then
                hi = pos
            else
                lo = pos + 1
            endif
        enddo

        ! Return the converged value.
        pos = hi

    end function


    subroutine stats_out_int32 (state, mc_col, num, title)

        type(write_state_t), intent(inout) :: state
        logical, intent(in) :: mc_col
        integer(kind=int32), intent(in) :: num
        character(*), intent(in) :: title
        character(5) :: colno
        character(20) :: title_str

        ! Increment the column counter
        state%cols = state%cols + 1
        if (mc_col) state%cols_mc = state%cols_mc + 1

        ! Add column spacing
        write(state%funit, '(" ")', advance='no')
        if (state%mc_out .and. mc_col) write(stdout, '(" ")', advance='no')

        if (state%init) then
            write(colno, '(i5)') state%cols
            title_str = trim(adjustl(colno)) // ". " // trim(adjustl(title))
            write(state%funit, '(a)', advance='no') trim(title_str)

            if (state%mc_out .and. mc_col) then
                write(stdout, '(a20)', advance='no') title_str
            end if
        else
            write(state%funit, '(i15)', advance='no') num
            if (state%mc_out .and. mc_col) then
                write(stdout, '(i15)', advance='no') num
            end if
        end if

    end subroutine

end module




module util_mod_numerical_int64

    ! This module deals with templated things that must be numbers.
    !
    ! Particularly, +, -, <, >, <=, >= must be defined

    use constants
    implicit none
    interface binary_search_first_ge
        module procedure binary_search_first_ge_int64
    end interface

    interface stats_out
        module procedure stats_out_int64
    end interface

contains
    pure function binary_search_first_ge_int64 (arr, val, ind_change, offset) result(pos)

        ! Find the first element in an array which is >= val, in an array
        ! which has been sorted. Note that we only use the < comparison.
        !
        ! If there is no such element, the function returns -1.
        !
        ! If ind_change is specified, then the values from this index onwards
        ! are modified by the value in offset. This allows us to essentially
        ! remove an element from an underlying list that this is the cumulative
        ! sum of.

        integer(kind=int64), intent(in) :: arr(:)
        integer(kind=int64), intent(in) :: val
        integer, intent(in), optional :: ind_change
        integer(kind=int64), intent(in), optional :: offset
        integer :: pos

        integer :: hi, lo, ind_change_
        integer(kind=int64) :: tmp

        ! The search range
        lo = lbound(arr, 1)
        hi = ubound(arr, 1)

        ! The change/offset values
        if (present(ind_change)) then
            ind_change_ = ind_change
        else
            ind_change_ = hi + 1
        end if

        ! Test if such an element exists
        tmp = arr(hi)
        if (hi >= ind_change_) tmp = tmp + offset
        if (tmp < val) then
            pos = -1
            return
        endif

        do while (hi /= lo)
            pos = int(real(hi + lo,sp) / 2_sp)

            !if (arr(pos) >= val) then
            tmp = arr(pos)
            if (pos >= ind_change_) tmp = tmp + offset
            if (.not.(tmp < val)) then
                hi = pos
            else
                lo = pos + 1
            endif
        enddo

        ! Return the converged value.
        pos = hi

    end function


    subroutine stats_out_int64 (state, mc_col, num, title)

        type(write_state_t), intent(inout) :: state
        logical, intent(in) :: mc_col
        integer(kind=int64), intent(in) :: num
        character(*), intent(in) :: title
        character(5) :: colno
        character(20) :: title_str

        ! Increment the column counter
        state%cols = state%cols + 1
        if (mc_col) state%cols_mc = state%cols_mc + 1

        ! Add column spacing
        write(state%funit, '(" ")', advance='no')
        if (state%mc_out .and. mc_col) write(stdout, '(" ")', advance='no')

        if (state%init) then
            write(colno, '(i5)') state%cols
            title_str = trim(adjustl(colno)) // ". " // trim(adjustl(title))
            write(state%funit, '(a)', advance='no') trim(title_str)

            if (state%mc_out .and. mc_col) then
                write(stdout, '(a20)', advance='no') title_str
            end if
        else
            write(state%funit, '(i16)', advance='no') num
            if (state%mc_out .and. mc_col) then
                write(stdout, '(i16)', advance='no') num
            end if
        end if

    end subroutine

end module



#if !defined(SX)

module util_mod_numerical_real

    ! This module deals with templated things that must be numbers.
    !
    ! Particularly, +, -, <, >, <=, >= must be defined

    use constants
    implicit none
    interface binary_search_first_ge
        module procedure binary_search_first_ge_real
    end interface

    interface stats_out
        module procedure stats_out_real
    end interface

contains
    pure function binary_search_first_ge_real (arr, val, ind_change, offset) result(pos)

        ! Find the first element in an array which is >= val, in an array
        ! which has been sorted. Note that we only use the < comparison.
        !
        ! If there is no such element, the function returns -1.
        !
        ! If ind_change is specified, then the values from this index onwards
        ! are modified by the value in offset. This allows us to essentially
        ! remove an element from an underlying list that this is the cumulative
        ! sum of.

        real(kind=sp), intent(in) :: arr(:)
        real(kind=sp), intent(in) :: val
        integer, intent(in), optional :: ind_change
        real(kind=sp), intent(in), optional :: offset
        integer :: pos

        integer :: hi, lo, ind_change_
        real(kind=sp) :: tmp

        ! The search range
        lo = lbound(arr, 1)
        hi = ubound(arr, 1)

        ! The change/offset values
        if (present(ind_change)) then
            ind_change_ = ind_change
        else
            ind_change_ = hi + 1
        end if

        ! Test if such an element exists
        tmp = arr(hi)
        if (hi >= ind_change_) tmp = tmp + offset
        if (tmp < val) then
            pos = -1
            return
        endif

        do while (hi /= lo)
            pos = int(real(hi + lo,sp) / 2_sp)

            !if (arr(pos) >= val) then
            tmp = arr(pos)
            if (pos >= ind_change_) tmp = tmp + offset
            if (.not.(tmp < val)) then
                hi = pos
            else
                lo = pos + 1
            endif
        enddo

        ! Return the converged value.
        pos = hi

    end function


    subroutine stats_out_real (state, mc_col, num, title)

        type(write_state_t), intent(inout) :: state
        logical, intent(in) :: mc_col
        real(kind=sp), intent(in) :: num
        character(*), intent(in) :: title
        character(5) :: colno
        character(20) :: title_str

        ! Increment the column counter
        state%cols = state%cols + 1
        if (mc_col) state%cols_mc = state%cols_mc + 1

        ! Add column spacing
        write(state%funit, '(" ")', advance='no')
        if (state%mc_out .and. mc_col) write(stdout, '(" ")', advance='no')

        if (state%init) then
            write(colno, '(i5)') state%cols
            title_str = trim(adjustl(colno)) // ". " // trim(adjustl(title))
            write(state%funit, '(a)', advance='no') trim(title_str)

            if (state%mc_out .and. mc_col) then
                write(stdout, '(a20)', advance='no') title_str
            end if
        else
            write(state%funit, '(g20.12e3)', advance='no') num
            if (state%mc_out .and. mc_col) then
                write(stdout, '(g20.12e3)', advance='no') num
            end if
        end if

    end subroutine

end module


#endif



module util_mod_numerical_doub

    ! This module deals with templated things that must be numbers.
    !
    ! Particularly, +, -, <, >, <=, >= must be defined

    use constants
    implicit none
    interface binary_search_first_ge
        module procedure binary_search_first_ge_doub
    end interface

    interface stats_out
        module procedure stats_out_doub
    end interface

contains
    pure function binary_search_first_ge_doub (arr, val, ind_change, offset) result(pos)

        ! Find the first element in an array which is >= val, in an array
        ! which has been sorted. Note that we only use the < comparison.
        !
        ! If there is no such element, the function returns -1.
        !
        ! If ind_change is specified, then the values from this index onwards
        ! are modified by the value in offset. This allows us to essentially
        ! remove an element from an underlying list that this is the cumulative
        ! sum of.

        real(kind=dp), intent(in) :: arr(:)
        real(kind=dp), intent(in) :: val
        integer, intent(in), optional :: ind_change
        real(kind=dp), intent(in), optional :: offset
        integer :: pos

        integer :: hi, lo, ind_change_
        real(kind=dp) :: tmp

        ! The search range
        lo = lbound(arr, 1)
        hi = ubound(arr, 1)

        ! The change/offset values
        if (present(ind_change)) then
            ind_change_ = ind_change
        else
            ind_change_ = hi + 1
        end if

        ! Test if such an element exists
        tmp = arr(hi)
        if (hi >= ind_change_) tmp = tmp + offset
        if (tmp < val) then
            pos = -1
            return
        endif

        do while (hi /= lo)
            pos = int(real(hi + lo,sp) / 2_sp)

            !if (arr(pos) >= val) then
            tmp = arr(pos)
            if (pos >= ind_change_) tmp = tmp + offset
            if (.not.(tmp < val)) then
                hi = pos
            else
                lo = pos + 1
            endif
        enddo

        ! Return the converged value.
        pos = hi

    end function


    subroutine stats_out_doub (state, mc_col, num, title)

        type(write_state_t), intent(inout) :: state
        logical, intent(in) :: mc_col
        real(kind=dp), intent(in) :: num
        character(*), intent(in) :: title
        character(5) :: colno
        character(40) :: title_str

        ! Increment the column counter
        state%cols = state%cols + 1
        if (mc_col) state%cols_mc = state%cols_mc + 1

        ! Add column spacing
        write(state%funit, '(" ")', advance='no')
        if (state%mc_out .and. mc_col) write(stdout, '(" ")', advance='no')

        if (state%init) then
            write(colno, '(i5)') state%cols
            title_str = trim(adjustl(colno)) // ". " // trim(adjustl(title))
            write(state%funit, '(a)', advance='no') trim(title_str)

            if (state%mc_out .and. mc_col) then
                write(stdout, '(a40)', advance='no') title_str
            end if
        else
            write(state%funit, '(g20.12e3)', advance='no') num
            if (state%mc_out .and. mc_col) then
                write(stdout, '(g20.12e3)', advance='no') num
            end if
        end if

    end subroutine

end module




module util_mod_numerical
    use util_mod_numerical_int32
    use util_mod_numerical_int64
#if !defined(SX)
    use util_mod_numerical_real
#endif
    use util_mod_numerical_doub
end module