Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for optimizing load balance between processes in HFX calculations
10 : !> \par History
11 : !> 04.2008 created [Manuel Guidon]
12 : !> \author Manuel Guidon
13 : ! **************************************************************************************************
14 : MODULE hfx_load_balance_methods
15 : USE cell_types, ONLY: cell_type
16 : USE cp_files, ONLY: close_file, &
17 : open_file
18 : USE message_passing, ONLY: mp_para_env_type
19 : USE hfx_pair_list_methods, ONLY: build_atomic_pair_list, &
20 : build_pair_list
21 : USE hfx_types, ONLY: &
22 : hfx_basis_type, hfx_block_range_type, hfx_distribution, hfx_load_balance_type, hfx_p_kind, &
23 : hfx_screen_coeff_type, hfx_set_distr_energy, hfx_set_distr_forces, hfx_type, &
24 : pair_list_type, pair_set_list_type
25 : USE input_constants, ONLY: hfx_do_eval_energy, &
26 : hfx_do_eval_forces
27 : USE kinds, ONLY: dp, &
28 : int_8
29 : USE message_passing, ONLY: mp_waitall, mp_request_type
30 : USE parallel_rng_types, ONLY: UNIFORM, &
31 : rng_stream_type
32 : USE particle_types, ONLY: particle_type
33 : USE util, ONLY: sort
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 : PRIVATE
38 :
39 : PUBLIC :: hfx_load_balance, &
40 : hfx_update_load_balance, &
41 : collect_load_balance_info, cost_model, p1_energy, p2_energy, p3_energy
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_load_balance_methods'
44 :
45 : REAL(KIND=dp), PARAMETER :: p1_energy(12) = [2.9461408209700424_dp, 1.0624718662999657_dp, &
46 : -1.91570128356921242E-002_dp, -1.6668495454436603_dp, &
47 : 1.7512639006523709_dp, -9.76074323945336081E-002_dp, &
48 : 2.6230786127311889_dp, -0.31870737623014189_dp, &
49 : 7.9588203912690973_dp, 1.8331423413134813_dp, &
50 : -0.15427618665346299_dp, 0.19749436090711650_dp]
51 : REAL(KIND=dp), PARAMETER :: p2_energy(12) = [2.3104682960662593_dp, 1.8744052737304417_dp, &
52 : -9.36564055598656797E-002_dp, 0.64284973765086939_dp, &
53 : 1.0137565430060556_dp, -6.80088178288954567E-003_dp, &
54 : 1.1692629207374552_dp, -2.6314710080507573_dp, &
55 : 19.237814781880786_dp, 1.0505934173661349_dp, &
56 : 0.80382371955699250_dp, 0.49903401991818103_dp]
57 : REAL(KIND=dp), PARAMETER :: p3_energy(2) = [7.82336287670072350E-002_dp, 0.38073304105744837_dp]
58 : REAL(KIND=dp), PARAMETER :: p1_forces(12) = [2.5746279948798874_dp, 1.3420575378609276_dp, &
59 : -9.41673106447732111E-002_dp, 0.94568006899317825_dp, &
60 : -1.4511897117448544_dp, 0.59178934677316952_dp, &
61 : 2.7291149361757236_dp, -0.50555512044800210_dp, &
62 : 8.3508180969609871_dp, 1.6829982496141809_dp, &
63 : -0.74895370472152600_dp, 0.43801726744197500_dp]
64 : REAL(KIND=dp), PARAMETER :: p2_forces(12) = [2.6398568961569020_dp, 2.3024918834564101_dp, &
65 : 5.33216585432061581E-003_dp, 0.45572145697283628_dp, &
66 : 1.8119743851500618_dp, -0.12533918548421166_dp, &
67 : -1.4040312084552751_dp, -4.5331650463917859_dp, &
68 : 12.593431549069477_dp, 1.1311978374487595_dp, &
69 : 1.4245996087624646_dp, 1.1425350529853495_dp]
70 : REAL(KIND=dp), PARAMETER :: p3_forces(2) = [0.12051930516830946_dp, 1.3828051586144336_dp]
71 :
72 : !***
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief Distributes the computation of eri's to all available processes.
78 : !> \param x_data Object that stores the indices array
79 : !> \param eps_schwarz screening parameter
80 : !> \param particle_set , atomic_kind_set, para_env ...
81 : !> \param max_set Maximum number of set to be considered
82 : !> \param para_env para_env
83 : !> \param coeffs_set screening functions
84 : !> \param coeffs_kind screening functions
85 : !> \param is_assoc_atomic_block_global KS-matrix sparsity
86 : !> \param do_periodic flag for periodicity
87 : !> \param load_balance_parameter Parameters for Monte-Carlo routines
88 : !> \param kind_of helper array for mapping
89 : !> \param basis_parameter Basis set parameters
90 : !> \param pmax_set Initial screening matrix
91 : !> \param pmax_atom ...
92 : !> \param i_thread Process ID of current Thread
93 : !> \param n_threads Total Number of Threads
94 : !> \param cell cell
95 : !> \param do_p_screening Flag for initial p screening
96 : !> \param map_atom_to_kind_atom ...
97 : !> \param nkind ...
98 : !> \param eval_type ...
99 : !> \param pmax_block ...
100 : !> \param use_virial ...
101 : !> \par History
102 : !> 06.2007 created [Manuel Guidon]
103 : !> 08.2007 new parallel scheme [Manuel Guidon]
104 : !> 09.2007 new 'modulo' parellel scheme and Monte Carlo step [Manuel Guidon]
105 : !> 11.2007 parallelize load balance on box_idx1 [Manuel Guidon]
106 : !> 02.2009 completely refactored [Manuel Guidon]
107 : !> \author Manuel Guidon
108 : !> \note
109 : !> The optimization is done via a binning procedure followed by simple
110 : !> Monte Carlo procedure:
111 : !> In a first step the total amount of integrals in the system is calculated,
112 : !> taking into account the sparsity of the KS-matrix , the screening based
113 : !> on near/farfield approximations and if desired the screening on an initial
114 : !> density matrix.
115 : !> In a second step, bins are generate that contain approximately the same number
116 : !> of integrals, and a cost for these bins is estimated (currently the number of integrals)
117 : !> In a third step, a Monte Carlo procedure optimizes the assignment
118 : !> of the different loads to each process
119 : !> At the end each process owns an unique array of *atomic* indices-ranges
120 : !> that are used to decide whether a process has to calculate a certain
121 : !> bunch of integrals or not
122 : ! **************************************************************************************************
123 2128 : SUBROUTINE hfx_load_balance(x_data, eps_schwarz, particle_set, max_set, para_env, &
124 : coeffs_set, coeffs_kind, &
125 2128 : is_assoc_atomic_block_global, do_periodic, &
126 : load_balance_parameter, kind_of, basis_parameter, pmax_set, &
127 : pmax_atom, i_thread, n_threads, cell, &
128 : do_p_screening, map_atom_to_kind_atom, nkind, eval_type, &
129 : pmax_block, use_virial)
130 : TYPE(hfx_type), POINTER :: x_data
131 : REAL(dp), INTENT(IN) :: eps_schwarz
132 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
133 : INTEGER, INTENT(IN) :: max_set
134 : TYPE(mp_para_env_type), POINTER :: para_env
135 : TYPE(hfx_screen_coeff_type), &
136 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
137 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
138 : POINTER :: coeffs_kind
139 : INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block_global
140 : LOGICAL :: do_periodic
141 : TYPE(hfx_load_balance_type), POINTER :: load_balance_parameter
142 : INTEGER :: kind_of(*)
143 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
144 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set
145 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom
146 : INTEGER, INTENT(IN) :: i_thread, n_threads
147 : TYPE(cell_type), POINTER :: cell
148 : LOGICAL, INTENT(IN) :: do_p_screening
149 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom
150 : INTEGER, INTENT(IN) :: nkind, eval_type
151 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_block
152 : LOGICAL, INTENT(IN) :: use_virial
153 :
154 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_load_balance'
155 : ! Cap per-thread cost cache at ~16 MB (16 MB * n_threads per rank);
156 : ! larger systems recompute estimate_block_cost on the second pass.
157 : INTEGER(int_8), PARAMETER :: max_cached_quartet_costs = 2000000_int_8
158 :
159 : CHARACTER(LEN=512) :: error_msg
160 : INTEGER :: block_size, current_block_id, data_from, dest, handle, handle_inner, &
161 : handle_range, i, iatom_block, iatom_end, iatom_start, ibin, icpu, j, jatom_block, &
162 : jatom_end, jatom_start, katom_block, katom_end, katom_start, latom_block, latom_end, &
163 : latom_start, mepos, my_process_id, n_processes, natom, nbins, nblocks, ncpu, &
164 : new_iatom_end, new_iatom_start, new_jatom_end, new_jatom_start, non_empty_blocks, &
165 : objective_block_size, objective_nblocks, source, total_blocks
166 6384 : TYPE(mp_request_type), DIMENSION(2) :: req
167 : INTEGER(int_8) :: atom_block, block_cost_val, cost_per_bin, cost_per_core, current_cost, &
168 : distribution_counter_end, distribution_counter_start, global_quartet_counter, &
169 : local_quartet_counter, n_my_quartets, q_idx, &
170 : self_cost_per_block, tmp_block, total_block_self_cost
171 2128 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: buffer_in, buffer_out, cached_block_cost
172 2128 : INTEGER(int_8), DIMENSION(:), POINTER :: local_cost_matrix, recbuffer, &
173 2128 : sendbuffer, swapbuffer
174 : INTEGER(int_8), DIMENSION(:), POINTER, SAVE :: cost_matrix
175 : INTEGER(int_8), SAVE :: shm_global_quartet_counter, &
176 : shm_local_quartet_counter
177 4256 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, tmp_index, tmp_pos, &
178 2128 : to_be_sorted
179 : INTEGER, DIMENSION(:), POINTER, SAVE :: shm_distribution_vector
180 : INTEGER, SAVE :: shm_nblocks
181 : LOGICAL :: changed, last_bin_needs_to_be_filled, &
182 : optimized
183 : LOGICAL, DIMENSION(:, :), POINTER, SAVE :: atomic_pair_list
184 : REAL(dp) :: coeffs_kind_max0, log10_eps_schwarz, &
185 : log_2, pmax_blocks
186 4256 : TYPE(hfx_block_range_type), DIMENSION(:), POINTER :: blocks_guess, tmp_blocks, tmp_blocks2
187 : TYPE(hfx_block_range_type), DIMENSION(:), &
188 : POINTER, SAVE :: shm_blocks
189 6384 : TYPE(hfx_distribution), DIMENSION(:), POINTER :: binned_dist, ptr_to_tmp_dist, tmp_dist
190 : TYPE(hfx_distribution), DIMENSION(:, :), POINTER, &
191 : SAVE :: full_dist
192 : TYPE(pair_list_type) :: list_ij, list_kl
193 : TYPE(pair_set_list_type), ALLOCATABLE, &
194 2128 : DIMENSION(:) :: set_list_ij, set_list_kl
195 :
196 2128 : !$OMP BARRIER
197 4256 : !$OMP MASTER
198 2128 : CALL timeset(routineN, handle)
199 : !$OMP END MASTER
200 2128 : !$OMP BARRIER
201 :
202 2128 : log10_eps_schwarz = LOG10(eps_schwarz)
203 2128 : log_2 = LOG10(2.0_dp)
204 14372 : coeffs_kind_max0 = MAXVAL(coeffs_kind(:, :)%x(2))
205 2128 : ncpu = para_env%num_pe
206 2128 : n_processes = ncpu*n_threads
207 2128 : natom = SIZE(particle_set)
208 :
209 2128 : block_size = load_balance_parameter%block_size
210 41068 : ALLOCATE (set_list_ij((max_set*block_size)**2))
211 38940 : ALLOCATE (set_list_kl((max_set*block_size)**2))
212 :
213 2128 : IF (.NOT. load_balance_parameter%blocks_initialized) THEN
214 1490 : !$OMP BARRIER
215 1490 : !$OMP MASTER
216 1490 : CALL timeset(routineN//"_range", handle_range)
217 :
218 1490 : nblocks = MAX((natom + block_size - 1)/block_size, 1)
219 9186 : ALLOCATE (blocks_guess(nblocks))
220 9186 : ALLOCATE (tmp_blocks(natom))
221 7696 : ALLOCATE (tmp_blocks2(natom))
222 :
223 1490 : pmax_blocks = 0.0_dp
224 2980 : SELECT CASE (eval_type)
225 : CASE (hfx_do_eval_energy)
226 1490 : atomic_pair_list => x_data%atomic_pair_list
227 : CASE (hfx_do_eval_forces)
228 1490 : atomic_pair_list => x_data%atomic_pair_list_forces
229 : END SELECT
230 25846 : atomic_pair_list = .TRUE.
231 : CALL init_blocks(nkind, para_env, natom, block_size, nblocks, blocks_guess, &
232 : list_ij, list_kl, set_list_ij, set_list_kl, &
233 : particle_set, &
234 : coeffs_set, coeffs_kind, &
235 : is_assoc_atomic_block_global, do_periodic, &
236 : kind_of, basis_parameter, pmax_set, pmax_atom, &
237 : pmax_blocks, cell, &
238 : do_p_screening, map_atom_to_kind_atom, eval_type, &
239 1490 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
240 :
241 1490 : total_block_self_cost = 0
242 :
243 6206 : DO i = 1, nblocks
244 6206 : total_block_self_cost = total_block_self_cost + blocks_guess(i)%cost
245 : END DO
246 :
247 1490 : CALL para_env%sum(total_block_self_cost)
248 :
249 1490 : objective_block_size = load_balance_parameter%block_size
250 1490 : objective_nblocks = MAX((natom + objective_block_size - 1)/objective_block_size, 1)
251 :
252 1490 : self_cost_per_block = (total_block_self_cost + objective_nblocks - 1)/(objective_nblocks)
253 :
254 6206 : DO i = 1, nblocks
255 6206 : tmp_blocks2(i) = blocks_guess(i)
256 : END DO
257 :
258 : optimized = .FALSE.
259 : i = 0
260 2980 : DO WHILE (.NOT. optimized)
261 1490 : i = i + 1
262 1490 : current_block_id = 0
263 1490 : changed = .FALSE.
264 6206 : DO atom_block = 1, nblocks
265 4716 : current_block_id = current_block_id + 1
266 4716 : iatom_start = tmp_blocks2(atom_block)%istart
267 4716 : iatom_end = tmp_blocks2(atom_block)%iend
268 6206 : IF (tmp_blocks2(atom_block)%cost > 1.5_dp*self_cost_per_block .AND. iatom_end - iatom_start > 0) THEN
269 0 : changed = .TRUE.
270 0 : new_iatom_start = iatom_start
271 0 : new_iatom_end = (iatom_end - iatom_start + 1)/2 + iatom_start - 1
272 0 : new_jatom_start = new_iatom_end + 1
273 0 : new_jatom_end = iatom_end
274 0 : tmp_blocks(current_block_id)%istart = new_iatom_start
275 0 : tmp_blocks(current_block_id)%iend = new_iatom_end
276 : tmp_blocks(current_block_id)%cost = estimate_block_cost( &
277 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
278 : new_iatom_start, new_iatom_end, new_iatom_start, new_iatom_end, &
279 : new_iatom_start, new_iatom_end, new_iatom_start, new_iatom_end, &
280 : particle_set, &
281 : coeffs_set, coeffs_kind, &
282 : is_assoc_atomic_block_global, do_periodic, &
283 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
284 : cell, &
285 : do_p_screening, map_atom_to_kind_atom, eval_type, &
286 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
287 0 : current_block_id = current_block_id + 1
288 0 : tmp_blocks(current_block_id)%istart = new_jatom_start
289 0 : tmp_blocks(current_block_id)%iend = new_jatom_end
290 : tmp_blocks(current_block_id)%cost = estimate_block_cost( &
291 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
292 : new_jatom_start, new_jatom_end, new_jatom_start, new_jatom_end, &
293 : new_jatom_start, new_jatom_end, new_jatom_start, new_jatom_end, &
294 : particle_set, &
295 : coeffs_set, coeffs_kind, &
296 : is_assoc_atomic_block_global, do_periodic, &
297 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
298 : cell, &
299 : do_p_screening, map_atom_to_kind_atom, eval_type, &
300 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
301 : ELSE
302 4716 : tmp_blocks(current_block_id)%istart = iatom_start
303 4716 : tmp_blocks(current_block_id)%iend = iatom_end
304 4716 : tmp_blocks(current_block_id)%cost = tmp_blocks2(atom_block)%cost
305 : END IF
306 : END DO
307 1490 : IF (.NOT. changed) optimized = .TRUE.
308 1490 : IF (i > 20) optimized = .TRUE.
309 1490 : nblocks = current_block_id
310 7696 : DO atom_block = 1, nblocks
311 6206 : tmp_blocks2(atom_block) = tmp_blocks(atom_block)
312 : END DO
313 : END DO
314 :
315 1490 : DEALLOCATE (tmp_blocks2)
316 :
317 : ! ** count number of non empty blocks on each node
318 1490 : non_empty_blocks = 0
319 6206 : DO atom_block = 1, nblocks
320 4716 : IF (tmp_blocks(atom_block)%istart == 0) CYCLE
321 6206 : non_empty_blocks = non_empty_blocks + 1
322 : END DO
323 :
324 4470 : ALLOCATE (rcount(ncpu))
325 1490 : rcount = 0
326 1490 : rcount(para_env%mepos + 1) = non_empty_blocks
327 1490 : CALL para_env%sum(rcount)
328 :
329 : ! ** sum all non_empty_blocks
330 1490 : total_blocks = 0
331 4468 : DO i = 1, ncpu
332 4468 : total_blocks = total_blocks + rcount(i)
333 : END DO
334 :
335 : ! ** calculate offsets
336 2980 : ALLOCATE (rdispl(ncpu))
337 4468 : rcount(:) = rcount(:)*3
338 1490 : rdispl(1) = 0
339 2978 : DO i = 2, ncpu
340 2978 : rdispl(i) = rdispl(i - 1) + rcount(i - 1)
341 : END DO
342 :
343 4440 : ALLOCATE (buffer_in(3*non_empty_blocks))
344 :
345 1490 : non_empty_blocks = 0
346 6206 : DO atom_block = 1, nblocks
347 4716 : IF (tmp_blocks(atom_block)%istart == 0) CYCLE
348 2361 : buffer_in(non_empty_blocks*3 + 1) = tmp_blocks(atom_block)%istart
349 2361 : buffer_in(non_empty_blocks*3 + 2) = tmp_blocks(atom_block)%iend
350 2361 : buffer_in(non_empty_blocks*3 + 3) = tmp_blocks(atom_block)%cost
351 6206 : non_empty_blocks = non_empty_blocks + 1
352 : END DO
353 :
354 1490 : nblocks = total_blocks
355 :
356 9186 : ALLOCATE (tmp_blocks2(nblocks))
357 :
358 4470 : ALLOCATE (buffer_out(3*nblocks))
359 :
360 : ! ** Gather all three arrays
361 1490 : CALL para_env%allgatherv(buffer_in, buffer_out, rcount, rdispl)
362 :
363 6206 : DO i = 1, nblocks
364 4716 : tmp_blocks2(i)%istart = INT(buffer_out((i - 1)*3 + 1))
365 4716 : tmp_blocks2(i)%iend = INT(buffer_out((i - 1)*3 + 2))
366 6206 : tmp_blocks2(i)%cost = buffer_out((i - 1)*3 + 3)
367 : END DO
368 :
369 : ! ** Now we sort the blocks
370 4470 : ALLOCATE (to_be_sorted(nblocks))
371 2980 : ALLOCATE (tmp_index(nblocks))
372 :
373 6206 : DO atom_block = 1, nblocks
374 6206 : to_be_sorted(atom_block) = tmp_blocks2(atom_block)%istart
375 : END DO
376 :
377 1490 : CALL sort(to_be_sorted, nblocks, tmp_index)
378 :
379 7696 : ALLOCATE (x_data%blocks(nblocks))
380 :
381 6206 : DO atom_block = 1, nblocks
382 6206 : x_data%blocks(atom_block) = tmp_blocks2(tmp_index(atom_block))
383 : END DO
384 :
385 1490 : shm_blocks => x_data%blocks
386 1490 : shm_nblocks = nblocks
387 :
388 : ! ** Set nblocks in structure
389 1490 : load_balance_parameter%nblocks = nblocks
390 :
391 1490 : DEALLOCATE (blocks_guess, tmp_blocks, tmp_blocks2)
392 :
393 1490 : DEALLOCATE (rcount, rdispl, buffer_in, buffer_out, to_be_sorted, tmp_index)
394 :
395 1490 : load_balance_parameter%blocks_initialized = .TRUE.
396 :
397 10922 : x_data%blocks = shm_blocks
398 1490 : load_balance_parameter%nblocks = shm_nblocks
399 1490 : load_balance_parameter%blocks_initialized = .TRUE.
400 :
401 5960 : ALLOCATE (x_data%pmax_block(shm_nblocks, shm_nblocks))
402 25846 : x_data%pmax_block = 0.0_dp
403 1490 : pmax_block => x_data%pmax_block
404 2980 : CALL timestop(handle_range)
405 : !$OMP END MASTER
406 1490 : !$OMP BARRIER
407 :
408 1490 : IF (.NOT. load_balance_parameter%blocks_initialized) THEN
409 0 : ALLOCATE (x_data%blocks(shm_nblocks))
410 0 : x_data%blocks = shm_blocks
411 0 : load_balance_parameter%nblocks = shm_nblocks
412 0 : load_balance_parameter%blocks_initialized = .TRUE.
413 : END IF
414 : !! ** precalculate maximum density matrix elements in blocks
415 1490 : !$OMP BARRIER
416 : END IF
417 :
418 2128 : !$OMP BARRIER
419 2128 : !$OMP MASTER
420 2128 : pmax_block => x_data%pmax_block
421 34544 : pmax_block = 0.0_dp
422 2128 : IF (do_p_screening) THEN
423 1736 : DO iatom_block = 1, shm_nblocks
424 1324 : iatom_start = x_data%blocks(iatom_block)%istart
425 1324 : iatom_end = x_data%blocks(iatom_block)%iend
426 6512 : DO jatom_block = 1, shm_nblocks
427 4776 : jatom_start = x_data%blocks(jatom_block)%istart
428 4776 : jatom_end = x_data%blocks(jatom_block)%iend
429 15652 : pmax_block(iatom_block, jatom_block) = MAXVAL(pmax_atom(iatom_start:iatom_end, jatom_start:jatom_end))
430 : END DO
431 : END DO
432 : END IF
433 :
434 3618 : SELECT CASE (eval_type)
435 : CASE (hfx_do_eval_energy)
436 1490 : atomic_pair_list => x_data%atomic_pair_list
437 : CASE (hfx_do_eval_forces)
438 2128 : atomic_pair_list => x_data%atomic_pair_list_forces
439 : END SELECT
440 : CALL build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
441 : do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
442 2128 : x_data%blocks)
443 :
444 : !$OMP END MASTER
445 2128 : !$OMP BARRIER
446 :
447 : !! If there is only 1 cpu skip the binning
448 2128 : IF (n_processes == 1) THEN
449 4 : ALLOCATE (tmp_dist(1))
450 2 : tmp_dist(1)%number_of_atom_quartets = HUGE(tmp_dist(1)%number_of_atom_quartets)
451 2 : tmp_dist(1)%istart = 0_int_8
452 2 : ptr_to_tmp_dist => tmp_dist(:)
453 4 : SELECT CASE (eval_type)
454 : CASE (hfx_do_eval_energy)
455 2 : CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
456 : CASE (hfx_do_eval_forces)
457 2 : CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
458 : END SELECT
459 2 : DEALLOCATE (tmp_dist)
460 : ELSE
461 : !! Calculate total numbers of integrals that have to be calculated (wrt screening and symmetry)
462 2126 : !$OMP BARRIER
463 2126 : !$OMP MASTER
464 2126 : CALL timeset(routineN//"_count", handle_inner)
465 : !$OMP END MASTER
466 2126 : !$OMP BARRIER
467 :
468 2126 : cost_per_core = 0_int_8
469 2126 : my_process_id = para_env%mepos*n_threads + i_thread
470 2126 : nblocks = load_balance_parameter%nblocks
471 2126 : n_my_quartets = 0_int_8
472 595542 : DO atom_block = my_process_id, INT(nblocks, KIND=int_8)**4 - 1_int_8, n_processes
473 593416 : latom_block = INT(MODULO(atom_block, INT(nblocks, KIND=int_8))) + 1
474 593416 : tmp_block = atom_block/nblocks
475 593416 : katom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
476 593416 : IF (latom_block < katom_block) CYCLE
477 332701 : tmp_block = tmp_block/nblocks
478 332701 : jatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
479 332701 : tmp_block = tmp_block/nblocks
480 332701 : iatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
481 332701 : IF (jatom_block < iatom_block) CYCLE
482 187570 : n_my_quartets = n_my_quartets + 1_int_8
483 189696 : IF (n_my_quartets == max_cached_quartet_costs) EXIT
484 : END DO
485 2096 : IF (n_my_quartets > 0_int_8) THEN
486 6288 : ALLOCATE (cached_block_cost(n_my_quartets))
487 : END IF
488 2126 : q_idx = 0_int_8
489 :
490 595542 : DO atom_block = my_process_id, INT(nblocks, KIND=int_8)**4 - 1, n_processes
491 :
492 593416 : latom_block = INT(MODULO(atom_block, INT(nblocks, KIND=int_8))) + 1
493 593416 : tmp_block = atom_block/nblocks
494 593416 : katom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
495 593416 : IF (latom_block < katom_block) CYCLE
496 332701 : tmp_block = tmp_block/nblocks
497 332701 : jatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
498 332701 : tmp_block = tmp_block/nblocks
499 332701 : iatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
500 332701 : IF (jatom_block < iatom_block) CYCLE
501 :
502 187570 : q_idx = q_idx + 1_int_8
503 :
504 187570 : iatom_start = x_data%blocks(iatom_block)%istart
505 187570 : iatom_end = x_data%blocks(iatom_block)%iend
506 187570 : jatom_start = x_data%blocks(jatom_block)%istart
507 187570 : jatom_end = x_data%blocks(jatom_block)%iend
508 187570 : katom_start = x_data%blocks(katom_block)%istart
509 187570 : katom_end = x_data%blocks(katom_block)%iend
510 187570 : latom_start = x_data%blocks(latom_block)%istart
511 187570 : latom_end = x_data%blocks(latom_block)%iend
512 :
513 357153 : SELECT CASE (eval_type)
514 : CASE (hfx_do_eval_energy)
515 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block), &
516 : pmax_block(latom_block, jatom_block), &
517 : pmax_block(latom_block, iatom_block), &
518 169583 : pmax_block(katom_block, jatom_block))
519 : CASE (hfx_do_eval_forces)
520 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block) + &
521 : pmax_block(latom_block, jatom_block), &
522 : pmax_block(latom_block, iatom_block) + &
523 187570 : pmax_block(katom_block, jatom_block))
524 : END SELECT
525 :
526 187570 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) THEN
527 652 : IF (q_idx <= n_my_quartets) cached_block_cost(q_idx) = -1_int_8
528 : CYCLE
529 : END IF
530 :
531 : block_cost_val = estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
532 : iatom_start, iatom_end, jatom_start, jatom_end, &
533 : katom_start, katom_end, latom_start, latom_end, &
534 : particle_set, &
535 : coeffs_set, coeffs_kind, &
536 : is_assoc_atomic_block_global, do_periodic, &
537 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
538 : cell, &
539 : do_p_screening, map_atom_to_kind_atom, eval_type, &
540 186918 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
541 186918 : IF (q_idx <= n_my_quartets) cached_block_cost(q_idx) = block_cost_val
542 594890 : cost_per_core = cost_per_core + block_cost_val
543 :
544 : END DO ! atom_block
545 :
546 2126 : nbins = load_balance_parameter%nbins
547 2126 : cost_per_bin = (cost_per_core + nbins - 1)/(nbins)
548 :
549 2126 : !$OMP BARRIER
550 2126 : !$OMP MASTER
551 2126 : CALL timestop(handle_inner)
552 : !$OMP END MASTER
553 2126 : !$OMP BARRIER
554 :
555 : ! new load balancing test
556 : IF (.FALSE.) THEN
557 : CALL hfx_recursive_load_balance(n_processes, my_process_id, nblocks, &
558 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
559 : particle_set, &
560 : coeffs_set, coeffs_kind, &
561 : is_assoc_atomic_block_global, do_periodic, &
562 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
563 : cell, x_data, para_env, pmax_block, &
564 : do_p_screening, map_atom_to_kind_atom, eval_type, &
565 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
566 : END IF
567 :
568 2126 : !$OMP BARRIER
569 2126 : !$OMP MASTER
570 2126 : CALL timeset(routineN//"_bin", handle_inner)
571 : !$OMP END MASTER
572 2126 : !$OMP BARRIER
573 :
574 142442 : ALLOCATE (binned_dist(nbins))
575 138190 : binned_dist(:)%istart = -1_int_8
576 138190 : binned_dist(:)%number_of_atom_quartets = 0_int_8
577 138190 : binned_dist(:)%cost = 0_int_8
578 138190 : binned_dist(:)%time_first_scf = 0.0_dp
579 138190 : binned_dist(:)%time_other_scf = 0.0_dp
580 138190 : binned_dist(:)%time_forces = 0.0_dp
581 :
582 2126 : current_cost = 0
583 2126 : mepos = 1
584 2126 : distribution_counter_start = 1
585 2126 : distribution_counter_end = 0
586 2126 : ibin = 1
587 :
588 2126 : global_quartet_counter = 0
589 2126 : local_quartet_counter = 0
590 2126 : last_bin_needs_to_be_filled = .FALSE.
591 595542 : DO atom_block = my_process_id, INT(nblocks, KIND=int_8)**4 - 1, n_processes
592 593416 : latom_block = INT(MODULO(atom_block, INT(nblocks, KIND=int_8))) + 1
593 593416 : tmp_block = atom_block/nblocks
594 593416 : katom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
595 593416 : IF (latom_block < katom_block) CYCLE
596 332701 : tmp_block = tmp_block/nblocks
597 332701 : jatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
598 332701 : tmp_block = tmp_block/nblocks
599 332701 : iatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
600 332701 : IF (jatom_block < iatom_block) CYCLE
601 :
602 187570 : distribution_counter_end = distribution_counter_end + 1
603 187570 : global_quartet_counter = global_quartet_counter + 1
604 187570 : last_bin_needs_to_be_filled = .TRUE.
605 :
606 187570 : IF (binned_dist(ibin)%istart == -1_int_8) binned_dist(ibin)%istart = atom_block
607 :
608 187570 : IF (distribution_counter_end <= n_my_quartets) THEN
609 187570 : block_cost_val = cached_block_cost(distribution_counter_end)
610 187570 : IF (block_cost_val < 0_int_8) CYCLE
611 186918 : current_cost = current_cost + block_cost_val
612 : ELSE
613 0 : iatom_start = x_data%blocks(iatom_block)%istart
614 0 : iatom_end = x_data%blocks(iatom_block)%iend
615 0 : jatom_start = x_data%blocks(jatom_block)%istart
616 0 : jatom_end = x_data%blocks(jatom_block)%iend
617 0 : katom_start = x_data%blocks(katom_block)%istart
618 0 : katom_end = x_data%blocks(katom_block)%iend
619 0 : latom_start = x_data%blocks(latom_block)%istart
620 0 : latom_end = x_data%blocks(latom_block)%iend
621 :
622 0 : SELECT CASE (eval_type)
623 : CASE (hfx_do_eval_energy)
624 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block), &
625 : pmax_block(latom_block, jatom_block), &
626 : pmax_block(latom_block, iatom_block), &
627 0 : pmax_block(katom_block, jatom_block))
628 : CASE (hfx_do_eval_forces)
629 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block) + &
630 : pmax_block(latom_block, jatom_block), &
631 : pmax_block(latom_block, iatom_block) + &
632 0 : pmax_block(katom_block, jatom_block))
633 : END SELECT
634 :
635 0 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
636 :
637 : current_cost = current_cost &
638 : + estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
639 : iatom_start, iatom_end, jatom_start, jatom_end, &
640 : katom_start, katom_end, latom_start, latom_end, &
641 : particle_set, &
642 : coeffs_set, coeffs_kind, &
643 : is_assoc_atomic_block_global, do_periodic, &
644 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
645 : cell, &
646 : do_p_screening, map_atom_to_kind_atom, eval_type, &
647 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
648 : END IF
649 :
650 189044 : IF (current_cost >= cost_per_bin) THEN
651 23009 : IF (ibin == nbins) THEN
652 : binned_dist(ibin)%number_of_atom_quartets = binned_dist(ibin)%number_of_atom_quartets + &
653 0 : distribution_counter_end - distribution_counter_start + 1
654 : ELSE
655 23009 : binned_dist(ibin)%number_of_atom_quartets = distribution_counter_end - distribution_counter_start + 1
656 : END IF
657 23009 : binned_dist(ibin)%cost = binned_dist(ibin)%cost + current_cost
658 23009 : ibin = MIN(ibin + 1, nbins)
659 23009 : distribution_counter_start = distribution_counter_end + 1
660 23009 : current_cost = 0
661 23009 : last_bin_needs_to_be_filled = .FALSE.
662 : END IF
663 : END DO
664 2126 : IF (ALLOCATED(cached_block_cost)) DEALLOCATE (cached_block_cost)
665 :
666 2126 : !$OMP BARRIER
667 2126 : !$OMP MASTER
668 2126 : CALL timestop(handle_inner)
669 2126 : CALL timeset(routineN//"_dist", handle_inner)
670 : !$OMP END MASTER
671 2126 : !$OMP BARRIER
672 : !! Fill the last bin if necessary
673 2126 : IF (last_bin_needs_to_be_filled) THEN
674 618 : binned_dist(ibin)%cost = binned_dist(ibin)%cost + current_cost
675 618 : IF (ibin == nbins) THEN
676 : binned_dist(ibin)%number_of_atom_quartets = binned_dist(ibin)%number_of_atom_quartets + &
677 0 : distribution_counter_end - distribution_counter_start + 1
678 : ELSE
679 618 : binned_dist(ibin)%number_of_atom_quartets = distribution_counter_end - distribution_counter_start + 1
680 : END IF
681 : END IF
682 :
683 : !! Sanity-Check
684 138190 : DO ibin = 1, nbins
685 138190 : local_quartet_counter = local_quartet_counter + binned_dist(ibin)%number_of_atom_quartets
686 : END DO
687 2126 : !$OMP BARRIER
688 2126 : !$OMP MASTER
689 2126 : shm_local_quartet_counter = 0
690 2126 : shm_global_quartet_counter = 0
691 : !$OMP END MASTER
692 2126 : !$OMP BARRIER
693 2126 : !$OMP ATOMIC
694 : shm_local_quartet_counter = shm_local_quartet_counter + local_quartet_counter
695 2126 : !$OMP ATOMIC
696 : shm_global_quartet_counter = shm_global_quartet_counter + global_quartet_counter
697 :
698 2126 : !$OMP BARRIER
699 2126 : !$OMP MASTER
700 2126 : CALL para_env%sum(shm_local_quartet_counter)
701 2126 : CALL para_env%sum(shm_global_quartet_counter)
702 2126 : IF (para_env%is_source()) THEN
703 1063 : IF (shm_local_quartet_counter /= shm_global_quartet_counter) THEN
704 : WRITE (error_msg, '(A,I0,A,I0,A)') "HFX Sanity check for parallel distribution failed. "// &
705 0 : "Number of local quartets (", shm_local_quartet_counter, &
706 0 : ") and number of global quartets (", shm_global_quartet_counter, &
707 0 : ") are different. Please send in a bug report."
708 0 : CPABORT(error_msg)
709 : END IF
710 : END IF
711 : !$OMP END MASTER
712 :
713 2126 : !$OMP BARRIER
714 2126 : !$OMP MASTER
715 6378 : ALLOCATE (cost_matrix(ncpu*nbins*n_threads))
716 274254 : cost_matrix = 0
717 : !$OMP END MASTER
718 2126 : !$OMP BARRIER
719 2126 : icpu = para_env%mepos + 1
720 138190 : DO i = 1, nbins
721 138190 : cost_matrix((icpu - 1)*nbins*n_threads + i_thread*nbins + i) = binned_dist(i)%cost
722 : END DO
723 2126 : mepos = para_env%mepos
724 2126 : !$OMP BARRIER
725 :
726 2126 : !$OMP MASTER
727 : ! sync before/after ring of isendrecv
728 2126 : CALL para_env%sync()
729 :
730 6378 : ALLOCATE (sendbuffer(nbins*n_threads))
731 4252 : ALLOCATE (recbuffer(nbins*n_threads))
732 :
733 274254 : sendbuffer = cost_matrix(mepos*nbins*n_threads + 1:mepos*nbins*n_threads + nbins*n_threads)
734 :
735 2126 : dest = MODULO(mepos + 1, ncpu)
736 2126 : source = MODULO(mepos - 1, ncpu)
737 6378 : DO icpu = 0, ncpu - 1
738 4252 : IF (icpu /= ncpu - 1) THEN
739 : CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
740 2126 : req(1), req(2), 13)
741 : END IF
742 4252 : data_from = MODULO(mepos - icpu, ncpu)
743 548508 : cost_matrix(data_from*nbins*n_threads + 1:data_from*nbins*n_threads + nbins*n_threads) = sendbuffer
744 4252 : IF (icpu /= ncpu - 1) THEN
745 2126 : CALL mp_waitall(req)
746 : END IF
747 4252 : swapbuffer => sendbuffer
748 4252 : sendbuffer => recbuffer
749 6378 : recbuffer => swapbuffer
750 : END DO
751 2126 : DEALLOCATE (recbuffer, sendbuffer)
752 : !$OMP END MASTER
753 2126 : !$OMP BARRIER
754 :
755 2126 : !$OMP BARRIER
756 2126 : !$OMP MASTER
757 2126 : CALL timestop(handle_inner)
758 2126 : CALL timeset(routineN//"_opt", handle_inner)
759 : !$OMP END MASTER
760 2126 : !$OMP BARRIER
761 :
762 : !! Find an optimal distribution i.e. assign each element of the cost matrix to a certain process
763 2126 : !$OMP BARRIER
764 6378 : ALLOCATE (local_cost_matrix(SIZE(cost_matrix, 1)))
765 546382 : local_cost_matrix = cost_matrix
766 2126 : !$OMP MASTER
767 6378 : ALLOCATE (shm_distribution_vector(ncpu*nbins*n_threads))
768 :
769 : CALL optimize_distribution(ncpu*nbins*n_threads, ncpu*n_threads, local_cost_matrix, &
770 2126 : shm_distribution_vector, x_data%load_balance_parameter%do_randomize)
771 :
772 2126 : CALL timestop(handle_inner)
773 2126 : CALL timeset(routineN//"_redist", handle_inner)
774 : !! Collect local data to global array
775 416696 : ALLOCATE (full_dist(ncpu*n_threads, nbins))
776 :
777 410318 : full_dist(:, :)%istart = 0_int_8
778 410318 : full_dist(:, :)%number_of_atom_quartets = 0_int_8
779 410318 : full_dist(:, :)%cost = 0_int_8
780 410318 : full_dist(:, :)%time_first_scf = 0.0_dp
781 410318 : full_dist(:, :)%time_other_scf = 0.0_dp
782 412444 : full_dist(:, :)%time_forces = 0.0_dp
783 : !$OMP END MASTER
784 2126 : !$OMP BARRIER
785 2126 : mepos = para_env%mepos + 1
786 274254 : full_dist((mepos - 1)*n_threads + i_thread + 1, :) = binned_dist(:)
787 :
788 2126 : !$OMP BARRIER
789 2126 : !$OMP MASTER
790 6378 : ALLOCATE (sendbuffer(3*nbins*n_threads))
791 4252 : ALLOCATE (recbuffer(3*nbins*n_threads))
792 2126 : mepos = para_env%mepos
793 4252 : DO j = 1, n_threads
794 140316 : DO i = 1, nbins
795 136064 : sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 1) = full_dist(mepos*n_threads + j, i)%istart
796 136064 : sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 2) = full_dist(mepos*n_threads + j, i)%number_of_atom_quartets
797 138190 : sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 3) = full_dist(mepos*n_threads + j, i)%cost
798 : END DO
799 : END DO
800 :
801 : ! sync before/after ring of isendrecv
802 2126 : CALL para_env%sync()
803 2126 : dest = MODULO(mepos + 1, ncpu)
804 2126 : source = MODULO(mepos - 1, ncpu)
805 6378 : DO icpu = 0, ncpu - 1
806 4252 : IF (icpu /= ncpu - 1) THEN
807 : CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
808 2126 : req(1), req(2), 13)
809 : END IF
810 4252 : data_from = MODULO(mepos - icpu, ncpu)
811 8504 : DO j = 1, n_threads
812 280632 : DO i = 1, nbins
813 272128 : full_dist(data_from*n_threads + j, i)%istart = sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 1)
814 272128 : full_dist(data_from*n_threads + j, i)%number_of_atom_quartets = sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 2)
815 276380 : full_dist(data_from*n_threads + j, i)%cost = sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 3)
816 : END DO
817 : END DO
818 :
819 4252 : IF (icpu /= ncpu - 1) THEN
820 2126 : CALL mp_waitall(req)
821 : END IF
822 4252 : swapbuffer => sendbuffer
823 4252 : sendbuffer => recbuffer
824 6378 : recbuffer => swapbuffer
825 : END DO
826 2126 : DEALLOCATE (recbuffer, sendbuffer)
827 :
828 : ! sync before/after ring of isendrecv
829 2126 : CALL para_env%sync()
830 : !$OMP END MASTER
831 2126 : !$OMP BARRIER
832 : !! reorder the distribution according to the distribution vector
833 6378 : ALLOCATE (tmp_pos(ncpu*n_threads))
834 6378 : tmp_pos = 1
835 278506 : ALLOCATE (tmp_dist(nbins*ncpu*n_threads))
836 :
837 274254 : tmp_dist(:)%istart = 0_int_8
838 274254 : tmp_dist(:)%number_of_atom_quartets = 0_int_8
839 274254 : tmp_dist(:)%cost = 0_int_8
840 274254 : tmp_dist(:)%time_first_scf = 0.0_dp
841 274254 : tmp_dist(:)%time_other_scf = 0.0_dp
842 274254 : tmp_dist(:)%time_forces = 0.0_dp
843 :
844 6378 : DO icpu = 1, n_processes
845 278506 : DO i = 1, nbins
846 272128 : mepos = my_process_id + 1
847 276380 : IF (shm_distribution_vector((icpu - 1)*nbins + i) == mepos) THEN
848 136064 : tmp_dist(tmp_pos(mepos)) = full_dist(icpu, i)
849 136064 : tmp_pos(mepos) = tmp_pos(mepos) + 1
850 : END IF
851 : END DO
852 : END DO
853 :
854 : !! Assign the load to each process
855 : NULLIFY (ptr_to_tmp_dist)
856 2126 : mepos = my_process_id + 1
857 2126 : ptr_to_tmp_dist => tmp_dist(1:tmp_pos(mepos) - 1)
858 3614 : SELECT CASE (eval_type)
859 : CASE (hfx_do_eval_energy)
860 1488 : CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
861 : CASE (hfx_do_eval_forces)
862 2126 : CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
863 : END SELECT
864 :
865 2126 : !$OMP BARRIER
866 2126 : !$OMP MASTER
867 2126 : DEALLOCATE (full_dist, cost_matrix, shm_distribution_vector)
868 : !$OMP END MASTER
869 2126 : !$OMP BARRIER
870 2126 : DEALLOCATE (tmp_dist, tmp_pos)
871 2126 : DEALLOCATE (binned_dist, local_cost_matrix)
872 2126 : DEALLOCATE (set_list_ij, set_list_kl)
873 :
874 2126 : !$OMP BARRIER
875 2126 : !$OMP MASTER
876 2126 : CALL timestop(handle_inner)
877 : !$OMP END MASTER
878 2126 : !$OMP BARRIER
879 : END IF
880 2128 : !$OMP BARRIER
881 2128 : !$OMP MASTER
882 2128 : CALL timestop(handle)
883 : !$OMP END MASTER
884 2128 : !$OMP BARRIER
885 4415600 : END SUBROUTINE hfx_load_balance
886 :
887 : ! **************************************************************************************************
888 : !> \brief Reference implementation of new recursive load balancing routine
889 : !> Computes a local list of atom_blocks (p_atom_blocks,q_atom_blocks) for
890 : !> each process in a P-Q grid such that every process has more or less the
891 : !> same amount of work. Has no output at the moment (not used) but writes
892 : !> its computed load balance values into a file. Possible output is ready
893 : !> to use in the two arrays p_atom_blocks & q_atom_blocks
894 : !> \param n_processes ...
895 : !> \param my_process_id ...
896 : !> \param nblocks ...
897 : !> \param natom ...
898 : !> \param nkind ...
899 : !> \param list_ij ...
900 : !> \param list_kl ...
901 : !> \param set_list_ij ...
902 : !> \param set_list_kl ...
903 : !> \param particle_set ...
904 : !> \param coeffs_set ...
905 : !> \param coeffs_kind ...
906 : !> \param is_assoc_atomic_block_global ...
907 : !> \param do_periodic ...
908 : !> \param kind_of ...
909 : !> \param basis_parameter ...
910 : !> \param pmax_set ...
911 : !> \param pmax_atom ...
912 : !> \param pmax_blocks ...
913 : !> \param cell ...
914 : !> \param x_data ...
915 : !> \param para_env ...
916 : !> \param pmax_block ...
917 : !> \param do_p_screening ...
918 : !> \param map_atom_to_kind_atom ...
919 : !> \param eval_type ...
920 : !> \param log10_eps_schwarz ...
921 : !> \param log_2 ...
922 : !> \param coeffs_kind_max0 ...
923 : !> \param use_virial ...
924 : !> \param atomic_pair_list ...
925 : !> \par History
926 : !> 03.2011 created [Michael Steinlechner]
927 : !> \author Michael Steinlechner
928 : ! **************************************************************************************************
929 :
930 0 : SUBROUTINE hfx_recursive_load_balance(n_processes, my_process_id, nblocks, &
931 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
932 : particle_set, &
933 : coeffs_set, coeffs_kind, &
934 0 : is_assoc_atomic_block_global, do_periodic, &
935 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
936 : cell, x_data, para_env, pmax_block, &
937 : do_p_screening, map_atom_to_kind_atom, eval_type, &
938 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
939 :
940 : ! input variables:
941 : INTEGER, INTENT(IN) :: n_processes, my_process_id, nblocks, &
942 : natom, nkind
943 : TYPE(pair_list_type), INTENT(IN) :: list_ij, list_kl
944 : TYPE(pair_set_list_type), ALLOCATABLE, &
945 : DIMENSION(:), INTENT(IN) :: set_list_ij, set_list_kl
946 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
947 : POINTER :: particle_set
948 : TYPE(hfx_screen_coeff_type), &
949 : DIMENSION(:, :, :, :), INTENT(IN), POINTER :: coeffs_set
950 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
951 : INTENT(IN), POINTER :: coeffs_kind
952 : INTEGER, DIMENSION(:, :), INTENT(IN) :: is_assoc_atomic_block_global
953 : LOGICAL, INTENT(IN) :: do_periodic
954 : INTEGER, INTENT(IN) :: kind_of(*)
955 : TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
956 : POINTER :: basis_parameter
957 : TYPE(hfx_p_kind), DIMENSION(:), INTENT(IN), &
958 : POINTER :: pmax_set
959 : REAL(dp), DIMENSION(:, :), INTENT(IN), POINTER :: pmax_atom
960 : REAL(dp) :: pmax_blocks
961 : TYPE(cell_type), INTENT(IN), POINTER :: cell
962 : TYPE(hfx_type), INTENT(IN), POINTER :: x_data
963 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
964 : REAL(dp), DIMENSION(:, :), INTENT(IN), POINTER :: pmax_block
965 : LOGICAL, INTENT(IN) :: do_p_screening
966 : INTEGER, DIMENSION(:), INTENT(IN), POINTER :: map_atom_to_kind_atom
967 : INTEGER, INTENT(IN) :: eval_type
968 : REAL(dp), INTENT(IN) :: log10_eps_schwarz, log_2, &
969 : coeffs_kind_max0
970 : LOGICAL, INTENT(IN) :: use_virial
971 : LOGICAL, DIMENSION(:, :), INTENT(IN), POINTER :: atomic_pair_list
972 :
973 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_recursive_load_balance'
974 :
975 : INTEGER :: handle, i, iatom_block, iatom_end, iatom_start, j, jatom_block, jatom_end, &
976 : jatom_start, katom_block, katom_end, katom_start, latom_block, latom_end, latom_start, &
977 : nP, nQ, numBins, p, q, sizeP, sizeQ, unit_nr
978 : INTEGER(int_8) :: local_cost, pidx, qidx, sumP, sumQ
979 0 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: local_cost_vector
980 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blocksize, p_atom_blocks, permute, &
981 0 : q_atom_blocks
982 : REAL(dp) :: maximum, mean
983 :
984 : ! internal variables:
985 :
986 0 : !$OMP BARRIER
987 0 : !$OMP MASTER
988 0 : CALL timeset(routineN, handle)
989 : !$OMP END MASTER
990 0 : !$OMP BARRIER
991 :
992 : ! calculate best p/q distribution grid for the n_processes
993 0 : CALL hfx_calculate_PQ(p, q, numBins, n_processes)
994 :
995 0 : ALLOCATE (blocksize(numBins))
996 0 : ALLOCATE (permute(nblocks**2))
997 0 : DO i = 1, nblocks**2
998 0 : permute(i) = i
999 : END DO
1000 :
1001 : ! call the main recursive permutation routine.
1002 : ! Output:
1003 : ! blocksize :: vector (size numBins) with the sizes for each column/row block
1004 : ! permute :: permutation vector
1005 : CALL hfx_recursive_permute(blocksize, 1, nblocks**2, numBins, &
1006 : permute, 1, &
1007 : my_process_id, n_processes, nblocks, &
1008 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
1009 : particle_set, &
1010 : coeffs_set, coeffs_kind, &
1011 : is_assoc_atomic_block_global, do_periodic, &
1012 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1013 : cell, x_data, para_env, pmax_block, &
1014 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1015 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1016 :
1017 : ! number of blocks per processor in p-direction (vertical)
1018 0 : nP = numBins/p
1019 : ! number of blocks per processor in q-direction (horizontal)
1020 0 : nQ = numBins/q
1021 :
1022 : ! calc own position in P-Q-processor grid (PQ-grid is column-major)
1023 0 : pidx = MODULO(INT(my_process_id), INT(p)) + 1
1024 0 : qidx = my_process_id/p + 1
1025 :
1026 0 : sizeP = SUM(blocksize((nP*(pidx - 1) + 1):(nP*pidx)))
1027 0 : sizeQ = SUM(blocksize((nQ*(qidx - 1) + 1):(nQ*qidx)))
1028 :
1029 0 : sumP = SUM(blocksize(1:(nP*(pidx - 1))))
1030 0 : sumQ = SUM(blocksize(1:(nQ*(qidx - 1))))
1031 :
1032 0 : ALLOCATE (p_atom_blocks(sizeP))
1033 0 : ALLOCATE (q_atom_blocks(sizeQ))
1034 :
1035 0 : p_atom_blocks(:) = permute((sumP + 1):(sumP + sizeP))
1036 0 : q_atom_blocks(:) = permute((sumQ + 1):(sumQ + sizeQ))
1037 :
1038 : ! from here on, we are actually finished, each process has been
1039 : ! assigned a (p_atom_blocks,q_atom_blocks) pair list.
1040 : ! what follows is just a small routine to calculate the local cost
1041 : ! for each processor which is then written to a file.
1042 :
1043 : ! calculate local cost for each processor!
1044 : ! ****************************************
1045 : local_cost = 0
1046 0 : DO i = 1, sizeP
1047 0 : DO j = 1, sizeQ
1048 :
1049 : ! get corresponding 4D block indices out of our own P-Q-block
1050 0 : latom_block = MODULO(q_atom_blocks(j), nblocks)
1051 0 : iatom_block = q_atom_blocks(j)/nblocks + 1
1052 0 : jatom_block = MODULO(p_atom_blocks(i), nblocks)
1053 0 : katom_block = p_atom_blocks(i)/nblocks + 1
1054 :
1055 : ! symmetry checks.
1056 0 : IF (latom_block < katom_block) CYCLE
1057 0 : IF (jatom_block < iatom_block) CYCLE
1058 :
1059 0 : iatom_start = x_data%blocks(iatom_block)%istart
1060 0 : iatom_end = x_data%blocks(iatom_block)%iend
1061 0 : jatom_start = x_data%blocks(jatom_block)%istart
1062 0 : jatom_end = x_data%blocks(jatom_block)%iend
1063 0 : katom_start = x_data%blocks(katom_block)%istart
1064 0 : katom_end = x_data%blocks(katom_block)%iend
1065 0 : latom_start = x_data%blocks(latom_block)%istart
1066 0 : latom_end = x_data%blocks(latom_block)%iend
1067 :
1068 : ! whatever.
1069 0 : SELECT CASE (eval_type)
1070 : CASE (hfx_do_eval_energy)
1071 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block), &
1072 : pmax_block(latom_block, jatom_block), &
1073 : pmax_block(latom_block, iatom_block), &
1074 0 : pmax_block(katom_block, jatom_block))
1075 : CASE (hfx_do_eval_forces)
1076 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block) + &
1077 : pmax_block(latom_block, jatom_block), &
1078 : pmax_block(latom_block, iatom_block) + &
1079 0 : pmax_block(katom_block, jatom_block))
1080 : END SELECT
1081 :
1082 : ! screening.
1083 0 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
1084 :
1085 : ! estimate the cost of this atom_block.
1086 : local_cost = local_cost + estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, &
1087 : set_list_kl, &
1088 : iatom_start, iatom_end, jatom_start, jatom_end, &
1089 : katom_start, katom_end, latom_start, latom_end, &
1090 : particle_set, &
1091 : coeffs_set, coeffs_kind, &
1092 : is_assoc_atomic_block_global, do_periodic, &
1093 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1094 : cell, &
1095 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1096 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1097 : END DO
1098 : END DO
1099 :
1100 0 : ALLOCATE (local_cost_vector(n_processes))
1101 0 : local_cost_vector = 0
1102 0 : local_cost_vector(my_process_id + 1) = local_cost
1103 0 : CALL para_env%sum(local_cost_vector)
1104 :
1105 0 : mean = SUM(local_cost_vector)/n_processes
1106 0 : maximum = MAXVAL(local_cost_vector)
1107 :
1108 0 : !$OMP BARRIER
1109 0 : !$OMP MASTER
1110 : ! only output once
1111 0 : IF (my_process_id == 0) THEN
1112 0 : CALL open_file(unit_number=unit_nr, file_name="loads.dat")
1113 0 : WRITE (unit_nr, *) 'maximum cost:', maximum
1114 0 : WRITE (unit_nr, *) 'mean cost:', mean
1115 0 : WRITE (unit_nr, *) 'load balance ratio max/mean: ', maximum/mean
1116 0 : WRITE (unit_nr, *) '-------- detailed per-process costs ---------'
1117 0 : DO i = 1, n_processes
1118 0 : WRITE (unit_nr, *) local_cost_vector(i)
1119 : END DO
1120 0 : CALL close_file(unit_nr)
1121 : END IF
1122 : !$OMP END MASTER
1123 0 : !$OMP BARRIER
1124 :
1125 0 : DEALLOCATE (local_cost_vector)
1126 0 : DEALLOCATE (p_atom_blocks, q_atom_blocks)
1127 0 : DEALLOCATE (blocksize, permute)
1128 :
1129 0 : !$OMP BARRIER
1130 0 : !$OMP MASTER
1131 0 : CALL timestop(handle)
1132 : !$OMP END MASTER
1133 0 : !$OMP BARRIER
1134 :
1135 0 : END SUBROUTINE hfx_recursive_load_balance
1136 :
1137 : ! **************************************************************************************************
1138 : !> \brief Small routine to calculate the optimal P-Q-processor grid distribution
1139 : !> for a given number of processors N
1140 : !> and the corresponding number of Bins for the load balancing routine
1141 : !> \param p number of rows on P-Q process grid (output)
1142 : !> \param q number of columns on P-Q process grid (output)
1143 : !> \param nBins number of Bins (output)
1144 : !> \param N number of processes (input)
1145 : !> \par History
1146 : !> 03.2011 created [Michael Steinlechner]
1147 : !> \author Michael Steinlechner
1148 : ! **************************************************************************************************
1149 0 : SUBROUTINE hfx_calculate_PQ(p, q, nBins, N)
1150 :
1151 : INTEGER, INTENT(OUT) :: p, q, nBins
1152 : INTEGER, INTENT(IN) :: N
1153 :
1154 : INTEGER :: a, b, k
1155 : REAL(dp) :: sqN
1156 :
1157 0 : k = 2
1158 0 : sqN = SQRT(REAL(N, KIND=dp))
1159 0 : p = 1
1160 :
1161 0 : DO WHILE (REAL(k, KIND=dp) <= sqN)
1162 0 : IF (MODULO(N, k) == 0) THEN
1163 0 : p = k
1164 : END IF
1165 0 : k = k + 1
1166 : END DO
1167 0 : q = N/p
1168 :
1169 : ! now compute the least common multiple of p & q to get the number of necessary bins
1170 : ! compute using the relation LCM(p,q) = abs(p*q) / GCD(p,q)
1171 : ! and use euclid's algorithm for GCD computation.
1172 0 : a = p
1173 0 : b = q
1174 :
1175 0 : DO WHILE (b /= 0)
1176 0 : IF (a > b) THEN
1177 0 : a = a - b
1178 : ELSE
1179 0 : b = b - a
1180 : END IF
1181 : END DO
1182 : ! gcd(p,q) is now saved in a
1183 :
1184 0 : nBins = p*q/a
1185 :
1186 0 : END SUBROUTINE hfx_calculate_PQ
1187 :
1188 : ! **************************************************************************************************
1189 : !> \brief Recursive permutation routine for the load balancing of the integral
1190 : !> computation
1191 : !> \param blocksize vector of blocksizes, size(nProc), which contains for
1192 : !> each process the local blocksize (OUTPUT)
1193 : !> \param blockstart starting row/column idx of the block which is to be examined
1194 : !> at this point (INPUT)
1195 : !> \param blockend ending row/column idx of the block which is to be examined
1196 : !> (INPUT)
1197 : !> \param nProc_in number of bins into which the current block has to be divided
1198 : !> (INPUT)
1199 : !> \param permute permutation vector which balances column/row cost
1200 : !> size(nblocks^2). (OUTPUT)
1201 : !> \param step ...
1202 : !> \param my_process_id ...
1203 : !> \param n_processes ...
1204 : !> \param nblocks ...
1205 : !> \param natom ...
1206 : !> \param nkind ...
1207 : !> \param list_ij ...
1208 : !> \param list_kl ...
1209 : !> \param set_list_ij ...
1210 : !> \param set_list_kl ...
1211 : !> \param particle_set ...
1212 : !> \param coeffs_set ...
1213 : !> \param coeffs_kind ...
1214 : !> \param is_assoc_atomic_block_global ...
1215 : !> \param do_periodic ...
1216 : !> \param kind_of ...
1217 : !> \param basis_parameter ...
1218 : !> \param pmax_set ...
1219 : !> \param pmax_atom ...
1220 : !> \param pmax_blocks ...
1221 : !> \param cell ...
1222 : !> \param x_data ...
1223 : !> \param para_env ...
1224 : !> \param pmax_block ...
1225 : !> \param do_p_screening ...
1226 : !> \param map_atom_to_kind_atom ...
1227 : !> \param eval_type ...
1228 : !> \param log10_eps_schwarz ...
1229 : !> \param log_2 ...
1230 : !> \param coeffs_kind_max0 ...
1231 : !> \param use_virial ...
1232 : !> \param atomic_pair_list ...
1233 : !> \par History
1234 : !> 03.2011 created [Michael Steinlechner]
1235 : !> \author Michael Steinlechner
1236 : ! **************************************************************************************************
1237 0 : RECURSIVE SUBROUTINE hfx_recursive_permute(blocksize, blockstart, blockend, nProc_in, &
1238 0 : permute, step, &
1239 : my_process_id, n_processes, nblocks, &
1240 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
1241 : particle_set, &
1242 : coeffs_set, coeffs_kind, &
1243 0 : is_assoc_atomic_block_global, do_periodic, &
1244 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1245 : cell, x_data, para_env, pmax_block, &
1246 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1247 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1248 :
1249 : INTEGER :: nProc_in, blockend, blockstart
1250 : INTEGER, DIMENSION(nProc_in) :: blocksize
1251 : INTEGER :: nblocks, n_processes, my_process_id
1252 : INTEGER, INTENT(IN) :: step
1253 : INTEGER, DIMENSION(nblocks*nblocks) :: permute
1254 : INTEGER :: natom
1255 : INTEGER, INTENT(IN) :: nkind
1256 : TYPE(pair_list_type) :: list_ij, list_kl
1257 : TYPE(pair_set_list_type), ALLOCATABLE, &
1258 : DIMENSION(:) :: set_list_ij, set_list_kl
1259 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1260 : TYPE(hfx_screen_coeff_type), &
1261 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
1262 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
1263 : POINTER :: coeffs_kind
1264 : INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block_global
1265 : LOGICAL :: do_periodic
1266 : INTEGER :: kind_of(*)
1267 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1268 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set
1269 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom
1270 : REAL(dp) :: pmax_blocks
1271 : TYPE(cell_type), POINTER :: cell
1272 : TYPE(hfx_type), POINTER :: x_data
1273 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1274 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_block
1275 : LOGICAL, INTENT(IN) :: do_p_screening
1276 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom
1277 : INTEGER, INTENT(IN) :: eval_type
1278 : REAL(dp) :: log10_eps_schwarz, log_2, &
1279 : coeffs_kind_max0
1280 : LOGICAL, INTENT(IN) :: use_virial
1281 : LOGICAL, DIMENSION(:, :), POINTER :: atomic_pair_list
1282 :
1283 : INTEGER :: col, endoffset, i, iatom_block, iatom_end, iatom_start, idx, inv_perm, &
1284 : jatom_block, jatom_end, jatom_start, katom_block, katom_end, katom_start, latom_block, &
1285 : latom_end, latom_start, nbins, nProc, row, startoffset
1286 : INTEGER(int_8) :: atom_block, tmp_block
1287 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ithblocksize, localblocksize
1288 0 : INTEGER, DIMENSION(blockend - blockstart + 1) :: bin_perm, tmp_perm
1289 : REAL(dp) :: partialcost
1290 0 : REAL(dp), DIMENSION(nblocks*nblocks) :: cost_vector
1291 :
1292 : nProc = nProc_in
1293 0 : cost_vector = 0.0_dp
1294 :
1295 : ! loop over local atom_blocks.
1296 0 : DO atom_block = my_process_id, INT(nblocks, KIND=int_8)**4 - 1, n_processes
1297 :
1298 : ! get corresponding 4D block indices
1299 0 : latom_block = INT(MODULO(atom_block, INT(nblocks, KIND=int_8))) + 1
1300 0 : tmp_block = atom_block/nblocks
1301 0 : katom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
1302 0 : IF (latom_block < katom_block) CYCLE
1303 0 : tmp_block = tmp_block/nblocks
1304 0 : jatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
1305 0 : tmp_block = tmp_block/nblocks
1306 0 : iatom_block = INT(MODULO(tmp_block, INT(nblocks, KIND=int_8))) + 1
1307 0 : IF (jatom_block < iatom_block) CYCLE
1308 :
1309 : ! get 2D indices of this atom_block (with permutation applied)
1310 : ! for this, we need to invert the permutation, this means
1311 : ! find position in permutation vector where value==idx
1312 :
1313 0 : row = (katom_block - 1)*nblocks + jatom_block
1314 0 : inv_perm = 1
1315 0 : DO WHILE (permute(inv_perm) /= row)
1316 0 : inv_perm = inv_perm + 1
1317 : END DO
1318 0 : row = inv_perm
1319 :
1320 0 : col = (iatom_block - 1)*nblocks + latom_block
1321 0 : inv_perm = 1
1322 0 : DO WHILE (permute(inv_perm) /= col)
1323 0 : inv_perm = inv_perm + 1
1324 : END DO
1325 0 : col = inv_perm
1326 :
1327 : ! if row/col outside our current diagonal block, skip calculation.
1328 0 : IF (col < blockstart .OR. col > blockend) CYCLE
1329 0 : IF (row < blockstart .OR. row > blockend) CYCLE
1330 :
1331 0 : iatom_start = x_data%blocks(iatom_block)%istart
1332 0 : iatom_end = x_data%blocks(iatom_block)%iend
1333 0 : jatom_start = x_data%blocks(jatom_block)%istart
1334 0 : jatom_end = x_data%blocks(jatom_block)%iend
1335 0 : katom_start = x_data%blocks(katom_block)%istart
1336 0 : katom_end = x_data%blocks(katom_block)%iend
1337 0 : latom_start = x_data%blocks(latom_block)%istart
1338 0 : latom_end = x_data%blocks(latom_block)%iend
1339 :
1340 : ! whatever.
1341 0 : SELECT CASE (eval_type)
1342 : CASE (hfx_do_eval_energy)
1343 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block), &
1344 : pmax_block(latom_block, jatom_block), &
1345 : pmax_block(latom_block, iatom_block), &
1346 0 : pmax_block(katom_block, jatom_block))
1347 : CASE (hfx_do_eval_forces)
1348 : pmax_blocks = MAX(pmax_block(katom_block, iatom_block) + &
1349 : pmax_block(latom_block, jatom_block), &
1350 : pmax_block(latom_block, iatom_block) + &
1351 0 : pmax_block(katom_block, jatom_block))
1352 : END SELECT
1353 :
1354 : ! screening.
1355 0 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
1356 :
1357 : ! every second recursion step, compute row sum instead of column sum
1358 :
1359 0 : IF (MODULO(step, 2) == 0) THEN
1360 : idx = row
1361 : ELSE
1362 0 : idx = col
1363 : END IF
1364 :
1365 : ! estimate the cost of this atom_block.
1366 : partialcost = estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, &
1367 : set_list_kl, &
1368 : iatom_start, iatom_end, jatom_start, jatom_end, &
1369 : katom_start, katom_end, latom_start, latom_end, &
1370 : particle_set, &
1371 : coeffs_set, coeffs_kind, &
1372 : is_assoc_atomic_block_global, do_periodic, &
1373 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1374 : cell, &
1375 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1376 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1377 :
1378 0 : cost_vector(idx) = cost_vector(idx) + partialcost
1379 : END DO ! atom_block
1380 :
1381 : ! sum costvector over all processes
1382 0 : CALL para_env%sum(cost_vector)
1383 :
1384 : ! calculate next prime factor of nProc
1385 0 : nBins = 2
1386 0 : DO WHILE (MODULO(INT(nProc), INT(nBins)) /= 0)
1387 0 : nBins = nBins + 1
1388 : END DO
1389 :
1390 0 : nProc = nProc/nBins
1391 :
1392 : ! ... do the binning...
1393 :
1394 0 : ALLOCATE (localblocksize(nBins))
1395 0 : CALL hfx_permute_binning(nBins, cost_vector(blockstart:blockend), blockend - blockstart + 1, bin_perm, localblocksize)
1396 :
1397 : !... and update the permutation vector
1398 :
1399 0 : tmp_perm = permute(blockstart:blockend)
1400 0 : permute(blockstart:blockend) = tmp_perm(bin_perm)
1401 :
1402 : ! split recursion into the nBins Bins
1403 0 : IF (nProc > 1) THEN
1404 0 : ALLOCATE (ithblocksize(nProc))
1405 0 : DO i = 1, nBins
1406 0 : startoffset = SUM(localblocksize(1:(i - 1)))
1407 0 : endoffset = SUM(localblocksize(1:i)) - 1
1408 :
1409 : CALL hfx_recursive_permute(ithblocksize, blockstart + startoffset, blockstart + endoffset, nProc, &
1410 : permute, step + 1, &
1411 : my_process_id, n_processes, nblocks, &
1412 : natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
1413 : particle_set, &
1414 : coeffs_set, coeffs_kind, &
1415 : is_assoc_atomic_block_global, do_periodic, &
1416 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1417 : cell, x_data, para_env, pmax_block, &
1418 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1419 0 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1420 0 : blocksize(((i - 1)*nProc + 1):(i*nProc)) = ithblocksize
1421 : END DO
1422 0 : DEALLOCATE (ithblocksize)
1423 : ELSE
1424 0 : DO i = 1, nBins
1425 0 : blocksize(i) = localblocksize(i)
1426 : END DO
1427 : END IF
1428 :
1429 0 : DEALLOCATE (localblocksize)
1430 :
1431 0 : END SUBROUTINE hfx_recursive_permute
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief small binning routine for the recursive load balancing
1435 : !>
1436 : !> \param nBins number of Bins (INPUT)
1437 : !> \param costvector vector of current row/column costs which have to be binned (INPUT)
1438 : !> \param maxbinsize upper bound for bin size (INPUT)
1439 : !> \param perm resulting permutation due to be binning routine (OUTPUT)
1440 : !> \param block_count vector of size(nbins) which contains the size of each bin (OUTPUT)
1441 : !> \par History
1442 : !> 03.2011 created [Michael Steinlechner]
1443 : !> \author Michael Steinlechner
1444 : ! **************************************************************************************************
1445 0 : SUBROUTINE hfx_permute_binning(nBins, costvector, maxbinsize, perm, block_count)
1446 :
1447 : INTEGER, INTENT(IN) :: nBins, maxbinsize
1448 : REAL(dp), DIMENSION(maxbinsize), INTENT(IN) :: costvector
1449 : INTEGER, DIMENSION(maxbinsize), INTENT(OUT) :: perm
1450 : INTEGER, DIMENSION(nBins), INTENT(OUT) :: block_count
1451 :
1452 : INTEGER :: i, j, mod_idx, offset
1453 0 : INTEGER, DIMENSION(nBins, maxbinsize) :: bin
1454 0 : INTEGER, DIMENSION(nBins) :: bin_idx
1455 0 : INTEGER, DIMENSION(maxbinsize) :: idx
1456 0 : REAL(dp), DIMENSION(maxbinsize) :: vec
1457 0 : REAL(dp), DIMENSION(nBins) :: bincosts
1458 :
1459 : ! be careful not to change costvector (copy it!)
1460 :
1461 0 : vec = costvector
1462 0 : block_count = 0
1463 0 : bincosts = 0
1464 :
1465 : !sort the array (ascending)
1466 0 : CALL sort(vec, maxbinsize, idx)
1467 :
1468 : ! count the loop down to distribute the largest cols/rows first
1469 0 : DO i = maxbinsize, 1, -1
1470 0 : IF (vec(i) == 0) THEN
1471 : ! spread zero-cost col/rows evenly among procs
1472 0 : mod_idx = MODULO(i, nBins) + 1 !(note the fortran offset by one!)
1473 0 : block_count(mod_idx) = block_count(mod_idx) + 1
1474 0 : bin(mod_idx, block_count(mod_idx)) = idx(i)
1475 : ELSE
1476 : ! sort the bins so that the one with the lowest cost is at the
1477 : ! first place, where we then assign the current col/row
1478 0 : CALL sort(bincosts, nBins, bin_idx)
1479 0 : block_count = block_count(bin_idx)
1480 0 : bin = bin(bin_idx, :)
1481 :
1482 0 : bincosts(1) = bincosts(1) + vec(i)
1483 0 : block_count(1) = block_count(1) + 1
1484 0 : bin(1, block_count(1)) = idx(i)
1485 : END IF
1486 : END DO
1487 :
1488 : ! construct permutation vector from the binning
1489 : offset = 0
1490 0 : DO i = 1, nBins
1491 0 : DO j = 1, block_count(i)
1492 0 : perm(offset + j) = bin(i, j)
1493 : END DO
1494 0 : offset = offset + block_count(i)
1495 : END DO
1496 :
1497 0 : END SUBROUTINE hfx_permute_binning
1498 :
1499 : ! **************************************************************************************************
1500 : !> \brief Cheap way of redistributing the eri's
1501 : !> \param x_data Object that stores the indices array
1502 : !> \param para_env para_env
1503 : !> \param load_balance_parameter contains parmameter for Monte-Carlo routines
1504 : !> \param i_thread current thread ID
1505 : !> \param n_threads Total Number of threads
1506 : !> \param eval_type ...
1507 : !> \par History
1508 : !> 12.2007 created [Manuel Guidon]
1509 : !> 02.2009 optimize Memory Usage [Manuel Guidon]
1510 : !> \author Manuel Guidon
1511 : !> \note
1512 : !> The cost matrix is given by the walltime for each bin that is measured
1513 : !> during the calculation
1514 : ! **************************************************************************************************
1515 1982 : SUBROUTINE hfx_update_load_balance(x_data, para_env, &
1516 : load_balance_parameter, &
1517 : i_thread, n_threads, eval_type)
1518 :
1519 : TYPE(hfx_type), POINTER :: x_data
1520 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1521 : TYPE(hfx_load_balance_type) :: load_balance_parameter
1522 : INTEGER, INTENT(IN) :: i_thread, n_threads, eval_type
1523 :
1524 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_update_load_balance'
1525 :
1526 : INTEGER :: data_from, dest, end_idx, handle, i, ibin, icpu, iprocess, j, mepos, my_bin_size, &
1527 : my_global_start_idx, my_process_id, n_processes, nbins, ncpu, source, start_idx
1528 5946 : TYPE(mp_request_type), DIMENSION(2) :: req
1529 1982 : INTEGER(int_8), DIMENSION(:), POINTER :: local_cost_matrix, recbuffer, &
1530 1982 : sendbuffer, swapbuffer
1531 : INTEGER(int_8), DIMENSION(:), POINTER, SAVE :: cost_matrix
1532 1982 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_pos
1533 : INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: bins_per_rank
1534 : INTEGER, ALLOCATABLE, DIMENSION(:, :), SAVE :: bin_histogram
1535 : INTEGER, DIMENSION(:), POINTER, SAVE :: shm_distribution_vector
1536 : INTEGER, SAVE :: max_bin_size
1537 3964 : TYPE(hfx_distribution), DIMENSION(:), POINTER :: binned_dist, ptr_to_tmp_dist, tmp_dist
1538 : TYPE(hfx_distribution), DIMENSION(:, :), POINTER, &
1539 : SAVE :: full_dist
1540 :
1541 1982 : !$OMP BARRIER
1542 3964 : !$OMP MASTER
1543 1982 : CALL timeset(routineN, handle)
1544 : !$OMP END MASTER
1545 1982 : !$OMP BARRIER
1546 :
1547 1982 : ncpu = para_env%num_pe
1548 1982 : n_processes = ncpu*n_threads
1549 : !! If there is only 1 cpu skip the binning
1550 1982 : IF (n_processes == 1) THEN
1551 0 : ALLOCATE (tmp_dist(1))
1552 0 : tmp_dist(1)%number_of_atom_quartets = HUGE(tmp_dist(1)%number_of_atom_quartets)
1553 0 : tmp_dist(1)%istart = 0_int_8
1554 0 : ptr_to_tmp_dist => tmp_dist(:)
1555 0 : SELECT CASE (eval_type)
1556 : CASE (hfx_do_eval_energy)
1557 0 : CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
1558 : CASE (hfx_do_eval_forces)
1559 0 : CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
1560 : END SELECT
1561 0 : DEALLOCATE (tmp_dist)
1562 : ELSE
1563 1982 : mepos = para_env%mepos
1564 1982 : my_process_id = para_env%mepos*n_threads + i_thread
1565 1982 : nbins = load_balance_parameter%nbins
1566 1982 : !$OMP MASTER
1567 5946 : ALLOCATE (bin_histogram(n_processes, 2))
1568 1982 : bin_histogram = 0
1569 : !$OMP END MASTER
1570 1982 : !$OMP BARRIER
1571 3088 : SELECT CASE (eval_type)
1572 : CASE (hfx_do_eval_energy)
1573 1106 : my_bin_size = SIZE(x_data%distribution_energy)
1574 : CASE (hfx_do_eval_forces)
1575 1982 : my_bin_size = SIZE(x_data%distribution_forces)
1576 : END SELECT
1577 1982 : bin_histogram(my_process_id + 1, 1) = my_bin_size
1578 1982 : !$OMP BARRIER
1579 1982 : !$OMP MASTER
1580 1982 : CALL para_env%sum(bin_histogram(:, 1))
1581 1982 : bin_histogram(1, 2) = bin_histogram(1, 1)
1582 3964 : DO iprocess = 2, n_processes
1583 3964 : bin_histogram(iprocess, 2) = bin_histogram(iprocess - 1, 2) + bin_histogram(iprocess, 1)
1584 : END DO
1585 :
1586 3964 : max_bin_size = MAXVAL(bin_histogram(para_env%mepos*n_threads + 1:para_env%mepos*n_threads + n_threads, 1))
1587 1982 : CALL para_env%max(max_bin_size)
1588 : !$OMP END MASTER
1589 1982 : !$OMP BARRIER
1590 132788 : ALLOCATE (binned_dist(my_bin_size))
1591 : !! Use old binned_dist, but with timings cost
1592 1106 : SELECT CASE (eval_type)
1593 : CASE (hfx_do_eval_energy)
1594 142674 : binned_dist = x_data%distribution_energy
1595 : CASE (hfx_do_eval_forces)
1596 114110 : binned_dist = x_data%distribution_forces
1597 : END SELECT
1598 :
1599 128830 : DO ibin = 1, my_bin_size
1600 128830 : IF (binned_dist(ibin)%number_of_atom_quartets == 0) THEN
1601 109449 : binned_dist(ibin)%cost = 0
1602 : ELSE
1603 10191 : SELECT CASE (eval_type)
1604 : CASE (hfx_do_eval_energy)
1605 10191 : IF (.NOT. load_balance_parameter%rtp_redistribute) THEN
1606 : binned_dist(ibin)%cost = INT((binned_dist(ibin)%time_first_scf + &
1607 10191 : binned_dist(ibin)%time_other_scf)*10000.0_dp, int_8)
1608 : ELSE
1609 0 : binned_dist(ibin)%cost = INT((binned_dist(ibin)%time_other_scf)*10000.0_dp, int_8)
1610 : END IF
1611 : CASE (hfx_do_eval_forces)
1612 17399 : binned_dist(ibin)%cost = INT((binned_dist(ibin)%time_forces)*10000.0_dp, int_8)
1613 : END SELECT
1614 : END IF
1615 : END DO
1616 1982 : !$OMP BARRIER
1617 1982 : !$OMP MASTER
1618 : !! store all local results in a big cost matrix
1619 5946 : ALLOCATE (cost_matrix(ncpu*nbins*n_threads))
1620 255678 : cost_matrix = 0
1621 5946 : ALLOCATE (sendbuffer(max_bin_size*n_threads))
1622 3964 : ALLOCATE (recbuffer(max_bin_size*n_threads))
1623 : !$OMP END MASTER
1624 1982 : !$OMP BARRIER
1625 1982 : my_global_start_idx = bin_histogram(my_process_id + 1, 2) - my_bin_size
1626 1982 : icpu = para_env%mepos + 1
1627 128830 : DO i = 1, my_bin_size
1628 128830 : cost_matrix(my_global_start_idx + i) = binned_dist(i)%cost
1629 : END DO
1630 :
1631 1982 : mepos = para_env%mepos
1632 1982 : !$OMP BARRIER
1633 1982 : !$OMP MASTER
1634 5946 : ALLOCATE (bins_per_rank(ncpu))
1635 1982 : bins_per_rank = 0
1636 5946 : DO icpu = 1, ncpu
1637 9910 : bins_per_rank(icpu) = SUM(bin_histogram((icpu - 1)*n_threads + 1:(icpu - 1)*n_threads + n_threads, 1))
1638 : END DO
1639 : sendbuffer(1:bins_per_rank(para_env%mepos + 1)) = &
1640 255678 : cost_matrix(my_global_start_idx + 1:my_global_start_idx + bins_per_rank(para_env%mepos + 1))
1641 :
1642 1982 : dest = MODULO(mepos + 1, ncpu)
1643 1982 : source = MODULO(mepos - 1, ncpu)
1644 : ! sync before/after ring of isendrecv
1645 1982 : CALL para_env%sync()
1646 5946 : DO icpu = 0, ncpu - 1
1647 3964 : IF (icpu /= ncpu - 1) THEN
1648 : CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
1649 1982 : req(1), req(2), 13)
1650 : END IF
1651 3964 : data_from = MODULO(mepos - icpu, ncpu)
1652 9910 : start_idx = SUM(bins_per_rank(1:data_from + 1)) - bins_per_rank(data_from + 1) + 1
1653 3964 : end_idx = start_idx + bins_per_rank(data_from + 1) - 1
1654 511356 : cost_matrix(start_idx:end_idx) = sendbuffer(1:end_idx - start_idx + 1)
1655 :
1656 3964 : IF (icpu /= ncpu - 1) THEN
1657 1982 : CALL mp_waitall(req)
1658 : END IF
1659 3964 : swapbuffer => sendbuffer
1660 3964 : sendbuffer => recbuffer
1661 5946 : recbuffer => swapbuffer
1662 : END DO
1663 1982 : DEALLOCATE (recbuffer, sendbuffer)
1664 : ! sync before/after ring of isendrecv
1665 1982 : CALL para_env%sync()
1666 : !$OMP END MASTER
1667 1982 : !$OMP BARRIER
1668 5946 : ALLOCATE (local_cost_matrix(SIZE(cost_matrix, 1)))
1669 509374 : local_cost_matrix = cost_matrix
1670 1982 : !$OMP MASTER
1671 5946 : ALLOCATE (shm_distribution_vector(ncpu*nbins*n_threads))
1672 : CALL optimize_distribution(ncpu*nbins*n_threads, ncpu*n_threads, local_cost_matrix, &
1673 1982 : shm_distribution_vector, x_data%load_balance_parameter%do_randomize)
1674 :
1675 718484 : ALLOCATE (full_dist(ncpu*n_threads, max_bin_size))
1676 :
1677 712538 : full_dist(:, :)%istart = 0_int_8
1678 712538 : full_dist(:, :)%number_of_atom_quartets = 0_int_8
1679 712538 : full_dist(:, :)%cost = 0_int_8
1680 712538 : full_dist(:, :)%time_first_scf = 0.0_dp
1681 712538 : full_dist(:, :)%time_other_scf = 0.0_dp
1682 712538 : full_dist(:, :)%time_forces = 0.0_dp
1683 : !$OMP END MASTER
1684 :
1685 1982 : !$OMP BARRIER
1686 1982 : mepos = para_env%mepos + 1
1687 255678 : full_dist((mepos - 1)*n_threads + i_thread + 1, 1:my_bin_size) = binned_dist(1:my_bin_size)
1688 1982 : !$OMP BARRIER
1689 1982 : !$OMP MASTER
1690 5946 : ALLOCATE (sendbuffer(3*max_bin_size*n_threads))
1691 3964 : ALLOCATE (recbuffer(3*max_bin_size*n_threads))
1692 1982 : mepos = para_env%mepos
1693 3964 : DO j = 1, n_threads
1694 240816 : DO i = 1, max_bin_size
1695 236852 : sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 1) = full_dist(mepos*n_threads + j, i)%istart
1696 236852 : sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 2) = full_dist(mepos*n_threads + j, i)%number_of_atom_quartets
1697 238834 : sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 3) = full_dist(mepos*n_threads + j, i)%cost
1698 : END DO
1699 : END DO
1700 1982 : dest = MODULO(mepos + 1, ncpu)
1701 1982 : source = MODULO(mepos - 1, ncpu)
1702 : ! sync before/after ring of isendrecv
1703 1982 : CALL para_env%sync()
1704 5946 : DO icpu = 0, ncpu - 1
1705 3964 : IF (icpu /= ncpu - 1) THEN
1706 : CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
1707 1982 : req(1), req(2), 13)
1708 : END IF
1709 3964 : data_from = MODULO(mepos - icpu, ncpu)
1710 7928 : DO j = 1, n_threads
1711 481632 : DO i = 1, max_bin_size
1712 473704 : full_dist(data_from*n_threads + j, i)%istart = sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 1)
1713 473704 : full_dist(data_from*n_threads + j, i)%number_of_atom_quartets = sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 2)
1714 477668 : full_dist(data_from*n_threads + j, i)%cost = sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 3)
1715 : END DO
1716 : END DO
1717 :
1718 3964 : IF (icpu /= ncpu - 1) THEN
1719 1982 : CALL mp_waitall(req)
1720 : END IF
1721 3964 : swapbuffer => sendbuffer
1722 3964 : sendbuffer => recbuffer
1723 5946 : recbuffer => swapbuffer
1724 : END DO
1725 : ! sync before/after ring of isendrecv
1726 1982 : DEALLOCATE (recbuffer, sendbuffer)
1727 1982 : CALL para_env%sync()
1728 : !$OMP END MASTER
1729 1982 : !$OMP BARRIER
1730 : !! reorder the distribution according to the distribution vector
1731 5946 : ALLOCATE (tmp_pos(ncpu*n_threads))
1732 5946 : tmp_pos = 1
1733 259642 : ALLOCATE (tmp_dist(nbins*ncpu*n_threads))
1734 :
1735 255678 : tmp_dist(:)%istart = 0_int_8
1736 255678 : tmp_dist(:)%number_of_atom_quartets = 0_int_8
1737 255678 : tmp_dist(:)%cost = 0_int_8
1738 255678 : tmp_dist(:)%time_first_scf = 0.0_dp
1739 255678 : tmp_dist(:)%time_other_scf = 0.0_dp
1740 255678 : tmp_dist(:)%time_forces = 0.0_dp
1741 :
1742 5946 : mepos = my_process_id + 1
1743 5946 : DO icpu = 1, n_processes
1744 259642 : DO i = 1, bin_histogram(icpu, 1)
1745 257660 : IF (shm_distribution_vector(bin_histogram(icpu, 2) - bin_histogram(icpu, 1) + i) == mepos) THEN
1746 126848 : tmp_dist(tmp_pos(mepos)) = full_dist(icpu, i)
1747 126848 : tmp_pos(mepos) = tmp_pos(mepos) + 1
1748 : END IF
1749 : END DO
1750 : END DO
1751 :
1752 : !! Assign the load to each process
1753 : NULLIFY (ptr_to_tmp_dist)
1754 1982 : mepos = my_process_id + 1
1755 1982 : ptr_to_tmp_dist => tmp_dist(1:tmp_pos(mepos) - 1)
1756 1106 : SELECT CASE (eval_type)
1757 : CASE (hfx_do_eval_energy)
1758 1106 : CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
1759 : CASE (hfx_do_eval_forces)
1760 1982 : CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
1761 : END SELECT
1762 :
1763 1982 : !$OMP BARRIER
1764 1982 : !$OMP MASTER
1765 1982 : DEALLOCATE (full_dist, cost_matrix, shm_distribution_vector)
1766 1982 : DEALLOCATE (bins_per_rank, bin_histogram)
1767 : !$OMP END MASTER
1768 1982 : !$OMP BARRIER
1769 1982 : DEALLOCATE (tmp_dist, tmp_pos)
1770 1982 : DEALLOCATE (binned_dist, local_cost_matrix)
1771 : END IF
1772 1982 : !$OMP BARRIER
1773 1982 : !$OMP MASTER
1774 1982 : CALL timestop(handle)
1775 : !$OMP END MASTER
1776 1982 : !$OMP BARRIER
1777 :
1778 3964 : END SUBROUTINE hfx_update_load_balance
1779 :
1780 : ! **************************************************************************************************
1781 : !> \brief estimates the cost of a set quartet with info available at load balance time
1782 : !> i.e. without much info on the primitives primitives
1783 : !> \param nsa ...
1784 : !> \param nsb ...
1785 : !> \param nsc ...
1786 : !> \param nsd ...
1787 : !> \param npgfa ...
1788 : !> \param npgfb ...
1789 : !> \param npgfc ...
1790 : !> \param npgfd ...
1791 : !> \param ratio ...
1792 : !> \param p1 ...
1793 : !> \param p2 ...
1794 : !> \param p3 ...
1795 : !> \return ...
1796 : !> \par History
1797 : !> 08.2009 created Joost VandeVondele
1798 : !> \author Joost VandeVondele
1799 : ! **************************************************************************************************
1800 5591994 : FUNCTION cost_model(nsa, nsb, nsc, nsd, npgfa, npgfb, npgfc, npgfd, ratio, p1, p2, p3) RESULT(res)
1801 : REAL(KIND=dp) :: estimate1, estimate2, estimate, ratio, switch, mu, sigma
1802 : INTEGER(KIND=int_8) :: res
1803 : REAL(KIND=dp), INTENT(IN) :: p1(12), p2(12), p3(2)
1804 :
1805 : INTEGER :: nsa, nsb, nsc, nsd, npgfa, npgfb, npgfc, npgfd
1806 :
1807 5591994 : estimate1 = estimate_basic(p1)
1808 5591994 : estimate2 = estimate_basic(p2)
1809 5591994 : mu = LOG(ABS(1.0E6_dp*p3(1)) + 1)
1810 5591994 : sigma = p3(2)*0.1_dp*mu
1811 5591994 : switch = 1.0_dp/(1.0_dp + EXP((LOG(estimate1) - mu)/sigma))
1812 5591994 : estimate = estimate1*(1.0_dp - switch) + estimate2*switch
1813 5591994 : res = INT(estimate*0.001_dp, KIND=int_8) + 1
1814 :
1815 : CONTAINS
1816 :
1817 : ! **************************************************************************************************
1818 : !> \brief ...
1819 : !> \param p ...
1820 : !> \return ...
1821 : ! **************************************************************************************************
1822 11183988 : REAL(KIND=dp) FUNCTION estimate_basic(p) RESULT(res)
1823 : REAL(KIND=dp) :: p(12)
1824 :
1825 : REAL(KIND=dp) :: p1, p10, p11, p12, p2, p3, p4, p5, p6, &
1826 : p7, p8, p9
1827 :
1828 11183988 : p1 = p(1); p2 = p(2); p3 = p(3); p4 = p(4)
1829 11183988 : p5 = p(5); p6 = p(6); p7 = p(7); p8 = p(8)
1830 11183988 : p9 = p(9); p10 = p(10); p11 = p(11); p12 = p(12)
1831 : res = poly2(nsa, p1, p2, p3)*poly2(nsb, p1, p2, p3)*poly2(nsc, p1, p2, p3)*poly2(nsd, p1, p2, p3)* &
1832 : poly2(npgfa, p4, p5, p6)*poly2(npgfb, p4, p5, p6)*poly2(npgfc, p4, p5, p6)* &
1833 : poly2(npgfd, p4, p5, p6)*EXP(-p7*ratio + p8*ratio**2) + &
1834 11183988 : 1000.0_dp*p9 + poly2(nsa, p10, p11, p12)*poly2(nsb, p10, p11, p12)*poly2(nsc, p10, p11, p12)*poly2(nsd, p10, p11, p12)
1835 11183988 : res = 1 + ABS(res)
1836 11183988 : END FUNCTION estimate_basic
1837 :
1838 : ! **************************************************************************************************
1839 : !> \brief ...
1840 : !> \param x ...
1841 : !> \param a0 ...
1842 : !> \param a1 ...
1843 : !> \param a2 ...
1844 : !> \return ...
1845 : ! **************************************************************************************************
1846 134207856 : REAL(KIND=dp) FUNCTION poly2(x, a0, a1, a2)
1847 : INTEGER, INTENT(IN) :: x
1848 : REAL(KIND=dp), INTENT(IN) :: a0, a1, a2
1849 : REAL(KIND=dp) :: r
1850 :
1851 134207856 : r = REAL(x, KIND=dp)
1852 134207856 : poly2 = a0 + (a1 + a2*r)*r
1853 134207856 : END FUNCTION poly2
1854 :
1855 : END FUNCTION cost_model
1856 : ! **************************************************************************************************
1857 : !> \brief Minimizes the maximum cost per cpu by shuffling around all bins
1858 : !> \param total_number_of_bins ...
1859 : !> \param number_of_processes ...
1860 : !> \param bin_costs costs per bin
1861 : !> \param distribution_vector will contain the final distribution
1862 : !> \param do_randomize ...
1863 : !> \par History
1864 : !> 03.2009 created from a hack by Joost [Manuel Guidon]
1865 : !> \author Manuel Guidon
1866 : ! **************************************************************************************************
1867 4108 : SUBROUTINE optimize_distribution(total_number_of_bins, number_of_processes, bin_costs, &
1868 : distribution_vector, do_randomize)
1869 : INTEGER :: total_number_of_bins, number_of_processes
1870 : INTEGER(int_8), DIMENSION(:), POINTER :: bin_costs
1871 : INTEGER, DIMENSION(:), POINTER :: distribution_vector
1872 : LOGICAL, INTENT(IN) :: do_randomize
1873 :
1874 : INTEGER :: i, itmp, j, nstep
1875 4108 : INTEGER(int_8), DIMENSION(:), POINTER :: my_cost_cpu, tmp_cost, tmp_cpu_cost
1876 : INTEGER, DIMENSION(:), POINTER :: tmp_cpu_index, tmp_index
1877 4108 : TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
1878 :
1879 4108 : nstep = MAX(1, INT(number_of_processes)/2)
1880 :
1881 12324 : ALLOCATE (tmp_cost(total_number_of_bins))
1882 12324 : ALLOCATE (tmp_index(total_number_of_bins))
1883 12324 : ALLOCATE (tmp_cpu_cost(number_of_processes))
1884 12324 : ALLOCATE (tmp_cpu_index(number_of_processes))
1885 8216 : ALLOCATE (my_cost_cpu(number_of_processes))
1886 1059864 : tmp_cost = bin_costs
1887 :
1888 4108 : CALL sort(tmp_cost, total_number_of_bins, tmp_index)
1889 12324 : my_cost_cpu = 0
1890 : !
1891 : ! assign the largest remaining bin to the CPU with the smallest load
1892 : ! gives near perfect distributions for a sufficient number of bins ...
1893 : ! doing this in chunks of nstep (where nstep ~ number_of_processes) makes this n log n and gives
1894 : ! each cpu a similar number of tasks.
1895 : ! it also avoids degenerate cases where thousands of zero sized tasks
1896 : ! are assigned to the same (least loaded) cpu
1897 : !
1898 4108 : IF (do_randomize) THEN
1899 : rng_stream = rng_stream_type(name="uniform_rng", &
1900 0 : distribution_type=UNIFORM)
1901 : END IF
1902 :
1903 529932 : DO i = total_number_of_bins, 1, -nstep
1904 3154944 : tmp_cpu_cost = my_cost_cpu
1905 525824 : CALL sort(tmp_cpu_cost, INT(number_of_processes), tmp_cpu_index)
1906 525824 : IF (do_randomize) THEN
1907 0 : CALL rng_stream%shuffle(tmp_cpu_index(1:MIN(i, nstep)))
1908 : END IF
1909 1055756 : DO j = 1, MIN(i, nstep)
1910 525824 : itmp = tmp_cpu_index(j)
1911 525824 : distribution_vector(tmp_index(i - j + 1)) = itmp
1912 1051648 : my_cost_cpu(itmp) = my_cost_cpu(itmp) + bin_costs(tmp_index(i - j + 1))
1913 : END DO
1914 : END DO
1915 :
1916 4108 : DEALLOCATE (tmp_cost, tmp_index, tmp_cpu_cost)
1917 4108 : DEALLOCATE (tmp_cpu_index, my_cost_cpu)
1918 8216 : END SUBROUTINE optimize_distribution
1919 :
1920 : ! **************************************************************************************************
1921 : !> \brief Given a 2d index pair, this function returns a 1d index pair for
1922 : !> a symmetric upper triangle NxN matrix
1923 : !> The compiler should inline this function, therefore it appears in
1924 : !> several modules
1925 : !> \param i 2d index
1926 : !> \param j 2d index
1927 : !> \param N matrix size
1928 : !> \return ...
1929 : !> \par History
1930 : !> 03.2009 created [Manuel Guidon]
1931 : !> \author Manuel Guidon
1932 : ! **************************************************************************************************
1933 34660 : PURE FUNCTION get_1D_idx(i, j, N)
1934 : INTEGER, INTENT(IN) :: i, j
1935 : INTEGER(int_8), INTENT(IN) :: N
1936 : INTEGER(int_8) :: get_1D_idx
1937 :
1938 : INTEGER(int_8) :: min_ij
1939 :
1940 34660 : min_ij = MIN(i, j)
1941 34660 : get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N
1942 :
1943 34660 : END FUNCTION get_1D_idx
1944 :
1945 : ! **************************************************************************************************
1946 : !> \brief ...
1947 : !> \param natom ...
1948 : !> \param nkind ...
1949 : !> \param list_ij ...
1950 : !> \param list_kl ...
1951 : !> \param set_list_ij ...
1952 : !> \param set_list_kl ...
1953 : !> \param iatom_start ...
1954 : !> \param iatom_end ...
1955 : !> \param jatom_start ...
1956 : !> \param jatom_end ...
1957 : !> \param katom_start ...
1958 : !> \param katom_end ...
1959 : !> \param latom_start ...
1960 : !> \param latom_end ...
1961 : !> \param particle_set ...
1962 : !> \param coeffs_set ...
1963 : !> \param coeffs_kind ...
1964 : !> \param is_assoc_atomic_block_global ...
1965 : !> \param do_periodic ...
1966 : !> \param kind_of ...
1967 : !> \param basis_parameter ...
1968 : !> \param pmax_set ...
1969 : !> \param pmax_atom ...
1970 : !> \param pmax_blocks ...
1971 : !> \param cell ...
1972 : !> \param do_p_screening ...
1973 : !> \param map_atom_to_kind_atom ...
1974 : !> \param eval_type ...
1975 : !> \param log10_eps_schwarz ...
1976 : !> \param log_2 ...
1977 : !> \param coeffs_kind_max0 ...
1978 : !> \param use_virial ...
1979 : !> \param atomic_pair_list ...
1980 : !> \return ...
1981 : ! **************************************************************************************************
1982 189279 : FUNCTION estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
1983 : iatom_start, iatom_end, jatom_start, jatom_end, &
1984 : katom_start, katom_end, latom_start, latom_end, &
1985 : particle_set, &
1986 189279 : coeffs_set, coeffs_kind, &
1987 189279 : is_assoc_atomic_block_global, do_periodic, &
1988 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
1989 : cell, &
1990 : do_p_screening, map_atom_to_kind_atom, eval_type, &
1991 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, &
1992 189279 : atomic_pair_list)
1993 :
1994 : INTEGER, INTENT(IN) :: natom, nkind
1995 : TYPE(pair_list_type) :: list_ij, list_kl
1996 : TYPE(pair_set_list_type), DIMENSION(:) :: set_list_ij, set_list_kl
1997 : INTEGER, INTENT(IN) :: iatom_start, iatom_end, jatom_start, &
1998 : jatom_end, katom_start, katom_end, &
1999 : latom_start, latom_end
2000 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2001 : TYPE(hfx_screen_coeff_type), &
2002 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
2003 : TYPE(hfx_screen_coeff_type), &
2004 : DIMENSION(nkind, nkind) :: coeffs_kind
2005 : INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block_global
2006 : LOGICAL :: do_periodic
2007 : INTEGER :: kind_of(*)
2008 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
2009 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set
2010 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom
2011 : REAL(dp) :: pmax_blocks
2012 : TYPE(cell_type), POINTER :: cell
2013 : LOGICAL, INTENT(IN) :: do_p_screening
2014 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom
2015 : INTEGER, INTENT(IN) :: eval_type
2016 : REAL(dp) :: log10_eps_schwarz, log_2, &
2017 : coeffs_kind_max0
2018 : LOGICAL, INTENT(IN) :: use_virial
2019 : LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
2020 : INTEGER(int_8) :: estimate_block_cost
2021 :
2022 : INTEGER :: i_list_ij, i_list_kl, i_set_list_ij, i_set_list_ij_start, i_set_list_ij_stop, &
2023 : i_set_list_kl, i_set_list_kl_start, i_set_list_kl_stop, iatom, ikind, iset, jatom, jkind, &
2024 : jset, katom, kind_kind_idx, kkind, kset, latom, lkind, lset, swap_id
2025 189279 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, npgfc, npgfd, nsgfa, &
2026 189279 : nsgfb, nsgfc, nsgfd
2027 : REAL(dp) :: actual_pmax_atom, cost_tmp, max_val1, &
2028 : max_val2, pmax_entry, rab2, rcd2, &
2029 : screen_kind_ij, screen_kind_kl
2030 189279 : REAL(dp), DIMENSION(:, :), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4
2031 :
2032 189279 : estimate_block_cost = 0_int_8
2033 :
2034 : CALL build_pair_list(natom, list_ij, set_list_ij, iatom_start, iatom_end, jatom_start, jatom_end, &
2035 : kind_of, basis_parameter, particle_set, &
2036 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, &
2037 189279 : log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
2038 :
2039 : CALL build_pair_list(natom, list_kl, set_list_kl, katom_start, katom_end, latom_start, latom_end, &
2040 : kind_of, basis_parameter, particle_set, &
2041 : do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, &
2042 189279 : log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
2043 :
2044 363427 : DO i_list_ij = 1, list_ij%n_element
2045 174148 : iatom = list_ij%elements(i_list_ij)%pair(1)
2046 174148 : jatom = list_ij%elements(i_list_ij)%pair(2)
2047 174148 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
2048 174148 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
2049 174148 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
2050 174148 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
2051 174148 : rab2 = list_ij%elements(i_list_ij)%dist2
2052 :
2053 174148 : nsgfa => basis_parameter(ikind)%nsgf
2054 174148 : nsgfb => basis_parameter(jkind)%nsgf
2055 174148 : npgfa => basis_parameter(ikind)%npgf
2056 174148 : npgfb => basis_parameter(jkind)%npgf
2057 :
2058 528416 : DO i_list_kl = 1, list_kl%n_element
2059 :
2060 164989 : katom = list_kl%elements(i_list_kl)%pair(1)
2061 164989 : latom = list_kl%elements(i_list_kl)%pair(2)
2062 :
2063 164989 : IF (.NOT. (katom + latom <= iatom + jatom)) CYCLE
2064 91295 : IF (((iatom + jatom) == (katom + latom)) .AND. (katom < iatom)) CYCLE
2065 :
2066 87629 : IF (eval_type == hfx_do_eval_forces) THEN
2067 8303 : IF (.NOT. use_virial) THEN
2068 7543 : IF ((iatom == jatom .AND. iatom == katom .AND. iatom == latom)) CYCLE
2069 : END IF
2070 : END IF
2071 :
2072 86836 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
2073 86836 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
2074 86836 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
2075 86836 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
2076 86836 : rcd2 = list_kl%elements(i_list_kl)%dist2
2077 :
2078 86836 : nsgfc => basis_parameter(kkind)%nsgf
2079 86836 : nsgfd => basis_parameter(lkind)%nsgf
2080 86836 : npgfc => basis_parameter(kkind)%npgf
2081 86836 : npgfd => basis_parameter(lkind)%npgf
2082 :
2083 86836 : IF (do_p_screening) THEN
2084 : actual_pmax_atom = MAX(pmax_atom(katom, iatom), &
2085 : pmax_atom(latom, jatom), &
2086 : pmax_atom(latom, iatom), &
2087 8845 : pmax_atom(katom, jatom))
2088 : ELSE
2089 : actual_pmax_atom = 0.0_dp
2090 : END IF
2091 :
2092 : screen_kind_ij = coeffs_kind(jkind, ikind)%x(1)*rab2 + &
2093 86836 : coeffs_kind(jkind, ikind)%x(2)
2094 : screen_kind_kl = coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
2095 86836 : coeffs_kind(lkind, kkind)%x(2)
2096 86836 : IF (screen_kind_ij + screen_kind_kl + actual_pmax_atom < log10_eps_schwarz) CYCLE
2097 :
2098 83247 : IF (.NOT. (is_assoc_atomic_block_global(latom, iatom) >= 1 .AND. &
2099 : is_assoc_atomic_block_global(katom, iatom) >= 1 .AND. &
2100 : is_assoc_atomic_block_global(katom, jatom) >= 1 .AND. &
2101 : is_assoc_atomic_block_global(latom, jatom) >= 1)) CYCLE
2102 :
2103 81655 : IF (do_p_screening) THEN
2104 3131 : SELECT CASE (eval_type)
2105 : CASE (hfx_do_eval_energy)
2106 3131 : swap_id = 0
2107 3131 : kind_kind_idx = INT(get_1D_idx(kkind, ikind, INT(nkind, int_8)))
2108 3131 : IF (ikind >= kkind) THEN
2109 : ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2110 : map_atom_to_kind_atom(katom), &
2111 2824 : map_atom_to_kind_atom(iatom))
2112 : ELSE
2113 : ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2114 : map_atom_to_kind_atom(iatom), &
2115 307 : map_atom_to_kind_atom(katom))
2116 307 : swap_id = swap_id + 1
2117 : END IF
2118 3131 : kind_kind_idx = INT(get_1D_idx(lkind, jkind, INT(nkind, int_8)))
2119 3131 : IF (jkind >= lkind) THEN
2120 : ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2121 : map_atom_to_kind_atom(latom), &
2122 2949 : map_atom_to_kind_atom(jatom))
2123 : ELSE
2124 : ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2125 : map_atom_to_kind_atom(jatom), &
2126 182 : map_atom_to_kind_atom(latom))
2127 182 : swap_id = swap_id + 2
2128 : END IF
2129 3131 : kind_kind_idx = INT(get_1D_idx(lkind, ikind, INT(nkind, int_8)))
2130 3131 : IF (ikind >= lkind) THEN
2131 : ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2132 : map_atom_to_kind_atom(latom), &
2133 1754 : map_atom_to_kind_atom(iatom))
2134 : ELSE
2135 : ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2136 : map_atom_to_kind_atom(iatom), &
2137 1377 : map_atom_to_kind_atom(latom))
2138 1377 : swap_id = swap_id + 4
2139 : END IF
2140 3131 : kind_kind_idx = INT(get_1D_idx(kkind, jkind, INT(nkind, int_8)))
2141 3131 : IF (jkind >= kkind) THEN
2142 : ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2143 : map_atom_to_kind_atom(katom), &
2144 3131 : map_atom_to_kind_atom(jatom))
2145 : ELSE
2146 : ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2147 : map_atom_to_kind_atom(jatom), &
2148 0 : map_atom_to_kind_atom(katom))
2149 0 : swap_id = swap_id + 8
2150 : END IF
2151 : CASE (hfx_do_eval_forces)
2152 5534 : swap_id = 16
2153 5534 : kind_kind_idx = INT(get_1D_idx(kkind, ikind, INT(nkind, int_8)))
2154 5534 : IF (ikind >= kkind) THEN
2155 : ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2156 : map_atom_to_kind_atom(katom), &
2157 5233 : map_atom_to_kind_atom(iatom))
2158 : ELSE
2159 : ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2160 : map_atom_to_kind_atom(iatom), &
2161 301 : map_atom_to_kind_atom(katom))
2162 301 : swap_id = swap_id + 1
2163 : END IF
2164 5534 : kind_kind_idx = INT(get_1D_idx(lkind, jkind, INT(nkind, int_8)))
2165 5534 : IF (jkind >= lkind) THEN
2166 : ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2167 : map_atom_to_kind_atom(latom), &
2168 5534 : map_atom_to_kind_atom(jatom))
2169 : ELSE
2170 : ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2171 : map_atom_to_kind_atom(jatom), &
2172 0 : map_atom_to_kind_atom(latom))
2173 0 : swap_id = swap_id + 2
2174 : END IF
2175 5534 : kind_kind_idx = INT(get_1D_idx(lkind, ikind, INT(nkind, int_8)))
2176 5534 : IF (ikind >= lkind) THEN
2177 : ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2178 : map_atom_to_kind_atom(latom), &
2179 4517 : map_atom_to_kind_atom(iatom))
2180 : ELSE
2181 : ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2182 : map_atom_to_kind_atom(iatom), &
2183 1017 : map_atom_to_kind_atom(latom))
2184 1017 : swap_id = swap_id + 4
2185 : END IF
2186 5534 : kind_kind_idx = INT(get_1D_idx(kkind, jkind, INT(nkind, int_8)))
2187 14199 : IF (jkind >= kkind) THEN
2188 : ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2189 : map_atom_to_kind_atom(katom), &
2190 5534 : map_atom_to_kind_atom(jatom))
2191 : ELSE
2192 : ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2193 : map_atom_to_kind_atom(jatom), &
2194 0 : map_atom_to_kind_atom(katom))
2195 0 : swap_id = swap_id + 8
2196 : END IF
2197 : END SELECT
2198 : END IF
2199 :
2200 652451 : DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
2201 396648 : iset = set_list_ij(i_set_list_ij)%pair(1)
2202 396648 : jset = set_list_ij(i_set_list_ij)%pair(2)
2203 :
2204 : max_val1 = coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + &
2205 396648 : coeffs_set(jset, iset, jkind, ikind)%x(2)
2206 :
2207 396648 : IF (max_val1 + screen_kind_kl + actual_pmax_atom < log10_eps_schwarz) CYCLE
2208 5734196 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
2209 5180743 : kset = set_list_kl(i_set_list_kl)%pair(1)
2210 5180743 : lset = set_list_kl(i_set_list_kl)%pair(2)
2211 :
2212 : max_val2 = max_val1 + (coeffs_set(lset, kset, lkind, kkind)%x(1)*rcd2 + &
2213 5180743 : coeffs_set(lset, kset, lkind, kkind)%x(2))
2214 :
2215 5180743 : IF (max_val2 + actual_pmax_atom < log10_eps_schwarz) CYCLE
2216 4678349 : IF (do_p_screening) THEN
2217 : CALL get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, &
2218 : iset, jset, kset, lset, &
2219 555612 : pmax_entry, swap_id)
2220 555612 : IF (eval_type == hfx_do_eval_forces) THEN
2221 349692 : pmax_entry = log_2 + pmax_entry
2222 : END IF
2223 : ELSE
2224 4122737 : pmax_entry = 0.0_dp
2225 : END IF
2226 4678349 : max_val2 = max_val2 + pmax_entry
2227 4678349 : IF (max_val2 < log10_eps_schwarz) CYCLE
2228 396648 : SELECT CASE (eval_type)
2229 : CASE (hfx_do_eval_energy)
2230 : cost_tmp = cost_model(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
2231 : npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
2232 : max_val2/log10_eps_schwarz, &
2233 4037248 : p1_energy, p2_energy, p3_energy)
2234 4037248 : estimate_block_cost = estimate_block_cost + INT(cost_tmp, KIND=int_8)
2235 : CASE (hfx_do_eval_forces)
2236 : cost_tmp = cost_model(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
2237 : npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
2238 : max_val2/log10_eps_schwarz, &
2239 542929 : p1_forces, p2_forces, p3_forces)
2240 5123106 : estimate_block_cost = estimate_block_cost + INT(cost_tmp, KIND=int_8)
2241 : END SELECT
2242 : END DO ! i_set_list_kl
2243 : END DO ! i_set_list_ij
2244 : END DO ! i_list_kl
2245 : END DO ! i_list_ij
2246 :
2247 189279 : END FUNCTION estimate_block_cost
2248 :
2249 : ! **************************************************************************************************
2250 : !> \brief ...
2251 : !> \param nkind ...
2252 : !> \param para_env ...
2253 : !> \param natom ...
2254 : !> \param block_size ...
2255 : !> \param nblock ...
2256 : !> \param blocks ...
2257 : !> \param list_ij ...
2258 : !> \param list_kl ...
2259 : !> \param set_list_ij ...
2260 : !> \param set_list_kl ...
2261 : !> \param particle_set ...
2262 : !> \param coeffs_set ...
2263 : !> \param coeffs_kind ...
2264 : !> \param is_assoc_atomic_block_global ...
2265 : !> \param do_periodic ...
2266 : !> \param kind_of ...
2267 : !> \param basis_parameter ...
2268 : !> \param pmax_set ...
2269 : !> \param pmax_atom ...
2270 : !> \param pmax_blocks ...
2271 : !> \param cell ...
2272 : !> \param do_p_screening ...
2273 : !> \param map_atom_to_kind_atom ...
2274 : !> \param eval_type ...
2275 : !> \param log10_eps_schwarz ...
2276 : !> \param log_2 ...
2277 : !> \param coeffs_kind_max0 ...
2278 : !> \param use_virial ...
2279 : !> \param atomic_pair_list ...
2280 : ! **************************************************************************************************
2281 1490 : SUBROUTINE init_blocks(nkind, para_env, natom, block_size, nblock, blocks, &
2282 1490 : list_ij, list_kl, set_list_ij, set_list_kl, &
2283 : particle_set, &
2284 : coeffs_set, coeffs_kind, &
2285 1490 : is_assoc_atomic_block_global, do_periodic, &
2286 : kind_of, basis_parameter, pmax_set, pmax_atom, &
2287 : pmax_blocks, cell, &
2288 : do_p_screening, map_atom_to_kind_atom, eval_type, &
2289 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, &
2290 1490 : atomic_pair_list)
2291 :
2292 : INTEGER, INTENT(IN) :: nkind
2293 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2294 : INTEGER :: natom, block_size, nblock
2295 : TYPE(hfx_block_range_type), DIMENSION(1:nblock) :: blocks
2296 : TYPE(pair_list_type) :: list_ij, list_kl
2297 : TYPE(pair_set_list_type), DIMENSION(:) :: set_list_ij, set_list_kl
2298 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2299 : TYPE(hfx_screen_coeff_type), &
2300 : DIMENSION(:, :, :, :), POINTER :: coeffs_set
2301 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
2302 : POINTER :: coeffs_kind
2303 : INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block_global
2304 : LOGICAL :: do_periodic
2305 : INTEGER :: kind_of(*)
2306 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
2307 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set
2308 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom
2309 : REAL(dp) :: pmax_blocks
2310 : TYPE(cell_type), POINTER :: cell
2311 : LOGICAL, INTENT(IN) :: do_p_screening
2312 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom
2313 : INTEGER, INTENT(IN) :: eval_type
2314 : REAL(dp) :: log10_eps_schwarz, log_2, &
2315 : coeffs_kind_max0
2316 : LOGICAL, INTENT(IN) :: use_virial
2317 : LOGICAL, DIMENSION(natom, natom) :: atomic_pair_list
2318 :
2319 : INTEGER :: atom_block, i, iatom_block, iatom_end, &
2320 : iatom_start, my_cpu_rank, ncpus
2321 :
2322 6206 : DO atom_block = 0, nblock - 1
2323 4716 : iatom_block = MODULO(atom_block, nblock) + 1
2324 4716 : iatom_start = (iatom_block - 1)*block_size + 1
2325 4716 : iatom_end = MIN(iatom_block*block_size, natom)
2326 4716 : blocks(atom_block + 1)%istart = iatom_start
2327 4716 : blocks(atom_block + 1)%iend = iatom_end
2328 6206 : blocks(atom_block + 1)%cost = 0_int_8
2329 : END DO
2330 :
2331 1490 : ncpus = para_env%num_pe
2332 1490 : my_cpu_rank = para_env%mepos
2333 6206 : DO i = 1, nblock
2334 4716 : IF (MODULO(i, ncpus) /= my_cpu_rank) THEN
2335 2355 : blocks(i)%istart = 0
2336 2355 : blocks(i)%iend = 0
2337 2355 : CYCLE
2338 : END IF
2339 2361 : iatom_start = blocks(i)%istart
2340 2361 : iatom_end = blocks(i)%iend
2341 : blocks(i)%cost = estimate_block_cost(natom, nkind, list_ij, list_kl, set_list_ij, set_list_kl, &
2342 : iatom_start, iatom_end, iatom_start, iatom_end, &
2343 : iatom_start, iatom_end, iatom_start, iatom_end, &
2344 : particle_set, &
2345 : coeffs_set, coeffs_kind, &
2346 : is_assoc_atomic_block_global, do_periodic, &
2347 : kind_of, basis_parameter, pmax_set, pmax_atom, pmax_blocks, &
2348 : cell, &
2349 : do_p_screening, map_atom_to_kind_atom, eval_type, &
2350 3851 : log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
2351 :
2352 : END DO
2353 1490 : END SUBROUTINE init_blocks
2354 :
2355 : ! **************************************************************************************************
2356 : !> \brief ...
2357 : !> \param para_env ...
2358 : !> \param x_data ...
2359 : !> \param iw ...
2360 : !> \param n_threads ...
2361 : !> \param i_thread ...
2362 : !> \param eval_type ...
2363 : ! **************************************************************************************************
2364 0 : SUBROUTINE collect_load_balance_info(para_env, x_data, iw, n_threads, i_thread, &
2365 : eval_type)
2366 :
2367 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2368 : TYPE(hfx_type), POINTER :: x_data
2369 : INTEGER, INTENT(IN) :: iw, n_threads, i_thread, eval_type
2370 :
2371 : INTEGER :: i, j, k, my_rank, nbins, nranks, &
2372 : total_bins
2373 : INTEGER(int_8) :: avg_bin, avg_rank, max_bin, max_rank, &
2374 : min_bin, min_rank, sum_bin, sum_rank
2375 0 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: buffer, buffer_in, buffer_out, summary
2376 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:), SAVE :: shm_cost_vector
2377 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bins_per_rank, rdispl, sort_idx
2378 : INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: shm_bins_per_rank, shm_displ
2379 :
2380 0 : SELECT CASE (eval_type)
2381 : CASE (hfx_do_eval_energy)
2382 0 : nbins = SIZE(x_data%distribution_energy)
2383 : CASE (hfx_do_eval_forces)
2384 0 : nbins = SIZE(x_data%distribution_forces)
2385 : END SELECT
2386 :
2387 0 : !$OMP MASTER
2388 0 : ALLOCATE (shm_bins_per_rank(n_threads))
2389 0 : ALLOCATE (shm_displ(n_threads + 1))
2390 : !$OMP END MASTER
2391 0 : !$OMP BARRIER
2392 :
2393 0 : shm_bins_per_rank(i_thread + 1) = nbins
2394 0 : !$OMP BARRIER
2395 0 : nbins = 0
2396 0 : DO i = 1, n_threads
2397 0 : nbins = nbins + shm_bins_per_rank(i)
2398 : END DO
2399 0 : my_rank = para_env%mepos
2400 0 : nranks = para_env%num_pe
2401 :
2402 0 : !$OMP BARRIER
2403 0 : !$OMP MASTER
2404 0 : ALLOCATE (bins_per_rank(nranks))
2405 0 : bins_per_rank = 0
2406 :
2407 0 : bins_per_rank(my_rank + 1) = nbins
2408 :
2409 0 : CALL para_env%sum(bins_per_rank)
2410 :
2411 0 : total_bins = 0
2412 0 : DO i = 1, nranks
2413 0 : total_bins = total_bins + bins_per_rank(i)
2414 : END DO
2415 :
2416 0 : ALLOCATE (shm_cost_vector(2*total_bins))
2417 0 : shm_cost_vector = -1_int_8
2418 0 : shm_displ(1) = 1
2419 0 : DO i = 2, n_threads
2420 0 : shm_displ(i) = shm_displ(i - 1) + shm_bins_per_rank(i - 1)
2421 : END DO
2422 0 : shm_displ(n_threads + 1) = nbins + 1
2423 : !$OMP END MASTER
2424 0 : !$OMP BARRIER
2425 0 : j = 0
2426 0 : SELECT CASE (eval_type)
2427 : CASE (hfx_do_eval_energy)
2428 0 : DO i = shm_displ(i_thread + 1), shm_displ(i_thread + 2) - 1
2429 0 : j = j + 1
2430 0 : shm_cost_vector(2*(i - 1) + 1) = x_data%distribution_energy(j)%cost
2431 0 : shm_cost_vector(2*i) = INT(x_data%distribution_energy(j)%time_first_scf*10000.0_dp, KIND=int_8)
2432 : END DO
2433 : CASE (hfx_do_eval_forces)
2434 0 : DO i = shm_displ(i_thread + 1), shm_displ(i_thread + 2) - 1
2435 0 : j = j + 1
2436 0 : shm_cost_vector(2*(i - 1) + 1) = x_data%distribution_forces(j)%cost
2437 0 : shm_cost_vector(2*i) = INT(x_data%distribution_forces(j)%time_forces*10000.0_dp, KIND=int_8)
2438 : END DO
2439 : END SELECT
2440 0 : !$OMP BARRIER
2441 0 : !$OMP MASTER
2442 : ! ** calculate offsets
2443 0 : ALLOCATE (rdispl(nranks))
2444 0 : bins_per_rank(:) = bins_per_rank(:)*2
2445 0 : rdispl(1) = 0
2446 0 : DO i = 2, nranks
2447 0 : rdispl(i) = rdispl(i - 1) + bins_per_rank(i - 1)
2448 : END DO
2449 :
2450 0 : ALLOCATE (buffer_in(2*nbins))
2451 0 : ALLOCATE (buffer_out(2*total_bins))
2452 :
2453 0 : DO i = 1, nbins
2454 0 : buffer_in(2*(i - 1) + 1) = shm_cost_vector(2*(i - 1) + 1)
2455 0 : buffer_in(2*i) = shm_cost_vector(2*i)
2456 : END DO
2457 :
2458 0 : CALL para_env%gatherv(buffer_in, buffer_out, bins_per_rank, rdispl)
2459 :
2460 0 : IF (iw > 0) THEN
2461 :
2462 0 : ALLOCATE (summary(2*nranks))
2463 0 : summary = 0_int_8
2464 :
2465 0 : WRITE (iw, '( /, 1X, 79("-") )')
2466 0 : WRITE (iw, '( " -", 77X, "-" )')
2467 0 : SELECT CASE (eval_type)
2468 : CASE (hfx_do_eval_energy)
2469 0 : WRITE (iw, '( " -", 20X, A, 19X, "-" )') ' HFX LOAD BALANCE INFORMATION - ENERGY '
2470 : CASE (hfx_do_eval_forces)
2471 0 : WRITE (iw, '( " -", 20X, A, 19X, "-" )') ' HFX LOAD BALANCE INFORMATION - FORCES '
2472 : END SELECT
2473 0 : WRITE (iw, '( " -", 77X, "-" )')
2474 0 : WRITE (iw, '( 1X, 79("-") )')
2475 :
2476 0 : WRITE (iw, FMT="(T3,A,T15,A,T35,A,T55,A)") "MPI RANK", "BIN #", "EST cost", "Processing time [s]"
2477 0 : WRITE (iw, '( 1X, 79("-"), / )')
2478 0 : k = 0
2479 0 : DO i = 1, nranks
2480 0 : DO j = 1, bins_per_rank(i)/2
2481 0 : k = k + 1
2482 : WRITE (iw, FMT="(T6,I5,T15,I5,T27,I16,T55,F19.8)") &
2483 0 : i - 1, j, buffer_out(2*(k - 1) + 1), REAL(buffer_out(2*k), dp)/10000.0_dp
2484 0 : summary(2*(i - 1) + 1) = summary(2*(i - 1) + 1) + buffer_out(2*(k - 1) + 1)
2485 0 : summary(2*i) = summary(2*i) + buffer_out(2*k)
2486 : END DO
2487 : END DO
2488 :
2489 : !** Summary
2490 : max_bin = 0_int_8
2491 : min_bin = HUGE(min_bin)
2492 : sum_bin = 0_int_8
2493 0 : DO i = 1, total_bins
2494 0 : sum_bin = sum_bin + buffer_out(2*i)
2495 0 : max_bin = MAX(max_bin, buffer_out(2*i))
2496 0 : min_bin = MIN(min_bin, buffer_out(2*i))
2497 : END DO
2498 0 : avg_bin = sum_bin/total_bins
2499 :
2500 0 : max_rank = 0_int_8
2501 0 : min_rank = HUGE(min_rank)
2502 0 : sum_rank = 0_int_8
2503 0 : DO i = 1, nranks
2504 0 : sum_rank = sum_rank + summary(2*i)
2505 0 : max_rank = MAX(max_rank, summary(2*i))
2506 0 : min_rank = MIN(min_rank, summary(2*i))
2507 : END DO
2508 0 : avg_rank = sum_rank/nranks
2509 :
2510 0 : WRITE (iw, FMT='(/,T3,A,/)') "SUMMARY:"
2511 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Max bin", REAL(max_bin, dp)/10000.0_dp
2512 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Min bin", REAL(min_bin, dp)/10000.0_dp
2513 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Sum bin", REAL(sum_bin, dp)/10000.0_dp
2514 0 : WRITE (iw, FMT="(T3,A,T35,F19.8,/)") "Avg bin", REAL(avg_bin, dp)/10000.0_dp
2515 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Max rank", REAL(max_rank, dp)/10000.0_dp
2516 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Min rank", REAL(min_rank, dp)/10000.0_dp
2517 0 : WRITE (iw, FMT="(T3,A,T35,F19.8)") "Sum rank", REAL(sum_rank, dp)/10000.0_dp
2518 0 : WRITE (iw, FMT="(T3,A,T35,F19.8,/)") "Avg rank", REAL(avg_rank, dp)/10000.0_dp
2519 :
2520 0 : ALLOCATE (buffer(nranks))
2521 0 : ALLOCATE (sort_idx(nranks))
2522 :
2523 0 : DO i = 1, nranks
2524 0 : buffer(i) = summary(2*i)
2525 : END DO
2526 :
2527 0 : CALL sort(buffer, nranks, sort_idx)
2528 :
2529 0 : WRITE (iw, FMT="(T3,A,T35,A,T55,A,/)") "MPI RANK", "EST cost", "Processing time [s]"
2530 0 : DO i = nranks, 1, -1
2531 0 : WRITE (iw, FMT="(T6,I5,T27,I16,T55,F19.8)") sort_idx(i) - 1, summary(2*(sort_idx(i) - 1) + 1), REAL(buffer(i), dp)/10000.0_dp
2532 : END DO
2533 :
2534 0 : DEALLOCATE (summary, buffer, sort_idx)
2535 :
2536 : END IF
2537 :
2538 0 : DEALLOCATE (buffer_in, buffer_out, rdispl)
2539 :
2540 0 : CALL para_env%sync()
2541 :
2542 0 : DEALLOCATE (shm_bins_per_rank, shm_displ, shm_cost_vector)
2543 : !$OMP END MASTER
2544 0 : !$OMP BARRIER
2545 :
2546 0 : END SUBROUTINE collect_load_balance_info
2547 :
2548 : #:include 'hfx_get_pmax_val.fypp'
2549 :
2550 : END MODULE hfx_load_balance_methods
|
