(git:db0abe2)
Loading...
Searching...
No Matches
task_list_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 generate the tasks lists used by collocate and integrate routines
10!> \par History
11!> 01.2008 [Joost VandeVondele] refactered out of qs_collocate / qs_integrate
12!> \author Joost VandeVondele
13! **************************************************************************************************
14MODULE task_list_methods
15 USE offload_api, ONLY: offload_create_buffer, offload_buffer_type
16 USE grid_api, ONLY: grid_create_basis_set, grid_create_task_list
17 USE ao_util, ONLY: exp_radius_very_extended
18 USE basis_set_types, ONLY: get_gto_basis_set, &
19 gto_basis_set_p_type, &
20 gto_basis_set_type
21 USE cell_types, ONLY: cell_type, &
22 pbc
23 USE cp_control_types, ONLY: dft_control_type
24 USE cube_utils, ONLY: compute_cube_center, &
25 cube_info_type, &
26 return_cube, &
27 return_cube_nonortho
28 USE cp_dbcsr_api, ONLY: dbcsr_convert_sizes_to_offsets, &
29 dbcsr_finalize, &
30 dbcsr_get_block_p, &
31 dbcsr_get_info, &
32 dbcsr_p_type, &
33 dbcsr_put_block, &
34 dbcsr_type, &
35 dbcsr_work_create
36 USE gaussian_gridlevels, ONLY: gaussian_gridlevel, &
37 gridlevel_info_type
38 USE kinds, ONLY: default_string_length, &
39 dp, &
40 int_8
41 USE kpoint_types, ONLY: get_kpoint_info, &
42 kpoint_type
43 USE memory_utilities, ONLY: reallocate
44 USE message_passing, ONLY: &
45 mp_comm_type
46 USE particle_types, ONLY: particle_type
47 USE particle_methods, ONLY: get_particle_set
48 USE pw_env_types, ONLY: pw_env_get, &
49 pw_env_type
50 USE qs_kind_types, ONLY: get_qs_kind, &
51 qs_kind_type
52 USE qs_ks_types, ONLY: get_ks_env, &
53 qs_ks_env_type
54 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
55 USE realspace_grid_types, ONLY: realspace_grid_desc_p_type, &
56 realspace_grid_desc_type, &
57 rs_grid_create, &
58 rs_grid_locate_rank, &
59 rs_grid_release, &
60 rs_grid_reorder_ranks, realspace_grid_type
61 USE task_list_types, ONLY: deserialize_task, &
62 reallocate_tasks, &
63 serialize_task, &
64 task_list_type, &
65 atom_pair_type, &
66 task_size_in_int8, &
67 task_type
68 USE util, ONLY: sort
69
70!$ USE OMP_LIB, ONLY: omp_destroy_lock, omp_get_num_threads, omp_init_lock, &
71!$ omp_lock_kind, omp_set_lock, omp_unset_lock
72#include "./base/base_uses.f90"
73
74
75 IMPLICIT NONE
76
77 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
78
79 PRIVATE
80
81 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'task_list_methods'
82
83 PUBLIC :: generate_qs_task_list, &
84 task_list_inner_loop
85 PUBLIC :: distribute_tasks, &
86 rs_distribute_matrix, &
87 rs_scatter_matrices, &
88 rs_gather_matrices, &
89 rs_copy_to_buffer, &
90 rs_copy_to_matrices
91
92CONTAINS
93
94! **************************************************************************************************
95!> \brief ...
96!> \param ks_env ...
97!> \param task_list ...
98!> \param basis_type ...
99!> \param reorder_rs_grid_ranks Flag that indicates if this routine should
100!> or should not overwrite the rs descriptor (see comment below)
101!> \param skip_load_balance_distributed ...
102!> \param pw_env_external ...
103!> \param sab_orb_external ...
104!> \param ext_kpoints ...
105!> \par History
106!> 01.2008 factored out of calculate_rho_elec [Joost VandeVondele]
107!> 04.2010 divides tasks into grid levels and atom pairs for integrate/collocate [Iain Bethune]
108!> (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
109!> 06.2015 adjusted to be used with multiple images (k-points) [JGH]
110!> \note If this routine is called several times with different task lists,
111!> the default behaviour is to re-optimize the grid ranks and overwrite
112!> the rs descriptor and grids. reorder_rs_grid_ranks = .FALSE. prevents the code
113!> of performing a new optimization by leaving the rank order in
114!> its current state.
115! **************************************************************************************************
116 SUBROUTINE generate_qs_task_list(ks_env, task_list, basis_type, &
117 reorder_rs_grid_ranks, skip_load_balance_distributed, &
118 pw_env_external, sab_orb_external, ext_kpoints)
119
120 TYPE(qs_ks_env_type), POINTER :: ks_env
121 TYPE(task_list_type), POINTER :: task_list
122 CHARACTER(LEN=*), INTENT(IN) :: basis_type
123 LOGICAL, INTENT(IN) :: reorder_rs_grid_ranks, &
124 skip_load_balance_distributed
125 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
126 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
127 OPTIONAL, POINTER :: sab_orb_external
128 TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
129
130 CHARACTER(LEN=*), PARAMETER :: routinen = 'generate_qs_task_list'
131 INTEGER, PARAMETER :: max_tasks = 2000
132
133 INTEGER :: cindex, curr_tasks, handle, i, iatom, iatom_old, igrid_level, igrid_level_old, &
134 ikind, ilevel, img, img_old, ipair, ipgf, iset, itask, jatom, jatom_old, jkind, jpgf, &
135 jset, maxpgf, maxset, natoms, nimages, nkind, nseta, nsetb, slot
136 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: blocks
137 INTEGER, DIMENSION(3) :: cellind
138 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
139 npgfb, nsgf
140 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
141 LOGICAL :: dokp
142 REAL(kind=dp) :: kind_radius_a, kind_radius_b
143 REAL(kind=dp), DIMENSION(3) :: ra, rab
144 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
145 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb
146 TYPE(cell_type), POINTER :: cell
147 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
148 TYPE(dft_control_type), POINTER :: dft_control
149 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
150 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
151 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, orb_basis_set
152 TYPE(kpoint_type), POINTER :: kpoints
153 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
154 POINTER :: sab_orb
155 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156 TYPE(pw_env_type), POINTER :: pw_env
157 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
158 TYPE(qs_kind_type), POINTER :: qs_kind
159 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
160 POINTER :: rs_descs
161 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
162 TYPE(task_type), DIMENSION(:), POINTER :: tasks
163
164 CALL timeset(routinen, handle)
165
166 CALL get_ks_env(ks_env, &
167 qs_kind_set=qs_kind_set, &
168 cell=cell, &
169 particle_set=particle_set, &
170 dft_control=dft_control)
171
172 CALL get_ks_env(ks_env, sab_orb=sab_orb)
173 IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
174
175 CALL get_ks_env(ks_env, pw_env=pw_env)
176 IF (PRESENT(pw_env_external)) pw_env => pw_env_external
177 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_grids)
178
179 ! *** assign from pw_env
180 gridlevel_info => pw_env%gridlevel_info
181 cube_info => pw_env%cube_info
182
183 ! find maximum numbers
184 nkind = SIZE(qs_kind_set)
185 natoms = SIZE(particle_set)
186 maxset = 0
187 maxpgf = 0
188 DO ikind = 1, nkind
189 qs_kind => qs_kind_set(ikind)
190 CALL get_qs_kind(qs_kind=qs_kind, &
191 basis_set=orb_basis_set, basis_type=basis_type)
192
193 IF (.NOT. ASSOCIATED(orb_basis_set)) cycle
194 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, npgf=npgfa, nset=nseta)
195
196 maxset = max(nseta, maxset)
197 maxpgf = max(maxval(npgfa), maxpgf)
198 END DO
199
200 ! kpoint related
201 nimages = dft_control%nimages
202 dokp = .false.
203 NULLIFY (cell_to_index)
204 NULLIFY (kpoints)
205 IF (PRESENT(ext_kpoints)) THEN
206 IF (ASSOCIATED(ext_kpoints)) THEN
207 dokp = .true.
208 CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
209 nimages = SIZE(ext_kpoints%index_to_cell, 2)
210 END IF
211 END IF
212 IF (.NOT. dokp .AND. nimages > 1) THEN
213 dokp = .true.
214 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
215 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
216 END IF
217
218 ! free the atom_pair lists if allocated
219 IF (ASSOCIATED(task_list%atom_pair_send)) DEALLOCATE (task_list%atom_pair_send)
220 IF (ASSOCIATED(task_list%atom_pair_recv)) DEALLOCATE (task_list%atom_pair_recv)
221
222 ! construct a list of all tasks
223 IF (.NOT. ASSOCIATED(task_list%tasks)) THEN
224 CALL reallocate_tasks(task_list%tasks, max_tasks)
225 END IF
226 task_list%ntasks = 0
227 curr_tasks = SIZE(task_list%tasks)
228
229 ALLOCATE (basis_set_list(nkind))
230 DO ikind = 1, nkind
231 qs_kind => qs_kind_set(ikind)
232 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, &
233 basis_type=basis_type)
234 IF (ASSOCIATED(basis_set_a)) THEN
235 basis_set_list(ikind)%gto_basis_set => basis_set_a
236 ELSE
237 NULLIFY (basis_set_list(ikind)%gto_basis_set)
238 END IF
239 END DO
240!!$OMP PARALLEL DEFAULT(NONE) &
241!!$OMP SHARED (sab_orb, dokp, basis_set_list, task_list, rs_descs, dft_control, cube_info, gridlevel_info, &
242!!$OMP curr_tasks, maxpgf, maxset, natoms, nimages, particle_set, cell_to_index, cell) &
243!!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, cellind, basis_set_a, basis_set_b, ra, &
244!!$OMP la_max, la_min, npgfa, nseta, rpgfa, set_radius_a, kind_radius_a, zeta, &
245!!$OMP lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, kind_radius_b, zetb, &
246!!$OMP cindex, slot)
247 ! Loop over neighbor list
248!!$OMP DO SCHEDULE(GUIDED)
249 DO slot = 1, sab_orb(1)%nl_size
250 ikind = sab_orb(1)%nlist_task(slot)%ikind
251 jkind = sab_orb(1)%nlist_task(slot)%jkind
252 iatom = sab_orb(1)%nlist_task(slot)%iatom
253 jatom = sab_orb(1)%nlist_task(slot)%jatom
254 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
255 cellind(1:3) = sab_orb(1)%nlist_task(slot)%cell(1:3)
256
257 basis_set_a => basis_set_list(ikind)%gto_basis_set
258 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
259 basis_set_b => basis_set_list(jkind)%gto_basis_set
260 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
261 ra(:) = pbc(particle_set(iatom)%r, cell)
262 ! basis ikind
263 la_max => basis_set_a%lmax
264 la_min => basis_set_a%lmin
265 npgfa => basis_set_a%npgf
266 nseta = basis_set_a%nset
267 rpgfa => basis_set_a%pgf_radius
268 set_radius_a => basis_set_a%set_radius
269 kind_radius_a = basis_set_a%kind_radius
270 zeta => basis_set_a%zet
271 ! basis jkind
272 lb_max => basis_set_b%lmax
273 lb_min => basis_set_b%lmin
274 npgfb => basis_set_b%npgf
275 nsetb = basis_set_b%nset
276 rpgfb => basis_set_b%pgf_radius
277 set_radius_b => basis_set_b%set_radius
278 kind_radius_b = basis_set_b%kind_radius
279 zetb => basis_set_b%zet
280
281 IF (dokp) THEN
282 cindex = cell_to_index(cellind(1), cellind(2), cellind(3))
283 ELSE
284 cindex = 1
285 END IF
286
287 CALL task_list_inner_loop(task_list%tasks, task_list%ntasks, curr_tasks, &
288 rs_descs, dft_control, cube_info, gridlevel_info, cindex, &
289 iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, &
290 set_radius_a, set_radius_b, ra, rab, &
291 la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
292
293 END DO
294!!$OMP END PARALLEL
295
296 ! redistribute the task list so that all tasks map on the local rs grids
297 CALL distribute_tasks( &
298 rs_descs=rs_descs, ntasks=task_list%ntasks, natoms=natoms, &
299 tasks=task_list%tasks, atom_pair_send=task_list%atom_pair_send, &
300 atom_pair_recv=task_list%atom_pair_recv, symmetric=.true., &
301 reorder_rs_grid_ranks=reorder_rs_grid_ranks, &
302 skip_load_balance_distributed=skip_load_balance_distributed)
303
304 ! compute offsets for rs_scatter_matrix / rs_copy_matrix
305 ALLOCATE (nsgf(natoms))
306 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list, nsgf=nsgf)
307 IF (ASSOCIATED(task_list%atom_pair_send)) THEN
308 ! Needed whenever tasks use matrix blocks from other ranks.
309 CALL rs_calc_offsets(pairs=task_list%atom_pair_send, &
310 nsgf=nsgf, &
311 group_size=rs_descs(1)%rs_desc%group_size, &
312 pair_offsets=task_list%pair_offsets_send, &
313 rank_offsets=task_list%rank_offsets_send, &
314 rank_sizes=task_list%rank_sizes_send, &
315 buffer_size=task_list%buffer_size_send)
316 ELSE
317 ! A rebuilt task list can return to owner-local execution.
318 IF (ASSOCIATED(task_list%pair_offsets_send)) DEALLOCATE (task_list%pair_offsets_send)
319 IF (ASSOCIATED(task_list%rank_offsets_send)) DEALLOCATE (task_list%rank_offsets_send)
320 IF (ASSOCIATED(task_list%rank_sizes_send)) DEALLOCATE (task_list%rank_sizes_send)
321 task_list%buffer_size_send = 0
322 END IF
323 CALL rs_calc_offsets(pairs=task_list%atom_pair_recv, &
324 nsgf=nsgf, &
325 group_size=rs_descs(1)%rs_desc%group_size, &
326 pair_offsets=task_list%pair_offsets_recv, &
327 rank_offsets=task_list%rank_offsets_recv, &
328 rank_sizes=task_list%rank_sizes_recv, &
329 buffer_size=task_list%buffer_size_recv)
330 IF (ASSOCIATED(task_list%atom_pair_send)) THEN
331 CALL rs_group_atom_pairs(task_list%atom_pair_send, &
332 task_list%unpack_block_starts, task_list%unpack_pair_indices)
333 ELSE
334 CALL rs_group_atom_pairs(task_list%atom_pair_recv, &
335 task_list%unpack_block_starts, task_list%unpack_pair_indices)
336 END IF
337 DEALLOCATE (basis_set_list, nsgf)
338
339 ! If the rank order has changed, reallocate any of the distributed rs_grids
340 IF (reorder_rs_grid_ranks) THEN
341 DO i = 1, gridlevel_info%ngrid_levels
342 IF (rs_descs(i)%rs_desc%distributed) THEN
343 CALL rs_grid_release(rs_grids(i))
344 CALL rs_grid_create(rs_grids(i), rs_descs(i)%rs_desc)
345 END IF
346 END DO
347 END IF
348
349 CALL create_grid_task_list(task_list=task_list, &
350 qs_kind_set=qs_kind_set, &
351 particle_set=particle_set, &
352 cell=cell, &
353 basis_type=basis_type, &
354 rs_grids=rs_grids)
355
356 ! Now we have the final list of tasks, setup the task_list with the
357 ! data needed for the loops in integrate_v/calculate_rho
358
359 IF (ASSOCIATED(task_list%taskstart)) THEN
360 DEALLOCATE (task_list%taskstart)
361 END IF
362 IF (ASSOCIATED(task_list%taskstop)) THEN
363 DEALLOCATE (task_list%taskstop)
364 END IF
365 IF (ASSOCIATED(task_list%npairs)) THEN
366 DEALLOCATE (task_list%npairs)
367 END IF
368
369 ! First, count the number of unique atom pairs per grid level
370
371 ALLOCATE (task_list%npairs(SIZE(rs_descs)))
372
373 iatom_old = -1; jatom_old = -1; igrid_level_old = -1; img_old = -1
374 ipair = 0
375 task_list%npairs = 0
376
377 DO i = 1, task_list%ntasks
378 igrid_level = task_list%tasks(i)%grid_level
379 img = task_list%tasks(i)%image
380 iatom = task_list%tasks(i)%iatom
381 jatom = task_list%tasks(i)%jatom
382 iset = task_list%tasks(i)%iset
383 jset = task_list%tasks(i)%jset
384 ipgf = task_list%tasks(i)%ipgf
385 jpgf = task_list%tasks(i)%jpgf
386 IF (igrid_level /= igrid_level_old) THEN
387 IF (igrid_level_old /= -1) THEN
388 task_list%npairs(igrid_level_old) = ipair
389 END IF
390 ipair = 1
391 igrid_level_old = igrid_level
392 iatom_old = iatom
393 jatom_old = jatom
394 img_old = img
395 ELSE IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
396 ipair = ipair + 1
397 iatom_old = iatom
398 jatom_old = jatom
399 img_old = img
400 END IF
401 END DO
402 ! Take care of the last iteration
403 IF (task_list%ntasks /= 0) THEN
404 task_list%npairs(igrid_level) = ipair
405 END IF
406
407 ! Second, for each atom pair, find the indices in the task list
408 ! of the first and last task
409
410 ! Array sized for worst case
411 ALLOCATE (task_list%taskstart(maxval(task_list%npairs), SIZE(rs_descs)))
412 ALLOCATE (task_list%taskstop(maxval(task_list%npairs), SIZE(rs_descs)))
413
414 iatom_old = -1; jatom_old = -1; igrid_level_old = -1; img_old = -1
415 ipair = 0
416 task_list%taskstart = 0
417 task_list%taskstop = 0
418
419 DO i = 1, task_list%ntasks
420 igrid_level = task_list%tasks(i)%grid_level
421 img = task_list%tasks(i)%image
422 iatom = task_list%tasks(i)%iatom
423 jatom = task_list%tasks(i)%jatom
424 iset = task_list%tasks(i)%iset
425 jset = task_list%tasks(i)%jset
426 ipgf = task_list%tasks(i)%ipgf
427 jpgf = task_list%tasks(i)%jpgf
428 IF (igrid_level /= igrid_level_old) THEN
429 IF (igrid_level_old /= -1) THEN
430 task_list%taskstop(ipair, igrid_level_old) = i - 1
431 END IF
432 ipair = 1
433 task_list%taskstart(ipair, igrid_level) = i
434 igrid_level_old = igrid_level
435 iatom_old = iatom
436 jatom_old = jatom
437 img_old = img
438 ELSE IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
439 ipair = ipair + 1
440 task_list%taskstart(ipair, igrid_level) = i
441 task_list%taskstop(ipair - 1, igrid_level) = i - 1
442 iatom_old = iatom
443 jatom_old = jatom
444 img_old = img
445 END IF
446 END DO
447 ! Take care of the last iteration
448 IF (task_list%ntasks /= 0) THEN
449 task_list%taskstop(ipair, igrid_level) = task_list%ntasks
450 END IF
451
452 ! Debug task destribution
453 IF (debug_this_module) THEN
454 tasks => task_list%tasks
455 WRITE (6, *)
456 WRITE (6, *) "Total number of tasks ", task_list%ntasks
457 DO igrid_level = 1, gridlevel_info%ngrid_levels
458 WRITE (6, *) "Total number of pairs(grid_level) ", &
459 igrid_level, task_list%npairs(igrid_level)
460 END DO
461 WRITE (6, *)
462
463 DO igrid_level = 1, gridlevel_info%ngrid_levels
464
465 ALLOCATE (blocks(natoms, natoms, nimages))
466 blocks = -1
467 DO ipair = 1, task_list%npairs(igrid_level)
468 itask = task_list%taskstart(ipair, igrid_level)
469 ilevel = task_list%tasks(itask)%grid_level
470 img = task_list%tasks(itask)%image
471 iatom = task_list%tasks(itask)%iatom
472 jatom = task_list%tasks(itask)%jatom
473 iset = task_list%tasks(itask)%iset
474 jset = task_list%tasks(itask)%jset
475 ipgf = task_list%tasks(itask)%ipgf
476 jpgf = task_list%tasks(itask)%jpgf
477 IF (blocks(iatom, jatom, img) == -1 .AND. blocks(jatom, iatom, img) == -1) THEN
478 blocks(iatom, jatom, img) = 1
479 blocks(jatom, iatom, img) = 1
480 ELSE
481 WRITE (6, *) "TASK LIST CONFLICT IN PAIR ", ipair
482 WRITE (6, *) "Reuse of iatom, jatom, image ", iatom, jatom, img
483 END IF
484
485 iatom_old = iatom
486 jatom_old = jatom
487 img_old = img
488 DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
489 ilevel = task_list%tasks(itask)%grid_level
490 img = task_list%tasks(itask)%image
491 iatom = task_list%tasks(itask)%iatom
492 jatom = task_list%tasks(itask)%jatom
493 iset = task_list%tasks(itask)%iset
494 jset = task_list%tasks(itask)%jset
495 ipgf = task_list%tasks(itask)%ipgf
496 jpgf = task_list%tasks(itask)%jpgf
497 IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
498 WRITE (6, *) "TASK LIST CONFLICT IN TASK ", itask
499 WRITE (6, *) "Inconsistent iatom, jatom, image ", iatom, jatom, img
500 WRITE (6, *) "Should be iatom, jatom, image ", iatom_old, jatom_old, img_old
501 END IF
502
503 END DO
504 END DO
505 DEALLOCATE (blocks)
506
507 END DO
508
509 END IF
510
511 CALL timestop(handle)
512
513 END SUBROUTINE generate_qs_task_list
514
515! **************************************************************************************************
516!> \brief Sends the task list data to the grid API.
517!> \author Ole Schuett
518! **************************************************************************************************
519 SUBROUTINE create_grid_task_list(task_list, qs_kind_set, particle_set, cell, basis_type, rs_grids)
520 TYPE(task_list_type), POINTER :: task_list
521 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
522 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
523 TYPE(cell_type), POINTER :: cell
524 CHARACTER(LEN=*) :: basis_type
525 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
526
527 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
528 INTEGER :: nset, natoms, nkinds, ntasks, &
529 ikind, iatom, itask, nsgf
530 INTEGER, DIMENSION(:), ALLOCATABLE :: atom_kinds, level_list, iatom_list, jatom_list, &
531 iset_list, jset_list, ipgf_list, jpgf_list, &
532 border_mask_list, block_num_list
533 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: radius_list
534 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: rab_list, atom_positions
535 TYPE(task_type), DIMENSION(:), POINTER :: tasks
536 INTEGER, DIMENSION(:, :), POINTER :: first_sgf
537 REAL(kind=dp), DIMENSION(:, :), POINTER :: sphi, zet
538 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
539
540 nkinds = SIZE(qs_kind_set)
541 natoms = SIZE(particle_set)
542 ntasks = task_list%ntasks
543 tasks => task_list%tasks
544
545 IF (.NOT. ASSOCIATED(task_list%grid_basis_sets)) THEN
546 ! Basis sets do not change during simulation - only need to create them once.
547 ALLOCATE (task_list%grid_basis_sets(nkinds))
548 DO ikind = 1, nkinds
549 CALL get_qs_kind(qs_kind_set(ikind), basis_type=basis_type, basis_set=orb_basis_set)
550 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
551 nset=nset, &
552 nsgf=nsgf, &
553 nsgf_set=nsgf_set, &
554 npgf=npgf, &
555 first_sgf=first_sgf, &
556 lmax=lmax, &
557 lmin=lmin, &
558 sphi=sphi, &
559 zet=zet)
560 CALL grid_create_basis_set(nset=nset, &
561 nsgf=nsgf, &
562 maxco=SIZE(sphi, 1), &
563 maxpgf=SIZE(zet, 1), &
564 lmin=lmin, &
565 lmax=lmax, &
566 npgf=npgf, &
567 nsgf_set=nsgf_set, &
568 first_sgf=first_sgf, &
569 sphi=sphi, &
570 zet=zet, &
571 basis_set=task_list%grid_basis_sets(ikind))
572 END DO
573 END IF
574
575 ! Pack task list infos
576 ALLOCATE (atom_kinds(natoms), atom_positions(3, natoms))
577 DO iatom = 1, natoms
578 atom_kinds(iatom) = particle_set(iatom)%atomic_kind%kind_number
579 atom_positions(:, iatom) = pbc(particle_set(iatom)%r, cell)
580 END DO
581
582 ALLOCATE (level_list(ntasks), iatom_list(ntasks), jatom_list(ntasks))
583 ALLOCATE (iset_list(ntasks), jset_list(ntasks), ipgf_list(ntasks), jpgf_list(ntasks))
584 ALLOCATE (border_mask_list(ntasks), block_num_list(ntasks))
585 ALLOCATE (radius_list(ntasks), rab_list(3, ntasks))
586
587 DO itask = 1, ntasks
588 level_list(itask) = tasks(itask)%grid_level
589 iatom_list(itask) = tasks(itask)%iatom
590 jatom_list(itask) = tasks(itask)%jatom
591 iset_list(itask) = tasks(itask)%iset
592 jset_list(itask) = tasks(itask)%jset
593 ipgf_list(itask) = tasks(itask)%ipgf
594 jpgf_list(itask) = tasks(itask)%jpgf
595 IF (tasks(itask)%dist_type == 2) THEN
596 border_mask_list(itask) = iand(63, not(tasks(itask)%subpatch_pattern)) ! invert last 6 bits
597 ELSE
598 border_mask_list(itask) = 0 ! no masking
599 END IF
600 block_num_list(itask) = tasks(itask)%pair_index ! change of nomenclature pair_index -> block_num
601 radius_list(itask) = tasks(itask)%radius
602 rab_list(:, itask) = tasks(itask)%rab(:)
603 END DO
604
605 CALL grid_create_task_list(ntasks=ntasks, &
606 natoms=natoms, &
607 nkinds=nkinds, &
608 nblocks=SIZE(task_list%pair_offsets_recv), &
609 block_offsets=task_list%pair_offsets_recv, &
610 atom_positions=atom_positions, &
611 atom_kinds=atom_kinds, &
612 basis_sets=task_list%grid_basis_sets, &
613 level_list=level_list, &
614 iatom_list=iatom_list, &
615 jatom_list=jatom_list, &
616 iset_list=iset_list, &
617 jset_list=jset_list, &
618 ipgf_list=ipgf_list, &
619 jpgf_list=jpgf_list, &
620 border_mask_list=border_mask_list, &
621 block_num_list=block_num_list, &
622 radius_list=radius_list, &
623 rab_list=rab_list, &
624 rs_grids=rs_grids, &
625 task_list=task_list%grid_task_list)
626
627 CALL offload_create_buffer(task_list%buffer_size_recv, task_list%pab_buffer)
628 CALL offload_create_buffer(task_list%buffer_size_recv, task_list%hab_buffer)
629
630 END SUBROUTINE create_grid_task_list
631
632! **************************************************************************************************
633!> \brief ...
634!> \param tasks ...
635!> \param ntasks ...
636!> \param curr_tasks ...
637!> \param rs_descs ...
638!> \param dft_control ...
639!> \param cube_info ...
640!> \param gridlevel_info ...
641!> \param cindex ...
642!> \param iatom ...
643!> \param jatom ...
644!> \param rpgfa ...
645!> \param rpgfb ...
646!> \param zeta ...
647!> \param zetb ...
648!> \param kind_radius_b ...
649!> \param set_radius_a ...
650!> \param set_radius_b ...
651!> \param ra ...
652!> \param rab ...
653!> \param la_max ...
654!> \param la_min ...
655!> \param lb_max ...
656!> \param lb_min ...
657!> \param npgfa ...
658!> \param npgfb ...
659!> \param nseta ...
660!> \param nsetb ...
661!> \par History
662!> Joost VandeVondele: 10.2008 refactored
663! **************************************************************************************************
664 SUBROUTINE task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, dft_control, &
665 cube_info, gridlevel_info, cindex, &
666 iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, set_radius_a, set_radius_b, ra, rab, &
667 la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
668
669 TYPE(task_type), DIMENSION(:), POINTER :: tasks
670 INTEGER :: ntasks, curr_tasks
671 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
672 POINTER :: rs_descs
673 TYPE(dft_control_type), POINTER :: dft_control
674 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
675 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
676 INTEGER :: cindex, iatom, jatom
677 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb
678 REAL(kind=dp) :: kind_radius_b
679 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
680 REAL(kind=dp), DIMENSION(3) :: ra, rab
681 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
682 npgfb
683 INTEGER :: nseta, nsetb
684
685 INTEGER :: cube_center(3), igrid_level, ipgf, iset, &
686 jpgf, jset, lb_cube(3), ub_cube(3)
687 REAL(kind=dp) :: dab, rab2, radius, zetp
688
689 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
690 dab = sqrt(rab2)
691
692 loop_iset: DO iset = 1, nseta
693
694 IF (set_radius_a(iset) + kind_radius_b < dab) cycle loop_iset
695
696 loop_jset: DO jset = 1, nsetb
697
698 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle loop_jset
699
700 loop_ipgf: DO ipgf = 1, npgfa(iset)
701
702 IF (rpgfa(ipgf, iset) + set_radius_b(jset) < dab) cycle loop_ipgf
703
704 loop_jpgf: DO jpgf = 1, npgfb(jset)
705
706 IF (rpgfa(ipgf, iset) + rpgfb(jpgf, jset) < dab) cycle loop_jpgf
707
708 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
709 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
710
711 CALL compute_pgf_properties(cube_center, lb_cube, ub_cube, radius, &
712 rs_descs(igrid_level)%rs_desc, cube_info(igrid_level), &
713 la_max(iset), zeta(ipgf, iset), la_min(iset), &
714 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
715 ra, rab, rab2, dft_control%qs_control%eps_rho_rspace)
716
717 CALL pgf_to_tasks(tasks, ntasks, curr_tasks, &
718 rab, cindex, iatom, jatom, iset, jset, ipgf, jpgf, &
719 la_max(iset), lb_max(jset), rs_descs(igrid_level)%rs_desc, &
720 igrid_level, gridlevel_info%ngrid_levels, cube_center, &
721 lb_cube, ub_cube, radius)
722
723 END DO loop_jpgf
724
725 END DO loop_ipgf
726
727 END DO loop_jset
728
729 END DO loop_iset
730
731 END SUBROUTINE task_list_inner_loop
732
733! **************************************************************************************************
734!> \brief combines the calculation of several basic properties of a given pgf:
735!> its center, the bounding cube, the radius, the cost,
736!> tries to predict the time needed for processing this task
737!> in this way an improved load balance might be obtained
738!> \param cube_center ...
739!> \param lb_cube ...
740!> \param ub_cube ...
741!> \param radius ...
742!> \param rs_desc ...
743!> \param cube_info ...
744!> \param la_max ...
745!> \param zeta ...
746!> \param la_min ...
747!> \param lb_max ...
748!> \param zetb ...
749!> \param lb_min ...
750!> \param ra ...
751!> \param rab ...
752!> \param rab2 ...
753!> \param eps ...
754!> \par History
755!> 10.2008 refactored [Joost VandeVondele]
756!> \note
757!> -) this requires the radius to be computed in the same way as
758!> collocate_pgf_product, we should factor that part into a subroutine
759!> -) we're assuming that integrate_pgf and collocate_pgf are the same cost for load balancing
760!> this is more or less true for map_consistent
761!> -) in principle, the computed radius could be recycled in integrate_pgf/collocate_pgf if it is certainly
762!> the same, this could lead to a small speedup
763!> -) the cost function is a fit through the median cost of mapping a pgf with a given l and a given radius (in grid points)
764!> fitting the measured data on an opteron/g95 using the expression
765!> a*(l+b)(r+c)**3+d which is based on the innerloop of the collocating routines
766! **************************************************************************************************
767 SUBROUTINE compute_pgf_properties(cube_center, lb_cube, ub_cube, radius, &
768 rs_desc, cube_info, la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, rab2, eps)
769
770 INTEGER, DIMENSION(3), INTENT(OUT) :: cube_center, lb_cube, ub_cube
771 REAL(kind=dp), INTENT(OUT) :: radius
772 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
773 TYPE(cube_info_type), INTENT(IN) :: cube_info
774 INTEGER, INTENT(IN) :: la_max
775 REAL(kind=dp), INTENT(IN) :: zeta
776 INTEGER, INTENT(IN) :: la_min, lb_max
777 REAL(kind=dp), INTENT(IN) :: zetb
778 INTEGER, INTENT(IN) :: lb_min
779 REAL(kind=dp), INTENT(IN) :: ra(3), rab(3), rab2, eps
780
781 INTEGER :: extent(3)
782 INTEGER, DIMENSION(:), POINTER :: sphere_bounds
783 REAL(kind=dp) :: cutoff, f, prefactor, rb(3), zetp
784 REAL(kind=dp), DIMENSION(3) :: rp
785
786! the radius for this task
787
788 zetp = zeta + zetb
789 rp(:) = ra(:) + zetb/zetp*rab(:)
790 rb(:) = ra(:) + rab(:)
791 cutoff = 1.0_dp
792 f = zetb/zetp
793 prefactor = exp(-zeta*f*rab2)
794 radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, &
795 zetp=zetp, eps=eps, prefactor=prefactor, cutoff=cutoff)
796
797 CALL compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab)
798 ! compute cube_center, the center of the gaussian product to map (folded to within the unit cell)
799 cube_center(:) = modulo(cube_center(:), rs_desc%npts(:))
800 cube_center(:) = cube_center(:) + rs_desc%lb(:)
801
802 IF (rs_desc%orthorhombic) THEN
803 CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
804 ELSE
805 CALL return_cube_nonortho(cube_info, radius, lb_cube, ub_cube, rp)
806 !? unclear if extent is computed correctly.
807 extent(:) = ub_cube(:) - lb_cube(:)
808 lb_cube(:) = -extent(:)/2 - 1
809 ub_cube(:) = extent(:)/2
810 END IF
811
812 END SUBROUTINE compute_pgf_properties
813! **************************************************************************************************
814!> \brief predicts the cost of a task in kcycles for a given task
815!> the model is based on a fit of actual data, and might need updating
816!> as collocate_pgf_product changes (or CPUs/compilers change)
817!> maybe some dynamic approach, improving the cost model on the fly could
818!> work as well
819!> the cost model does not yet take into account the fraction of space
820!> that is mapped locally for a given cube and rs_grid (generalised tasks)
821!> \param lb_cube ...
822!> \param ub_cube ...
823!> \param fraction ...
824!> \param lmax ...
825!> \param is_ortho ...
826!> \return ...
827! **************************************************************************************************
828 INTEGER FUNCTION cost_model(lb_cube, ub_cube, fraction, lmax, is_ortho)
829 INTEGER, DIMENSION(3), INTENT(IN) :: lb_cube, ub_cube
830 REAL(kind=dp), INTENT(IN) :: fraction
831 INTEGER :: lmax
832 LOGICAL :: is_ortho
833
834 INTEGER :: cmax
835 REAL(kind=dp) :: v1, v2, v3, v4, v5
836
837 cmax = maxval(((ub_cube - lb_cube) + 1)/2)
838
839 IF (is_ortho) THEN
840 v1 = 1.504760e+00_dp
841 v2 = 3.126770e+00_dp
842 v3 = 5.074106e+00_dp
843 v4 = 1.091568e+00_dp
844 v5 = 1.070187e+00_dp
845 ELSE
846 v1 = 7.831105e+00_dp
847 v2 = 2.675174e+00_dp
848 v3 = 7.546553e+00_dp
849 v4 = 6.122446e-01_dp
850 v5 = 3.886382e+00_dp
851 END IF
852 cost_model = ceiling(((lmax + v1)*(cmax + v2)**3*v3*fraction + v4 + v5*lmax**7)/1000.0_dp)
853
854 END FUNCTION cost_model
855! **************************************************************************************************
856!> \brief pgf_to_tasks converts a given pgf to one or more tasks, in particular
857!> this determines by which CPUs a given pgf gets collocated
858!> the format of the task array is as follows
859!> tasks(1,i) := destination
860!> tasks(2,i) := source
861!> tasks(3,i) := compressed type (iatom, jatom, ....)
862!> tasks(4,i) := type (0: replicated, 1: distributed local, 2: distributed generalised)
863!> tasks(5,i) := cost
864!> tasks(6,i) := alternate destination code (0 if none available)
865!>
866!> \param tasks ...
867!> \param ntasks ...
868!> \param curr_tasks ...
869!> \param rab ...
870!> \param cindex ...
871!> \param iatom ...
872!> \param jatom ...
873!> \param iset ...
874!> \param jset ...
875!> \param ipgf ...
876!> \param jpgf ...
877!> \param la_max ...
878!> \param lb_max ...
879!> \param rs_desc ...
880!> \param igrid_level ...
881!> \param n_levels ...
882!> \param cube_center ...
883!> \param lb_cube ...
884!> \param ub_cube ...
885!> \par History
886!> 10.2008 Refactored based on earlier routines by MattW [Joost VandeVondele]
887! **************************************************************************************************
888 SUBROUTINE pgf_to_tasks(tasks, ntasks, curr_tasks, &
889 rab, cindex, iatom, jatom, iset, jset, ipgf, jpgf, &
890 la_max, lb_max, rs_desc, igrid_level, n_levels, &
891 cube_center, lb_cube, ub_cube, radius)
892
893 TYPE(task_type), DIMENSION(:), POINTER :: tasks
894 INTEGER, INTENT(INOUT) :: ntasks, curr_tasks
895 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
896 INTEGER, INTENT(IN) :: cindex, iatom, jatom, iset, jset, ipgf, &
897 jpgf, la_max, lb_max
898 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
899 INTEGER, INTENT(IN) :: igrid_level, n_levels
900 INTEGER, DIMENSION(3), INTENT(IN) :: cube_center, lb_cube, ub_cube
901 REAL(kind=dp), INTENT(IN) :: radius
902
903 INTEGER, PARAMETER :: add_tasks = 1000
904 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp
905
906 INTEGER :: added_tasks, cost, j, lmax
907 LOGICAL :: is_ortho
908 REAL(kind=dp) :: tfraction
909
910!$OMP SINGLE
911 ntasks = ntasks + 1
912 IF (ntasks > curr_tasks) THEN
913 curr_tasks = int((curr_tasks + add_tasks)*mult_tasks)
914 CALL reallocate_tasks(tasks, curr_tasks)
915 END IF
916!$OMP END SINGLE
917
918 IF (rs_desc%distributed) THEN
919
920 ! finds the node(s) that need to process this task
921 ! on exit tasks(:)%dist_type is 1 for distributed tasks and 2 for generalised tasks
922 CALL rs_find_node(rs_desc, igrid_level, n_levels, cube_center, &
923 ntasks=ntasks, tasks=tasks, lb_cube=lb_cube, ub_cube=ub_cube, added_tasks=added_tasks)
924
925 ELSE
926 tasks(ntasks)%destination = encode_rank(rs_desc%my_pos, igrid_level, n_levels)
927 tasks(ntasks)%dist_type = 0
928 tasks(ntasks)%subpatch_pattern = 0
929 added_tasks = 1
930 END IF
931
932 lmax = la_max + lb_max
933 is_ortho = (tasks(ntasks)%dist_type == 0 .OR. tasks(ntasks)%dist_type == 1) .AND. rs_desc%orthorhombic
934 ! we assume the load is shared equally between processes dealing with a generalised Gaussian.
935 ! this could be refined in the future
936 tfraction = 1.0_dp/added_tasks
937
938 cost = cost_model(lb_cube, ub_cube, tfraction, lmax, is_ortho)
939
940 DO j = 1, added_tasks
941 tasks(ntasks - added_tasks + j)%source = encode_rank(rs_desc%my_pos, igrid_level, n_levels)
942 tasks(ntasks - added_tasks + j)%cost = cost
943 tasks(ntasks - added_tasks + j)%grid_level = igrid_level
944 tasks(ntasks - added_tasks + j)%image = cindex
945 tasks(ntasks - added_tasks + j)%iatom = iatom
946 tasks(ntasks - added_tasks + j)%jatom = jatom
947 tasks(ntasks - added_tasks + j)%iset = iset
948 tasks(ntasks - added_tasks + j)%jset = jset
949 tasks(ntasks - added_tasks + j)%ipgf = ipgf
950 tasks(ntasks - added_tasks + j)%jpgf = jpgf
951 tasks(ntasks - added_tasks + j)%rab = rab
952 tasks(ntasks - added_tasks + j)%radius = radius
953 END DO
954
955 END SUBROUTINE pgf_to_tasks
956
957! **************************************************************************************************
958!> \brief performs load balancing of the tasks on the distributed grids
959!> \param tasks ...
960!> \param ntasks ...
961!> \param rs_descs ...
962!> \param grid_level ...
963!> \param natoms ...
964!> \par History
965!> created 2008-10-03 [Joost VandeVondele]
966! **************************************************************************************************
967 SUBROUTINE load_balance_distributed(tasks, ntasks, rs_descs, grid_level, natoms)
968
969 TYPE(task_type), DIMENSION(:), POINTER :: tasks
970 INTEGER :: ntasks
971 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
972 POINTER :: rs_descs
973 INTEGER :: grid_level, natoms
974
975 CHARACTER(LEN=*), PARAMETER :: routinen = 'load_balance_distributed'
976
977 INTEGER :: handle
978 INTEGER, DIMENSION(:, :, :), POINTER :: list
979
980 CALL timeset(routinen, handle)
981
982 NULLIFY (list)
983 ! here we create for each cpu (0:ncpu-1) a list of possible destinations.
984 ! if a destination would not be in this list, it is a bug
985 CALL create_destination_list(list, rs_descs, grid_level)
986
987 ! now, walk over the tasks, filling in the loads of each destination
988 CALL compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list=.true.)
989
990 ! optimize loads & fluxes
991 CALL optimize_load_list(list, rs_descs(1)%rs_desc%group, rs_descs(1)%rs_desc%my_pos)
992
993 ! now, walk over the tasks, using the list to set the destinations
994 CALL compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list=.false.)
995
996 DEALLOCATE (list)
997
998 CALL timestop(handle)
999
1000 END SUBROUTINE load_balance_distributed
1001
1002! **************************************************************************************************
1003!> \brief this serial routine adjusts the fluxes in the global list
1004!>
1005!> \param list_global ...
1006!> \par History
1007!> created 2008-10-06 [Joost VandeVondele]
1008! **************************************************************************************************
1009 SUBROUTINE balance_global_list(list_global)
1010 INTEGER, DIMENSION(:, :, 0:) :: list_global
1011
1012 CHARACTER(LEN=*), PARAMETER :: routinen = 'balance_global_list'
1013 INTEGER, PARAMETER :: max_iter = 100
1014 REAL(kind=dp), PARAMETER :: tolerance_factor = 0.005_dp
1015
1016 INTEGER :: dest, handle, icpu, idest, iflux, &
1017 ilocal, k, maxdest, ncpu, nflux
1018 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: flux_connections
1019 LOGICAL :: solution_optimal
1020 REAL(kind=dp) :: average, load_shift, max_load_shift, &
1021 tolerance
1022 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: load, optimized_flux, optimized_load
1023 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: flux_limits
1024
1025 CALL timeset(routinen, handle)
1026
1027 ncpu = SIZE(list_global, 3)
1028 maxdest = SIZE(list_global, 2)
1029 ALLOCATE (load(0:ncpu - 1))
1030 load = 0.0_dp
1031 ALLOCATE (optimized_load(0:ncpu - 1))
1032
1033 ! figure out the number of fluxes
1034 ! we assume that the global_list is symmetric
1035 nflux = 0
1036 DO icpu = 0, ncpu - 1
1037 DO idest = 1, maxdest
1038 dest = list_global(1, idest, icpu)
1039 IF (dest < ncpu .AND. dest > icpu) nflux = nflux + 1
1040 END DO
1041 END DO
1042 ALLOCATE (optimized_flux(nflux))
1043 ALLOCATE (flux_limits(2, nflux))
1044 ALLOCATE (flux_connections(2, nflux))
1045
1046 ! reorder data
1047 flux_limits = 0
1048 nflux = 0
1049 DO icpu = 0, ncpu - 1
1050 load(icpu) = sum(list_global(2, :, icpu))
1051 DO idest = 1, maxdest
1052 dest = list_global(1, idest, icpu)
1053 IF (dest < ncpu) THEN
1054 IF (dest /= icpu) THEN
1055 IF (dest > icpu) THEN
1056 nflux = nflux + 1
1057 flux_limits(2, nflux) = list_global(2, idest, icpu)
1058 flux_connections(1, nflux) = icpu
1059 flux_connections(2, nflux) = dest
1060 ELSE
1061 DO iflux = 1, nflux
1062 IF (flux_connections(1, iflux) == dest .AND. flux_connections(2, iflux) == icpu) THEN
1063 flux_limits(1, iflux) = -list_global(2, idest, icpu)
1064 EXIT
1065 END IF
1066 END DO
1067 END IF
1068 END IF
1069 END IF
1070 END DO
1071 END DO
1072
1073 solution_optimal = .false.
1074 optimized_flux = 0.0_dp
1075
1076 ! an iterative solver, if iterated till convergence the maximum load is minimal
1077 ! we terminate before things are fully converged, since this does show up in the timings
1078 ! once the largest shift becomes less than a small fraction of the average load, we're done
1079 ! we're perfectly happy if the load balance is within 1 percent or so
1080 ! the maximum load normally converges even faster
1081 average = sum(load)/SIZE(load)
1082 tolerance = tolerance_factor*average
1083
1084 optimized_load(:) = load
1085 DO k = 1, max_iter
1086 max_load_shift = 0.0_dp
1087 DO iflux = 1, nflux
1088 load_shift = (optimized_load(flux_connections(1, iflux)) - optimized_load(flux_connections(2, iflux)))/2
1089 load_shift = max(flux_limits(1, iflux) - optimized_flux(iflux), load_shift)
1090 load_shift = min(flux_limits(2, iflux) - optimized_flux(iflux), load_shift)
1091 max_load_shift = max(abs(load_shift), max_load_shift)
1092 optimized_load(flux_connections(1, iflux)) = optimized_load(flux_connections(1, iflux)) - load_shift
1093 optimized_load(flux_connections(2, iflux)) = optimized_load(flux_connections(2, iflux)) + load_shift
1094 optimized_flux(iflux) = optimized_flux(iflux) + load_shift
1095 END DO
1096 IF (max_load_shift < tolerance) THEN
1097 solution_optimal = .true.
1098 EXIT
1099 END IF
1100 END DO
1101
1102 ! now adjust the load list to reflect the optimized fluxes
1103 ! reorder data
1104 nflux = 0
1105 DO icpu = 0, ncpu - 1
1106 DO idest = 1, maxdest
1107 IF (list_global(1, idest, icpu) == icpu) ilocal = idest
1108 END DO
1109 DO idest = 1, maxdest
1110 dest = list_global(1, idest, icpu)
1111 IF (dest < ncpu) THEN
1112 IF (dest /= icpu) THEN
1113 IF (dest > icpu) THEN
1114 nflux = nflux + 1
1115 IF (optimized_flux(nflux) > 0) THEN
1116 list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1117 list_global(2, idest, icpu) - nint(optimized_flux(nflux))
1118 list_global(2, idest, icpu) = nint(optimized_flux(nflux))
1119 ELSE
1120 list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1121 list_global(2, idest, icpu)
1122 list_global(2, idest, icpu) = 0
1123 END IF
1124 ELSE
1125 DO iflux = 1, nflux
1126 IF (flux_connections(1, iflux) == dest .AND. flux_connections(2, iflux) == icpu) THEN
1127 IF (optimized_flux(iflux) > 0) THEN
1128 list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1129 list_global(2, idest, icpu)
1130 list_global(2, idest, icpu) = 0
1131 ELSE
1132 list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1133 list_global(2, idest, icpu) + nint(optimized_flux(iflux))
1134 list_global(2, idest, icpu) = -nint(optimized_flux(iflux))
1135 END IF
1136 EXIT
1137 END IF
1138 END DO
1139 END IF
1140 END IF
1141 END IF
1142 END DO
1143 END DO
1144
1145 CALL timestop(handle)
1146
1147 END SUBROUTINE balance_global_list
1148
1149! **************************************************************************************************
1150!> \brief this routine gets back optimized loads for all destinations
1151!>
1152!> \param list ...
1153!> \param group ...
1154!> \param my_pos ...
1155!> \par History
1156!> created 2008-10-06 [Joost VandeVondele]
1157!> Modified 2016-01 [EPCC] Reduce memory requirements on P processes
1158!> from O(P^2) to O(P)
1159! **************************************************************************************************
1160 SUBROUTINE optimize_load_list(list, group, my_pos)
1161 INTEGER, DIMENSION(:, :, 0:) :: list
1162 TYPE(mp_comm_type), INTENT(IN) :: group
1163 INTEGER, INTENT(IN) :: my_pos
1164
1165 CHARACTER(LEN=*), PARAMETER :: routinen = 'optimize_load_list'
1166 INTEGER, PARAMETER :: rank_of_root = 0
1167
1168 INTEGER :: handle, icpu, idest, maxdest, ncpu
1169 INTEGER, ALLOCATABLE, DIMENSION(:) :: load_all
1170 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: load_partial
1171 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: list_global
1172
1173 CALL timeset(routinen, handle)
1174
1175 ncpu = SIZE(list, 3)
1176 maxdest = SIZE(list, 2)
1177
1178 !find total workload ...
1179 ALLOCATE (load_all(maxdest*ncpu))
1180 load_all(:) = reshape(list(2, :, :), [maxdest*ncpu])
1181 CALL group%sum(load_all(:), rank_of_root)
1182
1183 ! ... and optimise the work per process
1184 ALLOCATE (list_global(2, maxdest, ncpu))
1185 IF (rank_of_root == my_pos) THEN
1186 list_global(1, :, :) = list(1, :, :)
1187 list_global(2, :, :) = reshape(load_all, [maxdest, ncpu])
1188 CALL balance_global_list(list_global)
1189 END IF
1190 CALL group%bcast(list_global, rank_of_root)
1191
1192 !figure out how much can be sent to other processes
1193 ALLOCATE (load_partial(maxdest, ncpu))
1194 ! send 'load_all', which is a copy of 'list' (but without leading dimension/stride)
1195 CALL group%sum_partial(reshape(load_all, [maxdest, ncpu]), load_partial(:, :))
1196
1197 DO icpu = 1, ncpu
1198 DO idest = 1, maxdest
1199
1200 !need to deduct 1 because `list' was passed in to this routine as being indexed from zero
1201 IF (load_partial(idest, icpu) > list_global(2, idest, icpu)) THEN
1202 IF (load_partial(idest, icpu) - list(2, idest, icpu - 1) < list_global(2, idest, icpu)) THEN
1203 list(2, idest, icpu - 1) = list_global(2, idest, icpu) &
1204 - (load_partial(idest, icpu) - list(2, idest, icpu - 1))
1205 ELSE
1206 list(2, idest, icpu - 1) = 0
1207 END IF
1208 END IF
1209
1210 END DO
1211 END DO
1212
1213 !clean up before leaving
1214 DEALLOCATE (load_all)
1215 DEALLOCATE (list_global)
1216 DEALLOCATE (load_partial)
1217
1218 CALL timestop(handle)
1219 END SUBROUTINE optimize_load_list
1220
1221! **************************************************************************************************
1222!> \brief fill the load list with values derived from the tasks array
1223!> from the alternate locations, we select the alternate location that
1224!> can be used without increasing the number of matrix blocks needed to
1225!> distribute.
1226!> Replicated tasks are not yet considered
1227!>
1228!> \param list ...
1229!> \param rs_descs ...
1230!> \param grid_level ...
1231!> \param tasks ...
1232!> \param ntasks ...
1233!> \param natoms ...
1234!> \param create_list ...
1235!> \par History
1236!> created 2008-10-06 [Joost VandeVondele]
1237! **************************************************************************************************
1238 SUBROUTINE compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list)
1239 INTEGER, DIMENSION(:, :, 0:) :: list
1240 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1241 POINTER :: rs_descs
1242 INTEGER :: grid_level
1243 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1244 INTEGER :: ntasks, natoms
1245 LOGICAL :: create_list
1246
1247 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_load_list'
1248
1249 INTEGER :: cost, dest, handle, i, iatom, ilevel, img, img_old, iopt, ipgf, iset, itask, &
1250 itask_start, itask_stop, jatom, jpgf, jset, li, maxdest, ncpu, ndest_pair, nopt, nshort, &
1251 rank
1252 INTEGER(KIND=int_8) :: bit_pattern, ipair, ipair_old, natom8
1253 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: loads
1254 INTEGER, ALLOCATABLE, DIMENSION(:) :: all_dests, index
1255 INTEGER, DIMENSION(6) :: options
1256
1257 CALL timeset(routinen, handle)
1258
1259 ALLOCATE (loads(0:rs_descs(grid_level)%rs_desc%group_size - 1))
1260 CALL get_current_loads(loads, rs_descs, grid_level, ntasks, tasks, use_reordered_ranks=.false.)
1261
1262 maxdest = SIZE(list, 2)
1263 ncpu = SIZE(list, 3)
1264 natom8 = natoms
1265
1266 ! first find the tasks that deal with the same atom pair
1267 itask_stop = 0
1268 ipair_old = huge(ipair_old)
1269 img_old = -1
1270 ALLOCATE (all_dests(0))
1271 ALLOCATE (index(0))
1272
1273 DO
1274
1275 ! first find the range of tasks that deal with the same atom pair
1276 itask_start = itask_stop + 1
1277 itask_stop = itask_start
1278 IF (itask_stop > ntasks) EXIT
1279 ilevel = tasks(itask_stop)%grid_level
1280 img_old = tasks(itask_stop)%image
1281 iatom = tasks(itask_stop)%iatom
1282 jatom = tasks(itask_stop)%jatom
1283 iset = tasks(itask_stop)%iset
1284 jset = tasks(itask_stop)%jset
1285 ipgf = tasks(itask_stop)%ipgf
1286 jpgf = tasks(itask_stop)%jpgf
1287
1288 ipair_old = (iatom - 1)*natom8 + (jatom - 1)
1289 DO
1290 IF (itask_stop + 1 > ntasks) EXIT
1291 ilevel = tasks(itask_stop + 1)%grid_level
1292 img = tasks(itask_stop + 1)%image
1293 iatom = tasks(itask_stop + 1)%iatom
1294 jatom = tasks(itask_stop + 1)%jatom
1295 iset = tasks(itask_stop + 1)%iset
1296 jset = tasks(itask_stop + 1)%jset
1297 ipgf = tasks(itask_stop + 1)%ipgf
1298 jpgf = tasks(itask_stop + 1)%jpgf
1299
1300 ipair = (iatom - 1)*natom8 + (jatom - 1)
1301 IF (ipair == ipair_old .AND. img == img_old) THEN
1302 itask_stop = itask_stop + 1
1303 ELSE
1304 EXIT
1305 END IF
1306 END DO
1307 ipair = ipair_old
1308 nshort = itask_stop - itask_start + 1
1309
1310 ! find the unique list of destinations on this grid level only
1311 DEALLOCATE (all_dests)
1312 ALLOCATE (all_dests(nshort))
1313 DEALLOCATE (index)
1314 ALLOCATE (index(nshort))
1315 DO i = 1, nshort
1316 ilevel = tasks(itask_start + i - 1)%grid_level
1317 img = tasks(itask_start + i - 1)%image
1318 iatom = tasks(itask_start + i - 1)%iatom
1319 jatom = tasks(itask_start + i - 1)%jatom
1320 iset = tasks(itask_start + i - 1)%iset
1321 jset = tasks(itask_start + i - 1)%jset
1322 ipgf = tasks(itask_start + i - 1)%ipgf
1323 jpgf = tasks(itask_start + i - 1)%jpgf
1324
1325 IF (ilevel == grid_level) THEN
1326 all_dests(i) = decode_rank(tasks(itask_start + i - 1)%destination, SIZE(rs_descs))
1327 ELSE
1328 all_dests(i) = huge(all_dests(i))
1329 END IF
1330 END DO
1331 CALL sort(all_dests, nshort, index)
1332 ndest_pair = 1
1333 DO i = 2, nshort
1334 IF ((all_dests(ndest_pair) /= all_dests(i)) .AND. (all_dests(i) /= huge(all_dests(i)))) THEN
1335 ndest_pair = ndest_pair + 1
1336 all_dests(ndest_pair) = all_dests(i)
1337 END IF
1338 END DO
1339
1340 DO itask = itask_start, itask_stop
1341
1342 dest = decode_rank(tasks(itask)%destination, SIZE(rs_descs)) ! notice that dest can be changed
1343 ilevel = tasks(itask)%grid_level
1344 img = tasks(itask)%image
1345 iatom = tasks(itask)%iatom
1346 jatom = tasks(itask)%jatom
1347 iset = tasks(itask)%iset
1348 jset = tasks(itask)%jset
1349 ipgf = tasks(itask)%ipgf
1350 jpgf = tasks(itask)%jpgf
1351
1352 ! Only proceed with tasks which are on this grid level
1353 IF (ilevel /= grid_level) cycle
1354 ipair = (iatom - 1)*natom8 + (jatom - 1)
1355 cost = int(tasks(itask)%cost)
1356
1357 SELECT CASE (tasks(itask)%dist_type)
1358 CASE (1)
1359 bit_pattern = tasks(itask)%subpatch_pattern
1360 nopt = 0
1361 IF (btest(bit_pattern, 0)) THEN
1362 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [-1, 0, 0])
1363 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1364 nopt = nopt + 1
1365 options(nopt) = rank
1366 END IF
1367 END IF
1368 IF (btest(bit_pattern, 1)) THEN
1369 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [+1, 0, 0])
1370 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1371 nopt = nopt + 1
1372 options(nopt) = rank
1373 END IF
1374 END IF
1375 IF (btest(bit_pattern, 2)) THEN
1376 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, -1, 0])
1377 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1378 nopt = nopt + 1
1379 options(nopt) = rank
1380 END IF
1381 END IF
1382 IF (btest(bit_pattern, 3)) THEN
1383 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, +1, 0])
1384 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1385 nopt = nopt + 1
1386 options(nopt) = rank
1387 END IF
1388 END IF
1389 IF (btest(bit_pattern, 4)) THEN
1390 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, 0, -1])
1391 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1392 nopt = nopt + 1
1393 options(nopt) = rank
1394 END IF
1395 END IF
1396 IF (btest(bit_pattern, 5)) THEN
1397 rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, 0, +1])
1398 IF (any(all_dests(1:ndest_pair) == rank)) THEN
1399 nopt = nopt + 1
1400 options(nopt) = rank
1401 END IF
1402 END IF
1403 IF (nopt > 0) THEN
1404 ! set it to the rank with the lowest load
1405 rank = options(1)
1406 DO iopt = 2, nopt
1407 IF (loads(rank) > loads(options(iopt))) rank = options(iopt)
1408 END DO
1409 ELSE
1410 rank = dest
1411 END IF
1412 li = list_index(list, rank, dest)
1413 IF (create_list) THEN
1414 list(2, li, dest) = list(2, li, dest) + cost
1415 ELSE
1416 IF (list(1, li, dest) == dest) THEN
1417 tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1418 ELSE
1419 IF (list(2, li, dest) >= cost) THEN
1420 list(2, li, dest) = list(2, li, dest) - cost
1421 tasks(itask)%destination = encode_rank(list(1, li, dest), ilevel, SIZE(rs_descs))
1422 ELSE
1423 tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1424 END IF
1425 END IF
1426 END IF
1427 CASE (2) ! generalised
1428 li = list_index(list, dest, dest)
1429 IF (create_list) THEN
1430 list(2, li, dest) = list(2, li, dest) + cost
1431 ELSE
1432 IF (list(1, li, dest) == dest) THEN
1433 tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1434 ELSE
1435 IF (list(2, li, dest) >= cost) THEN
1436 list(2, li, dest) = list(2, li, dest) - cost
1437 tasks(itask)%destination = encode_rank(list(1, li, dest), ilevel, SIZE(rs_descs))
1438 ELSE
1439 tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1440 END IF
1441 END IF
1442 END IF
1443 CASE DEFAULT
1444 cpabort("Unknown task list distribution type")
1445 END SELECT
1446
1447 END DO
1448
1449 END DO
1450
1451 CALL timestop(handle)
1452
1453 END SUBROUTINE compute_load_list
1454! **************************************************************************************************
1455!> \brief small helper function to return the proper index in the list array
1456!>
1457!> \param list ...
1458!> \param rank ...
1459!> \param dest ...
1460!> \return ...
1461!> \par History
1462!> created 2008-10-06 [Joost VandeVondele]
1463! **************************************************************************************************
1464 INTEGER FUNCTION list_index(list, rank, dest)
1465 INTEGER, DIMENSION(:, :, 0:), INTENT(IN) :: list
1466 INTEGER, INTENT(IN) :: rank, dest
1467
1468 list_index = 1
1469 DO
1470 IF (list(1, list_index, dest) == rank) EXIT
1471 list_index = list_index + 1
1472 END DO
1473 END FUNCTION list_index
1474! **************************************************************************************************
1475!> \brief create a list with possible destinations (i.e. the central cpu and neighbors) for each cpu
1476!> note that we allocate it with an additional field to store the load of this destination
1477!>
1478!> \param list ...
1479!> \param rs_descs ...
1480!> \param grid_level ...
1481!> \par History
1482!> created 2008-10-06 [Joost VandeVondele]
1483! **************************************************************************************************
1484 SUBROUTINE create_destination_list(list, rs_descs, grid_level)
1485 INTEGER, DIMENSION(:, :, :), POINTER :: list
1486 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1487 POINTER :: rs_descs
1488 INTEGER, INTENT(IN) :: grid_level
1489
1490 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_destination_list'
1491
1492 INTEGER :: handle, i, icpu, j, maxcount, ncpu, &
1493 ultimate_max
1494 INTEGER, ALLOCATABLE, DIMENSION(:) :: index, sublist
1495
1496 CALL timeset(routinen, handle)
1497
1498 cpassert(.NOT. ASSOCIATED(list))
1499 ncpu = rs_descs(grid_level)%rs_desc%group_size
1500 ultimate_max = 7
1501
1502 ALLOCATE (list(2, ultimate_max, 0:ncpu - 1))
1503
1504 ALLOCATE (index(ultimate_max))
1505 ALLOCATE (sublist(ultimate_max))
1506 sublist = huge(sublist)
1507
1508 maxcount = 1
1509 DO icpu = 0, ncpu - 1
1510 sublist(1) = icpu
1511 sublist(2) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [-1, 0, 0])
1512 sublist(3) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [+1, 0, 0])
1513 sublist(4) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, -1, 0])
1514 sublist(5) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, +1, 0])
1515 sublist(6) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, 0, -1])
1516 sublist(7) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, 0, +1])
1517 ! only retain unique values of the destination
1518 CALL sort(sublist, ultimate_max, index)
1519 j = 1
1520 DO i = 2, 7
1521 IF (sublist(i) /= sublist(j)) THEN
1522 j = j + 1
1523 sublist(j) = sublist(i)
1524 END IF
1525 END DO
1526 maxcount = max(maxcount, j)
1527 sublist(j + 1:ultimate_max) = huge(sublist)
1528 list(1, :, icpu) = sublist
1529 list(2, :, icpu) = 0
1530 END DO
1531
1532 CALL reallocate(list, 1, 2, 1, maxcount, 0, ncpu - 1)
1533
1534 CALL timestop(handle)
1535
1536 END SUBROUTINE create_destination_list
1537
1538! **************************************************************************************************
1539!> \brief given a task list, compute the load of each process everywhere
1540!> giving this function the ability to loop over a (sub)set of rs_grids,
1541!> and do all the communication in one shot, would speed it up
1542!> \param loads ...
1543!> \param rs_descs ...
1544!> \param grid_level ...
1545!> \param ntasks ...
1546!> \param tasks ...
1547!> \param use_reordered_ranks ...
1548!> \par History
1549!> none
1550!> \author MattW 21/11/2007
1551! **************************************************************************************************
1552 SUBROUTINE get_current_loads(loads, rs_descs, grid_level, ntasks, tasks, use_reordered_ranks)
1553 INTEGER(KIND=int_8), DIMENSION(:) :: loads
1554 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1555 POINTER :: rs_descs
1556 INTEGER :: grid_level, ntasks
1557 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1558 LOGICAL, INTENT(IN) :: use_reordered_ranks
1559
1560 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_current_loads'
1561
1562 INTEGER :: handle, i, iatom, ilevel, img, ipgf, &
1563 iset, jatom, jpgf, jset
1564 INTEGER(KIND=int_8) :: total_cost_local
1565 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: recv_buf_i, send_buf_i
1566 TYPE(realspace_grid_desc_type), POINTER :: desc
1567
1568 CALL timeset(routinen, handle)
1569
1570 desc => rs_descs(grid_level)%rs_desc
1571
1572 ! allocate local arrays
1573 ALLOCATE (send_buf_i(desc%group_size))
1574 ALLOCATE (recv_buf_i(desc%group_size))
1575
1576 ! communication step 1 : compute the total local cost of the tasks
1577 ! each proc needs to know the amount of work he will receive
1578
1579 ! send buffer now contains for each target the cost of the tasks it will receive
1580 send_buf_i = 0
1581 DO i = 1, ntasks
1582 ilevel = tasks(i)%grid_level
1583 img = tasks(i)%image
1584 iatom = tasks(i)%iatom
1585 jatom = tasks(i)%jatom
1586 iset = tasks(i)%iset
1587 jset = tasks(i)%jset
1588 ipgf = tasks(i)%ipgf
1589 jpgf = tasks(i)%jpgf
1590 IF (ilevel /= grid_level) cycle
1591 IF (use_reordered_ranks) THEN
1592 send_buf_i(rs_descs(ilevel)%rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs))) + 1) = &
1593 send_buf_i(rs_descs(ilevel)%rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs))) + 1) &
1594 + tasks(i)%cost
1595 ELSE
1596 send_buf_i(decode_rank(tasks(i)%destination, SIZE(rs_descs)) + 1) = &
1597 send_buf_i(decode_rank(tasks(i)%destination, SIZE(rs_descs)) + 1) &
1598 + tasks(i)%cost
1599 END IF
1600 END DO
1601 CALL desc%group%alltoall(send_buf_i, recv_buf_i, 1)
1602
1603 ! communication step 2 : compute the global cost of the tasks
1604 total_cost_local = sum(recv_buf_i)
1605
1606 ! after this step, the recv buffer contains the local cost for each CPU
1607 CALL desc%group%allgather(total_cost_local, loads)
1608
1609 CALL timestop(handle)
1610
1611 END SUBROUTINE get_current_loads
1612! **************************************************************************************************
1613!> \brief performs load balancing shifting tasks on the replicated grids
1614!> this modifies the destination of some of the tasks on replicated
1615!> grids, and in this way balances the load
1616!> \param rs_descs ...
1617!> \param ntasks ...
1618!> \param tasks ...
1619!> \param redistributed true if any task changes rank; identical on all ranks
1620!> \par History
1621!> none
1622!> \author MattW 21/11/2007
1623! **************************************************************************************************
1624 SUBROUTINE load_balance_replicated(rs_descs, ntasks, tasks, redistributed)
1625
1626 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1627 POINTER :: rs_descs
1628 INTEGER :: ntasks
1629 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1630 LOGICAL, INTENT(OUT), OPTIONAL :: redistributed
1631
1632 CHARACTER(LEN=*), PARAMETER :: routinen = 'load_balance_replicated'
1633
1634 INTEGER :: handle, i, iatom, ilevel, img, ipgf, &
1635 iset, j, jatom, jpgf, jset, &
1636 no_overloaded, no_underloaded, &
1637 proc_receiving, moved
1638 INTEGER(KIND=int_8) :: average_cost, cost_task_rep, count, &
1639 max_cost, offset, total_cost_global
1640 LOGICAL :: fully_replicated
1641 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: load_imbalance, loads, recv_buf_i
1642 INTEGER, ALLOCATABLE, DIMENSION(:) :: index
1643 TYPE(realspace_grid_desc_type), POINTER :: desc
1644
1645 CALL timeset(routinen, handle)
1646
1647 desc => rs_descs(1)%rs_desc
1648
1649 IF (PRESENT(redistributed)) redistributed = .false.
1650 IF (desc%group_size == 1) THEN
1651 CALL timestop(handle)
1652 RETURN
1653 END IF
1654
1655 ! allocate local arrays
1656 ALLOCATE (recv_buf_i(desc%group_size))
1657 ALLOCATE (loads(desc%group_size))
1658
1659 fully_replicated = .true.
1660 DO i = 1, SIZE(rs_descs)
1661 fully_replicated = fully_replicated .AND. .NOT. rs_descs(i)%rs_desc%distributed
1662 END DO
1663
1664 IF (fully_replicated) THEN
1665 ! These tasks still reside on their matrix owners. One gather suffices for all levels.
1666 cost_task_rep = 0
1667 DO i = 1, ntasks
1668 cost_task_rep = cost_task_rep + tasks(i)%cost
1669 END DO
1670 CALL desc%group%allgather(cost_task_rep, recv_buf_i)
1671 ELSE
1672 recv_buf_i = 0
1673 DO i = 1, SIZE(rs_descs)
1674 CALL get_current_loads(loads, rs_descs, i, ntasks, tasks, use_reordered_ranks=.true.)
1675 recv_buf_i(:) = recv_buf_i + loads
1676 END DO
1677 END IF
1678
1679 total_cost_global = sum(recv_buf_i)
1680 average_cost = total_cost_global/desc%group_size
1681 max_cost = maxval(recv_buf_i)
1682
1683 IF (fully_replicated) THEN
1684 ! No assignment can improve on the rounded-up mean load.
1685 IF (max_cost <= (total_cost_global + desc%group_size - 1)/desc%group_size) THEN
1686 CALL timestop(handle)
1687 RETURN
1688 END IF
1689 END IF
1690
1691 !
1692 ! compute how to redistribute the replicated tasks so that the average cost is reached
1693 !
1694
1695 ! load imbalance measures the load of a given CPU relative
1696 ! to the optimal load distribution (load=average)
1697 ALLOCATE (load_imbalance(desc%group_size))
1698 ALLOCATE (index(desc%group_size))
1699
1700 load_imbalance(:) = recv_buf_i - average_cost
1701 no_overloaded = 0
1702 no_underloaded = 0
1703
1704 DO i = 1, desc%group_size
1705 IF (load_imbalance(i) > 0) no_overloaded = no_overloaded + 1
1706 IF (load_imbalance(i) < 0) no_underloaded = no_underloaded + 1
1707 END DO
1708
1709 ! sort the recv_buffer on number of tasks, gives us index which provides a
1710 ! mapping between processor ranks and how overloaded the processor
1711 CALL sort(recv_buf_i, SIZE(recv_buf_i), index)
1712
1713 ! find out the number of replicated tasks each proc has
1714 ! but only those tasks which have not yet been assigned
1715 cost_task_rep = 0
1716 DO i = 1, ntasks
1717 IF (tasks(i)%dist_type == 0 &
1718 .AND. decode_rank(tasks(i)%destination, SIZE(rs_descs)) == decode_rank(tasks(i)%source, SIZE(rs_descs))) THEN
1719 cost_task_rep = cost_task_rep + tasks(i)%cost
1720 END IF
1721 END DO
1722
1723 ! now, correct the load imbalance for the overloaded CPUs
1724 ! they will send away not more than the total load of replicated tasks
1725 CALL desc%group%allgather(cost_task_rep, recv_buf_i)
1726
1727 DO i = 1, desc%group_size
1728 ! At the moment we can only offload replicated tasks
1729 IF (load_imbalance(i) > 0) load_imbalance(i) = min(load_imbalance(i), recv_buf_i(i))
1730 END DO
1731
1732 ! simplest algorithm I can think of of is that the processor with the most
1733 ! excess tasks fills up the process needing most, then moves on to next most.
1734 ! At the moment if we've got less replicated tasks than we're overloaded then
1735 ! task balancing will be incomplete
1736
1737 ! only need to do anything if I've excess tasks
1738 IF (load_imbalance(desc%my_pos + 1) > 0) THEN
1739
1740 count = 0 ! weighted amount of tasks offloaded
1741 offset = 0 ! no of underloaded processes already filled by other more overloaded procs
1742
1743 ! calculate offset
1744 DO i = desc%group_size, desc%group_size - no_overloaded + 1, -1
1745 IF (index(i) == desc%my_pos + 1) THEN
1746 EXIT
1747 ELSE
1748 offset = offset + load_imbalance(index(i))
1749 END IF
1750 END DO
1751
1752 ! find my starting processor to send to
1753 proc_receiving = huge(proc_receiving)
1754 DO i = 1, no_underloaded
1755 offset = offset + load_imbalance(index(i))
1756 IF (offset < 0 .OR. (.NOT. fully_replicated .AND. offset == 0)) THEN
1757 proc_receiving = i
1758 EXIT
1759 END IF
1760 END DO
1761
1762 ! offset now contains minus the number of tasks proc_receiving requires
1763 ! we fill this up by adjusting the destination of tasks on the replicated grid,
1764 ! then move to next most underloaded proc
1765 DO j = 1, ntasks
1766 IF (tasks(j)%dist_type == 0 &
1767 .AND. decode_rank(tasks(j)%destination, SIZE(rs_descs)) == decode_rank(tasks(j)%source, SIZE(rs_descs))) THEN
1768 ! just avoid sending to non existing procs due to integer truncation
1769 ! in the computation of the average
1770 IF (proc_receiving > no_underloaded) EXIT
1771 ! set new destination
1772 ilevel = tasks(j)%grid_level
1773 img = tasks(j)%image
1774 iatom = tasks(j)%iatom
1775 jatom = tasks(j)%jatom
1776 iset = tasks(j)%iset
1777 jset = tasks(j)%jset
1778 ipgf = tasks(j)%ipgf
1779 jpgf = tasks(j)%jpgf
1780 tasks(j)%destination = encode_rank(index(proc_receiving) - 1, ilevel, SIZE(rs_descs))
1781 offset = offset + tasks(j)%cost
1782 count = count + tasks(j)%cost
1783 IF (count >= load_imbalance(desc%my_pos + 1)) EXIT
1784 IF (offset > 0 .OR. (fully_replicated .AND. offset == 0)) THEN
1785 proc_receiving = proc_receiving + 1
1786 ! just avoid sending to non existing procs due to integer truncation
1787 ! in the computation of the average
1788 IF (proc_receiving > no_underloaded) EXIT
1789 offset = load_imbalance(index(proc_receiving))
1790 END IF
1791 END IF
1792 END DO
1793 END IF
1794
1795 DEALLOCATE (index)
1796 DEALLOCATE (load_imbalance)
1797
1798 IF (fully_replicated) THEN
1799 ! Indivisible expensive tasks can defeat the greedy assignment. Keep only an improvement.
1800 loads = 0
1801 DO i = 1, ntasks
1802 j = decode_rank(tasks(i)%destination, SIZE(rs_descs)) + 1
1803 loads(j) = loads(j) + tasks(i)%cost
1804 END DO
1805 CALL desc%group%sum(loads)
1806 IF (maxval(loads) >= max_cost) THEN
1807 DO i = 1, ntasks
1808 tasks(i)%destination = tasks(i)%source
1809 END DO
1810 END IF
1811 END IF
1812
1813 IF (PRESENT(redistributed)) THEN
1814 moved = 0
1815 DO i = 1, ntasks
1816 IF (tasks(i)%destination /= tasks(i)%source) moved = 1
1817 END DO
1818 CALL desc%group%max(moved)
1819 redistributed = moved > 0
1820 END IF
1821
1822 CALL timestop(handle)
1823
1824 END SUBROUTINE load_balance_replicated
1825
1826! **************************************************************************************************
1827!> \brief given an input task list, redistribute so that all tasks can be processed locally,
1828!> i.e. dest equals rank
1829!> \param rs_descs ...
1830!> \param ntasks ...
1831!> \param tasks ...
1832!> \param ntasks_recv ...
1833!> \param tasks_recv ...
1834!> \par History
1835!> none
1836!> \author MattW 21/11/2007
1837! **************************************************************************************************
1838 SUBROUTINE create_local_tasks(rs_descs, ntasks, tasks, ntasks_recv, tasks_recv)
1839
1840 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1841 POINTER :: rs_descs
1842 INTEGER :: ntasks
1843 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1844 INTEGER :: ntasks_recv
1845 TYPE(task_type), DIMENSION(:), POINTER :: tasks_recv
1846
1847 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_local_tasks'
1848
1849 INTEGER :: handle, i, j, k, l, rank
1850 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: recv_buf, send_buf
1851 INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_disps, recv_sizes, send_disps, &
1852 send_sizes
1853 TYPE(realspace_grid_desc_type), POINTER :: desc
1854
1855 CALL timeset(routinen, handle)
1856
1857 desc => rs_descs(1)%rs_desc
1858
1859 ! allocate local arrays
1860 ALLOCATE (send_sizes(desc%group_size))
1861 ALLOCATE (recv_sizes(desc%group_size))
1862 ALLOCATE (send_disps(desc%group_size))
1863 ALLOCATE (recv_disps(desc%group_size))
1864 ALLOCATE (send_buf(desc%group_size))
1865 ALLOCATE (recv_buf(desc%group_size))
1866
1867 ! fill send buffer, now counting how many tasks will be send (stored in an int8 array for convenience only).
1868 send_buf = 0
1869 DO i = 1, ntasks
1870 rank = rs_descs(decode_level(tasks(i)%destination, SIZE(rs_descs))) &
1871 %rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs)))
1872 send_buf(rank + 1) = send_buf(rank + 1) + 1
1873 END DO
1874
1875 CALL desc%group%alltoall(send_buf, recv_buf, 1)
1876
1877 ! pack the tasks, and send them around
1878
1879 send_sizes = 0
1880 send_disps = 0
1881 recv_sizes = 0
1882 recv_disps = 0
1883
1884 send_sizes(1) = int(send_buf(1)*task_size_in_int8)
1885 recv_sizes(1) = int(recv_buf(1)*task_size_in_int8)
1886 DO i = 2, desc%group_size
1887 send_sizes(i) = int(send_buf(i)*task_size_in_int8)
1888 recv_sizes(i) = int(recv_buf(i)*task_size_in_int8)
1889 send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1890 recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1891 END DO
1892
1893 ! deallocate old send/recv buffers
1894 DEALLOCATE (send_buf)
1895 DEALLOCATE (recv_buf)
1896
1897 ! allocate them with new sizes
1898 ALLOCATE (send_buf(sum(send_sizes)))
1899 ALLOCATE (recv_buf(sum(recv_sizes)))
1900
1901 ! do packing
1902 send_buf = 0
1903 send_sizes = 0
1904 DO j = 1, ntasks
1905 i = rs_descs(decode_level(tasks(j)%destination, SIZE(rs_descs))) &
1906 %rs_desc%virtual2real(decode_rank(tasks(j)%destination, SIZE(rs_descs))) + 1
1907 l = send_disps(i) + send_sizes(i)
1908 CALL serialize_task(tasks(j), send_buf(l + 1:l + task_size_in_int8))
1909 send_sizes(i) = send_sizes(i) + task_size_in_int8
1910 END DO
1911
1912 ! do communication
1913 CALL desc%group%alltoall(send_buf, send_sizes, send_disps, recv_buf, recv_sizes, recv_disps)
1914
1915 DEALLOCATE (send_buf)
1916
1917 ntasks_recv = sum(recv_sizes)/task_size_in_int8
1918 ALLOCATE (tasks_recv(ntasks_recv))
1919
1920 ! do unpacking
1921 l = 0
1922 DO i = 1, desc%group_size
1923 DO j = 0, recv_sizes(i)/task_size_in_int8 - 1
1924 l = l + 1
1925 k = recv_disps(i) + j*task_size_in_int8
1926 CALL deserialize_task(tasks_recv(l), recv_buf(k + 1:k + task_size_in_int8))
1927 END DO
1928 END DO
1929
1930 DEALLOCATE (recv_buf)
1931 DEALLOCATE (send_sizes)
1932 DEALLOCATE (recv_sizes)
1933 DEALLOCATE (send_disps)
1934 DEALLOCATE (recv_disps)
1935
1936 CALL timestop(handle)
1937
1938 END SUBROUTINE create_local_tasks
1939
1940! **************************************************************************************************
1941!> \brief Assembles tasks to be performed on local grid
1942!> \param rs_descs the grids
1943!> \param ntasks Number of tasks for local processing
1944!> \param natoms ...
1945!> \param nimages ...
1946!> \param tasks the task set generated on this processor
1947!> \param rval ...
1948!> \param atom_pair_send ...
1949!> \param atom_pair_recv ...
1950!> \param symmetric ...
1951!> \param reorder_rs_grid_ranks ...
1952!> \param skip_load_balance_distributed ...
1953!> \param balance_replicated allow balancing when all grid levels are replicated (default true)
1954!> \par History
1955!> none
1956!> \author MattW 21/11/2007
1957! **************************************************************************************************
1958 SUBROUTINE distribute_tasks(rs_descs, ntasks, natoms, &
1959 tasks, atom_pair_send, atom_pair_recv, &
1960 symmetric, reorder_rs_grid_ranks, skip_load_balance_distributed, &
1961 balance_replicated)
1962
1963 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1964 POINTER :: rs_descs
1965 INTEGER :: ntasks, natoms
1966 TYPE(task_type), DIMENSION(:), POINTER :: tasks
1967 TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_send, atom_pair_recv
1968 LOGICAL, INTENT(IN) :: symmetric, reorder_rs_grid_ranks, &
1969 skip_load_balance_distributed
1970 LOGICAL, INTENT(IN), OPTIONAL :: balance_replicated
1971
1972 CHARACTER(LEN=*), PARAMETER :: routinen = 'distribute_tasks'
1973
1974 INTEGER :: handle, igrid_level, irank, ntasks_recv
1975 INTEGER(KIND=int_8) :: load_gap, max_load, replicated_load
1976 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: total_loads, total_loads_tmp, trial_loads
1977 INTEGER(KIND=int_8), DIMENSION(:, :), POINTER :: loads
1978 INTEGER, ALLOCATABLE, DIMENSION(:) :: indices, real2virtual, total_index
1979 LOGICAL :: distributed_grids, fixed_first_grid, &
1980 redistribute_matrices, allow_balance
1981 TYPE(realspace_grid_desc_type), POINTER :: desc
1982 TYPE(task_type), DIMENSION(:), POINTER :: tasks_recv
1983
1984 CALL timeset(routinen, handle)
1985
1986 cpassert(ASSOCIATED(tasks))
1987
1988 ! *** figure out if we have distributed grids
1989 distributed_grids = .false.
1990 DO igrid_level = 1, SIZE(rs_descs)
1991 IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
1992 distributed_grids = .true.
1993 END IF
1994 END DO
1995 desc => rs_descs(1)%rs_desc
1996
1997 IF (distributed_grids) THEN
1998
1999 ALLOCATE (loads(0:desc%group_size - 1, SIZE(rs_descs)))
2000 ALLOCATE (total_loads(0:desc%group_size - 1))
2001
2002 total_loads = 0
2003
2004 ! First round of balancing on the distributed grids
2005 ! we just balance each of the distributed grids independently
2006 DO igrid_level = 1, SIZE(rs_descs)
2007 IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
2008
2009 IF (.NOT. skip_load_balance_distributed) THEN
2010 CALL load_balance_distributed(tasks, ntasks, rs_descs, igrid_level, natoms)
2011 END IF
2012
2013 CALL get_current_loads(loads(:, igrid_level), rs_descs, igrid_level, ntasks, &
2014 tasks, use_reordered_ranks=.false.)
2015
2016 total_loads(:) = total_loads + loads(:, igrid_level)
2017
2018 END IF
2019 END DO
2020
2021 ! calculate the total load of replicated tasks, so we can decide if it is worth
2022 ! reordering the distributed grid levels
2023
2024 replicated_load = 0
2025 DO igrid_level = 1, SIZE(rs_descs)
2026 IF (.NOT. rs_descs(igrid_level)%rs_desc%distributed) THEN
2027 CALL get_current_loads(loads(:, igrid_level), rs_descs, igrid_level, ntasks, &
2028 tasks, use_reordered_ranks=.false.)
2029 replicated_load = replicated_load + sum(loads(:, igrid_level))
2030 END IF
2031 END DO
2032
2033 ! Now we adjust the rank ordering based on the current loads
2034 ! we leave the first distributed level and all the replicated levels in the default order
2035 IF (reorder_rs_grid_ranks) THEN
2036 fixed_first_grid = .false.
2037 DO igrid_level = 1, SIZE(rs_descs)
2038 IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
2039 IF (fixed_first_grid .EQV. .false.) THEN
2040 total_loads(:) = loads(:, igrid_level)
2041 fixed_first_grid = .true.
2042 ELSE
2043 ALLOCATE (trial_loads(0:desc%group_size - 1))
2044
2045 trial_loads(:) = total_loads + loads(:, igrid_level)
2046 max_load = maxval(trial_loads)
2047 load_gap = 0
2048 DO irank = 0, desc%group_size - 1
2049 load_gap = load_gap + max_load - trial_loads(irank)
2050 END DO
2051
2052 ! If there is not enough replicated load to load balance well enough
2053 ! then we will reorder this grid level
2054 IF (load_gap > replicated_load*1.05_dp) THEN
2055
2056 ALLOCATE (indices(0:desc%group_size - 1))
2057 ALLOCATE (total_index(0:desc%group_size - 1))
2058 ALLOCATE (total_loads_tmp(0:desc%group_size - 1))
2059 ALLOCATE (real2virtual(0:desc%group_size - 1))
2060
2061 total_loads_tmp(:) = total_loads
2062 CALL sort(total_loads_tmp, desc%group_size, total_index)
2063 CALL sort(loads(:, igrid_level), desc%group_size, indices)
2064
2065 ! Reorder so that the rank with smallest load on this grid level is paired with
2066 ! the highest load in total
2067 DO irank = 0, desc%group_size - 1
2068 total_loads(total_index(irank) - 1) = total_loads(total_index(irank) - 1) + &
2069 loads(desc%group_size - irank - 1, igrid_level)
2070 real2virtual(total_index(irank) - 1) = indices(desc%group_size - irank - 1) - 1
2071 END DO
2072
2073 CALL rs_grid_reorder_ranks(rs_descs(igrid_level)%rs_desc, real2virtual)
2074
2075 DEALLOCATE (indices)
2076 DEALLOCATE (total_index)
2077 DEALLOCATE (total_loads_tmp)
2078 DEALLOCATE (real2virtual)
2079 ELSE
2080 total_loads(:) = trial_loads
2081 END IF
2082
2083 DEALLOCATE (trial_loads)
2084
2085 END IF
2086 END IF
2087 END DO
2088 END IF
2089
2090 ! Now we use the replicated tasks to balance out the rest of the load
2091 CALL load_balance_replicated(rs_descs, ntasks, tasks)
2092
2093 DEALLOCATE (loads)
2094 DEALLOCATE (total_loads)
2095 END IF
2096
2097 redistribute_matrices = distributed_grids
2098 allow_balance = .true.
2099 IF (PRESENT(balance_replicated)) allow_balance = balance_replicated
2100 IF (.NOT. distributed_grids .AND. allow_balance) THEN
2101 CALL load_balance_replicated(rs_descs, ntasks, tasks, redistributed=redistribute_matrices)
2102 END IF
2103
2104 IF (redistribute_matrices) THEN
2105
2106 ! given a list of tasks, this will do the needed reshuffle so that all tasks will be local
2107 CALL create_local_tasks(rs_descs, ntasks, tasks, ntasks_recv, tasks_recv)
2108
2109 !
2110 ! tasks list are complete, we can compute the list of atomic blocks (atom pairs)
2111 ! we will be sending. These lists are needed for redistribute_matrix.
2112 !
2113 CALL get_atom_pair(atom_pair_send, tasks, ntasks=ntasks, send=.true., symmetric=symmetric, rs_descs=rs_descs)
2114
2115 CALL get_atom_pair(atom_pair_recv, tasks_recv, ntasks=ntasks_recv, send=.false., symmetric=symmetric, rs_descs=rs_descs)
2116
2117 ! cleanup, at this point we don't need the original tasks anymore
2118 DEALLOCATE (tasks)
2119
2120 ELSE
2121 tasks_recv => tasks
2122 ntasks_recv = ntasks
2123 CALL get_atom_pair(atom_pair_recv, tasks_recv, ntasks=ntasks_recv, send=.false., symmetric=symmetric, rs_descs=rs_descs)
2124 ! Tasks remain on matrix owners, so no matrix communication plan is needed.
2125 END IF
2126
2127 ! here we sort the task list we will process locally.
2128 ALLOCATE (indices(ntasks_recv))
2129 CALL tasks_sort(tasks_recv, ntasks_recv, indices)
2130 DEALLOCATE (indices)
2131
2132 !
2133 ! final lists are ready
2134 !
2135
2136 tasks => tasks_recv
2137 ntasks = ntasks_recv
2138
2139 CALL timestop(handle)
2140
2141 END SUBROUTINE distribute_tasks
2142
2143! **************************************************************************************************
2144!> \brief ...
2145!> \param atom_pair ...
2146!> \param my_tasks ...
2147!> \param send ...
2148!> \param symmetric ...
2149!> \param natoms ...
2150!> \param nimages ...
2151!> \param rs_descs ...
2152! **************************************************************************************************
2153 SUBROUTINE get_atom_pair(atom_pair, tasks, ntasks, send, symmetric, rs_descs)
2154
2155 TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair
2156 TYPE(task_type), DIMENSION(:), INTENT(INOUT) :: tasks
2157 INTEGER, INTENT(IN) :: ntasks
2158 LOGICAL, INTENT(IN) :: send, symmetric
2159 TYPE(realspace_grid_desc_p_type), DIMENSION(:), INTENT(IN) :: rs_descs
2160
2161 INTEGER :: i, ilevel, iatom, jatom, npairs, virt_rank
2162 INTEGER, DIMENSION(:), ALLOCATABLE :: indices
2163 TYPE(atom_pair_type), DIMENSION(:), ALLOCATABLE :: atom_pair_tmp
2164
2165 cpassert(.NOT. ASSOCIATED(atom_pair))
2166 IF (ntasks == 0) THEN
2167 ALLOCATE (atom_pair(0))
2168 RETURN
2169 END IF
2170
2171 ! calculate list of atom pairs
2172 ! fill pair list taking into account symmetry
2173 ALLOCATE (atom_pair_tmp(ntasks))
2174 DO i = 1, ntasks
2175 atom_pair_tmp(i)%image = tasks(i)%image
2176 iatom = tasks(i)%iatom
2177 jatom = tasks(i)%jatom
2178 IF (symmetric .AND. iatom > jatom) THEN
2179 ! iatom / jatom swapped
2180 atom_pair_tmp(i)%row = jatom
2181 atom_pair_tmp(i)%col = iatom
2182 ELSE
2183 atom_pair_tmp(i)%row = iatom
2184 atom_pair_tmp(i)%col = jatom
2185 END IF
2186
2187 IF (send) THEN
2188 ! If sending, we need to use the 'real rank' as the pair has to be sent to the process which
2189 ! actually has the correct part of the rs_grid to do the mapping
2190 ilevel = tasks(i)%grid_level
2191 virt_rank = decode_rank(tasks(i)%destination, SIZE(rs_descs))
2192 atom_pair_tmp(i)%rank = rs_descs(ilevel)%rs_desc%virtual2real(virt_rank)
2193 ELSE
2194 ! If we are receiving, then no conversion is needed as the rank is that of the process with the
2195 ! required matrix block, and the ordering of the rs grid is irrelevant
2196 atom_pair_tmp(i)%rank = decode_rank(tasks(i)%source, SIZE(rs_descs))
2197 END IF
2198 END DO
2199
2200 ! find unique atom pairs that I'm sending/receiving
2201 ALLOCATE (indices(ntasks))
2202 CALL atom_pair_sort(atom_pair_tmp, ntasks, indices)
2203 npairs = 1
2204 tasks(indices(1))%pair_index = 1
2205 DO i = 2, ntasks
2206 IF (atom_pair_less_than(atom_pair_tmp(i - 1), atom_pair_tmp(i))) THEN
2207 npairs = npairs + 1
2208 atom_pair_tmp(npairs) = atom_pair_tmp(i)
2209 END IF
2210 tasks(indices(i))%pair_index = npairs
2211 END DO
2212 DEALLOCATE (indices)
2213
2214 ! Copy unique pairs to final location.
2215 ALLOCATE (atom_pair(npairs))
2216 atom_pair(:) = atom_pair_tmp(:npairs)
2217 DEALLOCATE (atom_pair_tmp)
2218
2219 END SUBROUTINE get_atom_pair
2220
2221! **************************************************************************************************
2222!> \brief redistributes the matrix so that it can be used in realspace operations
2223!> i.e. according to the task lists for collocate and integrate.
2224!> This routine can become a bottleneck in large calculations.
2225!> \param rs_descs ...
2226!> \param pmats ...
2227!> \param atom_pair_send ...
2228!> \param atom_pair_recv ...
2229!> \param natoms ...
2230!> \param nimages ...
2231!> \param scatter ...
2232!> \param hmats ...
2233! **************************************************************************************************
2234 SUBROUTINE rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, &
2235 nimages, scatter, hmats)
2236
2237 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
2238 POINTER :: rs_descs
2239 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmats
2240 TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_send, atom_pair_recv
2241 INTEGER :: nimages
2242 LOGICAL :: scatter
2243 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
2244 POINTER :: hmats
2245
2246 CHARACTER(LEN=*), PARAMETER :: routinen = 'rs_distribute_matrix'
2247
2248 INTEGER :: acol, arow, handle, i, img, j, k, l, me, &
2249 nblkcols_total, nblkrows_total, ncol, &
2250 nrow, nthread, nthread_left
2251 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row, recv_disps, &
2252 recv_pair_count, recv_pair_disps, recv_sizes, send_disps, send_pair_count, &
2253 send_pair_disps, send_sizes
2254 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
2255 LOGICAL :: found
2256 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: recv_buf_r, send_buf_r
2257 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block
2258 TYPE(dbcsr_type), POINTER :: hmat, pmat
2259 TYPE(realspace_grid_desc_type), POINTER :: desc
2260 REAL(kind=dp), DIMENSION(:), POINTER :: vector
2261
2262!$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
2263
2264 CALL timeset(routinen, handle)
2265
2266 IF (.NOT. scatter) THEN
2267 cpassert(PRESENT(hmats))
2268 END IF
2269
2270 desc => rs_descs(1)%rs_desc
2271 me = desc%my_pos + 1
2272
2273 ! allocate local arrays
2274 ALLOCATE (send_sizes(desc%group_size))
2275 ALLOCATE (recv_sizes(desc%group_size))
2276 ALLOCATE (send_disps(desc%group_size))
2277 ALLOCATE (recv_disps(desc%group_size))
2278 ALLOCATE (send_pair_count(desc%group_size))
2279 ALLOCATE (recv_pair_count(desc%group_size))
2280 ALLOCATE (send_pair_disps(desc%group_size))
2281 ALLOCATE (recv_pair_disps(desc%group_size))
2282
2283 pmat => pmats(1)%matrix
2284 CALL dbcsr_get_info(pmat, &
2285 row_blk_size=row_blk_size, &
2286 col_blk_size=col_blk_size, &
2287 nblkrows_total=nblkrows_total, &
2288 nblkcols_total=nblkcols_total)
2289 ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total), &
2290 first_col(nblkcols_total), last_col(nblkcols_total))
2291 CALL dbcsr_convert_sizes_to_offsets(row_blk_size, first_row, last_row)
2292 CALL dbcsr_convert_sizes_to_offsets(col_blk_size, first_col, last_col)
2293
2294 ! set up send buffer sizes
2295 send_sizes = 0
2296 send_pair_count = 0
2297 DO i = 1, SIZE(atom_pair_send)
2298 k = atom_pair_send(i)%rank + 1 ! proc we're sending this block to
2299 arow = atom_pair_send(i)%row
2300 acol = atom_pair_send(i)%col
2301 nrow = last_row(arow) - first_row(arow) + 1
2302 ncol = last_col(acol) - first_col(acol) + 1
2303 send_sizes(k) = send_sizes(k) + nrow*ncol
2304 send_pair_count(k) = send_pair_count(k) + 1
2305 END DO
2306
2307 send_disps = 0
2308 send_pair_disps = 0
2309 DO i = 2, desc%group_size
2310 send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
2311 send_pair_disps(i) = send_pair_disps(i - 1) + send_pair_count(i - 1)
2312 END DO
2313
2314 ALLOCATE (send_buf_r(sum(send_sizes)))
2315
2316 ! set up recv buffer
2317
2318 recv_sizes = 0
2319 recv_pair_count = 0
2320 DO i = 1, SIZE(atom_pair_recv)
2321 k = atom_pair_recv(i)%rank + 1 ! proc we're receiving this data from
2322 arow = atom_pair_recv(i)%row
2323 acol = atom_pair_recv(i)%col
2324 nrow = last_row(arow) - first_row(arow) + 1
2325 ncol = last_col(acol) - first_col(acol) + 1
2326 recv_sizes(k) = recv_sizes(k) + nrow*ncol
2327 recv_pair_count(k) = recv_pair_count(k) + 1
2328 END DO
2329
2330 recv_disps = 0
2331 recv_pair_disps = 0
2332 DO i = 2, desc%group_size
2333 recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
2334 recv_pair_disps(i) = recv_pair_disps(i - 1) + recv_pair_count(i - 1)
2335 END DO
2336 ALLOCATE (recv_buf_r(sum(recv_sizes)))
2337
2338!$OMP PARALLEL DEFAULT(NONE) &
2339!$OMP SHARED(desc,send_pair_count,send_pair_disps,nimages),&
2340!$OMP SHARED(last_row,first_row,last_col,first_col),&
2341!$OMP SHARED(pmats,send_buf_r,send_disps,send_sizes),&
2342!$OMP SHARED(atom_pair_send,me,hmats,nblkrows_total),&
2343!$OMP SHARED(atom_pair_recv,recv_buf_r,scatter,recv_pair_disps), &
2344!$OMP SHARED(recv_sizes,recv_disps,recv_pair_count,locks), &
2345!$OMP PRIVATE(i,img,arow,acol,nrow,ncol,p_block,found,j,k,l),&
2346!$OMP PRIVATE(nthread,h_block,nthread_left,hmat,pmat,vector)
2347
2348 nthread = 1
2349!$ nthread = omp_get_num_threads()
2350 nthread_left = 1
2351!$ nthread_left = MAX(1, nthread - 1)
2352
2353 ! do packing
2354!$OMP DO schedule(guided)
2355 DO l = 1, desc%group_size
2356 IF (l == me) cycle
2357 send_sizes(l) = 0
2358 DO i = 1, send_pair_count(l)
2359 arow = atom_pair_send(send_pair_disps(l) + i)%row
2360 acol = atom_pair_send(send_pair_disps(l) + i)%col
2361 img = atom_pair_send(send_pair_disps(l) + i)%image
2362 nrow = last_row(arow) - first_row(arow) + 1
2363 ncol = last_col(acol) - first_col(acol) + 1
2364 pmat => pmats(img)%matrix
2365 CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, block=p_block, found=found)
2366 cpassert(found)
2367
2368 DO k = 1, ncol
2369 DO j = 1, nrow
2370 send_buf_r(send_disps(l) + send_sizes(l) + j + (k - 1)*nrow) = p_block(j, k)
2371 END DO
2372 END DO
2373 send_sizes(l) = send_sizes(l) + nrow*ncol
2374 END DO
2375 END DO
2376!$OMP END DO
2377
2378 IF (.NOT. scatter) THEN
2379 ! We need locks to protect concurrent summation into H
2380!$OMP SINGLE
2381!$ ALLOCATE (locks(nthread*10))
2382!$OMP END SINGLE
2383
2384!$OMP DO
2385!$ do i = 1, nthread*10
2386!$ call omp_init_lock(locks(i))
2387!$ end do
2388!$OMP END DO
2389 END IF
2390
2391!$OMP MASTER
2392 ! do communication
2393 CALL desc%group%alltoall(send_buf_r, send_sizes, send_disps, &
2394 recv_buf_r, recv_sizes, recv_disps)
2395!$OMP END MASTER
2396
2397 ! If this is a scatter, then no need to copy local blocks,
2398 ! If not, we sum them directly into H (bypassing the alltoall)
2399 IF (.NOT. scatter) THEN
2400
2401 ! Distribute work over remaining threads assuming one is still in the alltoall
2402!$OMP DO schedule(dynamic,MAX(1,send_pair_count(me)/nthread_left))
2403 DO i = 1, send_pair_count(me)
2404 arow = atom_pair_send(send_pair_disps(me) + i)%row
2405 acol = atom_pair_send(send_pair_disps(me) + i)%col
2406 img = atom_pair_send(send_pair_disps(me) + i)%image
2407 nrow = last_row(arow) - first_row(arow) + 1
2408 ncol = last_col(acol) - first_col(acol) + 1
2409 hmat => hmats(img)%matrix
2410 pmat => pmats(img)%matrix
2411 CALL dbcsr_get_block_p(matrix=hmat, row=arow, col=acol, block=h_block, found=found)
2412 cpassert(found)
2413 CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, block=p_block, found=found)
2414 cpassert(found)
2415
2416!$ call omp_set_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2417 DO k = 1, ncol
2418 DO j = 1, nrow
2419 h_block(j, k) = h_block(j, k) + p_block(j, k)
2420 END DO
2421 END DO
2422!$ call omp_unset_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2423 END DO
2424!$OMP END DO
2425 ELSE
2426 ! We will insert new blocks into P, so create mutable work matrices
2427 DO img = 1, nimages
2428 pmat => pmats(img)%matrix
2429 CALL dbcsr_work_create(pmat, work_mutable=.true., &
2430 nblks_guess=SIZE(atom_pair_recv)/nthread, sizedata_guess=SIZE(recv_buf_r)/nthread, &
2431 n=nthread)
2432 END DO
2433 END IF
2434
2435! wait for comm and setup to finish
2436!$OMP BARRIER
2437
2438 !do unpacking
2439!$OMP DO schedule(guided)
2440 DO l = 1, desc%group_size
2441 IF (l == me) cycle
2442 recv_sizes(l) = 0
2443 DO i = 1, recv_pair_count(l)
2444 arow = atom_pair_recv(recv_pair_disps(l) + i)%row
2445 acol = atom_pair_recv(recv_pair_disps(l) + i)%col
2446 img = atom_pair_recv(recv_pair_disps(l) + i)%image
2447 nrow = last_row(arow) - first_row(arow) + 1
2448 ncol = last_col(acol) - first_col(acol) + 1
2449 pmat => pmats(img)%matrix
2450 NULLIFY (p_block)
2451 CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, block=p_block, found=found)
2452
2453 IF (PRESENT(hmats)) THEN
2454 hmat => hmats(img)%matrix
2455 CALL dbcsr_get_block_p(matrix=hmat, row=arow, col=acol, block=h_block, found=found)
2456 cpassert(found)
2457 END IF
2458
2459 IF (scatter .AND. .NOT. ASSOCIATED(p_block)) THEN
2460 vector => recv_buf_r(recv_disps(l) + recv_sizes(l) + 1:recv_disps(l) + recv_sizes(l) + nrow*ncol)
2461 CALL dbcsr_put_block(pmat, arow, acol, block=reshape(vector, [nrow, ncol]))
2462 END IF
2463 IF (.NOT. scatter) THEN
2464!$ call omp_set_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2465 DO k = 1, ncol
2466 DO j = 1, nrow
2467 h_block(j, k) = h_block(j, k) + recv_buf_r(recv_disps(l) + recv_sizes(l) + j + (k - 1)*nrow)
2468 END DO
2469 END DO
2470!$ call omp_unset_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2471 END IF
2472 recv_sizes(l) = recv_sizes(l) + nrow*ncol
2473 END DO
2474 END DO
2475!$OMP END DO
2476
2477!$ IF (.not. scatter) THEN
2478!$OMP DO
2479!$ do i = 1, nthread*10
2480!$ call omp_destroy_lock(locks(i))
2481!$ end do
2482!$OMP END DO
2483!$ END IF
2484
2485!$OMP SINGLE
2486!$ IF (.not. scatter) THEN
2487!$ DEALLOCATE (locks)
2488!$ END IF
2489!$OMP END SINGLE NOWAIT
2490
2491 IF (scatter) THEN
2492 ! Blocks were added to P
2493 DO img = 1, nimages
2494 pmat => pmats(img)%matrix
2495 CALL dbcsr_finalize(pmat)
2496 END DO
2497 END IF
2498!$OMP END PARALLEL
2499
2500 DEALLOCATE (send_buf_r)
2501 DEALLOCATE (recv_buf_r)
2502
2503 DEALLOCATE (send_sizes)
2504 DEALLOCATE (recv_sizes)
2505 DEALLOCATE (send_disps)
2506 DEALLOCATE (recv_disps)
2507 DEALLOCATE (send_pair_count)
2508 DEALLOCATE (recv_pair_count)
2509 DEALLOCATE (send_pair_disps)
2510 DEALLOCATE (recv_pair_disps)
2511
2512 DEALLOCATE (first_row, last_row, first_col, last_col)
2513
2514 CALL timestop(handle)
2515
2516 END SUBROUTINE rs_distribute_matrix
2517
2518! **************************************************************************************************
2519!> \brief Calculates offsets and sizes for rs_scatter_matrix and rs_copy_matrix.
2520!> \author Ole Schuett
2521! **************************************************************************************************
2522 SUBROUTINE rs_calc_offsets(pairs, nsgf, group_size, &
2523 pair_offsets, rank_offsets, rank_sizes, buffer_size)
2524 TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: pairs
2525 INTEGER, DIMENSION(:), INTENT(IN) :: nsgf
2526 INTEGER, INTENT(IN) :: group_size
2527 INTEGER, DIMENSION(:), POINTER :: pair_offsets, rank_offsets, rank_sizes
2528 INTEGER, INTENT(INOUT) :: buffer_size
2529
2530 INTEGER :: acol, arow, i, block_size, total_size, k, prev_k
2531
2532 IF (ASSOCIATED(pair_offsets)) DEALLOCATE (pair_offsets)
2533 IF (ASSOCIATED(rank_offsets)) DEALLOCATE (rank_offsets)
2534 IF (ASSOCIATED(rank_sizes)) DEALLOCATE (rank_sizes)
2535
2536 ! calculate buffer_size and pair_offsets
2537 ALLOCATE (pair_offsets(SIZE(pairs)))
2538 total_size = 0
2539 DO i = 1, SIZE(pairs)
2540 pair_offsets(i) = total_size
2541 arow = pairs(i)%row
2542 acol = pairs(i)%col
2543 block_size = nsgf(arow)*nsgf(acol)
2544 total_size = total_size + block_size
2545 END DO
2546 buffer_size = total_size
2547
2548 ! calculate rank_offsets and rank_sizes
2549 ALLOCATE (rank_offsets(group_size))
2550 ALLOCATE (rank_sizes(group_size))
2551 rank_offsets = 0
2552 rank_sizes = 0
2553 IF (SIZE(pairs) > 0) THEN
2554 prev_k = pairs(1)%rank + 1
2555 DO i = 1, SIZE(pairs)
2556 k = pairs(i)%rank + 1
2557 cpassert(k >= prev_k) ! expecting the pairs to be ordered by rank
2558 IF (k > prev_k) THEN
2559 rank_offsets(k) = pair_offsets(i)
2560 rank_sizes(prev_k) = rank_offsets(k) - rank_offsets(prev_k)
2561 prev_k = k
2562 END IF
2563 END DO
2564 rank_sizes(k) = buffer_size - rank_offsets(k) ! complete last rank
2565 END IF
2566
2567 END SUBROUTINE rs_calc_offsets
2568
2569! **************************************************************************************************
2570!> \brief Group contributions to each matrix block without changing the communication order.
2571!> Rebuilt with the task list; contains indices only, independent of matrices and threads.
2572!> \param pairs atom pairs in communication order
2573!> \param block_starts first position of each block in pair_indices, plus an end sentinel
2574!> \param pair_indices original pair indices grouped by destination block
2575! **************************************************************************************************
2576 SUBROUTINE rs_group_atom_pairs(pairs, block_starts, pair_indices)
2577 TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: pairs
2578 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: block_starts, pair_indices
2579
2580 INTEGER :: i, nblocks, npairs
2581 INTEGER, ALLOCATABLE, DIMENSION(:) :: starts
2582 TYPE(atom_pair_type), ALLOCATABLE, DIMENSION(:) :: sorted_pairs
2583
2584 npairs = SIZE(pairs)
2585 ALLOCATE (sorted_pairs(npairs), pair_indices(npairs), starts(npairs + 1))
2586 sorted_pairs(:) = pairs(:)
2587 ! Ignore source rank for grouping. The stable sort preserves contribution order within each block.
2588 sorted_pairs(:)%rank = 0
2589 CALL atom_pair_sort(sorted_pairs, npairs, pair_indices)
2590
2591 nblocks = 0
2592 DO i = 1, npairs
2593 IF (i > 1) THEN
2594 IF (.NOT. atom_pair_less_than(sorted_pairs(i - 1), sorted_pairs(i))) cycle
2595 END IF
2596 nblocks = nblocks + 1
2597 starts(nblocks) = i
2598 END DO
2599 starts(nblocks + 1) = npairs + 1
2600 ALLOCATE (block_starts(nblocks + 1))
2601 block_starts(:) = starts(:nblocks + 1)
2602
2603 END SUBROUTINE rs_group_atom_pairs
2604
2605! **************************************************************************************************
2606!> \brief Scatters dbcsr matrix blocks and receives them into a buffer as needed before collocation.
2607!> \author Ole Schuett
2608! **************************************************************************************************
2609 SUBROUTINE rs_scatter_matrices(src_matrices, dest_buffer, task_list, group)
2610 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2611 TYPE(offload_buffer_type), INTENT(INOUT) :: dest_buffer
2612 TYPE(task_list_type), INTENT(IN) :: task_list
2613 TYPE(mp_comm_type), INTENT(IN) :: group
2614
2615 CHARACTER(LEN=*), PARAMETER :: routinen = 'rs_scatter_matrices'
2616
2617 INTEGER :: handle
2618 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: buffer_send
2619
2620 CALL timeset(routinen, handle)
2621
2622 IF (.NOT. ASSOCIATED(task_list%atom_pair_send)) THEN
2623 CALL rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
2624 CALL timestop(handle)
2625 RETURN
2626 END IF
2627 ALLOCATE (buffer_send(task_list%buffer_size_send))
2628
2629 ! pack dbcsr blocks into send buffer
2630 cpassert(ASSOCIATED(task_list%atom_pair_send))
2631 CALL rs_pack_buffer(src_matrices=src_matrices, &
2632 dest_buffer=buffer_send, &
2633 atom_pair=task_list%atom_pair_send, &
2634 pair_offsets=task_list%pair_offsets_send)
2635
2636 ! mpi all-to-all communication, receiving directly into blocks_recv%buffer.
2637 CALL group%alltoall(buffer_send, task_list%rank_sizes_send, task_list%rank_offsets_send, &
2638 dest_buffer%host_buffer, &
2639 task_list%rank_sizes_recv, task_list%rank_offsets_recv)
2640
2641 DEALLOCATE (buffer_send)
2642 CALL timestop(handle)
2643
2644 END SUBROUTINE rs_scatter_matrices
2645
2646! **************************************************************************************************
2647!> \brief Gather the dbcsr matrix blocks and receives them into a buffer as needed after integration.
2648!> \author Ole Schuett
2649! **************************************************************************************************
2650 SUBROUTINE rs_gather_matrices(src_buffer, dest_matrices, task_list, group)
2651 TYPE(offload_buffer_type), INTENT(IN) :: src_buffer
2652 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2653 TYPE(task_list_type), INTENT(IN) :: task_list
2654 TYPE(mp_comm_type), INTENT(IN) :: group
2655
2656 CHARACTER(LEN=*), PARAMETER :: routinen = 'rs_gather_matrices'
2657
2658 INTEGER :: handle
2659 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: buffer_send
2660
2661 CALL timeset(routinen, handle)
2662
2663 IF (.NOT. ASSOCIATED(task_list%atom_pair_send)) THEN
2664 CALL rs_copy_to_matrices(src_buffer, dest_matrices, task_list)
2665 CALL timestop(handle)
2666 RETURN
2667 END IF
2668
2669 ! Caution: The meaning of send and recv are reversed in this routine.
2670 ALLOCATE (buffer_send(task_list%buffer_size_send)) ! e.g. this is actually used for receiving
2671
2672 ! mpi all-to-all communication
2673 CALL group%alltoall(src_buffer%host_buffer, task_list%rank_sizes_recv, task_list%rank_offsets_recv, &
2674 buffer_send, task_list%rank_sizes_send, task_list%rank_offsets_send)
2675
2676 ! unpack dbcsr blocks from send buffer
2677 cpassert(ASSOCIATED(task_list%atom_pair_send))
2678 CALL rs_unpack_buffer(src_buffer=buffer_send, &
2679 dest_matrices=dest_matrices, &
2680 atom_pair=task_list%atom_pair_send, &
2681 pair_offsets=task_list%pair_offsets_send, &
2682 block_starts=task_list%unpack_block_starts, &
2683 pair_indices=task_list%unpack_pair_indices)
2684
2685 DEALLOCATE (buffer_send)
2686 CALL timestop(handle)
2687
2688 END SUBROUTINE rs_gather_matrices
2689
2690! **************************************************************************************************
2691!> \brief Copies the DBCSR blocks into buffer, replaces rs_scatter_matrix for non-distributed grids.
2692!> \author Ole Schuett
2693! **************************************************************************************************
2694 SUBROUTINE rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
2695 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2696 TYPE(offload_buffer_type), INTENT(INOUT) :: dest_buffer
2697 TYPE(task_list_type), INTENT(IN) :: task_list
2698
2699 CALL rs_pack_buffer(src_matrices=src_matrices, &
2700 dest_buffer=dest_buffer%host_buffer, &
2701 atom_pair=task_list%atom_pair_recv, &
2702 pair_offsets=task_list%pair_offsets_recv)
2703
2704 END SUBROUTINE rs_copy_to_buffer
2705
2706! **************************************************************************************************
2707!> \brief Copies from buffer into DBCSR matrics, replaces rs_gather_matrix for non-distributed grids.
2708!> \author Ole Schuett
2709! **************************************************************************************************
2710 SUBROUTINE rs_copy_to_matrices(src_buffer, dest_matrices, task_list)
2711 TYPE(offload_buffer_type), INTENT(IN) :: src_buffer
2712 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2713 TYPE(task_list_type), INTENT(IN) :: task_list
2714
2715 CALL rs_unpack_buffer(src_buffer=src_buffer%host_buffer, &
2716 dest_matrices=dest_matrices, &
2717 atom_pair=task_list%atom_pair_recv, &
2718 pair_offsets=task_list%pair_offsets_recv, &
2719 block_starts=task_list%unpack_block_starts, &
2720 pair_indices=task_list%unpack_pair_indices)
2721
2722 END SUBROUTINE rs_copy_to_matrices
2723
2724! **************************************************************************************************
2725!> \brief Helper routine for rs_scatter_matrix and rs_copy_to_buffer.
2726!> \author Ole Schuett
2727! **************************************************************************************************
2728 SUBROUTINE rs_pack_buffer(src_matrices, dest_buffer, atom_pair, pair_offsets)
2729 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2730 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: dest_buffer
2731 TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: atom_pair
2732 INTEGER, DIMENSION(:), INTENT(IN) :: pair_offsets
2733
2734 INTEGER :: acol, arow, img, i, offset, block_size
2735 LOGICAL :: found
2736 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
2737
2738!$OMP PARALLEL DEFAULT(NONE), &
2739!$OMP SHARED(src_matrices,atom_pair,pair_offsets,dest_buffer), &
2740!$OMP PRIVATE(acol,arow,img,i,offset,block_size,found,block)
2741!$OMP DO schedule(guided)
2742 DO i = 1, SIZE(atom_pair)
2743 arow = atom_pair(i)%row
2744 acol = atom_pair(i)%col
2745 img = atom_pair(i)%image
2746 CALL dbcsr_get_block_p(matrix=src_matrices(img)%matrix, row=arow, col=acol, &
2747 block=block, found=found)
2748 cpassert(found)
2749 block_size = SIZE(block)
2750 offset = pair_offsets(i)
2751 dest_buffer(offset + 1:offset + block_size) = reshape(block, shape=[block_size])
2752 END DO
2753!$OMP END DO
2754!$OMP END PARALLEL
2755
2756 END SUBROUTINE rs_pack_buffer
2757
2758! **************************************************************************************************
2759!> \brief Helper routine for rs_gather_matrix and rs_copy_to_matrices.
2760!> \author Ole Schuett
2761! **************************************************************************************************
2762 SUBROUTINE rs_unpack_buffer(src_buffer, dest_matrices, atom_pair, pair_offsets, block_starts, pair_indices)
2763 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: src_buffer
2764 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2765 TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: atom_pair
2766 INTEGER, DIMENSION(:), INTENT(IN) :: pair_offsets
2767 INTEGER, DIMENSION(:), INTENT(IN) :: block_starts, pair_indices
2768
2769 INTEGER :: acol, arow, i, iblock, img, ipair, &
2770 j, ncols, nrows, offset, pair
2771 LOGICAL :: found
2772 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
2773
2774 ! One thread owns each destination block and accumulates all ranks' contributions.
2775!$OMP PARALLEL DEFAULT(NONE), &
2776!$OMP SHARED(src_buffer,atom_pair,pair_offsets,dest_matrices,block_starts,pair_indices), &
2777!$OMP PRIVATE(acol,arow,i,iblock,img,ipair,j,ncols,nrows,offset,pair,found,block)
2778!$OMP DO schedule(guided)
2779 DO iblock = 1, SIZE(block_starts) - 1
2780 pair = pair_indices(block_starts(iblock))
2781 arow = atom_pair(pair)%row
2782 acol = atom_pair(pair)%col
2783 img = atom_pair(pair)%image
2784 CALL dbcsr_get_block_p(matrix=dest_matrices(img)%matrix, row=arow, col=acol, &
2785 block=block, found=found)
2786 cpassert(found)
2787 nrows = SIZE(block, 1)
2788 ncols = SIZE(block, 2)
2789 DO ipair = block_starts(iblock), block_starts(iblock + 1) - 1
2790 offset = pair_offsets(pair_indices(ipair))
2791 DO j = 1, ncols
2792 DO i = 1, nrows
2793 block(i, j) = block(i, j) + src_buffer(offset + i + (j - 1)*nrows)
2794 END DO
2795 END DO
2796 END DO
2797 END DO
2798!$OMP END DO
2799!$OMP END PARALLEL
2800
2801 END SUBROUTINE rs_unpack_buffer
2802
2803! **************************************************************************************************
2804!> \brief determines the rank of the processor who's real rs grid contains point
2805!> \param rs_desc ...
2806!> \param igrid_level ...
2807!> \param n_levels ...
2808!> \param cube_center ...
2809!> \param ntasks ...
2810!> \param tasks ...
2811!> \param lb_cube ...
2812!> \param ub_cube ...
2813!> \param added_tasks ...
2814!> \par History
2815!> 11.2007 created [MattW]
2816!> 10.2008 rewritten [Joost VandeVondele]
2817!> \author MattW
2818! **************************************************************************************************
2819 SUBROUTINE rs_find_node(rs_desc, igrid_level, n_levels, cube_center, ntasks, tasks, &
2820 lb_cube, ub_cube, added_tasks)
2821
2822 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
2823 INTEGER, INTENT(IN) :: igrid_level, n_levels
2824 INTEGER, DIMENSION(3), INTENT(IN) :: cube_center
2825 INTEGER, INTENT(INOUT) :: ntasks
2826 TYPE(task_type), DIMENSION(:), POINTER :: tasks
2827 INTEGER, DIMENSION(3), INTENT(IN) :: lb_cube, ub_cube
2828 INTEGER, INTENT(OUT) :: added_tasks
2829
2830 INTEGER, PARAMETER :: add_tasks = 1000
2831 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp
2832
2833 INTEGER :: bit_index, coord(3), curr_tasks, dest, i, icoord(3), idest, itask, ix, iy, iz, &
2834 lb_coord(3), lb_domain(3), lbc(3), ub_coord(3), ub_domain(3), ubc(3)
2835 INTEGER :: bit_pattern
2836 LOGICAL :: dir_periodic(3)
2837
2838 coord(1) = rs_desc%x2coord(cube_center(1))
2839 coord(2) = rs_desc%y2coord(cube_center(2))
2840 coord(3) = rs_desc%z2coord(cube_center(3))
2841 dest = rs_desc%coord2rank(coord(1), coord(2), coord(3))
2842
2843 ! the real cube coordinates
2844 lbc = lb_cube + cube_center
2845 ubc = ub_cube + cube_center
2846
2847 IF (all((rs_desc%lb_global(:, dest) - rs_desc%border) <= lbc) .AND. &
2848 all((rs_desc%ub_global(:, dest) + rs_desc%border) >= ubc)) THEN
2849 !standard distributed collocation/integration
2850 tasks(ntasks)%destination = encode_rank(dest, igrid_level, n_levels)
2851 tasks(ntasks)%dist_type = 1
2852 tasks(ntasks)%subpatch_pattern = 0
2853 added_tasks = 1
2854
2855 ! here we figure out if there is an alternate location for this task
2856 ! i.e. even though the cube_center is not in the real local domain,
2857 ! it might fully fit in the halo of the neighbor
2858 ! if its radius is smaller than the maximum radius
2859 ! the list of possible neighbors is stored here as a bit pattern
2860 ! to reconstruct what the rank of the neigbor is,
2861 ! we can use (note this requires the correct rs_grid)
2862 ! IF (BTEST(bit_pattern,0)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[-1,0,0])
2863 ! IF (BTEST(bit_pattern,1)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[+1,0,0])
2864 ! IF (BTEST(bit_pattern,2)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,-1,0])
2865 ! IF (BTEST(bit_pattern,3)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,+1,0])
2866 ! IF (BTEST(bit_pattern,4)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,0,-1])
2867 ! IF (BTEST(bit_pattern,5)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,0,+1])
2868 bit_index = 0
2869 bit_pattern = 0
2870 DO i = 1, 3
2871 IF (rs_desc%perd(i) == 1) THEN
2872 bit_pattern = ibclr(bit_pattern, bit_index)
2873 bit_index = bit_index + 1
2874 bit_pattern = ibclr(bit_pattern, bit_index)
2875 bit_index = bit_index + 1
2876 ELSE
2877 ! fits the left neighbor ?
2878 IF (ubc(i) <= rs_desc%lb_global(i, dest) - 1 + rs_desc%border) THEN
2879 bit_pattern = ibset(bit_pattern, bit_index)
2880 bit_index = bit_index + 1
2881 ELSE
2882 bit_pattern = ibclr(bit_pattern, bit_index)
2883 bit_index = bit_index + 1
2884 END IF
2885 ! fits the right neighbor ?
2886 IF (lbc(i) >= rs_desc%ub_global(i, dest) + 1 - rs_desc%border) THEN
2887 bit_pattern = ibset(bit_pattern, bit_index)
2888 bit_index = bit_index + 1
2889 ELSE
2890 bit_pattern = ibclr(bit_pattern, bit_index)
2891 bit_index = bit_index + 1
2892 END IF
2893 END IF
2894 END DO
2895 tasks(ntasks)%subpatch_pattern = bit_pattern
2896
2897 ELSE
2898 ! generalised collocation/integration needed
2899 ! first we figure out how many neighbors we have to add to include the lbc/ubc
2900 ! in the available domains (inclusive of halo points)
2901 ! first we 'ignore' periodic boundary conditions
2902 ! i.e. ub_coord-lb_coord+1 might be larger than group_dim
2903 lb_coord = coord
2904 ub_coord = coord
2905 lb_domain = rs_desc%lb_global(:, dest) - rs_desc%border
2906 ub_domain = rs_desc%ub_global(:, dest) + rs_desc%border
2907 DO i = 1, 3
2908 ! only if the grid is not periodic in this direction we need to take care of adding neighbors
2909 IF (rs_desc%perd(i) == 0) THEN
2910 ! if the domain lower bound is greater than the lbc we need to add the size of the neighbor domain
2911 DO
2912 IF (lb_domain(i) > lbc(i)) THEN
2913 lb_coord(i) = lb_coord(i) - 1
2914 icoord = modulo(lb_coord, rs_desc%group_dim)
2915 idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2916 lb_domain(i) = lb_domain(i) - (rs_desc%ub_global(i, idest) - rs_desc%lb_global(i, idest) + 1)
2917 ELSE
2918 EXIT
2919 END IF
2920 END DO
2921 ! same for the upper bound
2922 DO
2923 IF (ub_domain(i) < ubc(i)) THEN
2924 ub_coord(i) = ub_coord(i) + 1
2925 icoord = modulo(ub_coord, rs_desc%group_dim)
2926 idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2927 ub_domain(i) = ub_domain(i) + (rs_desc%ub_global(i, idest) - rs_desc%lb_global(i, idest) + 1)
2928 ELSE
2929 EXIT
2930 END IF
2931 END DO
2932 END IF
2933 END DO
2934
2935 ! some care is needed for the periodic boundaries
2936 DO i = 1, 3
2937 IF (ub_domain(i) - lb_domain(i) + 1 >= rs_desc%npts(i)) THEN
2938 dir_periodic(i) = .true.
2939 lb_coord(i) = 0
2940 ub_coord(i) = rs_desc%group_dim(i) - 1
2941 ELSE
2942 dir_periodic(i) = .false.
2943 END IF
2944 END DO
2945
2946 added_tasks = product(ub_coord - lb_coord + 1)
2947 itask = ntasks
2948 ntasks = ntasks + added_tasks - 1
2949 IF (ntasks > SIZE(tasks)) THEN
2950 curr_tasks = int((SIZE(tasks) + add_tasks)*mult_tasks)
2951 CALL reallocate_tasks(tasks, curr_tasks)
2952 END IF
2953 DO iz = lb_coord(3), ub_coord(3)
2954 DO iy = lb_coord(2), ub_coord(2)
2955 DO ix = lb_coord(1), ub_coord(1)
2956 icoord = modulo([ix, iy, iz], rs_desc%group_dim)
2957 idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2958 tasks(itask)%destination = encode_rank(idest, igrid_level, n_levels)
2959 tasks(itask)%dist_type = 2
2960 tasks(itask)%subpatch_pattern = 0
2961 ! encode the domain size for this task
2962 ! if the bit is set, we need to add the border in that direction
2963 IF (ix == lb_coord(1) .AND. .NOT. dir_periodic(1)) THEN
2964 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 0)
2965 END IF
2966 IF (ix == ub_coord(1) .AND. .NOT. dir_periodic(1)) THEN
2967 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 1)
2968 END IF
2969 IF (iy == lb_coord(2) .AND. .NOT. dir_periodic(2)) THEN
2970 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 2)
2971 END IF
2972 IF (iy == ub_coord(2) .AND. .NOT. dir_periodic(2)) THEN
2973 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 3)
2974 END IF
2975 IF (iz == lb_coord(3) .AND. .NOT. dir_periodic(3)) THEN
2976 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 4)
2977 END IF
2978 IF (iz == ub_coord(3) .AND. .NOT. dir_periodic(3)) THEN
2979 tasks(itask)%subpatch_pattern = ibset(tasks(itask)%subpatch_pattern, 5)
2980 END IF
2981 itask = itask + 1
2982 END DO
2983 END DO
2984 END DO
2985 END IF
2986
2987 END SUBROUTINE rs_find_node
2988
2989! **************************************************************************************************
2990!> \brief utility functions for encoding the grid level with a rank, allowing
2991!> different grid levels to maintain a different rank ordering without
2992!> losing information. These encoded_ints are stored in the tasks(1:2,:) array
2993!> \param rank ...
2994!> \param grid_level ...
2995!> \param n_levels ...
2996!> \return ...
2997!> \par History
2998!> 4.2009 created [Iain Bethune]
2999!> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
3000! **************************************************************************************************
3001 FUNCTION encode_rank(rank, grid_level, n_levels) RESULT(encoded_int)
3002
3003 INTEGER, INTENT(IN) :: rank, grid_level, n_levels
3004 INTEGER :: encoded_int
3005
3006! ordered so can still sort by rank
3007
3008 encoded_int = rank*n_levels + grid_level - 1
3009
3010 END FUNCTION encode_rank
3011
3012! **************************************************************************************************
3013!> \brief ...
3014!> \param encoded_int ...
3015!> \param n_levels ...
3016!> \return ...
3017! **************************************************************************************************
3018 FUNCTION decode_rank(encoded_int, n_levels) RESULT(rank)
3019
3020 INTEGER, INTENT(IN) :: encoded_int
3021 INTEGER, INTENT(IN) :: n_levels
3022 INTEGER :: rank
3023
3024 rank = int(encoded_int/n_levels)
3025
3026 END FUNCTION decode_rank
3027
3028! **************************************************************************************************
3029!> \brief ...
3030!> \param encoded_int ...
3031!> \param n_levels ...
3032!> \return ...
3033! **************************************************************************************************
3034 FUNCTION decode_level(encoded_int, n_levels) RESULT(grid_level)
3035
3036 INTEGER, INTENT(IN) :: encoded_int
3037 INTEGER, INTENT(IN) :: n_levels
3038 INTEGER :: grid_level
3039
3040 grid_level = int(modulo(encoded_int, n_levels)) + 1
3041
3042 END FUNCTION decode_level
3043
3044! **************************************************************************************************
3045!> \brief Sort pgf index pair (ipgf,iset,iatom),(jpgf,jset,jatom) for which all atom pairs are
3046!> grouped, and for each atom pair all set pairs are grouped and for each set pair,
3047!> all pgfs are grouped.
3048!> This yields the right order of the tasks for collocation after the sort
3049!> in distribute_tasks. E.g. for a atom pair, all sets and pgfs are computed in one go.
3050!> The exception is the gridlevel. Tasks are first ordered wrt to grid_level. This implies
3051!> that a given density matrix block will be decontracted several times,
3052!> but cache effects on the grid make up for this.
3053!> \param a ...
3054!> \param b ...
3055!> \return ...
3056!> \author Ole Schuett
3057! **************************************************************************************************
3058 PURE FUNCTION tasks_less_than(a, b) RESULT(res)
3059 TYPE(task_type), INTENT(IN) :: a, b
3060 LOGICAL :: res
3061
3062 IF (a%grid_level /= b%grid_level) THEN
3063 res = a%grid_level < b%grid_level
3064
3065 ELSE IF (a%image /= b%image) THEN
3066 res = a%image < b%image
3067
3068 ELSE IF (a%iatom /= b%iatom) THEN
3069 res = a%iatom < b%iatom
3070
3071 ELSE IF (a%jatom /= b%jatom) THEN
3072 res = a%jatom < b%jatom
3073
3074 ELSE IF (a%iset /= b%iset) THEN
3075 res = a%iset < b%iset
3076
3077 ELSE IF (a%jset /= b%jset) THEN
3078 res = a%jset < b%jset
3079
3080 ELSE IF (a%ipgf /= b%ipgf) THEN
3081 res = a%ipgf < b%ipgf
3082
3083 ELSE
3084 res = a%jpgf < b%jpgf
3085
3086 END IF
3087 END FUNCTION tasks_less_than
3088
3089
3090! **************************************************************************************************
3091!> \brief Sorts an array inplace using a combination of merge- and bubble-sort.
3092!> It also returns the indices, which the elements had before the sort.
3093!> \param arr the array to sort
3094!> \param n length of array
3095!> \param indices returns elements-indices before the sort
3096!> \par History
3097!> 12.2012 created [ole]
3098!> \author Ole Schuett
3099! **************************************************************************************************
3100 SUBROUTINE tasks_sort(arr, n, indices)
3101 INTEGER, INTENT(IN) :: n
3102 TYPE(task_type), DIMENSION(1:n), INTENT(INOUT) :: arr
3103 integer, DIMENSION(1:n), INTENT(INOUT) :: indices
3104
3105 INTEGER :: i
3106 TYPE(task_type), ALLOCATABLE :: tmp_arr(:)
3107 INTEGER, ALLOCATABLE :: tmp_idx(:)
3108
3109 IF (n > 1) THEN
3110 ! scratch space used during the merge step
3111 ALLOCATE (tmp_arr((n + 1)/2), tmp_idx((n + 1)/2))
3112
3113 indices = [(i, i=1, n)]
3114
3115 CALL tasks_sort_low(arr(1:n), indices, tmp_arr, tmp_idx)
3116
3117 DEALLOCATE (tmp_arr, tmp_idx)
3118 ELSE IF (n > 0) THEN ! can be a frequent case in cp2k
3119 indices(1) = 1
3120 END IF
3121
3122 END SUBROUTINE tasks_sort
3123
3124! **************************************************************************************************
3125!> \brief The actual sort routing. Only tasks_sort and itself should call this.
3126!> \param arr the array to sort
3127!> \param indices elements-indices before the sort
3128!> \param tmp_arr scratch space
3129!> \param tmp_idx scratch space
3130!> \par History
3131!> 12.2012 created [ole]
3132!> \author Ole Schuett
3133! **************************************************************************************************
3134 RECURSIVE SUBROUTINE tasks_sort_low(arr, indices, tmp_arr, tmp_idx)
3135 TYPE(task_type), DIMENSION(:), INTENT(INOUT) :: arr
3136 INTEGER, DIMENSION(size(arr)), INTENT(INOUT) :: indices
3137 TYPE(task_type), DIMENSION((size(arr) + 1)/2), INTENT(INOUT) :: tmp_arr
3138 INTEGER, DIMENSION((size(arr) + 1)/2), INTENT(INOUT) :: tmp_idx
3139 TYPE(task_type) :: a
3140 INTEGER :: t, m, i, j, k
3141 LOGICAL :: swapped
3142 ! a,t: used during swaping of elements in arr and indices
3143
3144 swapped = .true.
3145
3146 ! If only a few elements are left we switch to bubble-sort for efficiency.
3147 IF (size(arr) <= 7) THEN ! 7 seems to be a good choice for the moment
3148 DO j = size(arr) - 1, 1, -1
3149 swapped = .false.
3150 DO i = 1, j
3151 IF (tasks_less_than(arr(i + 1), arr(i))) THEN
3152 ! swap arr(i) with arr(i+1)
3153 a = arr(i)
3154 arr(i) = arr(i + 1)
3155 arr(i + 1) = a
3156 ! swap indices(i) with indices(i+1)
3157 t = indices(i)
3158 indices(i) = indices(i + 1)
3159 indices(i + 1) = t
3160 swapped = .true.
3161 END IF
3162 END DO
3163 IF (.NOT. swapped) EXIT
3164 END DO
3165 RETURN
3166 END IF
3167
3168 ! split list in half and recursively sort both sublists
3169 m = (size(arr) + 1)/2 ! index where we're going to divide the list in two
3170 CALL tasks_sort_low(arr(1:m), indices(1:m), tmp_arr, tmp_idx)
3171 CALL tasks_sort_low(arr(m + 1:), indices(m + 1:), tmp_arr, tmp_idx)
3172
3173 ! Check for a special case: Can we just concate the two sorted sublists?
3174 ! This leads to O(n) scaling if the input is already sorted.
3175 IF (tasks_less_than(arr(m + 1), arr(m))) THEN
3176 ! ...no - let's merge the two sorted sublists arr(:m) and arr(m+1:)
3177 ! Merge will be performed directly in arr. Need backup of first sublist.
3178 tmp_arr(1:m) = arr(1:m)
3179 tmp_idx(1:m) = indices(1:m)
3180 i = 1 ! number of elements consumed from 1st sublist
3181 j = 1 ! number of elements consumed from 2nd sublist
3182 k = 1 ! number of elements already merged
3183
3184 DO WHILE (i <= m .and. j <= size(arr) - m)
3185 IF (tasks_less_than(arr(m + j), tmp_arr(i))) THEN
3186 arr(k) = arr(m + j)
3187 indices(k) = indices(m + j)
3188 j = j + 1
3189 ELSE
3190 arr(k) = tmp_arr(i)
3191 indices(k) = tmp_idx(i)
3192 i = i + 1
3193 END IF
3194 k = k + 1
3195 END DO
3196
3197 ! One of the two sublist is now empty.
3198 ! Copy possibly remaining tail of 1st sublist
3199 DO WHILE (i <= m)
3200 arr(k) = tmp_arr(i)
3201 indices(k) = tmp_idx(i)
3202 i = i + 1
3203 k = k + 1
3204 END DO
3205 ! The possibly remaining tail of 2nd sublist is already at the right spot.
3206 END IF
3207
3208 END SUBROUTINE tasks_sort_low
3209
3210
3211! **************************************************************************************************
3212!> \brief Order atom pairs to find duplicates.
3213!> \param a ...
3214!> \param b ...
3215!> \return ...
3216!> \author Ole Schuett
3217! **************************************************************************************************
3218 PURE FUNCTION atom_pair_less_than(a, b) RESULT(res)
3219 TYPE(atom_pair_type), INTENT(IN) :: a, b
3220 LOGICAL :: res
3221
3222 IF (a%rank /= b%rank) THEN
3223 res = a%rank < b%rank
3224
3225 ELSE IF (a%row /= b%row) THEN
3226 res = a%row < b%row
3227
3228 ELSE IF (a%col /= b%col) THEN
3229 res = a%col < b%col
3230
3231 ELSE
3232 res = a%image < b%image
3233
3234 END IF
3235 END FUNCTION atom_pair_less_than
3236
3237
3238! **************************************************************************************************
3239!> \brief Sorts an array inplace using a combination of merge- and bubble-sort.
3240!> It also returns the indices, which the elements had before the sort.
3241!> \param arr the array to sort
3242!> \param n length of array
3243!> \param indices returns elements-indices before the sort
3244!> \par History
3245!> 12.2012 created [ole]
3246!> \author Ole Schuett
3247! **************************************************************************************************
3248 SUBROUTINE atom_pair_sort(arr, n, indices)
3249 INTEGER, INTENT(IN) :: n
3250 TYPE(atom_pair_type), DIMENSION(1:n), INTENT(INOUT) :: arr
3251 integer, DIMENSION(1:n), INTENT(INOUT) :: indices
3252
3253 INTEGER :: i
3254 TYPE(atom_pair_type), ALLOCATABLE :: tmp_arr(:)
3255 INTEGER, ALLOCATABLE :: tmp_idx(:)
3256
3257 IF (n > 1) THEN
3258 ! scratch space used during the merge step
3259 ALLOCATE (tmp_arr((n + 1)/2), tmp_idx((n + 1)/2))
3260
3261 indices = [(i, i=1, n)]
3262
3263 CALL atom_pair_sort_low(arr(1:n), indices, tmp_arr, tmp_idx)
3264
3265 DEALLOCATE (tmp_arr, tmp_idx)
3266 ELSE IF (n > 0) THEN ! can be a frequent case in cp2k
3267 indices(1) = 1
3268 END IF
3269
3270 END SUBROUTINE atom_pair_sort
3271
3272! **************************************************************************************************
3273!> \brief The actual sort routing. Only atom_pair_sort and itself should call this.
3274!> \param arr the array to sort
3275!> \param indices elements-indices before the sort
3276!> \param tmp_arr scratch space
3277!> \param tmp_idx scratch space
3278!> \par History
3279!> 12.2012 created [ole]
3280!> \author Ole Schuett
3281! **************************************************************************************************
3282 RECURSIVE SUBROUTINE atom_pair_sort_low(arr, indices, tmp_arr, tmp_idx)
3283 TYPE(atom_pair_type), DIMENSION(:), INTENT(INOUT) :: arr
3284 INTEGER, DIMENSION(size(arr)), INTENT(INOUT) :: indices
3285 TYPE(atom_pair_type), DIMENSION((size(arr) + 1)/2), INTENT(INOUT) :: tmp_arr
3286 INTEGER, DIMENSION((size(arr) + 1)/2), INTENT(INOUT) :: tmp_idx
3287 TYPE(atom_pair_type) :: a
3288 INTEGER :: t, m, i, j, k
3289 LOGICAL :: swapped
3290 ! a,t: used during swaping of elements in arr and indices
3291
3292 swapped = .true.
3293
3294 ! If only a few elements are left we switch to bubble-sort for efficiency.
3295 IF (size(arr) <= 7) THEN ! 7 seems to be a good choice for the moment
3296 DO j = size(arr) - 1, 1, -1
3297 swapped = .false.
3298 DO i = 1, j
3299 IF (atom_pair_less_than(arr(i + 1), arr(i))) THEN
3300 ! swap arr(i) with arr(i+1)
3301 a = arr(i)
3302 arr(i) = arr(i + 1)
3303 arr(i + 1) = a
3304 ! swap indices(i) with indices(i+1)
3305 t = indices(i)
3306 indices(i) = indices(i + 1)
3307 indices(i + 1) = t
3308 swapped = .true.
3309 END IF
3310 END DO
3311 IF (.NOT. swapped) EXIT
3312 END DO
3313 RETURN
3314 END IF
3315
3316 ! split list in half and recursively sort both sublists
3317 m = (size(arr) + 1)/2 ! index where we're going to divide the list in two
3318 CALL atom_pair_sort_low(arr(1:m), indices(1:m), tmp_arr, tmp_idx)
3319 CALL atom_pair_sort_low(arr(m + 1:), indices(m + 1:), tmp_arr, tmp_idx)
3320
3321 ! Check for a special case: Can we just concate the two sorted sublists?
3322 ! This leads to O(n) scaling if the input is already sorted.
3323 IF (atom_pair_less_than(arr(m + 1), arr(m))) THEN
3324 ! ...no - let's merge the two sorted sublists arr(:m) and arr(m+1:)
3325 ! Merge will be performed directly in arr. Need backup of first sublist.
3326 tmp_arr(1:m) = arr(1:m)
3327 tmp_idx(1:m) = indices(1:m)
3328 i = 1 ! number of elements consumed from 1st sublist
3329 j = 1 ! number of elements consumed from 2nd sublist
3330 k = 1 ! number of elements already merged
3331
3332 DO WHILE (i <= m .and. j <= size(arr) - m)
3333 IF (atom_pair_less_than(arr(m + j), tmp_arr(i))) THEN
3334 arr(k) = arr(m + j)
3335 indices(k) = indices(m + j)
3336 j = j + 1
3337 ELSE
3338 arr(k) = tmp_arr(i)
3339 indices(k) = tmp_idx(i)
3340 i = i + 1
3341 END IF
3342 k = k + 1
3343 END DO
3344
3345 ! One of the two sublist is now empty.
3346 ! Copy possibly remaining tail of 1st sublist
3347 DO WHILE (i <= m)
3348 arr(k) = tmp_arr(i)
3349 indices(k) = tmp_idx(i)
3350 i = i + 1
3351 k = k + 1
3352 END DO
3353 ! The possibly remaining tail of 2nd sublist is already at the right spot.
3354 END IF
3355
3356 END SUBROUTINE atom_pair_sort_low
3357
3358
3359END MODULE task_list_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
cell_types::pbc
Definition cell_types.F:105
memory_utilities::reallocate
Definition memory_utilities.F:27
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:208
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_work_create
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
Definition cp_dbcsr_api.F:1307
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition cube_utils.F:18
cube_utils::compute_cube_center
subroutine, public compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab)
unifies the computation of the cube center, so that differences in implementation,...
Definition cube_utils.F:68
cube_utils::return_cube
subroutine, public return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
...
Definition cube_utils.F:142
cube_utils::return_cube_nonortho
subroutine, public return_cube_nonortho(info, radius, lb, ub, rp)
...
Definition cube_utils.F:91
gaussian_gridlevels
Definition gaussian_gridlevels.F:13
gaussian_gridlevels::gaussian_gridlevel
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Definition gaussian_gridlevels.F:148
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::grid_create_task_list
subroutine, public grid_create_task_list(ntasks, natoms, nkinds, nblocks, block_offsets, atom_positions, atom_kinds, basis_sets, level_list, iatom_list, jatom_list, iset_list, jset_list, ipgf_list, jpgf_list, border_mask_list, block_num_list, radius_list, rab_list, rs_grids, task_list)
Allocates a task list which can be passed to grid_collocate_task_list.
Definition grid_api.F:721
grid_api::grid_create_basis_set
subroutine, public grid_create_basis_set(nset, nsgf, maxco, maxpgf, lmin, lmax, npgf, nsgf_set, first_sgf, sphi, zet, basis_set)
Allocates a basis set which can be passed to grid_create_task_list.
Definition grid_api.F:564
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
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:613
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
offload_api
Fortran API for the offload package, which is written in C.
Definition offload_api.F:12
offload_api::offload_create_buffer
subroutine, public offload_create_buffer(length, buffer)
Allocates a buffer of given length, ie. number of elements.
Definition offload_api.F:252
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition realspace_grid_types.F:478
realspace_grid_types::rs_grid_locate_rank
pure integer function, public rs_grid_locate_rank(rs_desc, rank_in, shift)
returns the 1D rank of the task which is a cartesian shift away from 1D rank rank_in only possible if...
Definition realspace_grid_types.F:173
realspace_grid_types::rs_grid_reorder_ranks
pure subroutine, public rs_grid_reorder_ranks(desc, real2virtual)
Defines a new ordering of ranks on this realspace grid, recalculating the data bounds and reallocatin...
Definition realspace_grid_types.F:552
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition realspace_grid_types.F:2020
task_list_methods
generate the tasks lists used by collocate and integrate routines
Definition task_list_methods.F:14
task_list_methods::rs_copy_to_matrices
subroutine, public rs_copy_to_matrices(src_buffer, dest_matrices, task_list)
Copies from buffer into DBCSR matrics, replaces rs_gather_matrix for non-distributed grids.
Definition task_list_methods.F:2711
task_list_methods::rs_scatter_matrices
subroutine, public rs_scatter_matrices(src_matrices, dest_buffer, task_list, group)
Scatters dbcsr matrix blocks and receives them into a buffer as needed before collocation.
Definition task_list_methods.F:2610
task_list_methods::rs_distribute_matrix
subroutine, public rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, nimages, scatter, hmats)
redistributes the matrix so that it can be used in realspace operations i.e. according to the task li...
Definition task_list_methods.F:2236
task_list_methods::distribute_tasks
subroutine, public distribute_tasks(rs_descs, ntasks, natoms, tasks, atom_pair_send, atom_pair_recv, symmetric, reorder_rs_grid_ranks, skip_load_balance_distributed, balance_replicated)
Assembles tasks to be performed on local grid.
Definition task_list_methods.F:1962
task_list_methods::rs_gather_matrices
subroutine, public rs_gather_matrices(src_buffer, dest_matrices, task_list, group)
Gather the dbcsr matrix blocks and receives them into a buffer as needed after integration.
Definition task_list_methods.F:2651
task_list_methods::task_list_inner_loop
subroutine, public task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, dft_control, cube_info, gridlevel_info, cindex, iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, set_radius_a, set_radius_b, ra, rab, la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
...
Definition task_list_methods.F:668
task_list_methods::generate_qs_task_list
subroutine, public generate_qs_task_list(ks_env, task_list, basis_type, reorder_rs_grid_ranks, skip_load_balance_distributed, pw_env_external, sab_orb_external, ext_kpoints)
...
Definition task_list_methods.F:119
task_list_methods::rs_copy_to_buffer
subroutine, public rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
Copies the DBCSR blocks into buffer, replaces rs_scatter_matrix for non-distributed grids.
Definition task_list_methods.F:2695
task_list_types
types for task lists
Definition task_list_types.F:14
task_list_types::serialize_task
subroutine, public serialize_task(task, serialized_task)
Serialize a task into an integer array. Used for MPI communication.
Definition task_list_types.F:210
task_list_types::deserialize_task
subroutine, public deserialize_task(task, serialized_task)
De-serialize a task from an integer array. Used for MPI communication.
Definition task_list_types.F:242
task_list_types::reallocate_tasks
subroutine, public reallocate_tasks(tasks, new_size)
Grow an array of tasks while preserving the existing entries.
Definition task_list_types.F:120
task_list_types::task_size_in_int8
integer, parameter, public task_size_in_int8
Definition task_list_types.F:85
util
All kind of helpful little routines.
Definition util.F:14
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cube_utils::cube_info_type
Definition cube_utils.F:38
gaussian_gridlevels::gridlevel_info_type
Definition gaussian_gridlevels.F:34
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_comm_type
Definition message_passing.F:179
offload_api::offload_buffer_type
Definition offload_api.F:38
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
realspace_grid_types::realspace_grid_desc_p_type
Definition realspace_grid_types.F:158
realspace_grid_types::realspace_grid_desc_type
Definition realspace_grid_types.F:90
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136
task_list_types::atom_pair_type
Definition task_list_types.F:51
task_list_types::task_list_type
Definition task_list_types.F:58
task_list_types::task_type
Definition task_list_types.F:31