From 06af75dd928eaf0cc4904dd61a7b2af149356a8f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 6 Oct 2021 10:48:18 -0400 Subject: [PATCH] Consolidated triangular encounter checks into a single subroutine that gets called by each of the variations. Still trying to track down the problem that is preventing parallel speedup. --- src/encounter/encounter_check.f90 | 255 +++++++++++++----------------- src/encounter/encounter_setup.f90 | 21 +-- src/modules/encounter_classes.f90 | 3 +- src/modules/swiftest_classes.f90 | 5 + src/util/util_index_array.f90 | 38 +++++ 5 files changed, 163 insertions(+), 159 deletions(-) create mode 100644 src/util/util_index_array.f90 diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 516fb84a9..d869e242e 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -375,8 +375,8 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, npl_last = npl end if - ! call timer%reset() - ! call timer%start() +! call timer%reset() +! call timer%start() !$omp parallel do default(private) schedule(static) & !$omp shared(x, v, renc, boundingbox) & !$omp firstprivate(dt, npl, n) @@ -392,8 +392,8 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt]) end do !$omp end parallel do - ! call timer%stop() - ! call timer%report(nsubsteps=nthreads, message="AABB sort plpl:") +! call timer%stop() +! call timer%report(nsubsteps=nthreads, message="AABB sort plpl:") call boundingbox%sweep(npl, nenc, index1, index2) @@ -454,8 +454,8 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt ntot_last = ntot end if - ! call timer%reset() - ! call timer%start() + call timer%reset() + call timer%start() !$omp parallel do default(private) schedule(static) & !$omp shared(xplm, xplt, vplm, vplt, rencm, renct, boundingbox) & !$omp firstprivate(dt, nplm, nplt, ntot) @@ -596,6 +596,48 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, end subroutine encounter_check_all_sort_and_sweep_pltp + subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, ind_arr, lenci) + !$omp declare simd(encounter_check_all_triangular_one) + implicit none + ! Arguments + integer(I4B), intent(in) :: i + integer(I4B), intent(in) :: n + real(DP), intent(in) :: xi, yi, zi + real(DP), intent(in) :: vxi, vyi, vzi + real(DP), dimension(:), intent(in) :: x, y, z + real(DP), dimension(:), intent(in) :: vx, vy, vz + real(DP), intent(in) :: renci + real(DP), dimension(:), intent(in) :: renc + real(DP), intent(in) :: dt + integer(I4B), dimension(:), intent(in) :: ind_arr + type(encounter_list), intent(out) :: lenci + ! Internals + integer(I4B) :: j, nenci + real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 + logical, dimension(n) :: lencounteri, lvdotri + + do concurrent(j = i+1:n) + xr = x(j) - xi + yr = y(j) - yi + zr = z(j) - zi + vxr = vx(j) - vxi + vyr = vy(j) - vyi + vzr = vz(j) - vzi + renc12 = renci + renc(j) + call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) + end do + nenci = count(lencounteri(i+1:n)) + if (nenci > 0) then + allocate(lenci%lvdotr(nenci), lenci%index2(nenci)) + lenci%nenc = nenci + lenci%lvdotr(:) = pack(lvdotri(i+1:n), lencounteri(i+1:n)) + lenci%index2(:) = pack(ind_arr(i+1:n), lencounteri(i+1:n)) + end if + + return + end subroutine encounter_check_all_triangular_one + + subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! @@ -613,55 +655,27 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter integer(I4B), intent(out) :: nenc !! Total number of encounters ! Internals - integer(I4B) :: i, j, nenci, j0, j1 + integer(I4B) :: i, j, k, nenci, j0, j1 real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(npl) :: lencounteri, lvdotri - integer(I4B), dimension(npl) :: ind_arr + integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(npl) :: lenc + + call util_index_array(ind_arr, npl) - ind_arr(:) = [(i, i = 1, npl)] !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, x, v, renc, dt, lenc, ind_arr) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, nenci, lencounteri, lvdotri) - do i = 1, npl - do concurrent(j = i+1:npl) - xr = x(1, j) - x(1, i) - yr = x(2, j) - x(2, i) - zr = x(3, j) - x(3, i) - vxr = v(1, j) - v(1, i) - vyr = v(2, j) - v(2, i) - vzr = v(3, j) - v(3, i) - renc12 = renc(i) + renc(j) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) - end do - nenci = count(lencounteri(i+1:npl)) - if (nenci > 0) then - allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) - lenc(i)%nenc = nenci - lenc(i)%lvdotr(:) = pack(lvdotri(i+1:npl), lencounteri(i+1:npl)) - lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl)) - end if + !$omp shared(x, v, renc, lenc, ind_arr) & + !$omp firstprivate(npl, dt) + do i = 1,npl + call encounter_check_all_triangular_one(i, npl, x(1,i), x(2,i), x(3,i), & + v(1,i), v(2,i), v(3,i), & + x(1,:), x(2,:), x(3,:), & + v(1,:), v(2,:), v(3,:), & + renc(i), renc(:), dt, ind_arr(:), lenc(i)) end do !$omp end parallel do - associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:npl)) - end associate - if (nenc > 0) then - allocate(lvdotr(nenc)) - allocate(index1(nenc)) - allocate(index2(nenc)) - j0 = 1 - do i = 1, npl - if (lenc(i)%nenc > 0) then - j1 = j0 + lenc(i)%nenc - 1 - lvdotr(j0:j1) = lenc(i)%lvdotr(:) - index1(j0:j1) = i - index2(j0:j1) = lenc(i)%index2(:) - j0 = j1 + 1 - end if - end do - end if + call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_triangular_plpl @@ -691,52 +705,24 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp integer(I4B) :: i, j, nenci, j0, j1 real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(nplt) :: lencounteri, lvdotri - integer(I4B), dimension(nplt) :: ind_arr + integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(nplm) :: lenc - ind_arr(:) = [(i, i = 1, nplt)] + call util_index_array(ind_arr, nplt) + !$omp parallel do default(private) schedule(static)& - !$omp shared(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lenc, ind_arr) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, nenci, lencounteri, lvdotri) + !$omp shared(xplm, vplm, xplt, vplt, rencm, renct, lenc, ind_arr) & + !$omp firstprivate(nplm, nplt, dt) do i = 1, nplm - do concurrent(j = 1:nplt) - xr = xplt(1, j) - xplm(1, i) - yr = xplt(2, j) - xplm(2, i) - zr = xplt(3, j) - xplm(3, i) - vxr = vplt(1, j) - vplm(1, i) - vyr = vplt(2, j) - vplm(2, i) - vzr = vplt(3, j) - vplm(3, i) - renc12 = rencm(i) + renct(j) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) - end do - nenci = count(lencounteri(1:nplt)) - if (nenci > 0) then - allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) - lenc(i)%nenc = nenci - lenc(i)%lvdotr(:) = pack(lvdotri(1:nplt), lencounteri(1:nplt)) - lenc(i)%index2(:) = pack(ind_arr(1:nplt), lencounteri(1:nplt)) - end if + call encounter_check_all_triangular_one(0, nplt, xplm(1,i), xplm(2,i), xplm(3,i), & + vplm(1,i), vplm(2,i), vplm(3,i), & + xplt(1,:), xplt(2,:), xplt(3,:), & + vplt(1,:), vplt(2,:), vplt(3,:), & + rencm(i), renct(:), dt, ind_arr(:), lenc(i)) end do !$omp end parallel do - associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:nplm)) - end associate - if (nenc > 0) then - allocate(lvdotr(nenc)) - allocate(index1(nenc)) - allocate(index2(nenc)) - j0 = 1 - do i = 1, nplm - if (lenc(i)%nenc > 0) then - j1 = j0 + lenc(i)%nenc - 1 - lvdotr(j0:j1) = lenc(i)%lvdotr(:) - index1(j0:j1) = i - index2(j0:j1) = lenc(i)%index2(:) + nplm - j0 = j1 + 1 - end if - end do - end if + call encounter_check_collapse_ragged_list(lenc, nplm, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_triangular_plplm @@ -765,51 +751,26 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren integer(I4B) :: i, j, nenci, j0, j1 real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc1, renc2 logical, dimension(ntp) :: lencounteri, lvdotri - integer(I4B), dimension(ntp) :: ind_arr + integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(npl) :: lenc + real(DP), dimension(ntp) :: renct + + call util_index_array(ind_arr, ntp) + renct(:) = 0.0_DP - ind_arr(:) = [(i, i = 1, ntp)] !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lenc, ind_arr) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, nenci, lencounteri, lvdotri) + !$omp shared(xpl, vpl, xtp, vtp, renc, renct, lenc, ind_arr) & + !$omp firstprivate(npl, ntp, dt) do i = 1, npl - do concurrent(j = 1:ntp) - xr = xtp(1, j) - xpl(1, i) - yr = xtp(2, j) - xpl(2, i) - zr = xtp(3, j) - xpl(3, i) - vxr = vtp(1, j) - vpl(1, i) - vyr = vtp(2, j) - vpl(2, i) - vzr = vtp(3, j) - vpl(3, i) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc(i), dt, lencounteri(j), lvdotri(j)) - end do - nenci = count(lencounteri(1:ntp)) - if (nenci > 0) then - allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) - lenc(i)%nenc = nenci - lenc(i)%lvdotr(:) = pack(lvdotri(1:ntp), lencounteri(1:ntp)) - lenc(i)%index2(:) = pack(ind_arr(1:ntp), lencounteri(1:ntp)) - end if + call encounter_check_all_triangular_one(0, ntp, xpl(1,i), xpl(2,i), xpl(3,i), & + vpl(1,i), vpl(2,i), vpl(3,i), & + xtp(1,:), xtp(2,:), xtp(3,:), & + vtp(1,:), vtp(2,:), vtp(3,:), & + renc(i), renct(:), dt, ind_arr(:), lenc(i)) end do !$omp end parallel do - associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:npl)) - end associate - if (nenc > 0) then - allocate(lvdotr(nenc)) - allocate(index1(nenc)) - allocate(index2(nenc)) - j0 = 1 - do i = 1, npl - if (lenc(i)%nenc > 0) then - j1 = j0 + lenc(i)%nenc - 1 - lvdotr(j0:j1) = lenc(i)%lvdotr(:) - index1(j0:j1) = i - index2(j0:j1) = lenc(i)%index2(:) - j0 = j1 + 1 - end if - end do - end if + call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_triangular_pltp @@ -836,22 +797,24 @@ module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, r2 = xr**2 + yr**2 + zr**2 r2crit = renc**2 - lencounter = (r2 < r2crit) - if (lencounter) return - - vdotr = vxr * xr + vyr * yr + vzr * zr - lvdotr = (vdotr < 0.0_DP) - if (.not.lvdotr) return - - v2 = vxr**2 + vyr**2 + vzr**2 - tmin = -vdotr / v2 - if (tmin < dt) then - r2min = r2 - vdotr**2 / v2 + if (r2 > r2crit) then + vdotr = vxr * xr + vyr * yr + vzr * zr + v2 = vxr**2 + vyr**2 + vzr**2 + tmin = -vdotr / v2 + + if (tmin < dt) then + r2min = r2 - vdotr**2 / v2 + else + r2min = r2 + 2 * vdotr * dt + v2 * dt**2 + end if else - r2min = r2 + 2 * vdotr * dt + v2 * dt**2 + vdotr = 0.0_DP + r2min = r2 end if - lencounter = (r2min <= r2crit) + + lencounter = (r2min <= r2crit) + lvdotr = (vdotr < 0.0_DP) return end subroutine encounter_check_one @@ -868,7 +831,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in integer(I4B), intent(out) :: nenc !! Total number of encountersj integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 - integer(I4B), dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching + logical, dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching ! Internals integer(I4B) :: i, j0, j1, nenci @@ -879,6 +842,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in allocate(index1(nenc)) allocate(index2(nenc)) + if (present(lvdotr)) allocate(lvdotr(nenc)) j0 = 1 do i = 1, n1 nenci = ragged_list(i)%nenc @@ -938,20 +902,22 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair !Internals - Integer(I4B) :: i, k, ntot + integer(I4B) :: i, k, ntot type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) type(walltimer) :: timer + integer(I4B), dimension(:), allocatable, save :: ind_arr ntot = n1 + n2 + call util_index_array(ind_arr, ntot) ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure ! call timer%reset() ! call timer%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(self, lenc) & + !$omp shared(self, lenc, ind_arr) & !$omp firstprivate(ntot, n1, n2) do i = 1, ntot - call encounter_check_sweep_aabb_one_double_list(i, n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), self%ind_arr(:), lenc(i)) + call encounter_check_sweep_aabb_one_double_list(i, n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(i)) end do !$omp end parallel do ! call timer%stop() @@ -986,16 +952,19 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1, Integer(I4B) :: i, k type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) type(walltimer) :: timer + integer(I4B), dimension(:), allocatable, save :: ind_arr + + call util_index_array(ind_arr, n) ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure ! call timer%reset() ! call timer%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(self, lenc) & + !$omp shared(self, lenc, ind_arr) & !$omp firstprivate(n) do i = 1, n - call encounter_check_sweep_aabb_one_single_list(i, n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), self%ind_arr(:), lenc(i)) + call encounter_check_sweep_aabb_one_single_list(i, n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(i)) end do !$omp end parallel do ! call timer%stop() diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index 9d747c565..cb2e2f3d8 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -13,30 +13,23 @@ module subroutine encounter_setup_aabb(self, n, n_last) integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called ! Internals integer(I4B) :: next, next_last, k, dim - integer(I4B), dimension(:), allocatable :: tmp + integer(I4B), dimension(:), allocatable :: itmp next = 2 * n next_last = 2 * n_last if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - allocate(tmp(n)) - if (n_last > 0) tmp(1:n_last) = self%ind_arr(1:n_last) - call move_alloc(tmp, self%ind_arr) - self%ind_arr(n_last+1:n) = [(k, k = n_last+1, n)] do dim = 1, SWEEPDIM - allocate(tmp(next)) - if (n_last > 0) tmp(1:next_last) = self%aabb(dim)%ind(1:next_last) - call move_alloc(tmp, self%aabb(dim)%ind) + allocate(itmp(next)) + if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) + call move_alloc(itmp, self%aabb(dim)%ind) self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] end do else ! The number of bodies has gone down. Resize and chop of the old indices - allocate(tmp(n)) - tmp(1:n) = self%ind_arr(1:n) - call move_alloc(tmp, self%ind_arr) do dim = 1, SWEEPDIM - allocate(tmp(next)) - tmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) - call move_alloc(tmp, self%aabb(dim)%ind) + allocate(itmp(next)) + itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) + call move_alloc(itmp, self%aabb(dim)%ind) end do end if diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index f87e35553..96dbb1e83 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -42,7 +42,6 @@ module encounter_classes type encounter_bounding_box type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb - integer(I4B), dimension(:), allocatable :: ind_arr !! Index array of bodies contains procedure :: setup => encounter_setup_aabb !! Setup a new axis-aligned bounding box structure procedure :: sweep_single => encounter_check_sweep_aabb_single_list !! Sweeps the sorted bounding box extents and returns the encounter candidates @@ -135,7 +134,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in integer(I4B), intent(out) :: nenc !! Total number of encountersj integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 - integer(I4B), dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching + logical, dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching end subroutine encounter_check_collapse_ragged_list module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index cd4484f8f..4d2b4f38f 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1408,6 +1408,11 @@ module subroutine util_flatten_eucl_pltp(self, pl, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine + module subroutine util_index_array(ind_arr, n) + implicit none + integer(I4B), dimension(:), allocatable, intent(inout) :: ind_arr !! Index array. Input is a pre-existing index array where n /= size(ind_arr). Output is a new index array ind_arr = [1, 2, ... n] + integer(I4B), intent(in) :: n !! The new size of the index array + end subroutine util_index_array module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1) use lambda_function diff --git a/src/util/util_index_array.f90 b/src/util/util_index_array.f90 new file mode 100644 index 000000000..d6b214267 --- /dev/null +++ b/src/util/util_index_array.f90 @@ -0,0 +1,38 @@ +submodule (swiftest_classes) s_util_index_array + use swiftest +contains + + module subroutine util_index_array(ind_arr, n) + !! author: David A. Minton + !! + !! Creates or resizes an index array of size n where ind_arr = [1, 2, ... n]. + !! This subroutine assumes that if ind_arr is already allocated, it is a pre-existing index array of a different size. + implicit none + ! Arguments + integer(I4B), dimension(:), allocatable, intent(inout) :: ind_arr !! Index array. Input is a pre-existing index array where n /= size(ind_arr). Output is a new index array ind_arr = [1, 2, ... n] + integer(I4B), intent(in) :: n !! The new size of the index array + ! Internals + integer(I4B) :: nold, i + integer(I4B), dimension(:), allocatable :: itmp + + if (allocated(ind_arr)) then + nold = size(ind_arr) + if (nold == n) return ! Nothing to do, so go home + else + nold = 0 + end if + + if (n >= nold) then + allocate(itmp(n)) + if (nold > 0) itmp(1:nold) = ind_arr(1:nold) + itmp(nold+1:n) = [(i, i = nold + 1, n)] + call move_alloc(itmp, ind_arr) + else + itmp(1:n) = ind_arr(1:n) + call move_alloc(itmp, ind_arr) + end if + + return + end subroutine util_index_array + +end submodule s_util_index_array \ No newline at end of file