From 4448d17ee1769c19873d6c96d2d28205d4348414 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 26 Aug 2016 13:32:36 -0500 Subject: [PATCH] Add procedures for sorting and finding --- src/algorithm.F90 | 275 +++++++++++++++++++++++++++++++++++ src/angle_distribution.F90 | 2 +- src/cmfd_execute.F90 | 2 +- src/cross_section.F90 | 2 +- src/eigenvalue.F90 | 2 +- src/endf_header.F90 | 2 +- src/energy_distribution.F90 | 2 +- src/mesh.F90 | 10 +- src/physics.F90 | 2 +- src/scattdata_header.F90 | 2 +- src/search.F90 | 143 ------------------ src/secondary_correlated.F90 | 2 +- src/secondary_kalbach.F90 | 2 +- src/source.F90 | 13 +- src/tally.F90 | 10 +- src/tally_filter.F90 | 2 +- 16 files changed, 302 insertions(+), 171 deletions(-) create mode 100644 src/algorithm.F90 delete mode 100644 src/search.F90 diff --git a/src/algorithm.F90 b/src/algorithm.F90 new file mode 100644 index 0000000000..101b11b9a7 --- /dev/null +++ b/src/algorithm.F90 @@ -0,0 +1,275 @@ +module algorithm + + use constants + use stl_vector, only: VectorInt, VectorReal + + implicit none + + integer, parameter :: MAX_ITERATION = 64 + + interface binary_search + module procedure binary_search_real, binary_search_int4, binary_search_int8 + end interface binary_search + + interface sort + module procedure sort_int, sort_real, sort_vector_int, sort_vector_real + end interface sort + + interface find + module procedure find_int, find_real, find_vector_int, find_vector_real + end interface find + +contains + +!=============================================================================== +! BINARY_SEARCH performs a binary search of an array to find where a specific +! value lies in the array. This is used extensively for energy grid searching +!=============================================================================== + + pure function binary_search_real(array, n, val) result(array_index) + + integer, intent(in) :: n + real(8), intent(in) :: array(n) + real(8), intent(in) :: val + integer :: array_index + + integer :: L + integer :: R + integer :: n_iteration + + L = 1 + R = n + + if (val < array(L) .or. val > array(R)) then + array_index = -1 + return + end if + + n_iteration = 0 + do while (R - L > 1) + ! Find values at midpoint + array_index = L + (R - L)/2 + if (val >= array(array_index)) then + L = array_index + else + R = array_index + end if + + ! check for large number of iterations + n_iteration = n_iteration + 1 + if (n_iteration == MAX_ITERATION) then + array_index = -2 + return + end if + end do + + array_index = L + + end function binary_search_real + + pure function binary_search_int4(array, n, val) result(array_index) + + integer, intent(in) :: n + integer, intent(in) :: array(n) + integer, intent(in) :: val + integer :: array_index + + integer :: L + integer :: R + integer :: n_iteration + + L = 1 + R = n + + if (val < array(L) .or. val > array(R)) then + array_index = -1 + return + end if + + n_iteration = 0 + do while (R - L > 1) + ! Find values at midpoint + array_index = L + (R - L)/2 + if (val >= array(array_index)) then + L = array_index + else + R = array_index + end if + + ! check for large number of iterations + n_iteration = n_iteration + 1 + if (n_iteration == MAX_ITERATION) then + array_index = -2 + return + end if + end do + + array_index = L + + end function binary_search_int4 + + pure function binary_search_int8(array, n, val) result(array_index) + + integer, intent(in) :: n + integer(8), intent(in) :: array(n) + integer(8), intent(in) :: val + integer :: array_index + + integer :: L + integer :: R + integer :: n_iteration + + L = 1 + R = n + + if (val < array(L) .or. val > array(R)) then + array_index = -1 + return + end if + + n_iteration = 0 + do while (R - L > 1) + ! Find values at midpoint + array_index = L + (R - L)/2 + if (val >= array(array_index)) then + L = array_index + else + R = array_index + end if + + ! check for large number of iterations + n_iteration = n_iteration + 1 + if (n_iteration == MAX_ITERATION) then + array_index = -2 + return + end if + end do + + array_index = L + + end function binary_search_int8 + +!=============================================================================== +! SORT sorts an array in place using an insertion sort. +!=============================================================================== + + pure subroutine sort_int(array) + integer, intent(inout) :: array(:) + + integer :: k, m + integer :: temp + + if (size(array) > 1) then + SORT: do k = 2, size(array) + ! Save value to move + m = k + temp = array(k) + + MOVE_OVER: do while (m > 1) + ! Check if insertion value is greater than (m-1)th value + if (temp >= array(m - 1)) exit + + ! Move values over until hitting one that's not larger + array(m) = array(m - 1) + m = m - 1 + end do MOVE_OVER + + ! Put the original value into its new position + array(m) = temp + end do SORT + end if + end subroutine sort_int + + pure subroutine sort_real(array) + real(8), intent(inout) :: array(:) + + integer :: k, m + real(8) :: temp + + if (size(array) > 1) then + SORT: do k = 2, size(array) + ! Save value to move + m = k + temp = array(k) + + MOVE_OVER: do while (m > 1) + ! Check if insertion value is greater than (m-1)th value + if (temp >= array(m - 1)) exit + + ! Move values over until hitting one that's not larger + array(m) = array(m - 1) + m = m - 1 + end do MOVE_OVER + + ! Put the original value into its new position + array(m) = temp + end do SORT + end if + end subroutine sort_real + + pure subroutine sort_vector_int(vec) + type(VectorInt), intent(inout) :: vec + + call sort_int(vec % data(1:vec%size())) + end subroutine sort_vector_int + + pure subroutine sort_vector_real(vec) + type(VectorReal), intent(inout) :: vec + + call sort_real(vec % data(1:vec%size())) + end subroutine sort_vector_real + +!=============================================================================== +! FIND determines the index of the first occurrence of a value in an array. If +! the value does not appear in the array, -1 is returned. +!=============================================================================== + + pure function find_int(array, val) result(index) + integer, intent(in) :: array(:) + integer, intent(in) :: val + integer :: index + + integer :: i + + index = -1 + do i = 1, size(array) + if (array(i) == val) then + index = i + exit + end if + end do + end function find_int + + pure function find_real(array, val) result(index) + real(8), intent(in) :: array(:) + real(8), intent(in) :: val + integer :: index + + integer :: i + + index = -1 + do i = 1, size(array) + if (array(i) == val) then + index = i + exit + end if + end do + end function find_real + + pure function find_vector_int(vec, val) result(index) + type(VectorInt), intent(in) :: vec + integer, intent(in) :: val + integer :: index + + index = find_int(vec % data(1:vec % size()), val) + end function find_vector_int + + pure function find_vector_real(vec, val) result(index) + type(VectorReal), intent(in) :: vec + real(8), intent(in) :: val + integer :: index + + index = find_real(vec % data(1:vec % size()), val) + end function find_vector_real + +end module algorithm diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 index a4fea6ff77..5d16f74242 100644 --- a/src/angle_distribution.F90 +++ b/src/angle_distribution.F90 @@ -2,12 +2,12 @@ module angle_distribution use hdf5, only: HID_T, HSIZE_T + use algorithm, only: binary_search use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR use distribution_univariate, only: DistributionContainer, Tabular use hdf5_interface, only: read_attribute, get_shape, read_dataset, & open_dataset, close_dataset use random_lcg, only: prn - use search, only: binary_search implicit none private diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index b7d0cc3879..d2631254b1 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -213,13 +213,13 @@ contains subroutine cmfd_reweight(new_weights) + use algorithm, only: binary_search use constants, only: ZERO, ONE use error, only: warning, fatal_error use global, only: meshes, source_bank, work, n_user_meshes, cmfd, & master use mesh_header, only: RegularMesh use mesh, only: count_bank_sites, get_mesh_indices - use search, only: binary_search use string, only: to_str #ifdef MPI diff --git a/src/cross_section.F90 b/src/cross_section.F90 index a8acb25a83..b67226a1b8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -1,5 +1,6 @@ module cross_section + use algorithm, only: binary_search use constants use energy_grid, only: grid_method, log_spacing use error, only: fatal_error @@ -14,7 +15,6 @@ module cross_section use particle_header, only: Particle use random_lcg, only: prn, future_prn, prn_set_stream use sab_header, only: SAlphaBeta - use search, only: binary_search implicit none diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 9befbe3c3d..713bbc351a 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -4,6 +4,7 @@ module eigenvalue use message_passing #endif + use algorithm, only: binary_search use constants, only: ZERO use error, only: fatal_error, warning use global @@ -11,7 +12,6 @@ module eigenvalue use mesh, only: count_bank_sites use mesh_header, only: RegularMesh use random_lcg, only: prn, set_particle_seed, advance_prn_seed - use search, only: binary_search use string, only: to_str implicit none diff --git a/src/endf_header.F90 b/src/endf_header.F90 index 8d8aefaa3c..e9e45ab751 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -2,10 +2,10 @@ module endf_header use hdf5, only: HID_T, HSIZE_T + use algorithm, only: binary_search use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, & LOG_LINEAR, LOG_LOG use hdf5_interface - use search, only: binary_search implicit none diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index c45762bb0a..770da617cc 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -2,12 +2,12 @@ module energy_distribution use hdf5 + use algorithm, only: binary_search use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR use endf_header, only: Tabulated1D use hdf5_interface use math, only: maxwell_spectrum, watt_spectrum use random_lcg, only: prn - use search, only: binary_search !=============================================================================== ! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a diff --git a/src/mesh.F90 b/src/mesh.F90 index cee29aa1ca..4a3def4282 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,14 +1,14 @@ module mesh - use constants - use global - use mesh_header - use search, only: binary_search - #ifdef MPI use message_passing #endif + use algorithm, only: binary_search + use constants + use global + use mesh_header + implicit none contains diff --git a/src/physics.F90 b/src/physics.F90 index a9f2263955..00bd2c20c9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,5 +1,6 @@ module physics + use algorithm, only: binary_search use constants use cross_section, only: elastic_xs_0K use endf, only: reaction_name @@ -15,7 +16,6 @@ module physics use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction - use search, only: binary_search use secondary_uncorrelated, only: UncorrelatedAngleEnergy use string, only: to_str diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 2066b36b92..d4643a0721 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -1,10 +1,10 @@ module scattdata_header + use algorithm, only: binary_search use constants use error, only: fatal_error use math use random_lcg, only: prn - use search, only: binary_search implicit none diff --git a/src/search.F90 b/src/search.F90 deleted file mode 100644 index f338105f4d..0000000000 --- a/src/search.F90 +++ /dev/null @@ -1,143 +0,0 @@ -module search - - use constants - - implicit none - - integer, parameter :: MAX_ITERATION = 64 - - interface binary_search - module procedure binary_search_real, binary_search_int4, binary_search_int8 - end interface binary_search - -contains - -!=============================================================================== -! BINARY_SEARCH performs a binary search of an array to find where a specific -! value lies in the array. This is used extensively for energy grid searching -!=============================================================================== - - pure function binary_search_real(array, n, val) result(array_index) - - integer, intent(in) :: n - real(8), intent(in) :: array(n) - real(8), intent(in) :: val - integer :: array_index - - integer :: L - integer :: R - integer :: n_iteration - - L = 1 - R = n - - if (val < array(L) .or. val > array(R)) then - array_index = -1 - return - end if - - n_iteration = 0 - do while (R - L > 1) - ! Find values at midpoint - array_index = L + (R - L)/2 - if (val >= array(array_index)) then - L = array_index - else - R = array_index - end if - - ! check for large number of iterations - n_iteration = n_iteration + 1 - if (n_iteration == MAX_ITERATION) then - array_index = -2 - return - end if - end do - - array_index = L - - end function binary_search_real - - pure function binary_search_int4(array, n, val) result(array_index) - - integer, intent(in) :: n - integer, intent(in) :: array(n) - integer, intent(in) :: val - integer :: array_index - - integer :: L - integer :: R - integer :: n_iteration - - L = 1 - R = n - - if (val < array(L) .or. val > array(R)) then - array_index = -1 - return - end if - - n_iteration = 0 - do while (R - L > 1) - ! Find values at midpoint - array_index = L + (R - L)/2 - if (val >= array(array_index)) then - L = array_index - else - R = array_index - end if - - ! check for large number of iterations - n_iteration = n_iteration + 1 - if (n_iteration == MAX_ITERATION) then - array_index = -2 - return - end if - end do - - array_index = L - - end function binary_search_int4 - - pure function binary_search_int8(array, n, val) result(array_index) - - integer, intent(in) :: n - integer(8), intent(in) :: array(n) - integer(8), intent(in) :: val - integer :: array_index - - integer :: L - integer :: R - integer :: n_iteration - - L = 1 - R = n - - if (val < array(L) .or. val > array(R)) then - array_index = -1 - return - end if - - n_iteration = 0 - do while (R - L > 1) - ! Find values at midpoint - array_index = L + (R - L)/2 - if (val >= array(array_index)) then - L = array_index - else - R = array_index - end if - - ! check for large number of iterations - n_iteration = n_iteration + 1 - if (n_iteration == MAX_ITERATION) then - array_index = -2 - return - end if - end do - - array_index = L - - end function binary_search_int8 - -end module search diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 index e163fdcc24..a0e203f33d 100644 --- a/src/secondary_correlated.F90 +++ b/src/secondary_correlated.F90 @@ -2,13 +2,13 @@ module secondary_correlated use hdf5, only: HID_T, HSIZE_T + use algorithm, only: binary_search use angleenergy_header, only: AngleEnergy use constants, only: ZERO, ONE, HALF, TWO, HISTOGRAM, LINEAR_LINEAR use distribution_univariate, only: DistributionContainer, Tabular use hdf5_interface, only: get_shape, read_attribute, open_dataset, & read_dataset, close_dataset use random_lcg, only: prn - use search, only: binary_search !=============================================================================== ! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 index 4b5e690b5d..f963cff3ff 100644 --- a/src/secondary_kalbach.F90 +++ b/src/secondary_kalbach.F90 @@ -2,12 +2,12 @@ module secondary_kalbach use hdf5, only: HID_T, HSIZE_T + use algorithm, only: binary_search use angleenergy_header, only: AngleEnergy use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR use hdf5_interface, only: read_attribute, read_dataset, open_dataset, & close_dataset, get_shape use random_lcg, only: prn - use search, only: binary_search !=============================================================================== ! KalbachMann represents a correlated angle-energy distribution with the angular diff --git a/src/source.F90 b/src/source.F90 index 452d8ddfce..9aeccde156 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,5 +1,11 @@ module source + use hdf5, only: HID_T +#ifdef MPI + use message_passing +#endif + + use algorithm, only: binary_search use bank_header, only: Bank use constants use distribution_univariate, only: Discrete @@ -12,17 +18,10 @@ module source use output, only: write_message use particle_header, only: Particle use random_lcg, only: prn, set_particle_seed, prn_set_stream - use search, only: binary_search use string, only: to_str use math use state_point, only: read_source_bank, write_source_bank -#ifdef MPI - use message_passing -#endif - - use hdf5, only: HID_T - implicit none contains diff --git a/src/tally.F90 b/src/tally.F90 index 3c3dc6a693..3ee4702e2f 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,5 +1,10 @@ module tally +#ifdef MPI + use message_passing +#endif + + use algorithm, only: binary_search use constants use error, only: fatal_error use geometry_header @@ -11,15 +16,10 @@ module tally use mesh_header, only: RegularMesh use output, only: header use particle_header, only: LocalCoord, Particle - use search, only: binary_search use string, only: to_str use tally_header, only: TallyResult use tally_filter -#ifdef MPI - use message_passing -#endif - implicit none integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 67d0452847..bd568e29ae 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -1,5 +1,6 @@ module tally_filter + use algorithm, only: binary_search use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION use dict_header, only: DictIntInt use geometry_header, only: BASE_UNIVERSE, RectLattice, HexLattice @@ -10,7 +11,6 @@ module tally_filter get_mesh_indices, mesh_indices_to_bin, & mesh_intersects_2d, mesh_intersects_3d use particle_header, only: Particle - use search, only: binary_search use string, only: to_str use tally_filter_header, only: TallyFilter, TallyFilterContainer