Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief All kind of helpful little routines
10 : !> \par History
11 : !> none
12 : !> \author CJM & JGH
13 : ! **************************************************************************************************
14 : MODULE util
15 : USE cp_array_sort, ONLY: cp_1d_i4_sort,&
16 : cp_1d_i8_sort,&
17 : cp_1d_r_sort,&
18 : cp_1d_s_sort
19 : USE kinds, ONLY: dp
20 :
21 : IMPLICIT NONE
22 :
23 : PRIVATE
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'util'
25 : PUBLIC :: sort, &
26 : get_limit, &
27 : locate, &
28 : find_boundary, &
29 : sort_unique
30 :
31 : INTERFACE sort
32 : MODULE PROCEDURE cp_1d_s_sort, cp_1d_r_sort, cp_1d_i4_sort, cp_1d_i8_sort, &
33 : sort_cv, sort_im, sort_cm
34 : END INTERFACE
35 :
36 : INTERFACE sort_unique
37 : MODULE PROCEDURE sort_unique1
38 : END INTERFACE
39 :
40 : INTERFACE find_boundary
41 : MODULE PROCEDURE find_boundary1, find_boundary2
42 : END INTERFACE
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Purpose: Given an array array(1:n), and given a value x, a value x_index
48 : !> is returned which is the index value of the array element equal
49 : !> to the value x: x = array(x_index)
50 : !> The array must be monotonic increasing.
51 : !> x_index = 0 is returned, if no array element equal to the value
52 : !> of x was found.
53 : !> \param array ...
54 : !> \param x ...
55 : !> \return ...
56 : !> \par History
57 : !> Derived from the locate function described in
58 : !> Numerical Recipes in Fortran 90 (09.01.2004,MK)
59 : ! **************************************************************************************************
60 10411243 : PURE FUNCTION locate(array, x) RESULT(x_index)
61 : INTEGER, DIMENSION(:), INTENT(IN) :: array
62 : INTEGER, INTENT(IN) :: x
63 : INTEGER :: x_index
64 :
65 : INTEGER :: jl, jm, ju, n
66 :
67 10411243 : x_index = 0
68 :
69 10411243 : IF (x < array(1)) RETURN
70 10411243 : n = SIZE(array)
71 10411243 : IF (x > array(n)) RETURN
72 10411169 : jl = 0
73 10411169 : ju = n + 1
74 36419265 : DO WHILE (ju - jl > 1)
75 26008096 : jm = (ju + jl)/2
76 36419265 : IF (x >= array(jm)) THEN
77 : jl = jm
78 : ELSE
79 19014277 : ju = jm
80 : END IF
81 : END DO
82 10411169 : IF (x == array(jl)) x_index = jl
83 : END FUNCTION locate
84 :
85 : ! **************************************************************************************************
86 : !> \brief Sorts and returns a logical that checks if all elements are unique
87 : !> \param arr ...
88 : !> \param unique ...
89 : !> \par History
90 : !> Teodoro Laino - Zurich University [tlaino] 04.2007
91 : ! **************************************************************************************************
92 2309 : SUBROUTINE sort_unique1(arr, unique)
93 : INTEGER, DIMENSION(:), INTENT(INOUT) :: arr
94 : LOGICAL, INTENT(OUT) :: unique
95 :
96 : INTEGER :: i, n
97 2309 : INTEGER, ALLOCATABLE, DIMENSION(:) :: wrk
98 :
99 2309 : n = SIZE(arr)
100 2309 : unique = .TRUE.
101 6277 : ALLOCATE (wrk(n))
102 2309 : CALL sort(arr, n, wrk)
103 3534 : DO i = 2, n
104 3534 : IF (arr(i) == arr(i - 1)) THEN
105 264 : unique = .FALSE.
106 264 : EXIT
107 : END IF
108 : END DO
109 2309 : DEALLOCATE (wrk)
110 2309 : END SUBROUTINE sort_unique1
111 :
112 : ! **************************************************************************************************
113 : !> \brief Sorts an array of strings
114 : !> \param arr ...
115 : !> \param n ...
116 : !> \param index ...
117 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
118 : ! **************************************************************************************************
119 23190 : SUBROUTINE sort_cv(arr, n, index)
120 : INTEGER, INTENT(IN) :: n
121 : CHARACTER(LEN=*), INTENT(INOUT) :: arr(1:n)
122 : INTEGER, INTENT(OUT) :: INDEX(1:n)
123 :
124 : INTEGER :: i, j, max_length
125 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: entries
126 :
127 23190 : max_length = 0
128 726928 : DO i = 1, n
129 726928 : max_length = MAX(max_length, LEN_TRIM(arr(i)))
130 : END DO
131 92760 : ALLOCATE (entries(max_length, SIZE(arr)))
132 726928 : DO i = 1, n
133 3203396 : DO j = 1, LEN_TRIM(arr(i))
134 3203396 : entries(j, i) = ICHAR(arr(i) (j:j))
135 : END DO
136 726928 : IF (j <= max_length) THEN
137 627454 : entries(j:max_length, i) = ICHAR(" ")
138 : END IF
139 : END DO
140 23190 : CALL sort_im(entries, istart=1, iend=n, j=1, jsize=max_length, INDEX=INDEX)
141 : ! Recover string once ordered
142 726928 : DO i = 1, n
143 3635566 : DO j = 1, max_length
144 3612376 : arr(i) (j:j) = CHAR(entries(j, i))
145 : END DO
146 : END DO
147 23190 : DEALLOCATE (entries)
148 23190 : END SUBROUTINE sort_cv
149 :
150 : #:for argtype, itemtype in [['INTEGER', 'INTEGER'], ['CHARACTER(LEN=*)', 'CHARACTER(LEN=LEN(matrix))']]
151 : ! **************************************************************************************************
152 : !> \brief Sorts a multiple arrays M(j,i), ordering iteratively over i with fixed j
153 : !> \param matrix ...
154 : !> \param istart ...
155 : !> \param iend ...
156 : !> \param j ...
157 : !> \param jsize ...
158 : !> \param index ...
159 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
160 : ! **************************************************************************************************
161 270110 : RECURSIVE SUBROUTINE sort_${argtype[0]}$m(matrix, istart, iend, j, jsize, index)
162 : ${argtype}$, DIMENSION(:, :), INTENT(INOUT) :: matrix
163 : INTEGER, INTENT(IN) :: istart, iend, j, jsize
164 : INTEGER, DIMENSION(:), INTENT(INOUT) :: index
165 :
166 :
167 23154 : ${itemtype}$ :: item
168 270110 : ${itemtype}$, ALLOCATABLE, DIMENSION(:) :: work, work2
169 : INTEGER :: i, ind, isize, k, kend, kstart
170 270110 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bck_index, tmp_index
171 :
172 270110 : isize = iend - istart + 1
173 : ! Initialize the INDEX array only for the first row..
174 270110 : IF (j == 1) THEN
175 981376 : DO i = 1, isize
176 981376 : INDEX(i) = i
177 : ENDDO
178 : END IF
179 :
180 : ! Allocate scratch arrays
181 1666968 : ALLOCATE (work(isize), work2(isize), tmp_index(isize), bck_index(isize))
182 4119790 : ind = 0
183 4119790 : DO i = istart, iend
184 3849680 : ind = ind + 1
185 3849680 : work(ind) = matrix(j, i)
186 4119790 : bck_index(ind) = INDEX(i)
187 : END DO
188 :
189 : ! Ordering row (j) interval istart..iend
190 270110 : CALL sort(work, isize, tmp_index)
191 :
192 : ! Copy into global INDEX array with a proper mapping
193 270110 : ind = 0
194 4119790 : DO i = istart, iend
195 3849680 : ind = ind + 1
196 3849680 : INDEX(i) = bck_index(tmp_index(ind))
197 4119790 : matrix(j, i) = work(ind)
198 : END DO
199 :
200 : ! Reorder the rest of the array according the present reordering
201 975296 : DO k = j + 1, jsize
202 : ind = 0
203 9143732 : DO i = istart, iend
204 8438546 : ind = ind + 1
205 9143732 : work2(ind) = matrix(k, i)
206 : END DO
207 : ind = 0
208 9413842 : DO i = istart, iend
209 8438546 : ind = ind + 1
210 9143732 : matrix(k, i) = work2(tmp_index(ind))
211 : END DO
212 : END DO
213 :
214 : ! There are more rows to order..
215 270110 : IF (j < jsize) THEN
216 210142 : kstart = istart
217 210142 : item = work(1)
218 195582 : ind = 0
219 3141600 : DO i = istart, iend
220 2931458 : ind = ind + 1
221 3141600 : IF (item /= work(ind)) THEN
222 76526 : kend = i - 1
223 76526 : IF (kstart /= kend) THEN
224 54220 : CALL sort(matrix, kstart, kend, j + 1, jsize, INDEX)
225 : END IF
226 76526 : item = work(ind)
227 76526 : kstart = i
228 : END IF
229 : END DO
230 210142 : kend = i - 1
231 210142 : IF (kstart /= kend) THEN
232 192592 : CALL sort(matrix, kstart, kend, j + 1, jsize, INDEX)
233 : END IF
234 : END IF
235 270110 : DEALLOCATE (work, work2, tmp_index, bck_index)
236 :
237 270110 : END SUBROUTINE sort_${argtype[0]}$m
238 : #:endfor
239 :
240 : ! **************************************************************************************************
241 : !> \brief divide m entries into n parts, return size of part me
242 : !> \param m ...
243 : !> \param n ...
244 : !> \param me ...
245 : !> \return ...
246 : ! **************************************************************************************************
247 10346348 : PURE FUNCTION get_limit(m, n, me) RESULT(nlim)
248 : INTEGER, INTENT(IN) :: m, n, me
249 : INTEGER :: nlim(2)
250 :
251 : INTEGER :: nl, nu
252 : REAL(KIND=dp) :: part
253 :
254 10346348 : part = REAL(m, KIND=dp)/REAL(n, KIND=dp)
255 10346348 : nl = NINT(REAL(me, KIND=dp)*part) + 1
256 10346348 : nu = NINT(REAL(me + 1, KIND=dp)*part)
257 10346348 : nlim(1) = MAX(1, nl)
258 10346348 : nlim(2) = MIN(m, nu)
259 :
260 10346348 : END FUNCTION get_limit
261 :
262 : ! **************************************************************************************************
263 : !> \brief finds boundary where element search starts and ends in a 1D array
264 : !> array1: XXXXXAAAAAAAAAXXDGFSFGWDDDDDDDAAAWE
265 : !> | |
266 : !> start end (searching for A)
267 : !> \param num_array ...
268 : !> \param ntot ...
269 : !> \param first ...
270 : !> \param last ...
271 : !> \param search ...
272 : ! **************************************************************************************************
273 2480 : PURE SUBROUTINE find_boundary1(num_array, ntot, first, last, search)
274 : INTEGER, POINTER :: num_array(:)
275 : INTEGER, INTENT(IN) :: ntot
276 : INTEGER, INTENT(OUT) :: first, last
277 : INTEGER, INTENT(IN) :: search
278 :
279 : INTEGER :: i
280 : LOGICAL :: found
281 :
282 2480 : found = .FALSE.
283 2480 : first = 0
284 2480 : last = ntot
285 :
286 1911983 : DO i = 1, ntot
287 1911983 : IF (num_array(i) == search) THEN
288 580855 : IF (.NOT. found) THEN
289 2480 : first = i
290 : END IF
291 : found = .TRUE.
292 : ELSE
293 1330432 : IF (found) THEN
294 1784 : last = i - 1
295 1784 : EXIT
296 : END IF
297 : found = .FALSE.
298 : END IF
299 : END DO
300 :
301 2480 : END SUBROUTINE find_boundary1
302 :
303 : ! **************************************************************************************************
304 : !> \brief finds boundary where element search1 starts and ends in array1 checking
305 : !> at the same time search2 in array2
306 : !> array1: XXXXXAAAAAAAAAXXDGFSFGWDDDDDDDAAAWE
307 : !> array2: XXXXASDEYYYYASDEFAAAARGASGASRGAWRRR
308 : !> | |
309 : !> start end (searching for A and Y)
310 : !> \param num_array1 ...
311 : !> \param num_array2 ...
312 : !> \param ntot ...
313 : !> \param first ...
314 : !> \param last ...
315 : !> \param search1 ...
316 : !> \param search2 ...
317 : ! **************************************************************************************************
318 1228 : PURE SUBROUTINE find_boundary2(num_array1, num_array2, ntot, first, last, search1, search2)
319 : INTEGER, POINTER :: num_array1(:), num_array2(:)
320 : INTEGER, INTENT(IN) :: ntot
321 : INTEGER, INTENT(OUT) :: first, last
322 : INTEGER, INTENT(IN) :: search1, search2
323 :
324 : INTEGER :: i, tfirst, tlast
325 : LOGICAL :: found
326 :
327 1228 : found = .FALSE.
328 1228 : first = 0
329 1228 : last = ntot
330 :
331 1228 : CALL find_boundary(num_array1, ntot, tfirst, tlast, search1)
332 1228 : last = tlast
333 41028 : DO i = tfirst, tlast
334 41028 : IF (num_array2(i) == search2) THEN
335 15986 : IF (.NOT. found) THEN
336 1228 : first = i
337 : END IF
338 : found = .TRUE.
339 : ELSE
340 24232 : IF (found) THEN
341 418 : last = i - 1
342 418 : EXIT
343 : END IF
344 : found = .FALSE.
345 : END IF
346 : END DO
347 :
348 1228 : END SUBROUTINE find_boundary2
349 :
350 : END MODULE util
|
