#include "macros.h" module util_mod use util_mod_comparisons, only: operator(.arrgt.), operator(.arrlt.), arr_gt, arr_lt use util_mod_numerical, only: binary_search_first_ge, stats_out use util_mod_cpts, only: arr_2d_ptr, arr_2d_dims, ptr_abuse_1d, & ptr_abuse_scalar, ptr_abuse_2d use basic_float_math, only: near_zero, operator(.isclose.), isclose use constants, only: sp, dp, int32, int64, n_int, inum_runs, lenof_sign, & sizeof_int, stdout use binomial_lookup, only: factrl => factorial, binomial_lookup_table_i64 #ifdef GFORTRAN_ use constants, only: int128 use binomial_lookup, only: binomial_lookup_table_i128 #endif use fmt_utils, only: int_fmt use dSFMT_interface, only: genrand_real2_dSFMT use DetBitOps, only: DetBitLt use, intrinsic :: iso_c_binding, only: c_char, c_int, c_double use mpi, only: MPI_WTIME use error_handling_neci, only: stop_all, neci_flush, warning_neci ! We want to use the builtin etime intrinsic with ifort to ! work around some broken behaviour. #if defined(IFORT_) || defined(INTELLLVM_) use ifport, only: etime #endif implicit none #if defined(IFORT_) || defined(INTELLLVM_) public #else private #endif public :: factrl, choose_i64, NECI_icopy, operator(.implies.), & abs_l1, abs_sign, near_zero, operator(.isclose.), isclose, operator(.div.), & stochastic_round, stochastic_round_r #ifdef GFORTRAN_ public :: choose_i128 #endif public :: error_function, error_function_c, stop_all public :: arr_2d_ptr, ptr_abuse_1d, ptr_abuse_2d, ptr_abuse_scalar public :: neci_etime, get_free_unit, int_fmt,& strlen_wrap, record_length, open_new_file, & append_ext, get_unique_filename, neci_flush, print_cstr_local, & stats_out public :: arr_lt, arr_gt, operator(.arrlt.), operator(.arrgt.), & find_next_comb, binary_search_ilut, binary_search_custom, binary_search_first_ge, & cumsum, pairswap, swap, lex_leq, lex_geq, & get_permutations, custom_findloc, addToIntArray, fuseIndex, linearIndex, & getSpinIndex, binary_search_int, binary_search_real, clamp public :: warning_neci public :: EnumBase_t interface binary_search_int module procedure binary_search_int_int32 module procedure binary_search_int_int64 end interface interface operator(.implies.) module procedure implies end interface interface choose_i64 module procedure choose_i64_int32 module procedure choose_i64_int64 end interface #ifdef GFORTRAN_ interface choose_i128 module procedure choose_i128_int32 module procedure choose_i128_int64 end interface #endif interface clamp !! If v compares less than lo, returns lo; !! otherwise if hi compares less than v, returns hi; otherwise returns v. !! Is also defined for lo > hi! module procedure clamp_integer_int32 module procedure clamp_integer_int64 module procedure clamp_real_sp module procedure clamp_real_dp end interface interface pure function strlen_wrap(str) result(len) bind(c) import :: c_char, c_int implicit none character(c_char), intent(in) :: str(*) integer(c_int) :: len end function pure function erf_local(x) result(e) bind(c, name='erf') import :: c_double implicit none real(c_double), intent(in) :: x real(c_double) :: e end function pure function erfc_local(x) result(ec) bind(c, name='erfc') import :: c_double implicit none real(c_double), intent(in) :: x real(c_double) :: ec end function end interface interface operator(.div.) module procedure div_int32, div_int64 #ifdef GFORTRAN_ module procedure div_int128 #endif end interface interface abs_sign module procedure abs_int4_sign module procedure abs_int8_sign module procedure abs_real_sign end interface interface abs_l1 module procedure abs_l1_dp module procedure abs_l1_sp module procedure abs_l1_cdp module procedure abs_l1_csp end interface abs_l1 interface fuseIndex module procedure fuseIndex_int32 module procedure fuseIndex_int64 end interface fuseIndex interface swap module procedure swap_int64 module procedure swap_int32 end interface swap interface custom_findloc module procedure custom_findloc_integer_int32 module procedure custom_findloc_integer_int64 module procedure custom_findloc_real_sp module procedure custom_findloc_real_dp module procedure custom_findloc_complex_sp module procedure custom_findloc_complex_dp module procedure custom_findloc_logical_ end interface custom_findloc interface cumsum module procedure cumsum_integer_int32 module procedure cumsum_integer_int64 module procedure cumsum_real_sp module procedure cumsum_real_dp module procedure cumsum_complex_sp module procedure cumsum_complex_dp end interface type, abstract :: EnumBase_t integer :: val contains private procedure :: eq_EnumBase_t procedure :: neq_EnumBase_t generic, public :: operator(==) => eq_EnumBase_t generic, public :: operator(/=) => neq_EnumBase_t end type contains function stochastic_round(r) result(i) ! Stochastically round the supplied real value to an integer. This is ! the primary method of introducing the monte-carlo nature of spawning ! or death into the algorithm. ! --> Probably nicer to use a centralised implementation than a bunch ! of hacked-in ones all over the place... ! ! Unfortunately, we cannot make this pure, as we would need to have ! a mutable variable in genrand_real2_dSFMT... real(dp), intent(in) :: r integer :: i real(dp) :: res i = int(r) res = r - real(i, dp) if(abs(res) >= 1.0e-12_dp) then if(abs(res) > genrand_real2_dSFMT()) & i = i + nint(sign(1.0_dp, r)) end if end function elemental function stochastic_round_r(num, r) result(i) ! Perform the stochastic rounding of the above function where the ! random number is already specified. real(dp), intent(in) :: num, r integer :: i real(dp) :: res i = int(num) res = num - real(i, dp) if(abs(res) >= 1.0e-12_dp) then if(abs(res) > r) & i = i + nint(sign(1.0_dp, num)) end if end function stochastic_round_r subroutine print_cstr_local(str, l) character(c_char), intent(in) :: str(*) integer, intent(in) :: l character(len=l) :: tmp_s tmp_s = transfer(str(1:l), tmp_s) write(stdout, '(a)', advance='no') tmp_s end subroutine ! routine to calculation the absolute magnitude of a complex integer ! variable (to nearest integer) pure real(dp) function abs_int4_sign(sgn) integer(int32), intent(in) :: sgn(lenof_sign / inum_runs) #ifdef CMPLX_ abs_int4_sign = real(int(sqrt(real(sgn(1), dp)**2 + real(sgn(2), dp)**2)), dp) ! The integerisation here is an approximation, but one that is ! used in the integer algorithm, so is retained in this real ! version of the algorithm #else abs_int4_sign = real(abs(sgn(1)), dp) #endif end function abs_int4_sign !routine to calculation the absolute magnitude of a complex integer(int64) variable (to nearest integer) pure integer(kind=int64) function abs_int8_sign(wsign) integer(kind=int64), dimension(lenof_sign/inum_runs), intent(in) :: wsign #ifdef CMPLX_ abs_int8_sign = nint(sqrt(real(wsign(1), dp)**2 + real(wsign(2), dp)**2), int64) #else abs_int8_sign = abs(wsign(1)) #endif end function abs_int8_sign pure real(dp) function abs_real_sign(sgn) real(dp), intent(in) :: sgn(lenof_sign / inum_runs) #ifdef CMPLX_ abs_real_sign = real(nint(sqrt(sum(sgn**2))), dp) #else abs_real_sign = abs(sgn(1)) #endif end function ! --------- L1 norm function ! These return the absolute L1 norm of the specified value ! ! --> for complex numbers this is not sqrt(r**2 + i**2), but is the sum ! of the absolute values of the terms !---------------------------- pure function abs_l1_sp(val) result(ret) real(sp), intent(in) :: val real(sp) :: ret ret = abs(val) end function pure function abs_l1_dp(val) result(ret) real(dp), intent(in) :: val real(dp) :: ret ret = abs(val) end function pure function abs_l1_csp(val) result(ret) complex(sp), intent(in) :: val real(sp) :: ret ret = abs(real(val, sp)) + abs(aimag(val)) end function pure function abs_l1_cdp(val) result(ret) complex(dp), intent(in) :: val real(dp) :: ret ret = abs(real(val, dp)) + abs(aimag(val)) end function !--- Array utilities --- SUBROUTINE NECI_ICOPY(N, A, IA, B, IB) ! Copy elements from integer array A to B. ! Simple version of BLAS routine ICOPY, which isn't always implemented ! in BLAS. ! Fortran 90 array features allow this to be done in one line of ! standard fortran, so this is just for legacy purposes. ! In: ! N: number of elements in A. ! A: vector to be copied. ! IA: increment between elements to be copied in A. ! IA=1 for continuous data blocks. ! IB: increment between elements to be copied to in B. ! IB=1 for continuous data blocks. ! Out: ! B: result vector. IMPLICIT NONE ! Arguments INTEGER, INTENT(IN) :: N, IA, IB INTEGER, INTENT(IN) :: A(IA * N) INTEGER, INTENT(OUT) :: B(IB * N) ! Variables INTEGER I, IAX, IBX DO I = 1, N IAX = (I - 1) * IA + 1 IBX = (I - 1) * IB + 1 B(IBX) = A(IAX) ENDDO RETURN END SUBROUTINE NECI_ICOPY subroutine addToIntArray(arr, ind, elem) integer, intent(inout), allocatable :: arr(:) integer, intent(in) :: ind, elem integer, allocatable :: tmp(:) integer :: nelems if(allocated(arr)) then nelems = size(arr) if(ind > nelems) then ! resize the array allocate(tmp(nelems)) tmp = arr deallocate(arr) allocate(arr(ind), source=0) arr(1:nelems) = tmp(1:nelems) endif else allocate(arr(ind), source=0) endif arr(ind) = elem end subroutine addToIntArray !------------------------------------------------------------------------------------------! !> Custom implementation of the findloc intrinsic (with somewhat reduced functionality) !! as it requires fortran2008 support and is thus not available for some relevant compilers pure function custom_findloc_integer_int32(arr, val, back) result(loc) integer(int32), intent(in) :: arr(:) integer(int32), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) == val) then loc = i return endif end do end function custom_findloc_integer_int32 pure function custom_findloc_integer_int64(arr, val, back) result(loc) integer(int64), intent(in) :: arr(:) integer(int64), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) == val) then loc = i return endif end do end function custom_findloc_integer_int64 pure function custom_findloc_real_sp(arr, val, back) result(loc) real(sp), intent(in) :: arr(:) real(sp), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) .isclose. val) then loc = i return endif end do end function custom_findloc_real_sp pure function custom_findloc_real_dp(arr, val, back) result(loc) real(dp), intent(in) :: arr(:) real(dp), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) .isclose. val) then loc = i return endif end do end function custom_findloc_real_dp pure function custom_findloc_complex_sp(arr, val, back) result(loc) complex(sp), intent(in) :: arr(:) complex(sp), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) .isclose. val) then loc = i return endif end do end function custom_findloc_complex_sp pure function custom_findloc_complex_dp(arr, val, back) result(loc) complex(dp), intent(in) :: arr(:) complex(dp), intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) .isclose. val) then loc = i return endif end do end function custom_findloc_complex_dp pure function custom_findloc_logical_(arr, val, back) result(loc) logical, intent(in) :: arr(:) logical, intent(in) :: val logical, intent(in), optional :: back integer :: loc integer :: i, first, last, step logical :: back_ def_default(back_, back, .false.) if(back_) then first = size(arr) last = 1 step = -1 else first = 1 last = size(arr) step = 1 end if loc = 0 do i = first, last, step if(arr(i) .eqv. val) then loc = i return endif end do end function custom_findloc_logical_ !--- Indexing utilities pure function fuseIndex_int32(q, p) result(ind) ! fuse p,q into one symmetric index ! the resulting index is not contigious in p or q ! Input: p,q - 2d-array indices ! Output: ind - 1d-array index assuming the array is symmetric w.r. p<->q integer(int32), intent(in) :: p, q integer(int32) :: ind ! qp and pq are considered to be the same index ! -> permutational symmetry ! implemented in terms of fuseIndex_int64 ind = int(fuseIndex_int64(int(q, int64), int(p, int64))) end function fuseIndex_int32 !------------------------------------------------------------------------------------------! pure function fuseIndex_int64(x, y) result(xy) ! create a composite index out of two indices, assuming they are unordered ! i.e. their ordering does not matter ! Input: p,q - 2d-array indices ! Output: ind - 1d-array index assuming the array is symmetric w.r. p<->q integer(int64), intent(in) :: x, y integer(int64) :: xy if(x < y) then xy = x + y * (y - 1) / 2 else xy = y + x * (x - 1) / 2 endif end function fuseIndex_int64 !------------------------------------------------------------------------------------------! elemental subroutine swap_int32(a, b) ! exchange the value of two integers a,b ! Input: a,b - integers to swapp (on return, a has the value of b on call and vice versa) integer(int32), intent(inout) :: a, b integer(int32) :: tmp tmp = a a = b b = tmp end subroutine swap_int32 !------------------------------------------------------------------------------------------! elemental subroutine swap_int64(a, b) ! exchange the value of two integers a,b ! Input: a,b - integers to swapp (on return, a has the value of b on call and vice versa) integer(int64), intent(inout) :: a, b integer(int64) :: tmp tmp = a a = b b = tmp end subroutine swap_int64 !------------------------------------------------------------------------------------------! pure subroutine pairSwap(a, i, b, j) ! exchange a pair of integers integer(int64), intent(inout) :: a, i, b, j call swap(a, b) call swap(i, j) end subroutine pairSwap !------------------------------------------------------------------------------------------! function linearIndex(p, q, dim) result(ind) ! fuse p,q into one contiguous index ! the resulting index is contiguous in q ! Input: p,q - 2d-array indices ! dim - dimension of the underlying array in q-direction ! Output: ind - contiguous 1d-array index integer, intent(in) :: p, q, dim integer :: ind ind = q + (p - 1) * dim end function linearIndex pure elemental function getSpinIndex(orb) result(ms) ! return a spin index of the orbital orb which can be used to address arrays ! Input: orb - spin orbital ! Output: ms - spin index of orb with the following values: ! 0 - alpha ! 1 - beta integer, intent(in) :: orb integer :: ms ms = mod(orb, 2) end function getSpinIndex !> @brief !> Calculate 1 + ... + n integer elemental function gauss_sum(n) integer, intent(in) :: n gauss_sum = (n * (n + 1)) .div. 2 end function !> @brief !> Get the index in the binomial_lookup_table integer elemental function get_index(n, k) integer, intent(in) :: n, k get_index = gauss_sum((n - 3) .div. 2) + gauss_sum((n - 4) .div. 2) + k - 1 end function ! Unfortunately there are no recursive elemental functions in Fortran. recursive pure function choose_i64_int32(n, r, signal_overflow) result(res) !! Return the binomail coefficient nCr(n, r) integer(int32), intent(in) :: n, r logical, intent(in), optional :: signal_overflow !! If true then the function returns -1 instead of aborting !! when overflow is encountered. integer(int64) :: res integer(int64) :: k character(*), parameter :: this_routine = "choose_i64" ! NOTE: This is highly optimized. If you change something, please time it! #ifdef DEBUG_ block if (.not. (n >= 0_int32)) then call stop_all (this_routine, "Assert fail: n >= 0_int32") end if end block #endif #ifdef DEBUG_ block if (.not. (r >= 0_int32)) then call stop_all (this_routine, "Assert fail: r >= 0_int32") end if end block #endif if(r > n) then res = 0_int64 return end if k = int(merge(r, n - r, r <= n - r), kind=int64) if (k == 0) then res = 1_int64 else if (k == 1) then res = int(n, int64) else if (n <= 66) then ! use lookup table res = binomial_lookup_table_i64(get_index(int(n), int(k))) else block integer(int64) :: prev prev = choose_i64_int32(n - 1, r - 1, signal_overflow) ! Note that the recursion stops at n = 66 res = (prev * n) .div. k check_for_overflow: if (prev < 0 .or. res < 0) then if (present(signal_overflow)) then if (signal_overflow) then res = -1 else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int64.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int64.') #endif end if else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int64.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int64.') #endif end if end if check_for_overflow end block end if end function #ifdef GFORTRAN_ recursive pure function choose_i128_int32(n, r, signal_overflow) result(res) !! Return the binomail coefficient nCr(n, r) integer(int32), intent(in) :: n, r logical, intent(in), optional :: signal_overflow !! If true then the function returns -1 instead of aborting !! when overflow is encountered. integer(int128) :: res integer(int128) :: k character(*), parameter :: this_routine = "choose_i128" ! NOTE: This is highly optimized. If you change something, please time it! #ifdef DEBUG_ block if (.not. (n >= 0_int32)) then call stop_all (this_routine, "Assert fail: n >= 0_int32") end if end block #endif #ifdef DEBUG_ block if (.not. (r >= 0_int32)) then call stop_all (this_routine, "Assert fail: r >= 0_int32") end if end block #endif if(r > n) then res = 0_int128 return end if k = int(merge(r, n - r, r <= n - r), kind=int128) if (k == 0) then res = 1_int128 else if (k == 1) then res = int(n, int128) else if (n <= 130) then ! use lookup table res = binomial_lookup_table_i128(get_index(int(n), int(k))) else ! Note that the recursion stops at n = 130 block integer(int128) :: prev prev = choose_i128_int32(n - 1, r - 1, signal_overflow) ! Note that the recursion stops at n = 66 res = (prev * n) .div. k check_for_overflow: if (prev < 0 .or. res < 0) then if (present(signal_overflow)) then if (signal_overflow) then res = -1 else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int128.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int128.') #endif end if else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int128.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int128.') #endif end if end if check_for_overflow end block end if end function #endif ! Unfortunately there are no recursive elemental functions in Fortran. recursive pure function choose_i64_int64(n, r, signal_overflow) result(res) !! Return the binomail coefficient nCr(n, r) integer(int64), intent(in) :: n, r logical, intent(in), optional :: signal_overflow !! If true then the function returns -1 instead of aborting !! when overflow is encountered. integer(int64) :: res integer(int64) :: k character(*), parameter :: this_routine = "choose_i64" ! NOTE: This is highly optimized. If you change something, please time it! #ifdef DEBUG_ block if (.not. (n >= 0_int64)) then call stop_all (this_routine, "Assert fail: n >= 0_int64") end if end block #endif #ifdef DEBUG_ block if (.not. (r >= 0_int64)) then call stop_all (this_routine, "Assert fail: r >= 0_int64") end if end block #endif if(r > n) then res = 0_int64 return end if k = int(merge(r, n - r, r <= n - r), kind=int64) if (k == 0) then res = 1_int64 else if (k == 1) then res = int(n, int64) else if (n <= 66) then ! use lookup table res = binomial_lookup_table_i64(get_index(int(n), int(k))) else block integer(int64) :: prev prev = choose_i64_int64(n - 1, r - 1, signal_overflow) ! Note that the recursion stops at n = 66 res = (prev * n) .div. k check_for_overflow: if (prev < 0 .or. res < 0) then if (present(signal_overflow)) then if (signal_overflow) then res = -1 else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int64.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int64.') #endif end if else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int64.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int64.') #endif end if end if check_for_overflow end block end if end function #ifdef GFORTRAN_ recursive pure function choose_i128_int64(n, r, signal_overflow) result(res) !! Return the binomail coefficient nCr(n, r) integer(int64), intent(in) :: n, r logical, intent(in), optional :: signal_overflow !! If true then the function returns -1 instead of aborting !! when overflow is encountered. integer(int128) :: res integer(int128) :: k character(*), parameter :: this_routine = "choose_i128" ! NOTE: This is highly optimized. If you change something, please time it! #ifdef DEBUG_ block if (.not. (n >= 0_int64)) then call stop_all (this_routine, "Assert fail: n >= 0_int64") end if end block #endif #ifdef DEBUG_ block if (.not. (r >= 0_int64)) then call stop_all (this_routine, "Assert fail: r >= 0_int64") end if end block #endif if(r > n) then res = 0_int128 return end if k = int(merge(r, n - r, r <= n - r), kind=int128) if (k == 0) then res = 1_int128 else if (k == 1) then res = int(n, int128) else if (n <= 130) then ! use lookup table res = binomial_lookup_table_i128(get_index(int(n), int(k))) else ! Note that the recursion stops at n = 130 block integer(int128) :: prev prev = choose_i128_int64(n - 1, r - 1, signal_overflow) ! Note that the recursion stops at n = 66 res = (prev * n) .div. k check_for_overflow: if (prev < 0 .or. res < 0) then if (present(signal_overflow)) then if (signal_overflow) then res = -1 else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int128.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int128.') #endif end if else #if defined(IFORT_) || defined(INTELLLVM_) error stop 'Binomial coefficient exceeds range of int128.' #else call stop_all(this_routine, 'Binomial coefficient exceeds range of int128.') #endif end if end if check_for_overflow end block end if end function #endif elemental integer(int32) function div_int32(a, b) integer(int32), intent(in) :: a, b #ifdef WARNING_WORKAROUND_ div_int32 = int(real(a, kind=sp) / real(b, kind=sp), kind=int32) #else div_int32 = a / b #endif end function elemental integer(int64) function div_int64(a, b) integer(int64), intent(in) :: a, b #ifdef WARNING_WORKAROUND_ div_int64 = int(real(a, kind=dp) / real(b, kind=dp), kind=int64) #else div_int64 = a / b #endif end function #ifdef GFORTRAN_ elemental integer(int128) function div_int128(a, b) integer(int128), intent(in) :: a, b #ifdef WARNING_WORKAROUND_ div_int128 = int(real(a, kind=dp) / real(b, kind=dp), kind=int128) #else div_int128 = a / b #endif end function #endif !--- Comparison of subarrays --- logical pure function det_int_arr_gt(a, b, len) ! Make a comparison we can sort determinant integer arrays by. Return true if the ! first differing integer of a, b is such that a(i) > b(i). ! ! In: a, b - The arrays to compare ! len - An optional argument to specify the size to consider. ! If not provided, then min(size(a), size(b)) is used. ! Ret: - If a > b !NOTE: These will sort by the bit-string integer length, n_int. !Therefore, these may be 32 or 64 bit integers and should !only be used as such. integer(kind=n_int), intent(in), dimension(:) :: a, b integer, intent(in), optional :: len integer llen, i if(present(len)) then llen = len else llen = min(size(a), size(b)) endif ! Sort by the first integer first ... i = 1 do i = 1, llen if(a(i) /= b(i)) exit enddo ! Make the comparison if(i > llen) then det_int_arr_gt = .false. else if(a(i) > b(i)) then det_int_arr_gt = .true. else det_int_arr_gt = .false. endif endif end function det_int_arr_gt logical pure function det_int_arr_eq(a, b, len) ! If two specified integer arrays are equal, return true. Otherwise ! return false. ! ! In: a, b - The arrays to consider ! len - The maximum length to consider. Otherwise will use whole ! length of array !NOTE: These will sort by the bit-string integer length, n_int. !Therefore, these may be 32 or 64 bit integers and should !only be used as such. integer(kind=n_int), intent(in), dimension(:) :: a, b integer, intent(in), optional :: len integer llen, i ! Obtain the lengths of the arrays if a bound is not specified. ! Return false if mismatched sizes and not specified. if(present(len)) then llen = len else if(size(a) /= size(b)) then det_int_arr_eq = .false. return endif llen = size(a) endif ! Run through the arrays. Return if they differ at any point. do i = 1, llen if(a(i) /= b(i)) then det_int_arr_eq = .false. return endif enddo ! If we get this far, they are equal det_int_arr_eq = .true. end function det_int_arr_eq !--- Searching --- ! NOTE: This can only be used for binary searching determinant bit ! strings now. We can template it if it wants to be more general ! in the future if needed. function binary_search_ilut(arr, val, cf_len) result(pos) integer(kind=n_int), intent(in) :: arr(:, :) integer(kind=n_int), intent(in) :: val(:) integer, intent(in), optional :: cf_len integer :: data_lo, data_hi, val_lo, val_hi integer :: pos integer :: hi, lo ! The search range lo = lbound(arr, 2) hi = ubound(arr, 2) ! Account for poor usage (i.e. array len == 0) if(hi < lo) then pos = -lo return endif ! Have we specified how much to look at? data_lo = lbound(arr, 1) val_lo = lbound(val, 1) if(present(cf_len)) then data_hi = data_lo + cf_len - 1 val_hi = val_lo + cf_len - 1 else data_hi = ubound(arr, 1) val_hi = ubound(val, 1) endif ! Narrow the search range down in steps. do while(hi /= lo) pos = int(real(hi + lo, sp) / 2_sp) if(all(arr(data_lo:data_hi, pos) == val(val_lo:val_hi))) then exit else if(arr_gt(val(val_lo:val_hi), arr(data_lo:data_hi, pos))) then ! val is "greater" than arr(:len,pos). ! The lowest position val can take is hence pos + 1 (i.e. if ! val is greater than pos by smaller than pos + 1). lo = pos + 1 else ! arr(:,pos) is "greater" than val. ! The highest position val can take is hence pos (i.e. if val is ! smaller than pos but greater than pos - 1). This is why ! we differ slightly from a standard binary search (where lo ! is set to be pos+1 and hi to be pos-1 accordingly), as ! a standard binary search assumes that the element you are ! searching for actually appears in the array being ! searched... hi = pos endif enddo ! If we have narrowed down to one position, and it is not the item, ! then return -pos to indicate that the item is not present, but that ! this is the location it should be in. if(hi == lo) then if(all(arr(data_lo:data_hi, hi) == val(val_lo:val_hi))) then pos = hi else if(arr_gt(val(val_lo:val_hi), arr(data_lo:data_hi, hi))) then pos = -hi - 1 else pos = -hi endif endif end function binary_search_ilut pure function binary_search_int_int32(arr, val) result(pos) integer(int32), intent(in) :: arr(:) integer(int32), intent(in) :: val integer(int64) :: pos integer(int64) :: hi, lo lo = 1 hi = size(arr) if(hi < lo) then pos = -lo return end if do while(hi /= lo) pos = int((hi + lo) / 2.0_dp, kind=int64) if(arr(pos) == val) then exit else if(val > arr(pos)) then lo = pos + 1 else hi = pos end if end do if(hi == lo) then if(arr(hi) == val) then pos = hi else pos = -1 end if end if end function binary_search_int_int32 pure function binary_search_int_int64(arr, val) result(pos) integer(int64), intent(in) :: arr(:) integer(int64), intent(in) :: val integer(int64) :: pos integer(int64) :: hi, lo lo = 1 hi = size(arr) if(hi < lo) then pos = -lo return end if do while(hi /= lo) pos = int((hi + lo) / 2.0_dp, kind=int64) if(arr(pos) == val) then exit else if(val > arr(pos)) then lo = pos + 1 else hi = pos end if end do if(hi == lo) then if(arr(hi) == val) then pos = hi else pos = -1 end if end if end function binary_search_int_int64 function binary_search_real(arr, val, thresh) & result(pos) real(dp), intent(in) :: arr(:) real(dp), intent(in) :: val real(dp), intent(in) :: thresh integer :: pos integer :: hi, lo ! The search range lo = lbound(arr, 1) hi = ubound(arr, 1) ! Account for poor usage (i.e. array len == 0) if(hi < lo) then pos = -lo return endif ! Narrow the search range down in steps. do while(hi /= lo) pos = int(real(hi + lo, dp) / 2_dp) if(abs(arr(pos) - val) < thresh) then exit else if(val > arr(pos)) then ! val is "greater" than arr(:len,pos). ! The lowest position val can take is hence pos + 1 (i.e. if ! val is greater than pos by smaller than pos + 1). lo = pos + 1 else ! arr(:,pos) is "greater" than val. ! The highest position val can take is hence pos (i.e. if val is ! smaller than pos but greater than pos - 1). This is why ! we differ slightly from a standard binary search (where lo ! is set to be pos+1 and hi to be pos-1 accordingly), as ! a standard binary search assumes that the element you are ! searching for actually appears in the array being ! searched... hi = pos endif enddo ! If we have narrowed down to one position, and it is not the item, ! then return -pos to indicate that the item is not present, but that ! this is the location it should be in. if(hi == lo) then if(abs(arr(hi) - val) < thresh) then pos = hi else if(val > arr(hi)) then pos = -hi - 1 else pos = -hi endif endif end function binary_search_real function binary_search_custom(arr, val, cf_len, custom_gt) & result(pos) interface pure function custom_gt(a, b) result(ret) import :: n_int implicit none logical :: ret integer(kind=n_int), intent(in) :: a(:), b(:) end function end interface integer(kind=n_int), intent(in) :: arr(:, :) integer(kind=n_int), intent(in) :: val(:) integer, intent(in), optional :: cf_len integer :: data_lo, data_hi, val_lo, val_hi integer :: pos integer :: hi, lo ! The search range lo = lbound(arr, 2) hi = ubound(arr, 2) ! Account for poor usage (i.e. array len == 0) if(hi < lo) then pos = -lo return endif ! Have we specified how much to look at? data_lo = lbound(arr, 1) val_lo = lbound(val, 1) if(present(cf_len)) then data_hi = data_lo + cf_len - 1 val_hi = val_lo + cf_len - 1 else data_hi = ubound(arr, 1) val_hi = ubound(val, 1) endif ! Narrow the search range down in steps. do while(hi /= lo) pos = int(real(hi + lo, sp) / 2_sp) if(DetBitLT(arr(data_lo:data_hi, pos), val(val_lo:val_hi)) == 0) then exit else if(custom_gt(val(val_lo:val_hi), arr(data_lo:data_hi, pos))) then ! val is "greater" than arr(:len,pos). ! The lowest position val can take is hence pos + 1 (i.e. if ! val is greater than pos by smaller than pos + 1). lo = pos + 1 else ! arr(:,pos) is "greater" than val. ! The highest position val can take is hence pos (i.e. if val is ! smaller than pos but greater than pos - 1). This is why ! we differ slightly from a standard binary search (where lo ! is set to be pos+1 and hi to be pos-1 accordingly), as ! a standard binary search assumes that the element you are ! searching for actually appears in the array being ! searched... hi = pos endif enddo ! If we have narrowed down to one position, and it is not the item, ! then return -pos to indicate that the item is not present, but that ! this is the location it should be in. if(hi == lo) then if(DetBitLT(arr(data_lo:data_hi, hi), val(val_lo:val_hi)) == 0) then pos = hi else if(custom_gt(val(val_lo:val_hi), arr(data_lo:data_hi, hi))) then pos = -hi - 1 else pos = -hi endif endif end function binary_search_custom !--- File utilities --- integer function record_length(bytes) ! Some compilers use record lengths in units of bytes. ! Some compilers use record lengths in units of words. ! This is an utter *pain* for reading unformatted files, ! where you must specify the record length. ! ! In: ! bytes: number of bytes in record type of interest (should ! be a multiple of 4). ! ! Returns: ! record_length: size of record in units of the compiler's ! choice. integer, intent(in) :: bytes integer :: record_length_loc inquire(iolength=record_length_loc) bytes record_length = (bytes / sizeof_int) * int(record_length_loc) ! 8 indicates 8-byte words I think end function record_length subroutine append_ext(stem, n, s) ! Returns stem.n in s. character(*), intent(in) :: stem integer, intent(in) :: n character(*), intent(out) :: s character(10) :: ext write(ext, '('//int_fmt(n, 0)//')') n s = stem//'.'//ext end subroutine append_ext subroutine get_unique_filename(stem, tincrement, tnext, istart, filename, & ext) ! Find a filename which is either the "newest" or the next to be used. ! The filename is assumed to be stem.x, where x is an integer. ! In: ! stem: stem of the filename. ! tincrement: the filename is given as stem.x if true, otherwise the ! filename is simply set to be equal to stem. ! tnext: the next unused filename is found if true, else the ! filename is set to be stem.x where stem.x exists and stem.x+1 ! doesn't and x is greater than istart ! istart: the integer of the first x value to check. ! If istart is negative, then the filename is set to be stem.x, ! where x = |istart+1|. This overrides everything else. ! ext: The file extension. Appended after the numbers. ! Out: ! filename. character(*), intent(in) :: stem logical, intent(in) :: tincrement, tnext integer, intent(in) :: istart character(*), intent(out) :: filename character(*), intent(in), optional :: ext integer :: i logical :: exists if(tincrement) then i = istart exists = .true. do while(exists) call append_ext(stem, i, filename) if(present(ext)) filename = trim(filename)//ext inquire(file=filename, exist=exists) i = i + 1 end do if(.not. tnext) then ! actually want the last file which existed. ! this will return stem.istart if stem.istart doesn't exist. i = max(istart, i - 2) call append_ext(stem, i, filename) if(present(ext)) filename = trim(filename)//ext end if else filename = stem if(present(ext)) filename = trim(filename)//ext end if if(.not. tnext) then inquire(file=filename, exist=exists) if(.not. exists) then inquire(file=stem, exist=exists) if(exists) then filename = stem if(present(ext)) filename = trim(filename)//ext endif end if end if if(istart < 0) then call append_ext(stem, abs(i + 1), filename) if(present(ext)) filename = trim(filename)//ext end if end subroutine get_unique_filename function get_free_unit() result(free_unit) ! Returns: ! The first free file unit above 10 and less than or equal to ! the paramater max_unit (currently set to 200). ! ! If max_unit is exceeded, the function returns -1 integer, parameter :: max_unit = 100 integer :: free_unit integer :: i logical :: t_open, t_exist free_unit = -1 do i = 10, max_unit inquire(unit=i, opened=t_open, exist=t_exist) if(.not. t_open .and. t_exist) then free_unit = i exit end if end do if(i == max_unit + 1) call stop_all('get_free_unit', 'Cannot find a free unit below max_unit.') end function get_free_unit function error_function_c(argument) result(res) real(dp), intent(in) :: argument real(dp) :: res res = erfc_local(real(argument, c_double)) end function error_function_c function error_function(argument) result(res) real(dp), intent(in) :: argument real(dp) :: res res = erf_local(real(argument, c_double)) end function error_function pure subroutine find_next_comb(comb, k, n, finish) integer, intent(in) :: k, n integer, intent(inout) :: comb(k) logical, intent(out) :: finish integer :: i if(k == 0 .or. n == 0) then finish = .true. return else if(comb(1) > n - k) then finish = .true. return else finish = .false. end if i = k comb(i) = comb(i) + 1 do if(i < 1 .or. comb(i) < n - k + i + 1) exit i = i - 1 comb(i) = comb(i) + 1 end do do i = i + 1, k comb(i) = comb(i - 1) + 1 end do end subroutine find_next_comb function neci_etime(time) result(ret) ! Return elapsed time for timing and calculation ending purposes. real(dp), intent(out) :: time(2) real(dp) :: ret #if defined(IFORT_) || defined(INTELLLVM_) ! intels etime takes a real(4) real(sp) :: ioTime(2) ! Ifort defines etime directly in its compatibility modules. ! Avoid timing inaccuracies from using cpu_time on cerebro. ret = real(etime(ioTime), dp) time = real(ioTime, dp) #else #ifdef BLUEGENE_HACKS time = 0.0_dp ret = 0.0_dp #else ! use MPI_WTIME - etime returns wall-clock time on multi-processor ! environments, so keep it consistent ret = MPI_WTIME() time(1) = ret time(2) = 0._dp #endif #endif end function neci_etime subroutine open_new_file(funit, filename) integer, intent(in) :: funit character(*), intent(in) :: filename logical :: exists integer :: ierr, i character(43) :: filename2 character(12) :: num character(*), parameter :: t_r = 'open_new_file' ! If we are doing a normal calculation, move existing fciqmc_stats ! files so that they are not overwritten, and then create a new one inquire(file=filename, exist=exists) if(exists) then ! Loop until we find an available spot to move the existing ! file to. i = 1 do while(exists) write(num, '(i12)') i filename2 = trim(adjustl(filename))//"."// & trim(adjustl(num)) inquire(file=filename2, exist=exists) if(i > 10000) & call stop_all(t_r, 'Error finding free fciqmc_stats.*') i = i + 1 end do ! Move the file call rename(filename, filename2) end if ! And finally open the file open(funit, file=filename, status='unknown', iostat=ierr) end subroutine open_new_file pure function cumsum_integer_int32(X) result(Y) integer(int32), intent(in) :: X(:) integer(int32) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function cumsum_integer_int64(X) result(Y) integer(int64), intent(in) :: X(:) integer(int64) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function cumsum_real_sp(X) result(Y) real(sp), intent(in) :: X(:) real(sp) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function cumsum_real_dp(X) result(Y) real(dp), intent(in) :: X(:) real(dp) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function cumsum_complex_sp(X) result(Y) complex(sp), intent(in) :: X(:) complex(sp) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function cumsum_complex_dp(X) result(Y) complex(dp), intent(in) :: X(:) complex(dp) :: Y(size(X)) integer :: i if(size(X) > 0) then Y(1) = X(1) do i = 2, size(Y) Y(i) = Y(i - 1) + X(i) end do end if end function pure function lex_leq(lhs, rhs) result(res) integer, intent(in) :: lhs(:), rhs(size(lhs)) logical :: res integer :: i res = .true. do i = 1, size(lhs) if (lhs(i) == rhs(i)) then cycle else if (lhs(i) < rhs(i)) then return else if (lhs(i) > rhs(i)) then res = .false. return end if end do end function pure function lex_geq(lhs, rhs) result(res) integer, intent(in) :: lhs(:), rhs(size(lhs)) logical :: res integer :: i res = .true. do i = 1, size(lhs) if (lhs(i) == rhs(i)) then cycle else if (lhs(i) > rhs(i)) then return else if (lhs(i) < rhs(i)) then res = .false. return end if end do end function !> @brief !> Create all possible permutations of [1, ..., n] pure function get_permutations(n) result(res) integer, intent(in) :: n integer :: res(n, factrl(n)) integer :: tmp(n), i, j, f tmp = [(i, i = 1, n)] res(:, 1) = tmp do f = 2, size(res, 2) i = 2 do while (tmp(i - 1) > tmp(i)) i = i + 1 end do j = 1 do while (tmp(j) > tmp(i)) j = j + 1 end do call swap(tmp(i), tmp(j)) i = i - 1 j = 1 do while (j < i) call swap(tmp(i), tmp(j)) i = i - 1 j = j + 1 end do res(:, f) = tmp end do end function !> @brief !> The logical operator P => Q !> !> @details !> P | Q | P => Q | ¬ P ∨ Q !> ------------------------------- !> T | T | T | T !> T | F | F | F !> F | T | T | T !> F | F | T | T logical elemental function implies(P, Q) logical, intent(in) :: P, Q implies = .not. P .or. Q end function logical elemental function eq_EnumBase_t(this, other) class(EnumBase_t), intent(in) :: this, other if (.not. SAME_TYPE_AS(this, other)) error stop 'Can only compare objects of same type' eq_EnumBase_t = this%val == other%val end function logical elemental function neq_EnumBase_t(this, other) class(EnumBase_t), intent(in) :: this, other if (.not. SAME_TYPE_AS(this, other)) error stop 'Can only compare objects of same type' neq_EnumBase_t = this%val /= other%val end function elemental function clamp_integer_int32(v, lo, hi) result(res) integer(int32), intent(in) :: v, lo, hi integer(int32) :: res res = merge(lo, merge(hi, v, v > hi), v < lo) end function elemental function clamp_integer_int64(v, lo, hi) result(res) integer(int64), intent(in) :: v, lo, hi integer(int64) :: res res = merge(lo, merge(hi, v, v > hi), v < lo) end function elemental function clamp_real_sp(v, lo, hi) result(res) real(sp), intent(in) :: v, lo, hi real(sp) :: res res = merge(lo, merge(hi, v, v > hi), v < lo) end function elemental function clamp_real_dp(v, lo, hi) result(res) real(dp), intent(in) :: v, lo, hi real(dp) :: res res = merge(lo, merge(hi, v, v > hi), v < lo) end function end module !Hacks for compiler specific system calls. integer function neci_iargc() implicit none integer :: command_argument_count neci_iargc = command_argument_count() end function subroutine neci_getarg(i, str) #ifdef NAGF95 use f90_unix_env, only: getarg #endif implicit none integer, intent(in) :: i character(len=*), intent(out) :: str #if defined(__OPEN64__) || defined(__PATHSCALE__) integer(int32) :: j #else integer :: j #endif #ifdef WARNING_WORKAROUND_ j = i #endif #if defined NAGF95 call getarg(i, str) #elif defined(BLUEGENE_HACKS) call getarg(int(i, 4), str) #elif defined(__OPEN64__) || defined(__PATHSCALE__) j = i call get_command_argument(j, str) #else call get_command_argument(i, str) #endif end subroutine neci_getarg ! Hacks for the IBM compilers on BlueGenes. ! --> The compiler intrinsics are provided as flush_, etime_, sleep_ etc. ! --> We need to either change the names used in the code, or provide wrappers #ifdef BLUEGENE_HACKS ! I presume that the function cpu_time will work here? ! If not, simply add BLUEGENE_HACKS to the neci_etime above. ! function etime (t) result(ret) ! implicit none ! real(4) :: t(2), etime_, ret ! ret = etime_(t) ! end function function hostnm(nm) result(ret) implicit none integer :: ret, hostnm_ character(255) :: nm ret = hostnm_(nm) end function #endif #ifdef CRAY_ETIME function etime(tarr) result(tret) implicit none real(4) :: tarr(2), tret, second tret = second() tarr = tret end function #endif