(git:3f0d10b)
Loading...
Searching...
No Matches
hfx_load_balance_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
74CONTAINS
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 SUBROUTINE hfx_load_balance(x_data, eps_schwarz, particle_set, max_set, para_env, &
124 coeffs_set, coeffs_kind, &
125 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 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 INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: buffer_in, buffer_out, cached_block_cost
172 INTEGER(int_8), DIMENSION(:), POINTER :: local_cost_matrix, recbuffer, &
173 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, tmp_index, tmp_pos, &
178 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 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 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 DIMENSION(:) :: set_list_ij, set_list_kl
195
196!$OMP BARRIER
197!$OMP MASTER
198 CALL timeset(routinen, handle)
199!$OMP END MASTER
200!$OMP BARRIER
201
202 log10_eps_schwarz = log10(eps_schwarz)
203 log_2 = log10(2.0_dp)
204 coeffs_kind_max0 = maxval(coeffs_kind(:, :)%x(2))
205 ncpu = para_env%num_pe
206 n_processes = ncpu*n_threads
207 natom = SIZE(particle_set)
208
209 block_size = load_balance_parameter%block_size
210 ALLOCATE (set_list_ij((max_set*block_size)**2))
211 ALLOCATE (set_list_kl((max_set*block_size)**2))
212
213 IF (.NOT. load_balance_parameter%blocks_initialized) THEN
214!$OMP BARRIER
215!$OMP MASTER
216 CALL timeset(routinen//"_range", handle_range)
217
218 nblocks = max((natom + block_size - 1)/block_size, 1)
219 ALLOCATE (blocks_guess(nblocks))
220 ALLOCATE (tmp_blocks(natom))
221 ALLOCATE (tmp_blocks2(natom))
222
223 pmax_blocks = 0.0_dp
224 SELECT CASE (eval_type)
225 CASE (hfx_do_eval_energy)
226 atomic_pair_list => x_data%atomic_pair_list
227 CASE (hfx_do_eval_forces)
228 atomic_pair_list => x_data%atomic_pair_list_forces
229 END SELECT
230 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
240
241 total_block_self_cost = 0
242
243 DO i = 1, nblocks
244 total_block_self_cost = total_block_self_cost + blocks_guess(i)%cost
245 END DO
246
247 CALL para_env%sum(total_block_self_cost)
248
249 objective_block_size = load_balance_parameter%block_size
250 objective_nblocks = max((natom + objective_block_size - 1)/objective_block_size, 1)
251
252 self_cost_per_block = (total_block_self_cost + objective_nblocks - 1)/(objective_nblocks)
253
254 DO i = 1, nblocks
255 tmp_blocks2(i) = blocks_guess(i)
256 END DO
257
258 optimized = .false.
259 i = 0
260 DO WHILE (.NOT. optimized)
261 i = i + 1
262 current_block_id = 0
263 changed = .false.
264 DO atom_block = 1, nblocks
265 current_block_id = current_block_id + 1
266 iatom_start = tmp_blocks2(atom_block)%istart
267 iatom_end = tmp_blocks2(atom_block)%iend
268 IF (tmp_blocks2(atom_block)%cost > 1.5_dp*self_cost_per_block .AND. iatom_end - iatom_start > 0) THEN
269 changed = .true.
270 new_iatom_start = iatom_start
271 new_iatom_end = (iatom_end - iatom_start + 1)/2 + iatom_start - 1
272 new_jatom_start = new_iatom_end + 1
273 new_jatom_end = iatom_end
274 tmp_blocks(current_block_id)%istart = new_iatom_start
275 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
287 current_block_id = current_block_id + 1
288 tmp_blocks(current_block_id)%istart = new_jatom_start
289 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
301 ELSE
302 tmp_blocks(current_block_id)%istart = iatom_start
303 tmp_blocks(current_block_id)%iend = iatom_end
304 tmp_blocks(current_block_id)%cost = tmp_blocks2(atom_block)%cost
305 END IF
306 END DO
307 IF (.NOT. changed) optimized = .true.
308 IF (i > 20) optimized = .true.
309 nblocks = current_block_id
310 DO atom_block = 1, nblocks
311 tmp_blocks2(atom_block) = tmp_blocks(atom_block)
312 END DO
313 END DO
314
315 DEALLOCATE (tmp_blocks2)
316
317 ! ** count number of non empty blocks on each node
318 non_empty_blocks = 0
319 DO atom_block = 1, nblocks
320 IF (tmp_blocks(atom_block)%istart == 0) cycle
321 non_empty_blocks = non_empty_blocks + 1
322 END DO
323
324 ALLOCATE (rcount(ncpu))
325 rcount = 0
326 rcount(para_env%mepos + 1) = non_empty_blocks
327 CALL para_env%sum(rcount)
328
329 ! ** sum all non_empty_blocks
330 total_blocks = 0
331 DO i = 1, ncpu
332 total_blocks = total_blocks + rcount(i)
333 END DO
334
335 ! ** calculate offsets
336 ALLOCATE (rdispl(ncpu))
337 rcount(:) = rcount(:)*3
338 rdispl(1) = 0
339 DO i = 2, ncpu
340 rdispl(i) = rdispl(i - 1) + rcount(i - 1)
341 END DO
342
343 ALLOCATE (buffer_in(3*non_empty_blocks))
344
345 non_empty_blocks = 0
346 DO atom_block = 1, nblocks
347 IF (tmp_blocks(atom_block)%istart == 0) cycle
348 buffer_in(non_empty_blocks*3 + 1) = tmp_blocks(atom_block)%istart
349 buffer_in(non_empty_blocks*3 + 2) = tmp_blocks(atom_block)%iend
350 buffer_in(non_empty_blocks*3 + 3) = tmp_blocks(atom_block)%cost
351 non_empty_blocks = non_empty_blocks + 1
352 END DO
353
354 nblocks = total_blocks
355
356 ALLOCATE (tmp_blocks2(nblocks))
357
358 ALLOCATE (buffer_out(3*nblocks))
359
360 ! ** Gather all three arrays
361 CALL para_env%allgatherv(buffer_in, buffer_out, rcount, rdispl)
362
363 DO i = 1, nblocks
364 tmp_blocks2(i)%istart = int(buffer_out((i - 1)*3 + 1))
365 tmp_blocks2(i)%iend = int(buffer_out((i - 1)*3 + 2))
366 tmp_blocks2(i)%cost = buffer_out((i - 1)*3 + 3)
367 END DO
368
369 ! ** Now we sort the blocks
370 ALLOCATE (to_be_sorted(nblocks))
371 ALLOCATE (tmp_index(nblocks))
372
373 DO atom_block = 1, nblocks
374 to_be_sorted(atom_block) = tmp_blocks2(atom_block)%istart
375 END DO
376
377 CALL sort(to_be_sorted, nblocks, tmp_index)
378
379 ALLOCATE (x_data%blocks(nblocks))
380
381 DO atom_block = 1, nblocks
382 x_data%blocks(atom_block) = tmp_blocks2(tmp_index(atom_block))
383 END DO
384
385 shm_blocks => x_data%blocks
386 shm_nblocks = nblocks
387
388 ! ** Set nblocks in structure
389 load_balance_parameter%nblocks = nblocks
390
391 DEALLOCATE (blocks_guess, tmp_blocks, tmp_blocks2)
392
393 DEALLOCATE (rcount, rdispl, buffer_in, buffer_out, to_be_sorted, tmp_index)
394
395 load_balance_parameter%blocks_initialized = .true.
396
397 x_data%blocks = shm_blocks
398 load_balance_parameter%nblocks = shm_nblocks
399 load_balance_parameter%blocks_initialized = .true.
400
401 ALLOCATE (x_data%pmax_block(shm_nblocks, shm_nblocks))
402 x_data%pmax_block = 0.0_dp
403 pmax_block => x_data%pmax_block
404 CALL timestop(handle_range)
405!$OMP END MASTER
406!$OMP BARRIER
407
408 IF (.NOT. load_balance_parameter%blocks_initialized) THEN
409 ALLOCATE (x_data%blocks(shm_nblocks))
410 x_data%blocks = shm_blocks
411 load_balance_parameter%nblocks = shm_nblocks
412 load_balance_parameter%blocks_initialized = .true.
413 END IF
414 !! ** precalculate maximum density matrix elements in blocks
415!$OMP BARRIER
416 END IF
417
418!$OMP BARRIER
419!$OMP MASTER
420 pmax_block => x_data%pmax_block
421 pmax_block = 0.0_dp
422 IF (do_p_screening) THEN
423 DO iatom_block = 1, shm_nblocks
424 iatom_start = x_data%blocks(iatom_block)%istart
425 iatom_end = x_data%blocks(iatom_block)%iend
426 DO jatom_block = 1, shm_nblocks
427 jatom_start = x_data%blocks(jatom_block)%istart
428 jatom_end = x_data%blocks(jatom_block)%iend
429 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 SELECT CASE (eval_type)
435 CASE (hfx_do_eval_energy)
436 atomic_pair_list => x_data%atomic_pair_list
437 CASE (hfx_do_eval_forces)
438 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 x_data%blocks)
443
444!$OMP END MASTER
445!$OMP BARRIER
446
447 !! If there is only 1 cpu skip the binning
448 IF (n_processes == 1) THEN
449 ALLOCATE (tmp_dist(1))
450 tmp_dist(1)%number_of_atom_quartets = huge(tmp_dist(1)%number_of_atom_quartets)
451 tmp_dist(1)%istart = 0_int_8
452 ptr_to_tmp_dist => tmp_dist(:)
453 SELECT CASE (eval_type)
454 CASE (hfx_do_eval_energy)
455 CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
456 CASE (hfx_do_eval_forces)
457 CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
458 END SELECT
459 DEALLOCATE (tmp_dist)
460 ELSE
461 !! Calculate total numbers of integrals that have to be calculated (wrt screening and symmetry)
462!$OMP BARRIER
463!$OMP MASTER
464 CALL timeset(routinen//"_count", handle_inner)
465!$OMP END MASTER
466!$OMP BARRIER
467
468 cost_per_core = 0_int_8
469 my_process_id = para_env%mepos*n_threads + i_thread
470 nblocks = load_balance_parameter%nblocks
471 n_my_quartets = 0_int_8
472 DO atom_block = my_process_id, int(nblocks, kind=int_8)**4 - 1_int_8, n_processes
473 latom_block = int(modulo(atom_block, int(nblocks, kind=int_8))) + 1
474 tmp_block = atom_block/nblocks
475 katom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
476 IF (latom_block < katom_block) cycle
477 tmp_block = tmp_block/nblocks
478 jatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
479 tmp_block = tmp_block/nblocks
480 iatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
481 IF (jatom_block < iatom_block) cycle
482 n_my_quartets = n_my_quartets + 1_int_8
483 IF (n_my_quartets == max_cached_quartet_costs) EXIT
484 END DO
485 IF (n_my_quartets > 0_int_8) THEN
486 ALLOCATE (cached_block_cost(n_my_quartets))
487 END IF
488 q_idx = 0_int_8
489
490 DO atom_block = my_process_id, int(nblocks, kind=int_8)**4 - 1, n_processes
491
492 latom_block = int(modulo(atom_block, int(nblocks, kind=int_8))) + 1
493 tmp_block = atom_block/nblocks
494 katom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
495 IF (latom_block < katom_block) cycle
496 tmp_block = tmp_block/nblocks
497 jatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
498 tmp_block = tmp_block/nblocks
499 iatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
500 IF (jatom_block < iatom_block) cycle
501
502 q_idx = q_idx + 1_int_8
503
504 iatom_start = x_data%blocks(iatom_block)%istart
505 iatom_end = x_data%blocks(iatom_block)%iend
506 jatom_start = x_data%blocks(jatom_block)%istart
507 jatom_end = x_data%blocks(jatom_block)%iend
508 katom_start = x_data%blocks(katom_block)%istart
509 katom_end = x_data%blocks(katom_block)%iend
510 latom_start = x_data%blocks(latom_block)%istart
511 latom_end = x_data%blocks(latom_block)%iend
512
513 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 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 pmax_block(katom_block, jatom_block))
524 END SELECT
525
526 IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) THEN
527 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
541 IF (q_idx <= n_my_quartets) cached_block_cost(q_idx) = block_cost_val
542 cost_per_core = cost_per_core + block_cost_val
543
544 END DO ! atom_block
545
546 nbins = load_balance_parameter%nbins
547 cost_per_bin = (cost_per_core + nbins - 1)/(nbins)
548
549!$OMP BARRIER
550!$OMP MASTER
551 CALL timestop(handle_inner)
552!$OMP END MASTER
553!$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!$OMP BARRIER
569!$OMP MASTER
570 CALL timeset(routinen//"_bin", handle_inner)
571!$OMP END MASTER
572!$OMP BARRIER
573
574 ALLOCATE (binned_dist(nbins))
575 binned_dist(:)%istart = -1_int_8
576 binned_dist(:)%number_of_atom_quartets = 0_int_8
577 binned_dist(:)%cost = 0_int_8
578 binned_dist(:)%time_first_scf = 0.0_dp
579 binned_dist(:)%time_other_scf = 0.0_dp
580 binned_dist(:)%time_forces = 0.0_dp
581
582 current_cost = 0
583 mepos = 1
584 distribution_counter_start = 1
585 distribution_counter_end = 0
586 ibin = 1
587
588 global_quartet_counter = 0
589 local_quartet_counter = 0
590 last_bin_needs_to_be_filled = .false.
591 DO atom_block = my_process_id, int(nblocks, kind=int_8)**4 - 1, n_processes
592 latom_block = int(modulo(atom_block, int(nblocks, kind=int_8))) + 1
593 tmp_block = atom_block/nblocks
594 katom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
595 IF (latom_block < katom_block) cycle
596 tmp_block = tmp_block/nblocks
597 jatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
598 tmp_block = tmp_block/nblocks
599 iatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
600 IF (jatom_block < iatom_block) cycle
601
602 distribution_counter_end = distribution_counter_end + 1
603 global_quartet_counter = global_quartet_counter + 1
604 last_bin_needs_to_be_filled = .true.
605
606 IF (binned_dist(ibin)%istart == -1_int_8) binned_dist(ibin)%istart = atom_block
607
608 IF (distribution_counter_end <= n_my_quartets) THEN
609 block_cost_val = cached_block_cost(distribution_counter_end)
610 IF (block_cost_val < 0_int_8) cycle
611 current_cost = current_cost + block_cost_val
612 ELSE
613 iatom_start = x_data%blocks(iatom_block)%istart
614 iatom_end = x_data%blocks(iatom_block)%iend
615 jatom_start = x_data%blocks(jatom_block)%istart
616 jatom_end = x_data%blocks(jatom_block)%iend
617 katom_start = x_data%blocks(katom_block)%istart
618 katom_end = x_data%blocks(katom_block)%iend
619 latom_start = x_data%blocks(latom_block)%istart
620 latom_end = x_data%blocks(latom_block)%iend
621
622 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 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 pmax_block(katom_block, jatom_block))
633 END SELECT
634
635 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
648 END IF
649
650 IF (current_cost >= cost_per_bin) THEN
651 IF (ibin == nbins) THEN
652 binned_dist(ibin)%number_of_atom_quartets = binned_dist(ibin)%number_of_atom_quartets + &
653 distribution_counter_end - distribution_counter_start + 1
654 ELSE
655 binned_dist(ibin)%number_of_atom_quartets = distribution_counter_end - distribution_counter_start + 1
656 END IF
657 binned_dist(ibin)%cost = binned_dist(ibin)%cost + current_cost
658 ibin = min(ibin + 1, nbins)
659 distribution_counter_start = distribution_counter_end + 1
660 current_cost = 0
661 last_bin_needs_to_be_filled = .false.
662 END IF
663 END DO
664 IF (ALLOCATED(cached_block_cost)) DEALLOCATE (cached_block_cost)
665
666!$OMP BARRIER
667!$OMP MASTER
668 CALL timestop(handle_inner)
669 CALL timeset(routinen//"_dist", handle_inner)
670!$OMP END MASTER
671!$OMP BARRIER
672 !! Fill the last bin if necessary
673 IF (last_bin_needs_to_be_filled) THEN
674 binned_dist(ibin)%cost = binned_dist(ibin)%cost + current_cost
675 IF (ibin == nbins) THEN
676 binned_dist(ibin)%number_of_atom_quartets = binned_dist(ibin)%number_of_atom_quartets + &
677 distribution_counter_end - distribution_counter_start + 1
678 ELSE
679 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 DO ibin = 1, nbins
685 local_quartet_counter = local_quartet_counter + binned_dist(ibin)%number_of_atom_quartets
686 END DO
687!$OMP BARRIER
688!$OMP MASTER
689 shm_local_quartet_counter = 0
690 shm_global_quartet_counter = 0
691!$OMP END MASTER
692!$OMP BARRIER
693!$OMP ATOMIC
694 shm_local_quartet_counter = shm_local_quartet_counter + local_quartet_counter
695!$OMP ATOMIC
696 shm_global_quartet_counter = shm_global_quartet_counter + global_quartet_counter
697
698!$OMP BARRIER
699!$OMP MASTER
700 CALL para_env%sum(shm_local_quartet_counter)
701 CALL para_env%sum(shm_global_quartet_counter)
702 IF (para_env%is_source()) THEN
703 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 "Number of local quartets (", shm_local_quartet_counter, &
706 ") and number of global quartets (", shm_global_quartet_counter, &
707 ") are different. Please send in a bug report."
708 cpabort(error_msg)
709 END IF
710 END IF
711!$OMP END MASTER
712
713!$OMP BARRIER
714!$OMP MASTER
715 ALLOCATE (cost_matrix(ncpu*nbins*n_threads))
716 cost_matrix = 0
717!$OMP END MASTER
718!$OMP BARRIER
719 icpu = para_env%mepos + 1
720 DO i = 1, nbins
721 cost_matrix((icpu - 1)*nbins*n_threads + i_thread*nbins + i) = binned_dist(i)%cost
722 END DO
723 mepos = para_env%mepos
724!$OMP BARRIER
725
726!$OMP MASTER
727 ! sync before/after ring of isendrecv
728 CALL para_env%sync()
729
730 ALLOCATE (sendbuffer(nbins*n_threads))
731 ALLOCATE (recbuffer(nbins*n_threads))
732
733 sendbuffer = cost_matrix(mepos*nbins*n_threads + 1:mepos*nbins*n_threads + nbins*n_threads)
734
735 dest = modulo(mepos + 1, ncpu)
736 source = modulo(mepos - 1, ncpu)
737 DO icpu = 0, ncpu - 1
738 IF (icpu /= ncpu - 1) THEN
739 CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
740 req(1), req(2), 13)
741 END IF
742 data_from = modulo(mepos - icpu, ncpu)
743 cost_matrix(data_from*nbins*n_threads + 1:data_from*nbins*n_threads + nbins*n_threads) = sendbuffer
744 IF (icpu /= ncpu - 1) THEN
745 CALL mp_waitall(req)
746 END IF
747 swapbuffer => sendbuffer
748 sendbuffer => recbuffer
749 recbuffer => swapbuffer
750 END DO
751 DEALLOCATE (recbuffer, sendbuffer)
752!$OMP END MASTER
753!$OMP BARRIER
754
755!$OMP BARRIER
756!$OMP MASTER
757 CALL timestop(handle_inner)
758 CALL timeset(routinen//"_opt", handle_inner)
759!$OMP END MASTER
760!$OMP BARRIER
761
762 !! Find an optimal distribution i.e. assign each element of the cost matrix to a certain process
763!$OMP BARRIER
764 ALLOCATE (local_cost_matrix(SIZE(cost_matrix, 1)))
765 local_cost_matrix = cost_matrix
766!$OMP MASTER
767 ALLOCATE (shm_distribution_vector(ncpu*nbins*n_threads))
768
769 CALL optimize_distribution(ncpu*nbins*n_threads, ncpu*n_threads, local_cost_matrix, &
770 shm_distribution_vector, x_data%load_balance_parameter%do_randomize)
771
772 CALL timestop(handle_inner)
773 CALL timeset(routinen//"_redist", handle_inner)
774 !! Collect local data to global array
775 ALLOCATE (full_dist(ncpu*n_threads, nbins))
776
777 full_dist(:, :)%istart = 0_int_8
778 full_dist(:, :)%number_of_atom_quartets = 0_int_8
779 full_dist(:, :)%cost = 0_int_8
780 full_dist(:, :)%time_first_scf = 0.0_dp
781 full_dist(:, :)%time_other_scf = 0.0_dp
782 full_dist(:, :)%time_forces = 0.0_dp
783!$OMP END MASTER
784!$OMP BARRIER
785 mepos = para_env%mepos + 1
786 full_dist((mepos - 1)*n_threads + i_thread + 1, :) = binned_dist(:)
787
788!$OMP BARRIER
789!$OMP MASTER
790 ALLOCATE (sendbuffer(3*nbins*n_threads))
791 ALLOCATE (recbuffer(3*nbins*n_threads))
792 mepos = para_env%mepos
793 DO j = 1, n_threads
794 DO i = 1, nbins
795 sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 1) = full_dist(mepos*n_threads + j, i)%istart
796 sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 2) = full_dist(mepos*n_threads + j, i)%number_of_atom_quartets
797 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 CALL para_env%sync()
803 dest = modulo(mepos + 1, ncpu)
804 source = modulo(mepos - 1, ncpu)
805 DO icpu = 0, ncpu - 1
806 IF (icpu /= ncpu - 1) THEN
807 CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
808 req(1), req(2), 13)
809 END IF
810 data_from = modulo(mepos - icpu, ncpu)
811 DO j = 1, n_threads
812 DO i = 1, nbins
813 full_dist(data_from*n_threads + j, i)%istart = sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 1)
814 full_dist(data_from*n_threads + j, i)%number_of_atom_quartets = sendbuffer((j - 1)*3*nbins + (i - 1)*3 + 2)
815 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 IF (icpu /= ncpu - 1) THEN
820 CALL mp_waitall(req)
821 END IF
822 swapbuffer => sendbuffer
823 sendbuffer => recbuffer
824 recbuffer => swapbuffer
825 END DO
826 DEALLOCATE (recbuffer, sendbuffer)
827
828 ! sync before/after ring of isendrecv
829 CALL para_env%sync()
830!$OMP END MASTER
831!$OMP BARRIER
832 !! reorder the distribution according to the distribution vector
833 ALLOCATE (tmp_pos(ncpu*n_threads))
834 tmp_pos = 1
835 ALLOCATE (tmp_dist(nbins*ncpu*n_threads))
836
837 tmp_dist(:)%istart = 0_int_8
838 tmp_dist(:)%number_of_atom_quartets = 0_int_8
839 tmp_dist(:)%cost = 0_int_8
840 tmp_dist(:)%time_first_scf = 0.0_dp
841 tmp_dist(:)%time_other_scf = 0.0_dp
842 tmp_dist(:)%time_forces = 0.0_dp
843
844 DO icpu = 1, n_processes
845 DO i = 1, nbins
846 mepos = my_process_id + 1
847 IF (shm_distribution_vector((icpu - 1)*nbins + i) == mepos) THEN
848 tmp_dist(tmp_pos(mepos)) = full_dist(icpu, i)
849 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 mepos = my_process_id + 1
857 ptr_to_tmp_dist => tmp_dist(1:tmp_pos(mepos) - 1)
858 SELECT CASE (eval_type)
859 CASE (hfx_do_eval_energy)
860 CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
861 CASE (hfx_do_eval_forces)
862 CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
863 END SELECT
864
865!$OMP BARRIER
866!$OMP MASTER
867 DEALLOCATE (full_dist, cost_matrix, shm_distribution_vector)
868!$OMP END MASTER
869!$OMP BARRIER
870 DEALLOCATE (tmp_dist, tmp_pos)
871 DEALLOCATE (binned_dist, local_cost_matrix)
872 DEALLOCATE (set_list_ij, set_list_kl)
873
874!$OMP BARRIER
875!$OMP MASTER
876 CALL timestop(handle_inner)
877!$OMP END MASTER
878!$OMP BARRIER
879 END IF
880!$OMP BARRIER
881!$OMP MASTER
882 CALL timestop(handle)
883!$OMP END MASTER
884!$OMP BARRIER
885 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 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 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 INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: local_cost_vector
980 INTEGER, ALLOCATABLE, DIMENSION(:) :: blocksize, p_atom_blocks, permute, &
981 q_atom_blocks
982 REAL(dp) :: maximum, mean
983
984! internal variables:
985
986!$OMP BARRIER
987!$OMP MASTER
988 CALL timeset(routinen, handle)
989!$OMP END MASTER
990!$OMP BARRIER
991
992 ! calculate best p/q distribution grid for the n_processes
993 CALL hfx_calculate_pq(p, q, numbins, n_processes)
994
995 ALLOCATE (blocksize(numbins))
996 ALLOCATE (permute(nblocks**2))
997 DO i = 1, nblocks**2
998 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 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 np = numbins/p
1019 ! number of blocks per processor in q-direction (horizontal)
1020 nq = numbins/q
1021
1022 ! calc own position in P-Q-processor grid (PQ-grid is column-major)
1023 pidx = modulo(int(my_process_id), int(p)) + 1
1024 qidx = my_process_id/p + 1
1025
1026 sizep = sum(blocksize((np*(pidx - 1) + 1):(np*pidx)))
1027 sizeq = sum(blocksize((nq*(qidx - 1) + 1):(nq*qidx)))
1028
1029 sump = sum(blocksize(1:(np*(pidx - 1))))
1030 sumq = sum(blocksize(1:(nq*(qidx - 1))))
1031
1032 ALLOCATE (p_atom_blocks(sizep))
1033 ALLOCATE (q_atom_blocks(sizeq))
1034
1035 p_atom_blocks(:) = permute((sump + 1):(sump + sizep))
1036 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 DO i = 1, sizep
1047 DO j = 1, sizeq
1048
1049 ! get corresponding 4D block indices out of our own P-Q-block
1050 latom_block = modulo(q_atom_blocks(j), nblocks)
1051 iatom_block = q_atom_blocks(j)/nblocks + 1
1052 jatom_block = modulo(p_atom_blocks(i), nblocks)
1053 katom_block = p_atom_blocks(i)/nblocks + 1
1054
1055 ! symmetry checks.
1056 IF (latom_block < katom_block) cycle
1057 IF (jatom_block < iatom_block) cycle
1058
1059 iatom_start = x_data%blocks(iatom_block)%istart
1060 iatom_end = x_data%blocks(iatom_block)%iend
1061 jatom_start = x_data%blocks(jatom_block)%istart
1062 jatom_end = x_data%blocks(jatom_block)%iend
1063 katom_start = x_data%blocks(katom_block)%istart
1064 katom_end = x_data%blocks(katom_block)%iend
1065 latom_start = x_data%blocks(latom_block)%istart
1066 latom_end = x_data%blocks(latom_block)%iend
1067
1068 ! whatever.
1069 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 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 pmax_block(katom_block, jatom_block))
1080 END SELECT
1081
1082 ! screening.
1083 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1097 END DO
1098 END DO
1099
1100 ALLOCATE (local_cost_vector(n_processes))
1101 local_cost_vector = 0
1102 local_cost_vector(my_process_id + 1) = local_cost
1103 CALL para_env%sum(local_cost_vector)
1104
1105 mean = sum(local_cost_vector)/n_processes
1106 maximum = maxval(local_cost_vector)
1107
1108!$OMP BARRIER
1109!$OMP MASTER
1110 ! only output once
1111 IF (my_process_id == 0) THEN
1112 CALL open_file(unit_number=unit_nr, file_name="loads.dat")
1113 WRITE (unit_nr, *) 'maximum cost:', maximum
1114 WRITE (unit_nr, *) 'mean cost:', mean
1115 WRITE (unit_nr, *) 'load balance ratio max/mean: ', maximum/mean
1116 WRITE (unit_nr, *) '-------- detailed per-process costs ---------'
1117 DO i = 1, n_processes
1118 WRITE (unit_nr, *) local_cost_vector(i)
1119 END DO
1120 CALL close_file(unit_nr)
1121 END IF
1122!$OMP END MASTER
1123!$OMP BARRIER
1124
1125 DEALLOCATE (local_cost_vector)
1126 DEALLOCATE (p_atom_blocks, q_atom_blocks)
1127 DEALLOCATE (blocksize, permute)
1128
1129!$OMP BARRIER
1130!$OMP MASTER
1131 CALL timestop(handle)
1132!$OMP END MASTER
1133!$OMP BARRIER
1134
1135 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 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 k = 2
1158 sqn = sqrt(real(n, kind=dp))
1159 p = 1
1160
1161 DO WHILE (real(k, kind=dp) <= sqn)
1162 IF (modulo(n, k) == 0) THEN
1163 p = k
1164 END IF
1165 k = k + 1
1166 END DO
1167 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 a = p
1173 b = q
1174
1175 DO WHILE (b /= 0)
1176 IF (a > b) THEN
1177 a = a - b
1178 ELSE
1179 b = b - a
1180 END IF
1181 END DO
1182 ! gcd(p,q) is now saved in a
1183
1184 nbins = p*q/a
1185
1186 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 RECURSIVE SUBROUTINE hfx_recursive_permute(blocksize, blockstart, blockend, nProc_in, &
1238 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: ithblocksize, localblocksize
1288 INTEGER, DIMENSION(blockend - blockstart + 1) :: bin_perm, tmp_perm
1289 REAL(dp) :: partialcost
1290 REAL(dp), DIMENSION(nblocks*nblocks) :: cost_vector
1291
1292 nproc = nproc_in
1293 cost_vector = 0.0_dp
1294
1295! loop over local atom_blocks.
1296 DO atom_block = my_process_id, int(nblocks, kind=int_8)**4 - 1, n_processes
1297
1298! get corresponding 4D block indices
1299 latom_block = int(modulo(atom_block, int(nblocks, kind=int_8))) + 1
1300 tmp_block = atom_block/nblocks
1301 katom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
1302 IF (latom_block < katom_block) cycle
1303 tmp_block = tmp_block/nblocks
1304 jatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
1305 tmp_block = tmp_block/nblocks
1306 iatom_block = int(modulo(tmp_block, int(nblocks, kind=int_8))) + 1
1307 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 row = (katom_block - 1)*nblocks + jatom_block
1314 inv_perm = 1
1315 DO WHILE (permute(inv_perm) /= row)
1316 inv_perm = inv_perm + 1
1317 END DO
1318 row = inv_perm
1319
1320 col = (iatom_block - 1)*nblocks + latom_block
1321 inv_perm = 1
1322 DO WHILE (permute(inv_perm) /= col)
1323 inv_perm = inv_perm + 1
1324 END DO
1325 col = inv_perm
1326
1327! if row/col outside our current diagonal block, skip calculation.
1328 IF (col < blockstart .OR. col > blockend) cycle
1329 IF (row < blockstart .OR. row > blockend) cycle
1330
1331 iatom_start = x_data%blocks(iatom_block)%istart
1332 iatom_end = x_data%blocks(iatom_block)%iend
1333 jatom_start = x_data%blocks(jatom_block)%istart
1334 jatom_end = x_data%blocks(jatom_block)%iend
1335 katom_start = x_data%blocks(katom_block)%istart
1336 katom_end = x_data%blocks(katom_block)%iend
1337 latom_start = x_data%blocks(latom_block)%istart
1338 latom_end = x_data%blocks(latom_block)%iend
1339
1340! whatever.
1341 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 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 pmax_block(katom_block, jatom_block))
1352 END SELECT
1353
1354! screening.
1355 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 IF (modulo(step, 2) == 0) THEN
1360 idx = row
1361 ELSE
1362 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1377
1378 cost_vector(idx) = cost_vector(idx) + partialcost
1379 END DO ! atom_block
1380
1381! sum costvector over all processes
1382 CALL para_env%sum(cost_vector)
1383
1384! calculate next prime factor of nProc
1385 nbins = 2
1386 DO WHILE (modulo(int(nproc), int(nbins)) /= 0)
1387 nbins = nbins + 1
1388 END DO
1389
1390 nproc = nproc/nbins
1391
1392! ... do the binning...
1393
1394 ALLOCATE (localblocksize(nbins))
1395 CALL hfx_permute_binning(nbins, cost_vector(blockstart:blockend), blockend - blockstart + 1, bin_perm, localblocksize)
1396
1397!... and update the permutation vector
1398
1399 tmp_perm = permute(blockstart:blockend)
1400 permute(blockstart:blockend) = tmp_perm(bin_perm)
1401
1402! split recursion into the nBins Bins
1403 IF (nproc > 1) THEN
1404 ALLOCATE (ithblocksize(nproc))
1405 DO i = 1, nbins
1406 startoffset = sum(localblocksize(1:(i - 1)))
1407 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
1420 blocksize(((i - 1)*nproc + 1):(i*nproc)) = ithblocksize
1421 END DO
1422 DEALLOCATE (ithblocksize)
1423 ELSE
1424 DO i = 1, nbins
1425 blocksize(i) = localblocksize(i)
1426 END DO
1427 END IF
1428
1429 DEALLOCATE (localblocksize)
1430
1431 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 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 INTEGER, DIMENSION(nBins, maxbinsize) :: bin
1454 INTEGER, DIMENSION(nBins) :: bin_idx
1455 INTEGER, DIMENSION(maxbinsize) :: idx
1456 REAL(dp), DIMENSION(maxbinsize) :: vec
1457 REAL(dp), DIMENSION(nBins) :: bincosts
1458
1459! be careful not to change costvector (copy it!)
1460
1461 vec = costvector
1462 block_count = 0
1463 bincosts = 0
1464
1465 !sort the array (ascending)
1466 CALL sort(vec, maxbinsize, idx)
1467
1468 ! count the loop down to distribute the largest cols/rows first
1469 DO i = maxbinsize, 1, -1
1470 IF (vec(i) == 0) THEN
1471 ! spread zero-cost col/rows evenly among procs
1472 mod_idx = modulo(i, nbins) + 1 !(note the fortran offset by one!)
1473 block_count(mod_idx) = block_count(mod_idx) + 1
1474 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 CALL sort(bincosts, nbins, bin_idx)
1479 block_count = block_count(bin_idx)
1480 bin = bin(bin_idx, :)
1481
1482 bincosts(1) = bincosts(1) + vec(i)
1483 block_count(1) = block_count(1) + 1
1484 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 DO i = 1, nbins
1491 DO j = 1, block_count(i)
1492 perm(offset + j) = bin(i, j)
1493 END DO
1494 offset = offset + block_count(i)
1495 END DO
1496
1497 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 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 TYPE(mp_request_type), DIMENSION(2) :: req
1529 INTEGER(int_8), DIMENSION(:), POINTER :: local_cost_matrix, recbuffer, &
1530 sendbuffer, swapbuffer
1531 INTEGER(int_8), DIMENSION(:), POINTER, SAVE :: cost_matrix
1532 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 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!$OMP BARRIER
1542!$OMP MASTER
1543 CALL timeset(routinen, handle)
1544!$OMP END MASTER
1545!$OMP BARRIER
1546
1547 ncpu = para_env%num_pe
1548 n_processes = ncpu*n_threads
1549 !! If there is only 1 cpu skip the binning
1550 IF (n_processes == 1) THEN
1551 ALLOCATE (tmp_dist(1))
1552 tmp_dist(1)%number_of_atom_quartets = huge(tmp_dist(1)%number_of_atom_quartets)
1553 tmp_dist(1)%istart = 0_int_8
1554 ptr_to_tmp_dist => tmp_dist(:)
1555 SELECT CASE (eval_type)
1556 CASE (hfx_do_eval_energy)
1557 CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
1558 CASE (hfx_do_eval_forces)
1559 CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
1560 END SELECT
1561 DEALLOCATE (tmp_dist)
1562 ELSE
1563 mepos = para_env%mepos
1564 my_process_id = para_env%mepos*n_threads + i_thread
1565 nbins = load_balance_parameter%nbins
1566!$OMP MASTER
1567 ALLOCATE (bin_histogram(n_processes, 2))
1568 bin_histogram = 0
1569!$OMP END MASTER
1570!$OMP BARRIER
1571 SELECT CASE (eval_type)
1572 CASE (hfx_do_eval_energy)
1573 my_bin_size = SIZE(x_data%distribution_energy)
1574 CASE (hfx_do_eval_forces)
1575 my_bin_size = SIZE(x_data%distribution_forces)
1576 END SELECT
1577 bin_histogram(my_process_id + 1, 1) = my_bin_size
1578!$OMP BARRIER
1579!$OMP MASTER
1580 CALL para_env%sum(bin_histogram(:, 1))
1581 bin_histogram(1, 2) = bin_histogram(1, 1)
1582 DO iprocess = 2, n_processes
1583 bin_histogram(iprocess, 2) = bin_histogram(iprocess - 1, 2) + bin_histogram(iprocess, 1)
1584 END DO
1585
1586 max_bin_size = maxval(bin_histogram(para_env%mepos*n_threads + 1:para_env%mepos*n_threads + n_threads, 1))
1587 CALL para_env%max(max_bin_size)
1588!$OMP END MASTER
1589!$OMP BARRIER
1590 ALLOCATE (binned_dist(my_bin_size))
1591 !! Use old binned_dist, but with timings cost
1592 SELECT CASE (eval_type)
1593 CASE (hfx_do_eval_energy)
1594 binned_dist = x_data%distribution_energy
1595 CASE (hfx_do_eval_forces)
1596 binned_dist = x_data%distribution_forces
1597 END SELECT
1598
1599 DO ibin = 1, my_bin_size
1600 IF (binned_dist(ibin)%number_of_atom_quartets == 0) THEN
1601 binned_dist(ibin)%cost = 0
1602 ELSE
1603 SELECT CASE (eval_type)
1604 CASE (hfx_do_eval_energy)
1605 IF (.NOT. load_balance_parameter%rtp_redistribute) THEN
1606 binned_dist(ibin)%cost = int((binned_dist(ibin)%time_first_scf + &
1607 binned_dist(ibin)%time_other_scf)*10000.0_dp, int_8)
1608 ELSE
1609 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 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!$OMP BARRIER
1617!$OMP MASTER
1618 !! store all local results in a big cost matrix
1619 ALLOCATE (cost_matrix(ncpu*nbins*n_threads))
1620 cost_matrix = 0
1621 ALLOCATE (sendbuffer(max_bin_size*n_threads))
1622 ALLOCATE (recbuffer(max_bin_size*n_threads))
1623!$OMP END MASTER
1624!$OMP BARRIER
1625 my_global_start_idx = bin_histogram(my_process_id + 1, 2) - my_bin_size
1626 icpu = para_env%mepos + 1
1627 DO i = 1, my_bin_size
1628 cost_matrix(my_global_start_idx + i) = binned_dist(i)%cost
1629 END DO
1630
1631 mepos = para_env%mepos
1632!$OMP BARRIER
1633!$OMP MASTER
1634 ALLOCATE (bins_per_rank(ncpu))
1635 bins_per_rank = 0
1636 DO icpu = 1, ncpu
1637 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 cost_matrix(my_global_start_idx + 1:my_global_start_idx + bins_per_rank(para_env%mepos + 1))
1641
1642 dest = modulo(mepos + 1, ncpu)
1643 source = modulo(mepos - 1, ncpu)
1644 ! sync before/after ring of isendrecv
1645 CALL para_env%sync()
1646 DO icpu = 0, ncpu - 1
1647 IF (icpu /= ncpu - 1) THEN
1648 CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
1649 req(1), req(2), 13)
1650 END IF
1651 data_from = modulo(mepos - icpu, ncpu)
1652 start_idx = sum(bins_per_rank(1:data_from + 1)) - bins_per_rank(data_from + 1) + 1
1653 end_idx = start_idx + bins_per_rank(data_from + 1) - 1
1654 cost_matrix(start_idx:end_idx) = sendbuffer(1:end_idx - start_idx + 1)
1655
1656 IF (icpu /= ncpu - 1) THEN
1657 CALL mp_waitall(req)
1658 END IF
1659 swapbuffer => sendbuffer
1660 sendbuffer => recbuffer
1661 recbuffer => swapbuffer
1662 END DO
1663 DEALLOCATE (recbuffer, sendbuffer)
1664 ! sync before/after ring of isendrecv
1665 CALL para_env%sync()
1666!$OMP END MASTER
1667!$OMP BARRIER
1668 ALLOCATE (local_cost_matrix(SIZE(cost_matrix, 1)))
1669 local_cost_matrix = cost_matrix
1670!$OMP MASTER
1671 ALLOCATE (shm_distribution_vector(ncpu*nbins*n_threads))
1672 CALL optimize_distribution(ncpu*nbins*n_threads, ncpu*n_threads, local_cost_matrix, &
1673 shm_distribution_vector, x_data%load_balance_parameter%do_randomize)
1674
1675 ALLOCATE (full_dist(ncpu*n_threads, max_bin_size))
1676
1677 full_dist(:, :)%istart = 0_int_8
1678 full_dist(:, :)%number_of_atom_quartets = 0_int_8
1679 full_dist(:, :)%cost = 0_int_8
1680 full_dist(:, :)%time_first_scf = 0.0_dp
1681 full_dist(:, :)%time_other_scf = 0.0_dp
1682 full_dist(:, :)%time_forces = 0.0_dp
1683!$OMP END MASTER
1684
1685!$OMP BARRIER
1686 mepos = para_env%mepos + 1
1687 full_dist((mepos - 1)*n_threads + i_thread + 1, 1:my_bin_size) = binned_dist(1:my_bin_size)
1688!$OMP BARRIER
1689!$OMP MASTER
1690 ALLOCATE (sendbuffer(3*max_bin_size*n_threads))
1691 ALLOCATE (recbuffer(3*max_bin_size*n_threads))
1692 mepos = para_env%mepos
1693 DO j = 1, n_threads
1694 DO i = 1, max_bin_size
1695 sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 1) = full_dist(mepos*n_threads + j, i)%istart
1696 sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 2) = full_dist(mepos*n_threads + j, i)%number_of_atom_quartets
1697 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 dest = modulo(mepos + 1, ncpu)
1701 source = modulo(mepos - 1, ncpu)
1702 ! sync before/after ring of isendrecv
1703 CALL para_env%sync()
1704 DO icpu = 0, ncpu - 1
1705 IF (icpu /= ncpu - 1) THEN
1706 CALL para_env%isendrecv(sendbuffer, dest, recbuffer, source, &
1707 req(1), req(2), 13)
1708 END IF
1709 data_from = modulo(mepos - icpu, ncpu)
1710 DO j = 1, n_threads
1711 DO i = 1, max_bin_size
1712 full_dist(data_from*n_threads + j, i)%istart = sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 1)
1713 full_dist(data_from*n_threads + j, i)%number_of_atom_quartets = sendbuffer((j - 1)*3*max_bin_size + (i - 1)*3 + 2)
1714 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 IF (icpu /= ncpu - 1) THEN
1719 CALL mp_waitall(req)
1720 END IF
1721 swapbuffer => sendbuffer
1722 sendbuffer => recbuffer
1723 recbuffer => swapbuffer
1724 END DO
1725 ! sync before/after ring of isendrecv
1726 DEALLOCATE (recbuffer, sendbuffer)
1727 CALL para_env%sync()
1728!$OMP END MASTER
1729!$OMP BARRIER
1730 !! reorder the distribution according to the distribution vector
1731 ALLOCATE (tmp_pos(ncpu*n_threads))
1732 tmp_pos = 1
1733 ALLOCATE (tmp_dist(nbins*ncpu*n_threads))
1734
1735 tmp_dist(:)%istart = 0_int_8
1736 tmp_dist(:)%number_of_atom_quartets = 0_int_8
1737 tmp_dist(:)%cost = 0_int_8
1738 tmp_dist(:)%time_first_scf = 0.0_dp
1739 tmp_dist(:)%time_other_scf = 0.0_dp
1740 tmp_dist(:)%time_forces = 0.0_dp
1741
1742 mepos = my_process_id + 1
1743 DO icpu = 1, n_processes
1744 DO i = 1, bin_histogram(icpu, 1)
1745 IF (shm_distribution_vector(bin_histogram(icpu, 2) - bin_histogram(icpu, 1) + i) == mepos) THEN
1746 tmp_dist(tmp_pos(mepos)) = full_dist(icpu, i)
1747 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 mepos = my_process_id + 1
1755 ptr_to_tmp_dist => tmp_dist(1:tmp_pos(mepos) - 1)
1756 SELECT CASE (eval_type)
1757 CASE (hfx_do_eval_energy)
1758 CALL hfx_set_distr_energy(ptr_to_tmp_dist, x_data)
1759 CASE (hfx_do_eval_forces)
1760 CALL hfx_set_distr_forces(ptr_to_tmp_dist, x_data)
1761 END SELECT
1762
1763!$OMP BARRIER
1764!$OMP MASTER
1765 DEALLOCATE (full_dist, cost_matrix, shm_distribution_vector)
1766 DEALLOCATE (bins_per_rank, bin_histogram)
1767!$OMP END MASTER
1768!$OMP BARRIER
1769 DEALLOCATE (tmp_dist, tmp_pos)
1770 DEALLOCATE (binned_dist, local_cost_matrix)
1771 END IF
1772!$OMP BARRIER
1773!$OMP MASTER
1774 CALL timestop(handle)
1775!$OMP END MASTER
1776!$OMP BARRIER
1777
1778 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 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 estimate1 = estimate_basic(p1)
1808 estimate2 = estimate_basic(p2)
1809 mu = log(abs(1.0e6_dp*p3(1)) + 1)
1810 sigma = p3(2)*0.1_dp*mu
1811 switch = 1.0_dp/(1.0_dp + exp((log(estimate1) - mu)/sigma))
1812 estimate = estimate1*(1.0_dp - switch) + estimate2*switch
1813 res = int(estimate*0.001_dp, kind=int_8) + 1
1814
1815 CONTAINS
1816
1817! **************************************************************************************************
1818!> \brief ...
1819!> \param p ...
1820!> \return ...
1821! **************************************************************************************************
1822 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 p1 = p(1); p2 = p(2); p3 = p(3); p4 = p(4)
1829 p5 = p(5); p6 = p(6); p7 = p(7); p8 = p(8)
1830 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 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 res = 1 + abs(res)
1836 END FUNCTION estimate_basic
1837
1838! **************************************************************************************************
1839!> \brief ...
1840!> \param x ...
1841!> \param a0 ...
1842!> \param a1 ...
1843!> \param a2 ...
1844!> \return ...
1845! **************************************************************************************************
1846 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 r = real(x, kind=dp)
1852 poly2 = a0 + (a1 + a2*r)*r
1853 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 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 INTEGER(int_8), DIMENSION(:), POINTER :: my_cost_cpu, tmp_cost, tmp_cpu_cost
1876 INTEGER, DIMENSION(:), POINTER :: tmp_cpu_index, tmp_index
1877 TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
1878
1879 nstep = max(1, int(number_of_processes)/2)
1880
1881 ALLOCATE (tmp_cost(total_number_of_bins))
1882 ALLOCATE (tmp_index(total_number_of_bins))
1883 ALLOCATE (tmp_cpu_cost(number_of_processes))
1884 ALLOCATE (tmp_cpu_index(number_of_processes))
1885 ALLOCATE (my_cost_cpu(number_of_processes))
1886 tmp_cost = bin_costs
1887
1888 CALL sort(tmp_cost, total_number_of_bins, tmp_index)
1889 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 IF (do_randomize) THEN
1899 rng_stream = rng_stream_type(name="uniform_rng", &
1900 distribution_type=uniform)
1901 END IF
1902
1903 DO i = total_number_of_bins, 1, -nstep
1904 tmp_cpu_cost = my_cost_cpu
1905 CALL sort(tmp_cpu_cost, int(number_of_processes), tmp_cpu_index)
1906 IF (do_randomize) THEN
1907 CALL rng_stream%shuffle(tmp_cpu_index(1:min(i, nstep)))
1908 END IF
1909 DO j = 1, min(i, nstep)
1910 itmp = tmp_cpu_index(j)
1911 distribution_vector(tmp_index(i - j + 1)) = itmp
1912 my_cost_cpu(itmp) = my_cost_cpu(itmp) + bin_costs(tmp_index(i - j + 1))
1913 END DO
1914 END DO
1915
1916 DEALLOCATE (tmp_cost, tmp_index, tmp_cpu_cost)
1917 DEALLOCATE (tmp_cpu_index, my_cost_cpu)
1918 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 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 min_ij = min(i, j)
1941 get_1d_idx = min_ij*n + max(i, j) - (min_ij - 1)*min_ij/2 - n
1942
1943 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 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 coeffs_set, coeffs_kind, &
1987 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 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 INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, npgfc, npgfd, nsgfa, &
2026 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 REAL(dp), DIMENSION(:, :), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4
2031
2032 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 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 log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
2043
2044 DO i_list_ij = 1, list_ij%n_element
2045 iatom = list_ij%elements(i_list_ij)%pair(1)
2046 jatom = list_ij%elements(i_list_ij)%pair(2)
2047 i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
2048 i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
2049 ikind = list_ij%elements(i_list_ij)%kind_pair(1)
2050 jkind = list_ij%elements(i_list_ij)%kind_pair(2)
2051 rab2 = list_ij%elements(i_list_ij)%dist2
2052
2053 nsgfa => basis_parameter(ikind)%nsgf
2054 nsgfb => basis_parameter(jkind)%nsgf
2055 npgfa => basis_parameter(ikind)%npgf
2056 npgfb => basis_parameter(jkind)%npgf
2057
2058 DO i_list_kl = 1, list_kl%n_element
2059
2060 katom = list_kl%elements(i_list_kl)%pair(1)
2061 latom = list_kl%elements(i_list_kl)%pair(2)
2062
2063 IF (.NOT. (katom + latom <= iatom + jatom)) cycle
2064 IF (((iatom + jatom) == (katom + latom)) .AND. (katom < iatom)) cycle
2065
2066 IF (eval_type == hfx_do_eval_forces) THEN
2067 IF (.NOT. use_virial) THEN
2068 IF ((iatom == jatom .AND. iatom == katom .AND. iatom == latom)) cycle
2069 END IF
2070 END IF
2071
2072 i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
2073 i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
2074 kkind = list_kl%elements(i_list_kl)%kind_pair(1)
2075 lkind = list_kl%elements(i_list_kl)%kind_pair(2)
2076 rcd2 = list_kl%elements(i_list_kl)%dist2
2077
2078 nsgfc => basis_parameter(kkind)%nsgf
2079 nsgfd => basis_parameter(lkind)%nsgf
2080 npgfc => basis_parameter(kkind)%npgf
2081 npgfd => basis_parameter(lkind)%npgf
2082
2083 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 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 coeffs_kind(jkind, ikind)%x(2)
2094 screen_kind_kl = coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
2095 coeffs_kind(lkind, kkind)%x(2)
2096 IF (screen_kind_ij + screen_kind_kl + actual_pmax_atom < log10_eps_schwarz) cycle
2097
2098 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 IF (do_p_screening) THEN
2104 SELECT CASE (eval_type)
2105 CASE (hfx_do_eval_energy)
2106 swap_id = 0
2107 kind_kind_idx = int(get_1d_idx(kkind, ikind, int(nkind, int_8)))
2108 IF (ikind >= kkind) THEN
2109 ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2110 map_atom_to_kind_atom(katom), &
2111 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 map_atom_to_kind_atom(katom))
2116 swap_id = swap_id + 1
2117 END IF
2118 kind_kind_idx = int(get_1d_idx(lkind, jkind, int(nkind, int_8)))
2119 IF (jkind >= lkind) THEN
2120 ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2121 map_atom_to_kind_atom(latom), &
2122 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 map_atom_to_kind_atom(latom))
2127 swap_id = swap_id + 2
2128 END IF
2129 kind_kind_idx = int(get_1d_idx(lkind, ikind, int(nkind, int_8)))
2130 IF (ikind >= lkind) THEN
2131 ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2132 map_atom_to_kind_atom(latom), &
2133 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 map_atom_to_kind_atom(latom))
2138 swap_id = swap_id + 4
2139 END IF
2140 kind_kind_idx = int(get_1d_idx(kkind, jkind, int(nkind, int_8)))
2141 IF (jkind >= kkind) THEN
2142 ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2143 map_atom_to_kind_atom(katom), &
2144 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 map_atom_to_kind_atom(katom))
2149 swap_id = swap_id + 8
2150 END IF
2151 CASE (hfx_do_eval_forces)
2152 swap_id = 16
2153 kind_kind_idx = int(get_1d_idx(kkind, ikind, int(nkind, int_8)))
2154 IF (ikind >= kkind) THEN
2155 ptr_p_1 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2156 map_atom_to_kind_atom(katom), &
2157 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 map_atom_to_kind_atom(katom))
2162 swap_id = swap_id + 1
2163 END IF
2164 kind_kind_idx = int(get_1d_idx(lkind, jkind, int(nkind, int_8)))
2165 IF (jkind >= lkind) THEN
2166 ptr_p_2 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2167 map_atom_to_kind_atom(latom), &
2168 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 map_atom_to_kind_atom(latom))
2173 swap_id = swap_id + 2
2174 END IF
2175 kind_kind_idx = int(get_1d_idx(lkind, ikind, int(nkind, int_8)))
2176 IF (ikind >= lkind) THEN
2177 ptr_p_3 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2178 map_atom_to_kind_atom(latom), &
2179 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 map_atom_to_kind_atom(latom))
2184 swap_id = swap_id + 4
2185 END IF
2186 kind_kind_idx = int(get_1d_idx(kkind, jkind, int(nkind, int_8)))
2187 IF (jkind >= kkind) THEN
2188 ptr_p_4 => pmax_set(kind_kind_idx)%p_kind(:, :, &
2189 map_atom_to_kind_atom(katom), &
2190 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 map_atom_to_kind_atom(katom))
2195 swap_id = swap_id + 8
2196 END IF
2197 END SELECT
2198 END IF
2199
2200 DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
2201 iset = set_list_ij(i_set_list_ij)%pair(1)
2202 jset = set_list_ij(i_set_list_ij)%pair(2)
2203
2204 max_val1 = coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + &
2205 coeffs_set(jset, iset, jkind, ikind)%x(2)
2206
2207 IF (max_val1 + screen_kind_kl + actual_pmax_atom < log10_eps_schwarz) cycle
2208 DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
2209 kset = set_list_kl(i_set_list_kl)%pair(1)
2210 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 coeffs_set(lset, kset, lkind, kkind)%x(2))
2214
2215 IF (max_val2 + actual_pmax_atom < log10_eps_schwarz) cycle
2216 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 pmax_entry, swap_id)
2220 IF (eval_type == hfx_do_eval_forces) THEN
2221 pmax_entry = log_2 + pmax_entry
2222 END IF
2223 ELSE
2224 pmax_entry = 0.0_dp
2225 END IF
2226 max_val2 = max_val2 + pmax_entry
2227 IF (max_val2 < log10_eps_schwarz) cycle
2228 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 p1_energy, p2_energy, p3_energy)
2234 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 p1_forces, p2_forces, p3_forces)
2240 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 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 SUBROUTINE init_blocks(nkind, para_env, natom, block_size, nblock, blocks, &
2282 list_ij, list_kl, set_list_ij, set_list_kl, &
2283 particle_set, &
2284 coeffs_set, coeffs_kind, &
2285 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 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 DO atom_block = 0, nblock - 1
2323 iatom_block = modulo(atom_block, nblock) + 1
2324 iatom_start = (iatom_block - 1)*block_size + 1
2325 iatom_end = min(iatom_block*block_size, natom)
2326 blocks(atom_block + 1)%istart = iatom_start
2327 blocks(atom_block + 1)%iend = iatom_end
2328 blocks(atom_block + 1)%cost = 0_int_8
2329 END DO
2330
2331 ncpus = para_env%num_pe
2332 my_cpu_rank = para_env%mepos
2333 DO i = 1, nblock
2334 IF (modulo(i, ncpus) /= my_cpu_rank) THEN
2335 blocks(i)%istart = 0
2336 blocks(i)%iend = 0
2337 cycle
2338 END IF
2339 iatom_start = blocks(i)%istart
2340 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 log10_eps_schwarz, log_2, coeffs_kind_max0, use_virial, atomic_pair_list)
2351
2352 END DO
2353 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 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 INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: buffer, buffer_in, buffer_out, summary
2376 INTEGER(int_8), ALLOCATABLE, DIMENSION(:), SAVE :: shm_cost_vector
2377 INTEGER, ALLOCATABLE, DIMENSION(:) :: bins_per_rank, rdispl, sort_idx
2378 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: shm_bins_per_rank, shm_displ
2379
2380 SELECT CASE (eval_type)
2381 CASE (hfx_do_eval_energy)
2382 nbins = SIZE(x_data%distribution_energy)
2383 CASE (hfx_do_eval_forces)
2384 nbins = SIZE(x_data%distribution_forces)
2385 END SELECT
2386
2387!$OMP MASTER
2388 ALLOCATE (shm_bins_per_rank(n_threads))
2389 ALLOCATE (shm_displ(n_threads + 1))
2390!$OMP END MASTER
2391!$OMP BARRIER
2392
2393 shm_bins_per_rank(i_thread + 1) = nbins
2394!$OMP BARRIER
2395 nbins = 0
2396 DO i = 1, n_threads
2397 nbins = nbins + shm_bins_per_rank(i)
2398 END DO
2399 my_rank = para_env%mepos
2400 nranks = para_env%num_pe
2401
2402!$OMP BARRIER
2403!$OMP MASTER
2404 ALLOCATE (bins_per_rank(nranks))
2405 bins_per_rank = 0
2406
2407 bins_per_rank(my_rank + 1) = nbins
2408
2409 CALL para_env%sum(bins_per_rank)
2410
2411 total_bins = 0
2412 DO i = 1, nranks
2413 total_bins = total_bins + bins_per_rank(i)
2414 END DO
2415
2416 ALLOCATE (shm_cost_vector(2*total_bins))
2417 shm_cost_vector = -1_int_8
2418 shm_displ(1) = 1
2419 DO i = 2, n_threads
2420 shm_displ(i) = shm_displ(i - 1) + shm_bins_per_rank(i - 1)
2421 END DO
2422 shm_displ(n_threads + 1) = nbins + 1
2423!$OMP END MASTER
2424!$OMP BARRIER
2425 j = 0
2426 SELECT CASE (eval_type)
2427 CASE (hfx_do_eval_energy)
2428 DO i = shm_displ(i_thread + 1), shm_displ(i_thread + 2) - 1
2429 j = j + 1
2430 shm_cost_vector(2*(i - 1) + 1) = x_data%distribution_energy(j)%cost
2431 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 DO i = shm_displ(i_thread + 1), shm_displ(i_thread + 2) - 1
2435 j = j + 1
2436 shm_cost_vector(2*(i - 1) + 1) = x_data%distribution_forces(j)%cost
2437 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!$OMP BARRIER
2441!$OMP MASTER
2442 ! ** calculate offsets
2443 ALLOCATE (rdispl(nranks))
2444 bins_per_rank(:) = bins_per_rank(:)*2
2445 rdispl(1) = 0
2446 DO i = 2, nranks
2447 rdispl(i) = rdispl(i - 1) + bins_per_rank(i - 1)
2448 END DO
2449
2450 ALLOCATE (buffer_in(2*nbins))
2451 ALLOCATE (buffer_out(2*total_bins))
2452
2453 DO i = 1, nbins
2454 buffer_in(2*(i - 1) + 1) = shm_cost_vector(2*(i - 1) + 1)
2455 buffer_in(2*i) = shm_cost_vector(2*i)
2456 END DO
2457
2458 CALL para_env%gatherv(buffer_in, buffer_out, bins_per_rank, rdispl)
2459
2460 IF (iw > 0) THEN
2461
2462 ALLOCATE (summary(2*nranks))
2463 summary = 0_int_8
2464
2465 WRITE (iw, '( /, 1X, 79("-") )')
2466 WRITE (iw, '( " -", 77X, "-" )')
2467 SELECT CASE (eval_type)
2468 CASE (hfx_do_eval_energy)
2469 WRITE (iw, '( " -", 20X, A, 19X, "-" )') ' HFX LOAD BALANCE INFORMATION - ENERGY '
2470 CASE (hfx_do_eval_forces)
2471 WRITE (iw, '( " -", 20X, A, 19X, "-" )') ' HFX LOAD BALANCE INFORMATION - FORCES '
2472 END SELECT
2473 WRITE (iw, '( " -", 77X, "-" )')
2474 WRITE (iw, '( 1X, 79("-") )')
2475
2476 WRITE (iw, fmt="(T3,A,T15,A,T35,A,T55,A)") "MPI RANK", "BIN #", "EST cost", "Processing time [s]"
2477 WRITE (iw, '( 1X, 79("-"), / )')
2478 k = 0
2479 DO i = 1, nranks
2480 DO j = 1, bins_per_rank(i)/2
2481 k = k + 1
2482 WRITE (iw, fmt="(T6,I5,T15,I5,T27,I16,T55,F19.8)") &
2483 i - 1, j, buffer_out(2*(k - 1) + 1), real(buffer_out(2*k), dp)/10000.0_dp
2484 summary(2*(i - 1) + 1) = summary(2*(i - 1) + 1) + buffer_out(2*(k - 1) + 1)
2485 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 DO i = 1, total_bins
2494 sum_bin = sum_bin + buffer_out(2*i)
2495 max_bin = max(max_bin, buffer_out(2*i))
2496 min_bin = min(min_bin, buffer_out(2*i))
2497 END DO
2498 avg_bin = sum_bin/total_bins
2499
2500 max_rank = 0_int_8
2501 min_rank = huge(min_rank)
2502 sum_rank = 0_int_8
2503 DO i = 1, nranks
2504 sum_rank = sum_rank + summary(2*i)
2505 max_rank = max(max_rank, summary(2*i))
2506 min_rank = min(min_rank, summary(2*i))
2507 END DO
2508 avg_rank = sum_rank/nranks
2509
2510 WRITE (iw, fmt='(/,T3,A,/)') "SUMMARY:"
2511 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Max bin", real(max_bin, dp)/10000.0_dp
2512 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Min bin", real(min_bin, dp)/10000.0_dp
2513 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Sum bin", real(sum_bin, dp)/10000.0_dp
2514 WRITE (iw, fmt="(T3,A,T35,F19.8,/)") "Avg bin", real(avg_bin, dp)/10000.0_dp
2515 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Max rank", real(max_rank, dp)/10000.0_dp
2516 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Min rank", real(min_rank, dp)/10000.0_dp
2517 WRITE (iw, fmt="(T3,A,T35,F19.8)") "Sum rank", real(sum_rank, dp)/10000.0_dp
2518 WRITE (iw, fmt="(T3,A,T35,F19.8,/)") "Avg rank", real(avg_rank, dp)/10000.0_dp
2519
2520 ALLOCATE (buffer(nranks))
2521 ALLOCATE (sort_idx(nranks))
2522
2523 DO i = 1, nranks
2524 buffer(i) = summary(2*i)
2525 END DO
2526
2527 CALL sort(buffer, nranks, sort_idx)
2528
2529 WRITE (iw, fmt="(T3,A,T35,A,T55,A,/)") "MPI RANK", "EST cost", "Processing time [s]"
2530 DO i = nranks, 1, -1
2531 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 DEALLOCATE (summary, buffer, sort_idx)
2535
2536 END IF
2537
2538 DEALLOCATE (buffer_in, buffer_out, rdispl)
2539
2540 CALL para_env%sync()
2541
2542 DEALLOCATE (shm_bins_per_rank, shm_displ, shm_cost_vector)
2543!$OMP END MASTER
2544!$OMP BARRIER
2545
2546 END SUBROUTINE collect_load_balance_info
2547
2548
2549! **************************************************************************************************
2550!> \brief This routine calculates the maximum density matrix element, when
2551!> screening on an initial density matrix is applied. Due to symmetry of
2552!> the ERI's, there are always 4 matrix elements to be considered.
2553!> CASE 0-15 belong to an energy calculation (linear screening)
2554!> CASE 16-31 belong to a force calculation (square screening)
2555!> \param ptr_p_1 Pointers to atomic density matrices
2556!> \param ptr_p_2 Pointers to atomic density matrices
2557!> \param ptr_p_3 Pointers to atomic density matrices
2558!> \param ptr_p_4 Pointers to atomic density matrices
2559!> \param iset Current set
2560!> \param jset Current set
2561!> \param kset Current set
2562!> \param lset Current set
2563!> \param pmax_val value to be calculated
2564!> \param swap_id Defines how the matrices are accessed
2565!> \par History
2566!> 06.2009 created [Manuel Guidon]
2567!> \author Manuel Guidon
2568! **************************************************************************************************
2569PURE SUBROUTINE get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, iset, jset, kset, lset, pmax_val, swap_id)
2570
2571 REAL(dp), DIMENSION(:, :), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4
2572 INTEGER, INTENT(IN) :: iset, jset, kset, lset
2573
2574 REAL(dp), INTENT(OUT) :: pmax_val
2575 INTEGER, INTENT(IN) :: swap_id
2576
2577 REAL(dp) :: pmax_1, pmax_2, pmax_3, pmax_4
2578
2579 SELECT CASE (swap_id)
2580 CASE (0)
2581 pmax_1 = ptr_p_1(kset, iset)
2582 pmax_2 = ptr_p_2(lset, jset)
2583 pmax_3 = ptr_p_3(lset, iset)
2584 pmax_4 = ptr_p_4(kset, jset)
2585 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2586 CASE (1)
2587 pmax_1 = ptr_p_1(iset, kset)
2588 pmax_2 = ptr_p_2(lset, jset)
2589 pmax_3 = ptr_p_3(lset, iset)
2590 pmax_4 = ptr_p_4(kset, jset)
2591 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2592 CASE (2)
2593 pmax_1 = ptr_p_1(kset, iset)
2594 pmax_2 = ptr_p_2(jset, lset)
2595 pmax_3 = ptr_p_3(lset, iset)
2596 pmax_4 = ptr_p_4(kset, jset)
2597 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2598 CASE (3)
2599 pmax_1 = ptr_p_1(iset, kset)
2600 pmax_2 = ptr_p_2(jset, lset)
2601 pmax_3 = ptr_p_3(lset, iset)
2602 pmax_4 = ptr_p_4(kset, jset)
2603 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2604 CASE (4)
2605 pmax_1 = ptr_p_1(kset, iset)
2606 pmax_2 = ptr_p_2(lset, jset)
2607 pmax_3 = ptr_p_3(iset, lset)
2608 pmax_4 = ptr_p_4(kset, jset)
2609 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2610 CASE (5)
2611 pmax_1 = ptr_p_1(iset, kset)
2612 pmax_2 = ptr_p_2(lset, jset)
2613 pmax_3 = ptr_p_3(iset, lset)
2614 pmax_4 = ptr_p_4(kset, jset)
2615 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2616 CASE (6)
2617 pmax_1 = ptr_p_1(kset, iset)
2618 pmax_2 = ptr_p_2(jset, lset)
2619 pmax_3 = ptr_p_3(iset, lset)
2620 pmax_4 = ptr_p_4(kset, jset)
2621 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2622 CASE (7)
2623 pmax_1 = ptr_p_1(iset, kset)
2624 pmax_2 = ptr_p_2(jset, lset)
2625 pmax_3 = ptr_p_3(iset, lset)
2626 pmax_4 = ptr_p_4(kset, jset)
2627 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2628 CASE (8)
2629 pmax_1 = ptr_p_1(kset, iset)
2630 pmax_2 = ptr_p_2(lset, jset)
2631 pmax_3 = ptr_p_3(lset, iset)
2632 pmax_4 = ptr_p_4(jset, kset)
2633 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2634 CASE (9)
2635 pmax_1 = ptr_p_1(iset, kset)
2636 pmax_2 = ptr_p_2(lset, jset)
2637 pmax_3 = ptr_p_3(lset, iset)
2638 pmax_4 = ptr_p_4(jset, kset)
2639 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2640 CASE (10)
2641 pmax_1 = ptr_p_1(kset, iset)
2642 pmax_2 = ptr_p_2(jset, lset)
2643 pmax_3 = ptr_p_3(lset, iset)
2644 pmax_4 = ptr_p_4(jset, kset)
2645 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2646 CASE (11)
2647 pmax_1 = ptr_p_1(iset, kset)
2648 pmax_2 = ptr_p_2(jset, lset)
2649 pmax_3 = ptr_p_3(lset, iset)
2650 pmax_4 = ptr_p_4(jset, kset)
2651 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2652 CASE (12)
2653 pmax_1 = ptr_p_1(kset, iset)
2654 pmax_2 = ptr_p_2(lset, jset)
2655 pmax_3 = ptr_p_3(iset, lset)
2656 pmax_4 = ptr_p_4(jset, kset)
2657 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2658 CASE (13)
2659 pmax_1 = ptr_p_1(iset, kset)
2660 pmax_2 = ptr_p_2(lset, jset)
2661 pmax_3 = ptr_p_3(iset, lset)
2662 pmax_4 = ptr_p_4(jset, kset)
2663 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2664 CASE (14)
2665 pmax_1 = ptr_p_1(kset, iset)
2666 pmax_2 = ptr_p_2(jset, lset)
2667 pmax_3 = ptr_p_3(iset, lset)
2668 pmax_4 = ptr_p_4(jset, kset)
2669 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2670 CASE (15)
2671 pmax_1 = ptr_p_1(iset, kset)
2672 pmax_2 = ptr_p_2(jset, lset)
2673 pmax_3 = ptr_p_3(iset, lset)
2674 pmax_4 = ptr_p_4(jset, kset)
2675 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
2676 CASE (16)
2677 pmax_1 = ptr_p_1(kset, iset)
2678 pmax_2 = ptr_p_2(lset, jset)
2679 pmax_3 = ptr_p_3(lset, iset)
2680 pmax_4 = ptr_p_4(kset, jset)
2681 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2682 CASE (17)
2683 pmax_1 = ptr_p_1(iset, kset)
2684 pmax_2 = ptr_p_2(lset, jset)
2685 pmax_3 = ptr_p_3(lset, iset)
2686 pmax_4 = ptr_p_4(kset, jset)
2687 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2688 CASE (18)
2689 pmax_1 = ptr_p_1(kset, iset)
2690 pmax_2 = ptr_p_2(jset, lset)
2691 pmax_3 = ptr_p_3(lset, iset)
2692 pmax_4 = ptr_p_4(kset, jset)
2693 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2694 CASE (19)
2695 pmax_1 = ptr_p_1(iset, kset)
2696 pmax_2 = ptr_p_2(jset, lset)
2697 pmax_3 = ptr_p_3(lset, iset)
2698 pmax_4 = ptr_p_4(kset, jset)
2699 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2700 CASE (20)
2701 pmax_1 = ptr_p_1(kset, iset)
2702 pmax_2 = ptr_p_2(lset, jset)
2703 pmax_3 = ptr_p_3(iset, lset)
2704 pmax_4 = ptr_p_4(kset, jset)
2705 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2706 CASE (21)
2707 pmax_1 = ptr_p_1(iset, kset)
2708 pmax_2 = ptr_p_2(lset, jset)
2709 pmax_3 = ptr_p_3(iset, lset)
2710 pmax_4 = ptr_p_4(kset, jset)
2711 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2712 CASE (22)
2713 pmax_1 = ptr_p_1(kset, iset)
2714 pmax_2 = ptr_p_2(jset, lset)
2715 pmax_3 = ptr_p_3(iset, lset)
2716 pmax_4 = ptr_p_4(kset, jset)
2717 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2718 CASE (23)
2719 pmax_1 = ptr_p_1(iset, kset)
2720 pmax_2 = ptr_p_2(jset, lset)
2721 pmax_3 = ptr_p_3(iset, lset)
2722 pmax_4 = ptr_p_4(kset, jset)
2723 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2724 CASE (24)
2725 pmax_1 = ptr_p_1(kset, iset)
2726 pmax_2 = ptr_p_2(lset, jset)
2727 pmax_3 = ptr_p_3(lset, iset)
2728 pmax_4 = ptr_p_4(jset, kset)
2729 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2730 CASE (25)
2731 pmax_1 = ptr_p_1(iset, kset)
2732 pmax_2 = ptr_p_2(lset, jset)
2733 pmax_3 = ptr_p_3(lset, iset)
2734 pmax_4 = ptr_p_4(jset, kset)
2735 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2736 CASE (26)
2737 pmax_1 = ptr_p_1(kset, iset)
2738 pmax_2 = ptr_p_2(jset, lset)
2739 pmax_3 = ptr_p_3(lset, iset)
2740 pmax_4 = ptr_p_4(jset, kset)
2741 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2742 CASE (27)
2743 pmax_1 = ptr_p_1(iset, kset)
2744 pmax_2 = ptr_p_2(jset, lset)
2745 pmax_3 = ptr_p_3(lset, iset)
2746 pmax_4 = ptr_p_4(jset, kset)
2747 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2748 CASE (28)
2749 pmax_1 = ptr_p_1(kset, iset)
2750 pmax_2 = ptr_p_2(lset, jset)
2751 pmax_3 = ptr_p_3(iset, lset)
2752 pmax_4 = ptr_p_4(jset, kset)
2753 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2754 CASE (29)
2755 pmax_1 = ptr_p_1(iset, kset)
2756 pmax_2 = ptr_p_2(lset, jset)
2757 pmax_3 = ptr_p_3(iset, lset)
2758 pmax_4 = ptr_p_4(jset, kset)
2759 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2760 CASE (30)
2761 pmax_1 = ptr_p_1(kset, iset)
2762 pmax_2 = ptr_p_2(jset, lset)
2763 pmax_3 = ptr_p_3(iset, lset)
2764 pmax_4 = ptr_p_4(jset, kset)
2765 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2766 CASE (31)
2767 pmax_1 = ptr_p_1(iset, kset)
2768 pmax_2 = ptr_p_2(jset, lset)
2769 pmax_3 = ptr_p_3(iset, lset)
2770 pmax_4 = ptr_p_4(jset, kset)
2771 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
2772 END SELECT
2773
2774END SUBROUTINE get_pmax_val
2775
2776
2777END MODULE hfx_load_balance_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:162
message_passing::mp_waitall
Definition message_passing.F:892
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
hfx_load_balance_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_load_balance_methods.F:14
hfx_load_balance_methods::p1_energy
real(kind=dp), dimension(12), parameter, public p1_energy
Definition hfx_load_balance_methods.F:45
hfx_load_balance_methods::p3_energy
real(kind=dp), dimension(2), parameter, public p3_energy
Definition hfx_load_balance_methods.F:57
hfx_load_balance_methods::p2_energy
real(kind=dp), dimension(12), parameter, public p2_energy
Definition hfx_load_balance_methods.F:51
hfx_load_balance_methods::collect_load_balance_info
subroutine, public collect_load_balance_info(para_env, x_data, iw, n_threads, i_thread, eval_type)
...
Definition hfx_load_balance_methods.F:2366
hfx_load_balance_methods::hfx_load_balance
subroutine, public hfx_load_balance(x_data, eps_schwarz, particle_set, max_set, para_env, coeffs_set, coeffs_kind, is_assoc_atomic_block_global, do_periodic, load_balance_parameter, kind_of, basis_parameter, pmax_set, pmax_atom, i_thread, n_threads, cell, do_p_screening, map_atom_to_kind_atom, nkind, eval_type, pmax_block, use_virial)
Distributes the computation of eri's to all available processes.
Definition hfx_load_balance_methods.F:130
hfx_load_balance_methods::cost_model
integer(kind=int_8) function, public cost_model(nsa, nsb, nsc, nsd, npgfa, npgfb, npgfc, npgfd, ratio, p1, p2, p3)
estimates the cost of a set quartet with info available at load balance time i.e. without much info o...
Definition hfx_load_balance_methods.F:1801
hfx_load_balance_methods::hfx_update_load_balance
subroutine, public hfx_update_load_balance(x_data, para_env, load_balance_parameter, i_thread, n_threads, eval_type)
Cheap way of redistributing the eri's.
Definition hfx_load_balance_methods.F:1518
hfx_pair_list_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_pair_list_methods.F:15
hfx_pair_list_methods::build_atomic_pair_list
subroutine, public build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, blocks)
...
Definition hfx_pair_list_methods.F:287
hfx_pair_list_methods::build_pair_list
subroutine, public build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
...
Definition hfx_pair_list_methods.F:192
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
hfx_types::hfx_set_distr_energy
subroutine, public hfx_set_distr_energy(ptr_to_distr, x_data)
This routine stores the data obtained from the load balance routine for the energy
Definition hfx_types.F:2588
hfx_types::hfx_set_distr_forces
subroutine, public hfx_set_distr_forces(ptr_to_distr, x_data)
This routine stores the data obtained from the load balance routine for the forces
Definition hfx_types.F:2608
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::hfx_do_eval_energy
integer, parameter, public hfx_do_eval_energy
Definition input_constants.F:890
input_constants::hfx_do_eval_forces
integer, parameter, public hfx_do_eval_forces
Definition input_constants.F:890
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::uniform
integer, parameter, public uniform
Definition parallel_rng_types.F:73
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
util
All kind of helpful little routines.
Definition util.F:14
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
hfx_types::hfx_basis_type
Definition hfx_types.F:256
hfx_types::hfx_block_range_type
Definition hfx_types.F:296
hfx_types::hfx_distribution
Definition hfx_types.F:197
hfx_types::hfx_load_balance_type
Definition hfx_types.F:180
hfx_types::hfx_p_kind
Definition hfx_types.F:285
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:474
hfx_types::pair_list_type
Definition hfx_types.F:222
hfx_types::pair_set_list_type
Definition hfx_types.F:217
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
message_passing::mp_request_type
Definition message_passing.F:634
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_types::particle_type
Definition particle_types.F:35