module list_operations_m !! Routines operating on lists, like sorting, removing duplicates etc... use precision_m, only : FP use screen_io_m, only : get_stderr implicit none save public :: unique_tuples public :: insertion_sort public :: sort_tuples public :: findindex_2d public :: sort_and_sum public :: findloc_i public :: rank_list private :: i_mrgref private :: fp_mrgref interface rank_list module procedure i_mrgref, fp_mrgref end interface contains subroutine unique_tuples(ndim, nz, arr) !! Removes duplicate entries from integer tuple arr integer, intent(in) :: ndim !! Dimension of tuple integer, intent(inout) :: nz !! Number of tuples, on return number of unique tuples integer, intent(inout), dimension(ndim, nz) :: arr !! Tuples, on return unique tuples (nz adapted) integer :: i, j, k, res(ndim, nz) if (nz <= 0) then return endif res = 0 k = 1 res(1:ndim, k) = arr(1:ndim, k) outer: do i = 2, nz do j = 1, k if ( all(arr(:, i) == res(:, j)) ) then cycle outer endif enddo k = k + 1 res(1:ndim, k) = arr(1:ndim, i) end do outer nz = k arr = res end subroutine subroutine insertion_sort(nz, ndim, decrow, arr) !! Sorts integer tuple arr in ascending order according to row decr integer, intent(in) :: ndim !! Dimension of tuples integer, intent(in) :: nz !! Number of tuples integer, intent(in) :: decrow !! Is sorted in ascending order according to row decrow integer, intent(inout), dimension(ndim, nz) :: arr !! Tuples, on return sorted integer :: i, j, val(ndim) if (nz>1) then do i = 2, nz val(:) = arr(:, i) j = i lp: do if (j == 1) exit lp if ( arr(decrow, j-1) <= val(decrow) ) exit lp arr(:,j) = arr(:,j-1) j = j - 1 arr(:, j) = val(:) enddo lp enddo endif end subroutine subroutine sort_tuples(nz, ndim, precedence, arr) !! Sorts integer tuples arr in ascending order, !! where pecedence of rows in sorting is prescribed integer, intent(in) :: ndim !! Dimension of tuples integer, intent(in) :: nz !! Number of tuples integer, intent(in), dimension(ndim) :: precedence !! Precedence of sorting, i.e. for ndim=2 !! precedence=[2,1] sorts firstly according to second row integer, intent(inout), dimension(ndim,nz) :: arr !! Tuples, on return sorted integer :: j, val, kinit, kend, nk call insertion_sort(nz, ndim, precedence(1), arr) if (ndim > 1) then do j = 2, ndim kinit = 1 kend = 0 do while (kend < nz) val = arr(precedence(j-1), kinit) kend = kinit lp: do if (kend == nz) exit lp if ( arr(precedence(j - 1), kend + 1) /= val ) exit lp kend = kend + 1 enddo lp nk = kend - kinit + 1 call insertion_sort(nk, ndim, precedence(j), arr(1, kinit) ) kinit = kend + 1 enddo enddo endif end subroutine subroutine findindex_2d(lpi, lpf, li, lj, i, j, lg) !! Searches in grid lists (li, lj) if it contains point (i, j) !! and returns its index in lg. !! If index cannot be found, returns lg=0 !! Note: (li, lj) shall be uniqued tuple ! jumps left right around initial guess with increasing step size ! --> fast if initial guess of index close to real index integer, intent(in) :: lpi !! Start index of lists li and lj integer, intent(in) :: lpf !! Final index of lists li and lj integer, intent(in) :: i !! Integer coordinate i that is looked for in li integer, intent(in) :: j !! Integer coordinate j that is looked for in lj integer, intent(in), dimension(lpi:lpf) :: li !! Grid list li integer, intent(in), dimension(lpi:lpf) :: lj !! Grid list lj integer, intent(out)::lg !! If present (i,j)=(li(lg), lj(lg)), if not present lg=0 is returned integer, save :: lguess = 0 !! Routine is accelarated by storing a guess lguess from previous call !$OMP THREADPRIVATE(lguess) ! Threadprivate to avoid race condition at call integer :: ig, jg, nlist, sgn, k, lnext logical :: forward, backward, found if (lpi <= 0) then write(get_stderr(),*) & 'error(findindex_2d)::lpi must be greater 0 but is', lpi stop endif nlist = lpf-lpi+1 ! return 0 if list is empty if (nlist <= 0) then lg = 0 lguess = lg return endif lg = lguess if (lg < lpi) then lg = lpi endif if (lg > lpf) then lg = lpf endif forward = .false. backward = .false. found = .false. searchloop: do k = 1, nlist ig = li(lg) jg = lj(lg) if ((i == ig) .and. (j == jg)) then lguess = lg found = .true. exit searchloop else if (forward) then lg = lg + 1 cycle searchloop endif if (backward) then lg = lg - 1 cycle searchloop endif if (mod(k,2)==1) then sgn = 1 else sgn = -1 endif lnext = lg+sgn*k if (lnext < lpi) then lg = lg+1 forward = .true. elseif (lnext > lpf) then lg = lg-1 backward = .true. else lg = lnext endif endif enddo searchloop if (.not.found) then lg = 0 endif end subroutine subroutine sort_row(size, col_ind, values) !! Simple bubble sort for small arrays. !! Different from the sort_and_sum function below !! as it does not include the summed-up part. implicit none integer, intent(in) :: size !! Dimension of the arrays integer, intent(inout) :: col_ind(size) !! Column index of the input and the output sorted array real(FP), intent(inout) :: values(size) !! Element value integer :: i, j, temp_ind real(FP) :: temp_value do i = 1, size-1 do j = i+1, size if (col_ind(i) > col_ind(j)) then temp_ind = col_ind(i) col_ind(i) = col_ind(j) col_ind(j) = temp_ind temp_value = values(i) values(i) = values(j) values(j) = temp_value end if end do end do end subroutine sort_row subroutine sort_and_sum(nz, j, val) !! Uniques and sorts values according to j-ascending ordering !! Same j-values are summed up in val integer, intent(inout) :: nz !! On input: Dimension of j and val !! On output: Number of unique integer, intent(inout), dimension(nz) :: j !! Integer list, on output uniqued and sorted in ascending order real(FP), intent(inout), dimension(nz) :: val !! Values, on output sorted and summed according to new j-ordering logical :: mask(nz) integer :: jmin, k, r, jnew(nz), nzn, cnt real(FP), dimension(nz) :: valnew mask = .true. jnew = 0 valnew = 0.0_FP cnt = 0 r = 0 nzn = nz do while (r < nzn) r = r + 1 valnew(r) = 0 jmin = minval(j, mask) cnt = 0 do k = 1, nz if (j(k) == jmin) then mask(k) = .false. valnew(r) = valnew(r) + val(k) cnt = cnt + 1 endif enddo nzn = nzn - cnt + 1 jnew(r) = jmin enddo nz = nzn j = jnew val = valnew end subroutine #if defined __GFORTRAN__ && __GNUC__ < 9 ! Implementation of findloc only available for gcc from version 9 integer function findloc_i(arr, val) !! Finds first index i where arr(i) = val. !! Returns 0 if val not found or size of array is 0 !! !! Explicit implementation integer, dimension(:), intent(in) :: arr !! Array integer, intent(in) :: val !! Value to be searched for integer :: ndim, i ndim = size(arr) findloc_i = 0 do i = 1, ndim if (arr(i)==val) then findloc_i = i return endif enddo end function #else integer function findloc_i(arr, val) !! Finds first index i where arr(i) = val. !! Returns 0 if val not found or size of array is 0 !! !! Based on Fortran intrinsic findloc implementation integer, dimension(:), intent(in) :: arr !! Array integer, intent(in) :: val !! Value to be searched for integer, dimension(1) :: res_tmp res_tmp = findloc(arr, val) findloc_i = res_tmp(1) end function #endif Subroutine i_mrgref (XVALT, IRNGT) !! Ranks array XVALT into index array IRNGT, using merge-sort ! This routine is taken from ORDERPACK 2.0 (see http://www.fortran-2000.com/rank/#3.2) ! and therefore left in the original formatting. ! __________________________________________________________ ! This version is not optimized for performance, and is thus ! not as difficult to read as some other ones. ! Michel Olagnon - April 2000 ! __________________________________________________________ ! __________________________________________________________ Integer, Dimension (:), Intent (In) :: XVALT !! Integer, Dimension (:), Intent (Out) :: IRNGT ! __________________________________________________________ ! Integer, Dimension (:), Allocatable :: JWRKT Integer :: LMTNA, LMTNC Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! NVAL = Min (SIZE(XVALT), SIZE(IRNGT)) If (NVAL <= 0) Then Return End If ! ! Fill-in the index array, creating ordered couples ! Do IIND = 2, NVAL, 2 If (XVALT(IIND-1) < XVALT(IIND)) Then IRNGT (IIND-1) = IIND - 1 IRNGT (IIND) = IIND Else IRNGT (IIND-1) = IIND IRNGT (IIND) = IIND - 1 End If End Do If (Modulo (NVAL, 2) /= 0) Then IRNGT (NVAL) = NVAL End If ! ! We will now have ordered subsets A - B - A - B - ... ! and merge A and B couples into C - C - ... ! Allocate (JWRKT(1:NVAL)) LMTNC = 2 LMTNA = 2 ! ! Iteration. Each time, the length of the ordered subsets ! is doubled. ! Do If (LMTNA >= NVAL) Exit IWRKF = 0 LMTNC = 2 * LMTNC IWRK = 0 ! ! Loop on merges of A and B into C ! Do IINDA = IWRKF IWRKD = IWRKF + 1 IWRKF = IINDA + LMTNC JINDA = IINDA + LMTNA If (IWRKF >= NVAL) Then If (JINDA >= NVAL) Exit IWRKF = NVAL End If IINDB = JINDA ! ! Shortcut for the case when the max of A is smaller ! than the min of B (no need to do anything) ! If (XVALT(IRNGT(JINDA)) <= XVALT(IRNGT(JINDA+1))) Then IWRK = IWRKF Cycle End If ! ! One steps in the C subset, that we create in the final rank array ! Do If (IWRK >= IWRKF) Then ! ! Make a copy of the rank array for next iteration ! IRNGT (IWRKD:IWRKF) = JWRKT (IWRKD:IWRKF) Exit End If ! IWRK = IWRK + 1 ! ! We still have unprocessed values in both A and B ! If (IINDA < JINDA) Then If (IINDB < IWRKF) Then If (XVALT(IRNGT(IINDA+1)) > XVALT(IRNGT(IINDB+1))) & & Then IINDB = IINDB + 1 JWRKT (IWRK) = IRNGT (IINDB) Else IINDA = IINDA + 1 JWRKT (IWRK) = IRNGT (IINDA) End If Else ! ! Only A still with unprocessed values ! IINDA = IINDA + 1 JWRKT (IWRK) = IRNGT (IINDA) End If Else ! ! Only B still with unprocessed values ! IRNGT (IWRKD:IINDB) = JWRKT (IWRKD:IINDB) IWRK = IWRKF Exit End If ! End Do End Do ! ! The Cs become As and Bs ! LMTNA = 2 * LMTNA End Do ! ! Clean up ! Deallocate (JWRKT) Return ! End Subroutine i_mrgref Subroutine fp_mrgref (XVALT, IRNGT) !! Ranks array XVALT into index array IRNGT, using merge-sort ! This routine is taken from ORDERPACK 2.0 (see http://www.fortran-2000.com/rank/#3.2) ! and therefore left in the original formatting. ! __________________________________________________________ ! This version is not optimized for performance, and is thus ! not as difficult to read as some other ones. ! Michel Olagnon - April 2000 ! __________________________________________________________ ! __________________________________________________________ real(kind = FP), Dimension (:), Intent (In) :: xvalt Integer, Dimension (:), Intent (Out) :: IRNGT ! __________________________________________________________ ! Integer, Dimension (:), Allocatable :: JWRKT Integer :: LMTNA, LMTNC Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! nval = min(SIZE(xvalt), SIZE(irngt)) If (nval <= 0) Then Return End If ! ! Fill-in the index array, creating ordered couples ! Do iind = 2, nval, 2 If (xvalt(iind - 1) <= xvalt(iind)) Then irngt(iind - 1) = iind - 1 irngt(iind) = iind Else irngt(iind - 1) = iind irngt(iind) = iind - 1 End If End Do If (modulo(nval, 2) /= 0) Then irngt(nval) = nval End If ! ! We will now have ordered subsets A - B - A - B - ... ! and merge A and B couples into C - C - ... ! Allocate (jwrkt(1 : nval)) lmtnc = 2 lmtna = 2 ! ! Iteration. Each time, the length of the ordered subsets ! is doubled. ! Do If (lmtna >= nval) Exit iwrkf = 0 lmtnc = 2 * lmtnc iwrk = 0 ! ! Loop on merges of A and B into C ! Do iinda = iwrkf iwrkd = iwrkf + 1 iwrkf = iinda + lmtnc jinda = iinda + lmtna If (iwrkf >= nval) Then If (jinda >= nval) Exit iwrkf = nval End If iindb = jinda ! ! Shortcut for the case when the max of A is smaller ! than the min of B (no need to do anything) ! If (xvalt(irngt(jinda)) <= xvalt(irngt(jinda + 1))) Then iwrk = iwrkf cycle End If ! ! One steps in the C subset, that we create in the final rank array ! Do If (iwrk >= iwrkf) Then ! ! Make a copy of the rank array for next iteration ! irngt(iwrkd : iwrkf) = jwrkt(iwrkd : iwrkf) Exit End If ! iwrk = iwrk + 1 ! ! We still have unprocessed values in both A and B ! If (iinda < jinda) Then If (iindb < iwrkf) Then If (xvalt(irngt(iinda + 1)) > xvalt(irngt(iindb + 1))) & & Then iindb = iindb + 1 jwrkt(iwrk) = irngt(iindb) Else iinda = iinda + 1 jwrkt(iwrk) = irngt(iinda) End If Else ! ! Only A still with unprocessed values ! iinda = iinda + 1 jwrkt(iwrk) = irngt(iinda) End If Else ! ! Only B still with unprocessed values ! irngt(iwrkd : iindb) = jwrkt(iwrkd : iindb) iwrk = iwrkf Exit End If ! End Do End Do ! ! The Cs become As and Bs ! lmtna = 2 * lmtna End Do ! ! Clean up ! Deallocate (jwrkt) Return ! End Subroutine fp_mrgref end module