Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
14 : MODULE 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 : #:include './common/array_sort.fypp'
75 :
76 : IMPLICIT NONE
77 :
78 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
79 :
80 : PRIVATE
81 :
82 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'task_list_methods'
83 :
84 : PUBLIC :: generate_qs_task_list, &
85 : task_list_inner_loop
86 : PUBLIC :: distribute_tasks, &
87 : rs_distribute_matrix, &
88 : rs_scatter_matrices, &
89 : rs_gather_matrices, &
90 : rs_copy_to_buffer, &
91 : rs_copy_to_matrices
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief ...
97 : !> \param ks_env ...
98 : !> \param task_list ...
99 : !> \param basis_type ...
100 : !> \param reorder_rs_grid_ranks Flag that indicates if this routine should
101 : !> or should not overwrite the rs descriptor (see comment below)
102 : !> \param skip_load_balance_distributed ...
103 : !> \param pw_env_external ...
104 : !> \param sab_orb_external ...
105 : !> \param ext_kpoints ...
106 : !> \par History
107 : !> 01.2008 factored out of calculate_rho_elec [Joost VandeVondele]
108 : !> 04.2010 divides tasks into grid levels and atom pairs for integrate/collocate [Iain Bethune]
109 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project
110 : !> 06.2015 adjusted to be used with multiple images (k-points) [JGH]
111 : !> \note If this routine is called several times with different task lists,
112 : !> the default behaviour is to re-optimize the grid ranks and overwrite
113 : !> the rs descriptor and grids. reorder_rs_grid_ranks = .FALSE. prevents the code
114 : !> of performing a new optimization by leaving the rank order in
115 : !> its current state.
116 : ! **************************************************************************************************
117 17655 : SUBROUTINE generate_qs_task_list(ks_env, task_list, basis_type, &
118 : reorder_rs_grid_ranks, skip_load_balance_distributed, &
119 : pw_env_external, sab_orb_external, ext_kpoints)
120 :
121 : TYPE(qs_ks_env_type), POINTER :: ks_env
122 : TYPE(task_list_type), POINTER :: task_list
123 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
124 : LOGICAL, INTENT(IN) :: reorder_rs_grid_ranks, &
125 : skip_load_balance_distributed
126 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
127 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
128 : OPTIONAL, POINTER :: sab_orb_external
129 : TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
130 :
131 : CHARACTER(LEN=*), PARAMETER :: routineN = 'generate_qs_task_list'
132 : INTEGER, PARAMETER :: max_tasks = 2000
133 :
134 : INTEGER :: cindex, curr_tasks, handle, i, iatom, iatom_old, igrid_level, igrid_level_old, &
135 : ikind, ilevel, img, img_old, ipair, ipgf, iset, itask, jatom, jatom_old, jkind, jpgf, &
136 : jset, maxpgf, maxset, natoms, nimages, nkind, nseta, nsetb, slot
137 17655 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: blocks
138 : INTEGER, DIMENSION(3) :: cellind
139 17655 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
140 17655 : npgfb, nsgf
141 17655 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
142 : LOGICAL :: dokp
143 : REAL(KIND=dp) :: kind_radius_a, kind_radius_b
144 : REAL(KIND=dp), DIMENSION(3) :: ra, rab
145 17655 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
146 17655 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb
147 : TYPE(cell_type), POINTER :: cell
148 17655 : TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
149 : TYPE(dft_control_type), POINTER :: dft_control
150 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
151 17655 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
152 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, orb_basis_set
153 : TYPE(kpoint_type), POINTER :: kpoints
154 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
155 17655 : POINTER :: sab_orb
156 17655 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
157 : TYPE(pw_env_type), POINTER :: pw_env
158 17655 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
159 : TYPE(qs_kind_type), POINTER :: qs_kind
160 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
161 17655 : POINTER :: rs_descs
162 17655 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
163 17655 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
164 :
165 17655 : CALL timeset(routineN, handle)
166 :
167 : CALL get_ks_env(ks_env, &
168 : qs_kind_set=qs_kind_set, &
169 : cell=cell, &
170 : particle_set=particle_set, &
171 17655 : dft_control=dft_control)
172 :
173 17655 : CALL get_ks_env(ks_env, sab_orb=sab_orb)
174 17655 : IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
175 :
176 17655 : CALL get_ks_env(ks_env, pw_env=pw_env)
177 17655 : IF (PRESENT(pw_env_external)) pw_env => pw_env_external
178 17655 : CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_grids)
179 :
180 : ! *** assign from pw_env
181 17655 : gridlevel_info => pw_env%gridlevel_info
182 17655 : cube_info => pw_env%cube_info
183 :
184 : ! find maximum numbers
185 17655 : nkind = SIZE(qs_kind_set)
186 17655 : natoms = SIZE(particle_set)
187 17655 : maxset = 0
188 17655 : maxpgf = 0
189 48580 : DO ikind = 1, nkind
190 30925 : qs_kind => qs_kind_set(ikind)
191 : CALL get_qs_kind(qs_kind=qs_kind, &
192 30925 : basis_set=orb_basis_set, basis_type=basis_type)
193 :
194 30925 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
195 30925 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, npgf=npgfa, nset=nseta)
196 :
197 30925 : maxset = MAX(nseta, maxset)
198 112223 : maxpgf = MAX(MAXVAL(npgfa), maxpgf)
199 : END DO
200 :
201 : ! kpoint related
202 17655 : nimages = dft_control%nimages
203 17655 : dokp = .FALSE.
204 17655 : NULLIFY (cell_to_index)
205 17655 : NULLIFY (kpoints)
206 17655 : IF (PRESENT(ext_kpoints)) THEN
207 864 : IF (ASSOCIATED(ext_kpoints)) THEN
208 22 : dokp = .TRUE.
209 22 : CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
210 22 : nimages = SIZE(ext_kpoints%index_to_cell, 2)
211 : END IF
212 : END IF
213 17655 : IF (.NOT. dokp .AND. nimages > 1) THEN
214 1174 : dokp = .TRUE.
215 1174 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
216 1174 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
217 : END IF
218 :
219 : ! free the atom_pair lists if allocated
220 17655 : IF (ASSOCIATED(task_list%atom_pair_send)) DEALLOCATE (task_list%atom_pair_send)
221 17655 : IF (ASSOCIATED(task_list%atom_pair_recv)) DEALLOCATE (task_list%atom_pair_recv)
222 :
223 : ! construct a list of all tasks
224 17655 : IF (.NOT. ASSOCIATED(task_list%tasks)) THEN
225 10459 : CALL reallocate_tasks(task_list%tasks, max_tasks)
226 : END IF
227 17655 : task_list%ntasks = 0
228 17655 : curr_tasks = SIZE(task_list%tasks)
229 :
230 83890 : ALLOCATE (basis_set_list(nkind))
231 48580 : DO ikind = 1, nkind
232 30925 : qs_kind => qs_kind_set(ikind)
233 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, &
234 30925 : basis_type=basis_type)
235 48580 : IF (ASSOCIATED(basis_set_a)) THEN
236 30925 : basis_set_list(ikind)%gto_basis_set => basis_set_a
237 : ELSE
238 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
239 : END IF
240 : END DO
241 : !!$OMP PARALLEL DEFAULT(NONE) &
242 : !!$OMP SHARED (sab_orb, dokp, basis_set_list, task_list, rs_descs, dft_control, cube_info, gridlevel_info, &
243 : !!$OMP curr_tasks, maxpgf, maxset, natoms, nimages, particle_set, cell_to_index, cell) &
244 : !!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, cellind, basis_set_a, basis_set_b, ra, &
245 : !!$OMP la_max, la_min, npgfa, nseta, rpgfa, set_radius_a, kind_radius_a, zeta, &
246 : !!$OMP lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, kind_radius_b, zetb, &
247 : !!$OMP cindex, slot)
248 : ! Loop over neighbor list
249 : !!$OMP DO SCHEDULE(GUIDED)
250 1935814 : DO slot = 1, sab_orb(1)%nl_size
251 1918159 : ikind = sab_orb(1)%nlist_task(slot)%ikind
252 1918159 : jkind = sab_orb(1)%nlist_task(slot)%jkind
253 1918159 : iatom = sab_orb(1)%nlist_task(slot)%iatom
254 1918159 : jatom = sab_orb(1)%nlist_task(slot)%jatom
255 7672636 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
256 7672636 : cellind(1:3) = sab_orb(1)%nlist_task(slot)%cell(1:3)
257 :
258 1918159 : basis_set_a => basis_set_list(ikind)%gto_basis_set
259 1918159 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
260 1918159 : basis_set_b => basis_set_list(jkind)%gto_basis_set
261 1918159 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
262 1918159 : ra(:) = pbc(particle_set(iatom)%r, cell)
263 : ! basis ikind
264 1918159 : la_max => basis_set_a%lmax
265 1918159 : la_min => basis_set_a%lmin
266 1918159 : npgfa => basis_set_a%npgf
267 1918159 : nseta = basis_set_a%nset
268 1918159 : rpgfa => basis_set_a%pgf_radius
269 1918159 : set_radius_a => basis_set_a%set_radius
270 1918159 : kind_radius_a = basis_set_a%kind_radius
271 1918159 : zeta => basis_set_a%zet
272 : ! basis jkind
273 1918159 : lb_max => basis_set_b%lmax
274 1918159 : lb_min => basis_set_b%lmin
275 1918159 : npgfb => basis_set_b%npgf
276 1918159 : nsetb = basis_set_b%nset
277 1918159 : rpgfb => basis_set_b%pgf_radius
278 1918159 : set_radius_b => basis_set_b%set_radius
279 1918159 : kind_radius_b = basis_set_b%kind_radius
280 1918159 : zetb => basis_set_b%zet
281 :
282 1918159 : IF (dokp) THEN
283 574869 : cindex = cell_to_index(cellind(1), cellind(2), cellind(3))
284 : ELSE
285 1343290 : cindex = 1
286 : END IF
287 :
288 : CALL task_list_inner_loop(task_list%tasks, task_list%ntasks, curr_tasks, &
289 : rs_descs, dft_control, cube_info, gridlevel_info, cindex, &
290 : iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, &
291 : set_radius_a, set_radius_b, ra, rab, &
292 1935814 : la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
293 :
294 : END DO
295 : !!$OMP END PARALLEL
296 :
297 : ! redistribute the task list so that all tasks map on the local rs grids
298 : CALL distribute_tasks( &
299 : rs_descs=rs_descs, ntasks=task_list%ntasks, natoms=natoms, &
300 : tasks=task_list%tasks, atom_pair_send=task_list%atom_pair_send, &
301 : atom_pair_recv=task_list%atom_pair_recv, symmetric=.TRUE., &
302 : reorder_rs_grid_ranks=reorder_rs_grid_ranks, &
303 17655 : skip_load_balance_distributed=skip_load_balance_distributed)
304 :
305 : ! compute offsets for rs_scatter_matrix / rs_copy_matrix
306 52965 : ALLOCATE (nsgf(natoms))
307 17655 : CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list, nsgf=nsgf)
308 17655 : IF (ASSOCIATED(task_list%atom_pair_send)) THEN
309 : ! only needed when there is a distributed grid
310 : CALL rs_calc_offsets(pairs=task_list%atom_pair_send, &
311 : nsgf=nsgf, &
312 : group_size=rs_descs(1)%rs_desc%group_size, &
313 : pair_offsets=task_list%pair_offsets_send, &
314 : rank_offsets=task_list%rank_offsets_send, &
315 : rank_sizes=task_list%rank_sizes_send, &
316 50 : buffer_size=task_list%buffer_size_send)
317 : END IF
318 : CALL rs_calc_offsets(pairs=task_list%atom_pair_recv, &
319 : nsgf=nsgf, &
320 : group_size=rs_descs(1)%rs_desc%group_size, &
321 : pair_offsets=task_list%pair_offsets_recv, &
322 : rank_offsets=task_list%rank_offsets_recv, &
323 : rank_sizes=task_list%rank_sizes_recv, &
324 17655 : buffer_size=task_list%buffer_size_recv)
325 17655 : IF (ASSOCIATED(task_list%atom_pair_send)) THEN
326 : CALL rs_group_atom_pairs(task_list%atom_pair_send, &
327 50 : task_list%unpack_block_starts, task_list%unpack_pair_indices)
328 : ELSE
329 : CALL rs_group_atom_pairs(task_list%atom_pair_recv, &
330 17605 : task_list%unpack_block_starts, task_list%unpack_pair_indices)
331 : END IF
332 17655 : DEALLOCATE (basis_set_list, nsgf)
333 :
334 : ! If the rank order has changed, reallocate any of the distributed rs_grids
335 17655 : IF (reorder_rs_grid_ranks) THEN
336 75096 : DO i = 1, gridlevel_info%ngrid_levels
337 75096 : IF (rs_descs(i)%rs_desc%distributed) THEN
338 56 : CALL rs_grid_release(rs_grids(i))
339 56 : CALL rs_grid_create(rs_grids(i), rs_descs(i)%rs_desc)
340 : END IF
341 : END DO
342 : END IF
343 :
344 : CALL create_grid_task_list(task_list=task_list, &
345 : qs_kind_set=qs_kind_set, &
346 : particle_set=particle_set, &
347 : cell=cell, &
348 : basis_type=basis_type, &
349 17655 : rs_grids=rs_grids)
350 :
351 : ! Now we have the final list of tasks, setup the task_list with the
352 : ! data needed for the loops in integrate_v/calculate_rho
353 :
354 17655 : IF (ASSOCIATED(task_list%taskstart)) THEN
355 7196 : DEALLOCATE (task_list%taskstart)
356 : END IF
357 17655 : IF (ASSOCIATED(task_list%taskstop)) THEN
358 7196 : DEALLOCATE (task_list%taskstop)
359 : END IF
360 17655 : IF (ASSOCIATED(task_list%npairs)) THEN
361 7196 : DEALLOCATE (task_list%npairs)
362 : END IF
363 :
364 : ! First, count the number of unique atom pairs per grid level
365 :
366 52965 : ALLOCATE (task_list%npairs(SIZE(rs_descs)))
367 :
368 17655 : iatom_old = -1; jatom_old = -1; igrid_level_old = -1; img_old = -1
369 17655 : ipair = 0
370 87546 : task_list%npairs = 0
371 :
372 10687311 : DO i = 1, task_list%ntasks
373 10669656 : igrid_level = task_list%tasks(i)%grid_level
374 10669656 : img = task_list%tasks(i)%image
375 10669656 : iatom = task_list%tasks(i)%iatom
376 10669656 : jatom = task_list%tasks(i)%jatom
377 10669656 : iset = task_list%tasks(i)%iset
378 10669656 : jset = task_list%tasks(i)%jset
379 10669656 : ipgf = task_list%tasks(i)%ipgf
380 10669656 : jpgf = task_list%tasks(i)%jpgf
381 10687311 : IF (igrid_level /= igrid_level_old) THEN
382 48170 : IF (igrid_level_old /= -1) THEN
383 30854 : task_list%npairs(igrid_level_old) = ipair
384 : END IF
385 : ipair = 1
386 : igrid_level_old = igrid_level
387 : iatom_old = iatom
388 : jatom_old = jatom
389 : img_old = img
390 10621486 : ELSE IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
391 1266191 : ipair = ipair + 1
392 1266191 : iatom_old = iatom
393 1266191 : jatom_old = jatom
394 1266191 : img_old = img
395 : END IF
396 : END DO
397 : ! Take care of the last iteration
398 17655 : IF (task_list%ntasks /= 0) THEN
399 17316 : task_list%npairs(igrid_level) = ipair
400 : END IF
401 :
402 : ! Second, for each atom pair, find the indices in the task list
403 : ! of the first and last task
404 :
405 : ! Array sized for worst case
406 140172 : ALLOCATE (task_list%taskstart(MAXVAL(task_list%npairs), SIZE(rs_descs)))
407 140172 : ALLOCATE (task_list%taskstop(MAXVAL(task_list%npairs), SIZE(rs_descs)))
408 :
409 17655 : iatom_old = -1; jatom_old = -1; igrid_level_old = -1; img_old = -1
410 17655 : ipair = 0
411 2903361 : task_list%taskstart = 0
412 2903361 : task_list%taskstop = 0
413 :
414 10687311 : DO i = 1, task_list%ntasks
415 10669656 : igrid_level = task_list%tasks(i)%grid_level
416 10669656 : img = task_list%tasks(i)%image
417 10669656 : iatom = task_list%tasks(i)%iatom
418 10669656 : jatom = task_list%tasks(i)%jatom
419 10669656 : iset = task_list%tasks(i)%iset
420 10669656 : jset = task_list%tasks(i)%jset
421 10669656 : ipgf = task_list%tasks(i)%ipgf
422 10669656 : jpgf = task_list%tasks(i)%jpgf
423 10687311 : IF (igrid_level /= igrid_level_old) THEN
424 48170 : IF (igrid_level_old /= -1) THEN
425 30854 : task_list%taskstop(ipair, igrid_level_old) = i - 1
426 : END IF
427 48170 : ipair = 1
428 48170 : task_list%taskstart(ipair, igrid_level) = i
429 48170 : igrid_level_old = igrid_level
430 48170 : iatom_old = iatom
431 48170 : jatom_old = jatom
432 48170 : img_old = img
433 10621486 : ELSE IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
434 1266191 : ipair = ipair + 1
435 1266191 : task_list%taskstart(ipair, igrid_level) = i
436 1266191 : task_list%taskstop(ipair - 1, igrid_level) = i - 1
437 1266191 : iatom_old = iatom
438 1266191 : jatom_old = jatom
439 1266191 : img_old = img
440 : END IF
441 : END DO
442 : ! Take care of the last iteration
443 17655 : IF (task_list%ntasks /= 0) THEN
444 17316 : task_list%taskstop(ipair, igrid_level) = task_list%ntasks
445 : END IF
446 :
447 : ! Debug task destribution
448 : IF (debug_this_module) THEN
449 : tasks => task_list%tasks
450 : WRITE (6, *)
451 : WRITE (6, *) "Total number of tasks ", task_list%ntasks
452 : DO igrid_level = 1, gridlevel_info%ngrid_levels
453 : WRITE (6, *) "Total number of pairs(grid_level) ", &
454 : igrid_level, task_list%npairs(igrid_level)
455 : END DO
456 : WRITE (6, *)
457 :
458 : DO igrid_level = 1, gridlevel_info%ngrid_levels
459 :
460 : ALLOCATE (blocks(natoms, natoms, nimages))
461 : blocks = -1
462 : DO ipair = 1, task_list%npairs(igrid_level)
463 : itask = task_list%taskstart(ipair, igrid_level)
464 : ilevel = task_list%tasks(itask)%grid_level
465 : img = task_list%tasks(itask)%image
466 : iatom = task_list%tasks(itask)%iatom
467 : jatom = task_list%tasks(itask)%jatom
468 : iset = task_list%tasks(itask)%iset
469 : jset = task_list%tasks(itask)%jset
470 : ipgf = task_list%tasks(itask)%ipgf
471 : jpgf = task_list%tasks(itask)%jpgf
472 : IF (blocks(iatom, jatom, img) == -1 .AND. blocks(jatom, iatom, img) == -1) THEN
473 : blocks(iatom, jatom, img) = 1
474 : blocks(jatom, iatom, img) = 1
475 : ELSE
476 : WRITE (6, *) "TASK LIST CONFLICT IN PAIR ", ipair
477 : WRITE (6, *) "Reuse of iatom, jatom, image ", iatom, jatom, img
478 : END IF
479 :
480 : iatom_old = iatom
481 : jatom_old = jatom
482 : img_old = img
483 : DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
484 : ilevel = task_list%tasks(itask)%grid_level
485 : img = task_list%tasks(itask)%image
486 : iatom = task_list%tasks(itask)%iatom
487 : jatom = task_list%tasks(itask)%jatom
488 : iset = task_list%tasks(itask)%iset
489 : jset = task_list%tasks(itask)%jset
490 : ipgf = task_list%tasks(itask)%ipgf
491 : jpgf = task_list%tasks(itask)%jpgf
492 : IF (iatom /= iatom_old .OR. jatom /= jatom_old .OR. img /= img_old) THEN
493 : WRITE (6, *) "TASK LIST CONFLICT IN TASK ", itask
494 : WRITE (6, *) "Inconsistent iatom, jatom, image ", iatom, jatom, img
495 : WRITE (6, *) "Should be iatom, jatom, image ", iatom_old, jatom_old, img_old
496 : END IF
497 :
498 : END DO
499 : END DO
500 : DEALLOCATE (blocks)
501 :
502 : END DO
503 :
504 : END IF
505 :
506 17655 : CALL timestop(handle)
507 :
508 17655 : END SUBROUTINE generate_qs_task_list
509 :
510 : ! **************************************************************************************************
511 : !> \brief Sends the task list data to the grid API.
512 : !> \author Ole Schuett
513 : ! **************************************************************************************************
514 17655 : SUBROUTINE create_grid_task_list(task_list, qs_kind_set, particle_set, cell, basis_type, rs_grids)
515 : TYPE(task_list_type), POINTER :: task_list
516 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
517 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
518 : TYPE(cell_type), POINTER :: cell
519 : CHARACTER(LEN=*) :: basis_type
520 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
521 :
522 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
523 : INTEGER :: nset, natoms, nkinds, ntasks, &
524 : ikind, iatom, itask, nsgf
525 17655 : INTEGER, DIMENSION(:), ALLOCATABLE :: atom_kinds, level_list, iatom_list, jatom_list, &
526 17655 : iset_list, jset_list, ipgf_list, jpgf_list, &
527 17655 : border_mask_list, block_num_list
528 17655 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: radius_list
529 17655 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: rab_list, atom_positions
530 17655 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
531 17655 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
532 17655 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi, zet
533 17655 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
534 :
535 17655 : nkinds = SIZE(qs_kind_set)
536 17655 : natoms = SIZE(particle_set)
537 17655 : ntasks = task_list%ntasks
538 17655 : tasks => task_list%tasks
539 :
540 17655 : IF (.NOT. ASSOCIATED(task_list%grid_basis_sets)) THEN
541 : ! Basis sets do not change during simulation - only need to create them once.
542 49540 : ALLOCATE (task_list%grid_basis_sets(nkinds))
543 28622 : DO ikind = 1, nkinds
544 18163 : CALL get_qs_kind(qs_kind_set(ikind), basis_type=basis_type, basis_set=orb_basis_set)
545 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
546 : nset=nset, &
547 : nsgf=nsgf, &
548 : nsgf_set=nsgf_set, &
549 : npgf=npgf, &
550 : first_sgf=first_sgf, &
551 : lmax=lmax, &
552 : lmin=lmin, &
553 : sphi=sphi, &
554 18163 : zet=zet)
555 : CALL grid_create_basis_set(nset=nset, &
556 : nsgf=nsgf, &
557 : maxco=SIZE(sphi, 1), &
558 : maxpgf=SIZE(zet, 1), &
559 : lmin=lmin, &
560 : lmax=lmax, &
561 : npgf=npgf, &
562 : nsgf_set=nsgf_set, &
563 : first_sgf=first_sgf, &
564 : sphi=sphi, &
565 : zet=zet, &
566 28622 : basis_set=task_list%grid_basis_sets(ikind))
567 : END DO
568 : END IF
569 :
570 : ! Pack task list infos
571 88275 : ALLOCATE (atom_kinds(natoms), atom_positions(3, natoms))
572 81060 : DO iatom = 1, natoms
573 63405 : atom_kinds(iatom) = particle_set(iatom)%atomic_kind%kind_number
574 81060 : atom_positions(:, iatom) = pbc(particle_set(iatom)%r, cell)
575 : END DO
576 :
577 87258 : ALLOCATE (level_list(ntasks), iatom_list(ntasks), jatom_list(ntasks))
578 86919 : ALLOCATE (iset_list(ntasks), jset_list(ntasks), ipgf_list(ntasks), jpgf_list(ntasks))
579 52287 : ALLOCATE (border_mask_list(ntasks), block_num_list(ntasks))
580 87597 : ALLOCATE (radius_list(ntasks), rab_list(3, ntasks))
581 :
582 10687311 : DO itask = 1, ntasks
583 10669656 : level_list(itask) = tasks(itask)%grid_level
584 10669656 : iatom_list(itask) = tasks(itask)%iatom
585 10669656 : jatom_list(itask) = tasks(itask)%jatom
586 10669656 : iset_list(itask) = tasks(itask)%iset
587 10669656 : jset_list(itask) = tasks(itask)%jset
588 10669656 : ipgf_list(itask) = tasks(itask)%ipgf
589 10669656 : jpgf_list(itask) = tasks(itask)%jpgf
590 10669656 : IF (tasks(itask)%dist_type == 2) THEN
591 0 : border_mask_list(itask) = IAND(63, NOT(tasks(itask)%subpatch_pattern)) ! invert last 6 bits
592 : ELSE
593 10669656 : border_mask_list(itask) = 0 ! no masking
594 : END IF
595 10669656 : block_num_list(itask) = tasks(itask)%pair_index ! change of nomenclature pair_index -> block_num
596 10669656 : radius_list(itask) = tasks(itask)%radius
597 42696279 : rab_list(:, itask) = tasks(itask)%rab(:)
598 : END DO
599 :
600 : CALL grid_create_task_list(ntasks=ntasks, &
601 : natoms=natoms, &
602 : nkinds=nkinds, &
603 : nblocks=SIZE(task_list%pair_offsets_recv), &
604 : block_offsets=task_list%pair_offsets_recv, &
605 : atom_positions=atom_positions, &
606 : atom_kinds=atom_kinds, &
607 : basis_sets=task_list%grid_basis_sets, &
608 : level_list=level_list, &
609 : iatom_list=iatom_list, &
610 : jatom_list=jatom_list, &
611 : iset_list=iset_list, &
612 : jset_list=jset_list, &
613 : ipgf_list=ipgf_list, &
614 : jpgf_list=jpgf_list, &
615 : border_mask_list=border_mask_list, &
616 : block_num_list=block_num_list, &
617 : radius_list=radius_list, &
618 : rab_list=rab_list, &
619 : rs_grids=rs_grids, &
620 17655 : task_list=task_list%grid_task_list)
621 :
622 17655 : CALL offload_create_buffer(task_list%buffer_size_recv, task_list%pab_buffer)
623 17655 : CALL offload_create_buffer(task_list%buffer_size_recv, task_list%hab_buffer)
624 :
625 35310 : END SUBROUTINE create_grid_task_list
626 :
627 : ! **************************************************************************************************
628 : !> \brief ...
629 : !> \param tasks ...
630 : !> \param ntasks ...
631 : !> \param curr_tasks ...
632 : !> \param rs_descs ...
633 : !> \param dft_control ...
634 : !> \param cube_info ...
635 : !> \param gridlevel_info ...
636 : !> \param cindex ...
637 : !> \param iatom ...
638 : !> \param jatom ...
639 : !> \param rpgfa ...
640 : !> \param rpgfb ...
641 : !> \param zeta ...
642 : !> \param zetb ...
643 : !> \param kind_radius_b ...
644 : !> \param set_radius_a ...
645 : !> \param set_radius_b ...
646 : !> \param ra ...
647 : !> \param rab ...
648 : !> \param la_max ...
649 : !> \param la_min ...
650 : !> \param lb_max ...
651 : !> \param lb_min ...
652 : !> \param npgfa ...
653 : !> \param npgfb ...
654 : !> \param nseta ...
655 : !> \param nsetb ...
656 : !> \par History
657 : !> Joost VandeVondele: 10.2008 refactored
658 : ! **************************************************************************************************
659 1929220 : SUBROUTINE task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, dft_control, &
660 : cube_info, gridlevel_info, cindex, &
661 : iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, set_radius_a, set_radius_b, ra, rab, &
662 : la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
663 :
664 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
665 : INTEGER :: ntasks, curr_tasks
666 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
667 : POINTER :: rs_descs
668 : TYPE(dft_control_type), POINTER :: dft_control
669 : TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
670 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
671 : INTEGER :: cindex, iatom, jatom
672 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, zeta, zetb
673 : REAL(KIND=dp) :: kind_radius_b
674 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
675 : REAL(KIND=dp), DIMENSION(3) :: ra, rab
676 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
677 : npgfb
678 : INTEGER :: nseta, nsetb
679 :
680 : INTEGER :: cube_center(3), igrid_level, ipgf, iset, &
681 : jpgf, jset, lb_cube(3), ub_cube(3)
682 : REAL(KIND=dp) :: dab, rab2, radius, zetp
683 :
684 1929220 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
685 1929220 : dab = SQRT(rab2)
686 :
687 5704238 : loop_iset: DO iset = 1, nseta
688 :
689 3775018 : IF (set_radius_a(iset) + kind_radius_b < dab) CYCLE loop_iset
690 :
691 11377116 : loop_jset: DO jset = 1, nsetb
692 :
693 6695318 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE loop_jset
694 :
695 19345993 : loop_ipgf: DO ipgf = 1, npgfa(iset)
696 :
697 11271771 : IF (rpgfa(ipgf, iset) + set_radius_b(jset) < dab) CYCLE loop_ipgf
698 :
699 33771488 : loop_jpgf: DO jpgf = 1, npgfb(jset)
700 :
701 20414466 : IF (rpgfa(ipgf, iset) + rpgfb(jpgf, jset) < dab) CYCLE loop_jpgf
702 :
703 10823718 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
704 10823718 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
705 :
706 : CALL compute_pgf_properties(cube_center, lb_cube, ub_cube, radius, &
707 : rs_descs(igrid_level)%rs_desc, cube_info(igrid_level), &
708 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
709 : lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
710 10823718 : ra, rab, rab2, dft_control%qs_control%eps_rho_rspace)
711 :
712 : CALL pgf_to_tasks(tasks, ntasks, curr_tasks, &
713 : rab, cindex, iatom, jatom, iset, jset, ipgf, jpgf, &
714 : la_max(iset), lb_max(jset), rs_descs(igrid_level)%rs_desc, &
715 : igrid_level, gridlevel_info%ngrid_levels, cube_center, &
716 31686237 : lb_cube, ub_cube, radius)
717 :
718 : END DO loop_jpgf
719 :
720 : END DO loop_ipgf
721 :
722 : END DO loop_jset
723 :
724 : END DO loop_iset
725 :
726 1929220 : END SUBROUTINE task_list_inner_loop
727 :
728 : ! **************************************************************************************************
729 : !> \brief combines the calculation of several basic properties of a given pgf:
730 : !> its center, the bounding cube, the radius, the cost,
731 : !> tries to predict the time needed for processing this task
732 : !> in this way an improved load balance might be obtained
733 : !> \param cube_center ...
734 : !> \param lb_cube ...
735 : !> \param ub_cube ...
736 : !> \param radius ...
737 : !> \param rs_desc ...
738 : !> \param cube_info ...
739 : !> \param la_max ...
740 : !> \param zeta ...
741 : !> \param la_min ...
742 : !> \param lb_max ...
743 : !> \param zetb ...
744 : !> \param lb_min ...
745 : !> \param ra ...
746 : !> \param rab ...
747 : !> \param rab2 ...
748 : !> \param eps ...
749 : !> \par History
750 : !> 10.2008 refactored [Joost VandeVondele]
751 : !> \note
752 : !> -) this requires the radius to be computed in the same way as
753 : !> collocate_pgf_product, we should factor that part into a subroutine
754 : !> -) we're assuming that integrate_pgf and collocate_pgf are the same cost for load balancing
755 : !> this is more or less true for map_consistent
756 : !> -) in principle, the computed radius could be recycled in integrate_pgf/collocate_pgf if it is certainly
757 : !> the same, this could lead to a small speedup
758 : !> -) 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)
759 : !> fitting the measured data on an opteron/g95 using the expression
760 : !> a*(l+b)(r+c)**3+d which is based on the innerloop of the collocating routines
761 : ! **************************************************************************************************
762 10823718 : SUBROUTINE compute_pgf_properties(cube_center, lb_cube, ub_cube, radius, &
763 : rs_desc, cube_info, la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, rab2, eps)
764 :
765 : INTEGER, DIMENSION(3), INTENT(OUT) :: cube_center, lb_cube, ub_cube
766 : REAL(KIND=dp), INTENT(OUT) :: radius
767 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
768 : TYPE(cube_info_type), INTENT(IN) :: cube_info
769 : INTEGER, INTENT(IN) :: la_max
770 : REAL(KIND=dp), INTENT(IN) :: zeta
771 : INTEGER, INTENT(IN) :: la_min, lb_max
772 : REAL(KIND=dp), INTENT(IN) :: zetb
773 : INTEGER, INTENT(IN) :: lb_min
774 : REAL(KIND=dp), INTENT(IN) :: ra(3), rab(3), rab2, eps
775 :
776 : INTEGER :: extent(3)
777 10823718 : INTEGER, DIMENSION(:), POINTER :: sphere_bounds
778 : REAL(KIND=dp) :: cutoff, f, prefactor, rb(3), zetp
779 : REAL(KIND=dp), DIMENSION(3) :: rp
780 :
781 : ! the radius for this task
782 :
783 10823718 : zetp = zeta + zetb
784 43294872 : rp(:) = ra(:) + zetb/zetp*rab(:)
785 43294872 : rb(:) = ra(:) + rab(:)
786 10823718 : cutoff = 1.0_dp
787 10823718 : f = zetb/zetp
788 10823718 : prefactor = EXP(-zeta*f*rab2)
789 : radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, &
790 10823718 : zetp=zetp, eps=eps, prefactor=prefactor, cutoff=cutoff)
791 :
792 10823718 : CALL compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab)
793 : ! compute cube_center, the center of the gaussian product to map (folded to within the unit cell)
794 43294872 : cube_center(:) = MODULO(cube_center(:), rs_desc%npts(:))
795 43294872 : cube_center(:) = cube_center(:) + rs_desc%lb(:)
796 :
797 10823718 : IF (rs_desc%orthorhombic) THEN
798 10092425 : CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
799 : ELSE
800 731293 : CALL return_cube_nonortho(cube_info, radius, lb_cube, ub_cube, rp)
801 : !? unclear if extent is computed correctly.
802 2925172 : extent(:) = ub_cube(:) - lb_cube(:)
803 2925172 : lb_cube(:) = -extent(:)/2 - 1
804 2925172 : ub_cube(:) = extent(:)/2
805 : END IF
806 :
807 10823718 : END SUBROUTINE compute_pgf_properties
808 : ! **************************************************************************************************
809 : !> \brief predicts the cost of a task in kcycles for a given task
810 : !> the model is based on a fit of actual data, and might need updating
811 : !> as collocate_pgf_product changes (or CPUs/compilers change)
812 : !> maybe some dynamic approach, improving the cost model on the fly could
813 : !> work as well
814 : !> the cost model does not yet take into account the fraction of space
815 : !> that is mapped locally for a given cube and rs_grid (generalised tasks)
816 : !> \param lb_cube ...
817 : !> \param ub_cube ...
818 : !> \param fraction ...
819 : !> \param lmax ...
820 : !> \param is_ortho ...
821 : !> \return ...
822 : ! **************************************************************************************************
823 10823718 : INTEGER FUNCTION cost_model(lb_cube, ub_cube, fraction, lmax, is_ortho)
824 : INTEGER, DIMENSION(3), INTENT(IN) :: lb_cube, ub_cube
825 : REAL(KIND=dp), INTENT(IN) :: fraction
826 : INTEGER :: lmax
827 : LOGICAL :: is_ortho
828 :
829 : INTEGER :: cmax
830 : REAL(KIND=dp) :: v1, v2, v3, v4, v5
831 :
832 43294872 : cmax = MAXVAL(((ub_cube - lb_cube) + 1)/2)
833 :
834 10823718 : IF (is_ortho) THEN
835 : v1 = 1.504760E+00_dp
836 : v2 = 3.126770E+00_dp
837 : v3 = 5.074106E+00_dp
838 : v4 = 1.091568E+00_dp
839 : v5 = 1.070187E+00_dp
840 : ELSE
841 731293 : v1 = 7.831105E+00_dp
842 731293 : v2 = 2.675174E+00_dp
843 731293 : v3 = 7.546553E+00_dp
844 731293 : v4 = 6.122446E-01_dp
845 731293 : v5 = 3.886382E+00_dp
846 : END IF
847 10823718 : cost_model = CEILING(((lmax + v1)*(cmax + v2)**3*v3*fraction + v4 + v5*lmax**7)/1000.0_dp)
848 :
849 10823718 : END FUNCTION cost_model
850 : ! **************************************************************************************************
851 : !> \brief pgf_to_tasks converts a given pgf to one or more tasks, in particular
852 : !> this determines by which CPUs a given pgf gets collocated
853 : !> the format of the task array is as follows
854 : !> tasks(1,i) := destination
855 : !> tasks(2,i) := source
856 : !> tasks(3,i) := compressed type (iatom, jatom, ....)
857 : !> tasks(4,i) := type (0: replicated, 1: distributed local, 2: distributed generalised)
858 : !> tasks(5,i) := cost
859 : !> tasks(6,i) := alternate destination code (0 if none available)
860 : !>
861 : !> \param tasks ...
862 : !> \param ntasks ...
863 : !> \param curr_tasks ...
864 : !> \param rab ...
865 : !> \param cindex ...
866 : !> \param iatom ...
867 : !> \param jatom ...
868 : !> \param iset ...
869 : !> \param jset ...
870 : !> \param ipgf ...
871 : !> \param jpgf ...
872 : !> \param la_max ...
873 : !> \param lb_max ...
874 : !> \param rs_desc ...
875 : !> \param igrid_level ...
876 : !> \param n_levels ...
877 : !> \param cube_center ...
878 : !> \param lb_cube ...
879 : !> \param ub_cube ...
880 : !> \par History
881 : !> 10.2008 Refactored based on earlier routines by MattW [Joost VandeVondele]
882 : ! **************************************************************************************************
883 10823718 : SUBROUTINE pgf_to_tasks(tasks, ntasks, curr_tasks, &
884 : rab, cindex, iatom, jatom, iset, jset, ipgf, jpgf, &
885 : la_max, lb_max, rs_desc, igrid_level, n_levels, &
886 : cube_center, lb_cube, ub_cube, radius)
887 :
888 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
889 : INTEGER, INTENT(INOUT) :: ntasks, curr_tasks
890 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
891 : INTEGER, INTENT(IN) :: cindex, iatom, jatom, iset, jset, ipgf, &
892 : jpgf, la_max, lb_max
893 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
894 : INTEGER, INTENT(IN) :: igrid_level, n_levels
895 : INTEGER, DIMENSION(3), INTENT(IN) :: cube_center, lb_cube, ub_cube
896 : REAL(KIND=dp), INTENT(IN) :: radius
897 :
898 : INTEGER, PARAMETER :: add_tasks = 1000
899 : REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp
900 :
901 : INTEGER :: added_tasks, cost, j, lmax
902 : LOGICAL :: is_ortho
903 : REAL(KIND=dp) :: tfraction
904 :
905 : !$OMP SINGLE
906 10823718 : ntasks = ntasks + 1
907 10823718 : IF (ntasks > curr_tasks) THEN
908 1346 : curr_tasks = INT((curr_tasks + add_tasks)*mult_tasks)
909 1346 : CALL reallocate_tasks(tasks, curr_tasks)
910 : END IF
911 : !$OMP END SINGLE
912 :
913 10823718 : IF (rs_desc%distributed) THEN
914 :
915 : ! finds the node(s) that need to process this task
916 : ! on exit tasks(:)%dist_type is 1 for distributed tasks and 2 for generalised tasks
917 : CALL rs_find_node(rs_desc, igrid_level, n_levels, cube_center, &
918 1476 : ntasks=ntasks, tasks=tasks, lb_cube=lb_cube, ub_cube=ub_cube, added_tasks=added_tasks)
919 :
920 : ELSE
921 10822242 : tasks(ntasks)%destination = encode_rank(rs_desc%my_pos, igrid_level, n_levels)
922 10822242 : tasks(ntasks)%dist_type = 0
923 10822242 : tasks(ntasks)%subpatch_pattern = 0
924 10822242 : added_tasks = 1
925 : END IF
926 :
927 10823718 : lmax = la_max + lb_max
928 10823718 : is_ortho = (tasks(ntasks)%dist_type == 0 .OR. tasks(ntasks)%dist_type == 1) .AND. rs_desc%orthorhombic
929 : ! we assume the load is shared equally between processes dealing with a generalised Gaussian.
930 : ! this could be refined in the future
931 10823718 : tfraction = 1.0_dp/added_tasks
932 :
933 10823718 : cost = cost_model(lb_cube, ub_cube, tfraction, lmax, is_ortho)
934 :
935 21647436 : DO j = 1, added_tasks
936 10823718 : tasks(ntasks - added_tasks + j)%source = encode_rank(rs_desc%my_pos, igrid_level, n_levels)
937 10823718 : tasks(ntasks - added_tasks + j)%cost = cost
938 10823718 : tasks(ntasks - added_tasks + j)%grid_level = igrid_level
939 10823718 : tasks(ntasks - added_tasks + j)%image = cindex
940 10823718 : tasks(ntasks - added_tasks + j)%iatom = iatom
941 10823718 : tasks(ntasks - added_tasks + j)%jatom = jatom
942 10823718 : tasks(ntasks - added_tasks + j)%iset = iset
943 10823718 : tasks(ntasks - added_tasks + j)%jset = jset
944 10823718 : tasks(ntasks - added_tasks + j)%ipgf = ipgf
945 10823718 : tasks(ntasks - added_tasks + j)%jpgf = jpgf
946 43294872 : tasks(ntasks - added_tasks + j)%rab = rab
947 21647436 : tasks(ntasks - added_tasks + j)%radius = radius
948 : END DO
949 :
950 10823718 : END SUBROUTINE pgf_to_tasks
951 :
952 : ! **************************************************************************************************
953 : !> \brief performs load balancing of the tasks on the distributed grids
954 : !> \param tasks ...
955 : !> \param ntasks ...
956 : !> \param rs_descs ...
957 : !> \param grid_level ...
958 : !> \param natoms ...
959 : !> \par History
960 : !> created 2008-10-03 [Joost VandeVondele]
961 : ! **************************************************************************************************
962 56 : SUBROUTINE load_balance_distributed(tasks, ntasks, rs_descs, grid_level, natoms)
963 :
964 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
965 : INTEGER :: ntasks
966 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
967 : POINTER :: rs_descs
968 : INTEGER :: grid_level, natoms
969 :
970 : CHARACTER(LEN=*), PARAMETER :: routineN = 'load_balance_distributed'
971 :
972 : INTEGER :: handle
973 56 : INTEGER, DIMENSION(:, :, :), POINTER :: list
974 :
975 56 : CALL timeset(routineN, handle)
976 :
977 56 : NULLIFY (list)
978 : ! here we create for each cpu (0:ncpu-1) a list of possible destinations.
979 : ! if a destination would not be in this list, it is a bug
980 56 : CALL create_destination_list(list, rs_descs, grid_level)
981 :
982 : ! now, walk over the tasks, filling in the loads of each destination
983 56 : CALL compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list=.TRUE.)
984 :
985 : ! optimize loads & fluxes
986 56 : CALL optimize_load_list(list, rs_descs(1)%rs_desc%group, rs_descs(1)%rs_desc%my_pos)
987 :
988 : ! now, walk over the tasks, using the list to set the destinations
989 56 : CALL compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list=.FALSE.)
990 :
991 56 : DEALLOCATE (list)
992 :
993 56 : CALL timestop(handle)
994 :
995 56 : END SUBROUTINE load_balance_distributed
996 :
997 : ! **************************************************************************************************
998 : !> \brief this serial routine adjusts the fluxes in the global list
999 : !>
1000 : !> \param list_global ...
1001 : !> \par History
1002 : !> created 2008-10-06 [Joost VandeVondele]
1003 : ! **************************************************************************************************
1004 28 : SUBROUTINE balance_global_list(list_global)
1005 : INTEGER, DIMENSION(:, :, 0:) :: list_global
1006 :
1007 : CHARACTER(LEN=*), PARAMETER :: routineN = 'balance_global_list'
1008 : INTEGER, PARAMETER :: Max_Iter = 100
1009 : REAL(KIND=dp), PARAMETER :: Tolerance_factor = 0.005_dp
1010 :
1011 : INTEGER :: dest, handle, icpu, idest, iflux, &
1012 : ilocal, k, maxdest, Ncpu, Nflux
1013 28 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: flux_connections
1014 : LOGICAL :: solution_optimal
1015 : REAL(KIND=dp) :: average, load_shift, max_load_shift, &
1016 : tolerance
1017 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: load, optimized_flux, optimized_load
1018 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: flux_limits
1019 :
1020 28 : CALL timeset(routineN, handle)
1021 :
1022 28 : Ncpu = SIZE(list_global, 3)
1023 28 : maxdest = SIZE(list_global, 2)
1024 84 : ALLOCATE (load(0:Ncpu - 1))
1025 28 : load = 0.0_dp
1026 56 : ALLOCATE (optimized_load(0:Ncpu - 1))
1027 :
1028 : ! figure out the number of fluxes
1029 : ! we assume that the global_list is symmetric
1030 84 : Nflux = 0
1031 84 : DO icpu = 0, ncpu - 1
1032 196 : DO idest = 1, maxdest
1033 112 : dest = list_global(1, idest, icpu)
1034 168 : IF (dest < ncpu .AND. dest > icpu) Nflux = Nflux + 1
1035 : END DO
1036 : END DO
1037 84 : ALLOCATE (optimized_flux(Nflux))
1038 84 : ALLOCATE (flux_limits(2, Nflux))
1039 84 : ALLOCATE (flux_connections(2, Nflux))
1040 :
1041 : ! reorder data
1042 28 : flux_limits = 0
1043 28 : Nflux = 0
1044 84 : DO icpu = 0, ncpu - 1
1045 168 : load(icpu) = SUM(list_global(2, :, icpu))
1046 196 : DO idest = 1, maxdest
1047 112 : dest = list_global(1, idest, icpu)
1048 168 : IF (dest < ncpu) THEN
1049 112 : IF (dest /= icpu) THEN
1050 56 : IF (dest > icpu) THEN
1051 28 : Nflux = Nflux + 1
1052 28 : flux_limits(2, Nflux) = list_global(2, idest, icpu)
1053 28 : flux_connections(1, Nflux) = icpu
1054 28 : flux_connections(2, Nflux) = dest
1055 : ELSE
1056 28 : DO iflux = 1, Nflux
1057 28 : IF (flux_connections(1, iflux) == dest .AND. flux_connections(2, iflux) == icpu) THEN
1058 28 : flux_limits(1, iflux) = -list_global(2, idest, icpu)
1059 28 : EXIT
1060 : END IF
1061 : END DO
1062 : END IF
1063 : END IF
1064 : END IF
1065 : END DO
1066 : END DO
1067 :
1068 28 : solution_optimal = .FALSE.
1069 28 : optimized_flux = 0.0_dp
1070 :
1071 : ! an iterative solver, if iterated till convergence the maximum load is minimal
1072 : ! we terminate before things are fully converged, since this does show up in the timings
1073 : ! once the largest shift becomes less than a small fraction of the average load, we're done
1074 : ! we're perfectly happy if the load balance is within 1 percent or so
1075 : ! the maximum load normally converges even faster
1076 84 : average = SUM(load)/SIZE(load)
1077 28 : tolerance = Tolerance_factor*average
1078 :
1079 84 : optimized_load(:) = load
1080 333 : DO k = 1, Max_iter
1081 : max_load_shift = 0.0_dp
1082 660 : DO iflux = 1, Nflux
1083 330 : load_shift = (optimized_load(flux_connections(1, iflux)) - optimized_load(flux_connections(2, iflux)))/2
1084 330 : load_shift = MAX(flux_limits(1, iflux) - optimized_flux(iflux), load_shift)
1085 330 : load_shift = MIN(flux_limits(2, iflux) - optimized_flux(iflux), load_shift)
1086 330 : max_load_shift = MAX(ABS(load_shift), max_load_shift)
1087 330 : optimized_load(flux_connections(1, iflux)) = optimized_load(flux_connections(1, iflux)) - load_shift
1088 330 : optimized_load(flux_connections(2, iflux)) = optimized_load(flux_connections(2, iflux)) + load_shift
1089 660 : optimized_flux(iflux) = optimized_flux(iflux) + load_shift
1090 : END DO
1091 333 : IF (max_load_shift < tolerance) THEN
1092 84 : solution_optimal = .TRUE.
1093 : EXIT
1094 : END IF
1095 : END DO
1096 :
1097 : ! now adjust the load list to reflect the optimized fluxes
1098 : ! reorder data
1099 : Nflux = 0
1100 84 : DO icpu = 0, ncpu - 1
1101 168 : DO idest = 1, maxdest
1102 168 : IF (list_global(1, idest, icpu) == icpu) ilocal = idest
1103 : END DO
1104 196 : DO idest = 1, maxdest
1105 112 : dest = list_global(1, idest, icpu)
1106 168 : IF (dest < ncpu) THEN
1107 112 : IF (dest /= icpu) THEN
1108 56 : IF (dest > icpu) THEN
1109 28 : Nflux = Nflux + 1
1110 28 : IF (optimized_flux(Nflux) > 0) THEN
1111 : list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1112 0 : list_global(2, idest, icpu) - NINT(optimized_flux(Nflux))
1113 0 : list_global(2, idest, icpu) = NINT(optimized_flux(Nflux))
1114 : ELSE
1115 : list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1116 28 : list_global(2, idest, icpu)
1117 28 : list_global(2, idest, icpu) = 0
1118 : END IF
1119 : ELSE
1120 28 : DO iflux = 1, Nflux
1121 28 : IF (flux_connections(1, iflux) == dest .AND. flux_connections(2, iflux) == icpu) THEN
1122 28 : IF (optimized_flux(iflux) > 0) THEN
1123 : list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1124 0 : list_global(2, idest, icpu)
1125 0 : list_global(2, idest, icpu) = 0
1126 : ELSE
1127 : list_global(2, ilocal, icpu) = list_global(2, ilocal, icpu) + &
1128 28 : list_global(2, idest, icpu) + NINT(optimized_flux(iflux))
1129 28 : list_global(2, idest, icpu) = -NINT(optimized_flux(iflux))
1130 : END IF
1131 : EXIT
1132 : END IF
1133 : END DO
1134 : END IF
1135 : END IF
1136 : END IF
1137 : END DO
1138 : END DO
1139 :
1140 28 : CALL timestop(handle)
1141 :
1142 56 : END SUBROUTINE balance_global_list
1143 :
1144 : ! **************************************************************************************************
1145 : !> \brief this routine gets back optimized loads for all destinations
1146 : !>
1147 : !> \param list ...
1148 : !> \param group ...
1149 : !> \param my_pos ...
1150 : !> \par History
1151 : !> created 2008-10-06 [Joost VandeVondele]
1152 : !> Modified 2016-01 [EPCC] Reduce memory requirements on P processes
1153 : !> from O(P^2) to O(P)
1154 : ! **************************************************************************************************
1155 56 : SUBROUTINE optimize_load_list(list, group, my_pos)
1156 : INTEGER, DIMENSION(:, :, 0:) :: list
1157 : TYPE(mp_comm_type), INTENT(IN) :: group
1158 : INTEGER, INTENT(IN) :: my_pos
1159 :
1160 : CHARACTER(LEN=*), PARAMETER :: routineN = 'optimize_load_list'
1161 : INTEGER, PARAMETER :: rank_of_root = 0
1162 :
1163 : INTEGER :: handle, icpu, idest, maxdest, ncpu
1164 : INTEGER, ALLOCATABLE, DIMENSION(:) :: load_all
1165 56 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: load_partial
1166 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: list_global
1167 :
1168 56 : CALL timeset(routineN, handle)
1169 :
1170 56 : ncpu = SIZE(list, 3)
1171 56 : maxdest = SIZE(list, 2)
1172 :
1173 : !find total workload ...
1174 168 : ALLOCATE (load_all(maxdest*ncpu))
1175 112 : load_all(:) = RESHAPE(list(2, :, :), [maxdest*ncpu])
1176 56 : CALL group%sum(load_all(:), rank_of_root)
1177 :
1178 : ! ... and optimise the work per process
1179 224 : ALLOCATE (list_global(2, maxdest, ncpu))
1180 56 : IF (rank_of_root == my_pos) THEN
1181 196 : list_global(1, :, :) = list(1, :, :)
1182 252 : list_global(2, :, :) = RESHAPE(load_all, [maxdest, ncpu])
1183 28 : CALL balance_global_list(list_global)
1184 : END IF
1185 56 : CALL group%bcast(list_global, rank_of_root)
1186 :
1187 : !figure out how much can be sent to other processes
1188 224 : ALLOCATE (load_partial(maxdest, ncpu))
1189 : ! send 'load_all', which is a copy of 'list' (but without leading dimension/stride)
1190 168 : CALL group%sum_partial(RESHAPE(load_all, [maxdest, ncpu]), load_partial(:, :))
1191 :
1192 168 : DO icpu = 1, ncpu
1193 392 : DO idest = 1, maxdest
1194 :
1195 : !need to deduct 1 because `list' was passed in to this routine as being indexed from zero
1196 336 : IF (load_partial(idest, icpu) > list_global(2, idest, icpu)) THEN
1197 38 : IF (load_partial(idest, icpu) - list(2, idest, icpu - 1) < list_global(2, idest, icpu)) THEN
1198 : list(2, idest, icpu - 1) = list_global(2, idest, icpu) &
1199 10 : - (load_partial(idest, icpu) - list(2, idest, icpu - 1))
1200 : ELSE
1201 28 : list(2, idest, icpu - 1) = 0
1202 : END IF
1203 : END IF
1204 :
1205 : END DO
1206 : END DO
1207 :
1208 : !clean up before leaving
1209 56 : DEALLOCATE (load_all)
1210 56 : DEALLOCATE (list_global)
1211 56 : DEALLOCATE (load_partial)
1212 :
1213 56 : CALL timestop(handle)
1214 56 : END SUBROUTINE optimize_load_list
1215 :
1216 : ! **************************************************************************************************
1217 : !> \brief fill the load list with values derived from the tasks array
1218 : !> from the alternate locations, we select the alternate location that
1219 : !> can be used without increasing the number of matrix blocks needed to
1220 : !> distribute.
1221 : !> Replicated tasks are not yet considered
1222 : !>
1223 : !> \param list ...
1224 : !> \param rs_descs ...
1225 : !> \param grid_level ...
1226 : !> \param tasks ...
1227 : !> \param ntasks ...
1228 : !> \param natoms ...
1229 : !> \param create_list ...
1230 : !> \par History
1231 : !> created 2008-10-06 [Joost VandeVondele]
1232 : ! **************************************************************************************************
1233 112 : SUBROUTINE compute_load_list(list, rs_descs, grid_level, tasks, ntasks, natoms, create_list)
1234 : INTEGER, DIMENSION(:, :, 0:) :: list
1235 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1236 : POINTER :: rs_descs
1237 : INTEGER :: grid_level
1238 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1239 : INTEGER :: ntasks, natoms
1240 : LOGICAL :: create_list
1241 :
1242 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_load_list'
1243 :
1244 : INTEGER :: cost, dest, handle, i, iatom, ilevel, img, img_old, iopt, ipgf, iset, itask, &
1245 : itask_start, itask_stop, jatom, jpgf, jset, li, maxdest, ncpu, ndest_pair, nopt, nshort, &
1246 : rank
1247 : INTEGER(KIND=int_8) :: bit_pattern, ipair, ipair_old, natom8
1248 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: loads
1249 112 : INTEGER, ALLOCATABLE, DIMENSION(:) :: all_dests, index
1250 : INTEGER, DIMENSION(6) :: options
1251 :
1252 112 : CALL timeset(routineN, handle)
1253 :
1254 336 : ALLOCATE (loads(0:rs_descs(grid_level)%rs_desc%group_size - 1))
1255 112 : CALL get_current_loads(loads, rs_descs, grid_level, ntasks, tasks, use_reordered_ranks=.FALSE.)
1256 :
1257 112 : maxdest = SIZE(list, 2)
1258 112 : ncpu = SIZE(list, 3)
1259 112 : natom8 = natoms
1260 :
1261 : ! first find the tasks that deal with the same atom pair
1262 112 : itask_stop = 0
1263 112 : ipair_old = HUGE(ipair_old)
1264 112 : img_old = -1
1265 112 : ALLOCATE (all_dests(0))
1266 112 : ALLOCATE (INDEX(0))
1267 :
1268 : DO
1269 :
1270 : ! first find the range of tasks that deal with the same atom pair
1271 306 : itask_start = itask_stop + 1
1272 306 : itask_stop = itask_start
1273 306 : IF (itask_stop > ntasks) EXIT
1274 194 : ilevel = tasks(itask_stop)%grid_level
1275 194 : img_old = tasks(itask_stop)%image
1276 194 : iatom = tasks(itask_stop)%iatom
1277 194 : jatom = tasks(itask_stop)%jatom
1278 194 : iset = tasks(itask_stop)%iset
1279 194 : jset = tasks(itask_stop)%jset
1280 194 : ipgf = tasks(itask_stop)%ipgf
1281 194 : jpgf = tasks(itask_stop)%jpgf
1282 :
1283 194 : ipair_old = (iatom - 1)*natom8 + (jatom - 1)
1284 : DO
1285 7842 : IF (itask_stop + 1 > ntasks) EXIT
1286 7750 : ilevel = tasks(itask_stop + 1)%grid_level
1287 7750 : img = tasks(itask_stop + 1)%image
1288 7750 : iatom = tasks(itask_stop + 1)%iatom
1289 7750 : jatom = tasks(itask_stop + 1)%jatom
1290 7750 : iset = tasks(itask_stop + 1)%iset
1291 7750 : jset = tasks(itask_stop + 1)%jset
1292 7750 : ipgf = tasks(itask_stop + 1)%ipgf
1293 7750 : jpgf = tasks(itask_stop + 1)%jpgf
1294 :
1295 7750 : ipair = (iatom - 1)*natom8 + (jatom - 1)
1296 7842 : IF (ipair == ipair_old .AND. img == img_old) THEN
1297 : itask_stop = itask_stop + 1
1298 : ELSE
1299 : EXIT
1300 : END IF
1301 : END DO
1302 194 : ipair = ipair_old
1303 194 : nshort = itask_stop - itask_start + 1
1304 :
1305 : ! find the unique list of destinations on this grid level only
1306 194 : DEALLOCATE (all_dests)
1307 582 : ALLOCATE (all_dests(nshort))
1308 194 : DEALLOCATE (index)
1309 388 : ALLOCATE (INDEX(nshort))
1310 8036 : DO i = 1, nshort
1311 7842 : ilevel = tasks(itask_start + i - 1)%grid_level
1312 7842 : img = tasks(itask_start + i - 1)%image
1313 7842 : iatom = tasks(itask_start + i - 1)%iatom
1314 7842 : jatom = tasks(itask_start + i - 1)%jatom
1315 7842 : iset = tasks(itask_start + i - 1)%iset
1316 7842 : jset = tasks(itask_start + i - 1)%jset
1317 7842 : ipgf = tasks(itask_start + i - 1)%ipgf
1318 7842 : jpgf = tasks(itask_start + i - 1)%jpgf
1319 :
1320 8036 : IF (ilevel == grid_level) THEN
1321 2952 : all_dests(i) = decode_rank(tasks(itask_start + i - 1)%destination, SIZE(rs_descs))
1322 : ELSE
1323 4890 : all_dests(i) = HUGE(all_dests(i))
1324 : END IF
1325 : END DO
1326 194 : CALL sort(all_dests, nshort, index)
1327 194 : ndest_pair = 1
1328 7842 : DO i = 2, nshort
1329 7842 : IF ((all_dests(ndest_pair) /= all_dests(i)) .AND. (all_dests(i) /= HUGE(all_dests(i)))) THEN
1330 10 : ndest_pair = ndest_pair + 1
1331 10 : all_dests(ndest_pair) = all_dests(i)
1332 : END IF
1333 : END DO
1334 :
1335 8148 : DO itask = itask_start, itask_stop
1336 :
1337 7842 : dest = decode_rank(tasks(itask)%destination, SIZE(rs_descs)) ! notice that dest can be changed
1338 7842 : ilevel = tasks(itask)%grid_level
1339 7842 : img = tasks(itask)%image
1340 7842 : iatom = tasks(itask)%iatom
1341 7842 : jatom = tasks(itask)%jatom
1342 7842 : iset = tasks(itask)%iset
1343 7842 : jset = tasks(itask)%jset
1344 7842 : ipgf = tasks(itask)%ipgf
1345 7842 : jpgf = tasks(itask)%jpgf
1346 :
1347 : ! Only proceed with tasks which are on this grid level
1348 7842 : IF (ilevel /= grid_level) CYCLE
1349 2952 : ipair = (iatom - 1)*natom8 + (jatom - 1)
1350 2952 : cost = INT(tasks(itask)%cost)
1351 :
1352 194 : SELECT CASE (tasks(itask)%dist_type)
1353 : CASE (1)
1354 2952 : bit_pattern = tasks(itask)%subpatch_pattern
1355 2952 : nopt = 0
1356 2952 : IF (BTEST(bit_pattern, 0)) THEN
1357 0 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [-1, 0, 0])
1358 0 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1359 0 : nopt = nopt + 1
1360 0 : options(nopt) = rank
1361 : END IF
1362 : END IF
1363 2952 : IF (BTEST(bit_pattern, 1)) THEN
1364 0 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [+1, 0, 0])
1365 0 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1366 0 : nopt = nopt + 1
1367 0 : options(nopt) = rank
1368 : END IF
1369 : END IF
1370 2952 : IF (BTEST(bit_pattern, 2)) THEN
1371 48 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, -1, 0])
1372 96 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1373 0 : nopt = nopt + 1
1374 0 : options(nopt) = rank
1375 : END IF
1376 : END IF
1377 2952 : IF (BTEST(bit_pattern, 3)) THEN
1378 0 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, +1, 0])
1379 0 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1380 0 : nopt = nopt + 1
1381 0 : options(nopt) = rank
1382 : END IF
1383 : END IF
1384 2952 : IF (BTEST(bit_pattern, 4)) THEN
1385 1342 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, 0, -1])
1386 2584 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1387 100 : nopt = nopt + 1
1388 100 : options(nopt) = rank
1389 : END IF
1390 : END IF
1391 2952 : IF (BTEST(bit_pattern, 5)) THEN
1392 500 : rank = rs_grid_locate_rank(rs_descs(ilevel)%rs_desc, dest, [0, 0, +1])
1393 1000 : IF (ANY(all_dests(1:ndest_pair) == rank)) THEN
1394 80 : nopt = nopt + 1
1395 80 : options(nopt) = rank
1396 : END IF
1397 : END IF
1398 2952 : IF (nopt > 0) THEN
1399 : ! set it to the rank with the lowest load
1400 180 : rank = options(1)
1401 180 : DO iopt = 2, nopt
1402 180 : IF (loads(rank) > loads(options(iopt))) rank = options(iopt)
1403 : END DO
1404 : ELSE
1405 2772 : rank = dest
1406 : END IF
1407 2952 : li = list_index(list, rank, dest)
1408 2952 : IF (create_list) THEN
1409 1476 : list(2, li, dest) = list(2, li, dest) + cost
1410 : ELSE
1411 1476 : IF (list(1, li, dest) == dest) THEN
1412 1386 : tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1413 : ELSE
1414 90 : IF (list(2, li, dest) >= cost) THEN
1415 10 : list(2, li, dest) = list(2, li, dest) - cost
1416 10 : tasks(itask)%destination = encode_rank(list(1, li, dest), ilevel, SIZE(rs_descs))
1417 : ELSE
1418 80 : tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1419 : END IF
1420 : END IF
1421 : END IF
1422 : CASE (2) ! generalised
1423 0 : li = list_index(list, dest, dest)
1424 0 : IF (create_list) THEN
1425 0 : list(2, li, dest) = list(2, li, dest) + cost
1426 : ELSE
1427 0 : IF (list(1, li, dest) == dest) THEN
1428 0 : tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1429 : ELSE
1430 0 : IF (list(2, li, dest) >= cost) THEN
1431 0 : list(2, li, dest) = list(2, li, dest) - cost
1432 0 : tasks(itask)%destination = encode_rank(list(1, li, dest), ilevel, SIZE(rs_descs))
1433 : ELSE
1434 0 : tasks(itask)%destination = encode_rank(dest, ilevel, SIZE(rs_descs))
1435 : END IF
1436 : END IF
1437 : END IF
1438 : CASE DEFAULT
1439 7842 : CPABORT("Unknown task list distribution type")
1440 : END SELECT
1441 :
1442 : END DO
1443 :
1444 : END DO
1445 :
1446 112 : CALL timestop(handle)
1447 :
1448 224 : END SUBROUTINE compute_load_list
1449 : ! **************************************************************************************************
1450 : !> \brief small helper function to return the proper index in the list array
1451 : !>
1452 : !> \param list ...
1453 : !> \param rank ...
1454 : !> \param dest ...
1455 : !> \return ...
1456 : !> \par History
1457 : !> created 2008-10-06 [Joost VandeVondele]
1458 : ! **************************************************************************************************
1459 2952 : INTEGER FUNCTION list_index(list, rank, dest)
1460 : INTEGER, DIMENSION(:, :, 0:), INTENT(IN) :: list
1461 : INTEGER, INTENT(IN) :: rank, dest
1462 :
1463 2952 : list_index = 1
1464 1424 : DO
1465 4376 : IF (list(1, list_index, dest) == rank) EXIT
1466 1424 : list_index = list_index + 1
1467 : END DO
1468 2952 : END FUNCTION list_index
1469 : ! **************************************************************************************************
1470 : !> \brief create a list with possible destinations (i.e. the central cpu and neighbors) for each cpu
1471 : !> note that we allocate it with an additional field to store the load of this destination
1472 : !>
1473 : !> \param list ...
1474 : !> \param rs_descs ...
1475 : !> \param grid_level ...
1476 : !> \par History
1477 : !> created 2008-10-06 [Joost VandeVondele]
1478 : ! **************************************************************************************************
1479 56 : SUBROUTINE create_destination_list(list, rs_descs, grid_level)
1480 : INTEGER, DIMENSION(:, :, :), POINTER :: list
1481 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1482 : POINTER :: rs_descs
1483 : INTEGER, INTENT(IN) :: grid_level
1484 :
1485 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_destination_list'
1486 :
1487 : INTEGER :: handle, i, icpu, j, maxcount, ncpu, &
1488 : ultimate_max
1489 56 : INTEGER, ALLOCATABLE, DIMENSION(:) :: index, sublist
1490 :
1491 56 : CALL timeset(routineN, handle)
1492 :
1493 56 : CPASSERT(.NOT. ASSOCIATED(list))
1494 56 : ncpu = rs_descs(grid_level)%rs_desc%group_size
1495 56 : ultimate_max = 7
1496 :
1497 168 : ALLOCATE (list(2, ultimate_max, 0:ncpu - 1))
1498 :
1499 56 : ALLOCATE (INDEX(ultimate_max))
1500 56 : ALLOCATE (sublist(ultimate_max))
1501 448 : sublist = HUGE(sublist)
1502 :
1503 56 : maxcount = 1
1504 168 : DO icpu = 0, ncpu - 1
1505 112 : sublist(1) = icpu
1506 112 : sublist(2) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [-1, 0, 0])
1507 112 : sublist(3) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [+1, 0, 0])
1508 112 : sublist(4) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, -1, 0])
1509 112 : sublist(5) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, +1, 0])
1510 112 : sublist(6) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, 0, -1])
1511 112 : sublist(7) = rs_grid_locate_rank(rs_descs(grid_level)%rs_desc, icpu, [0, 0, +1])
1512 : ! only retain unique values of the destination
1513 112 : CALL sort(sublist, ultimate_max, index)
1514 112 : j = 1
1515 784 : DO i = 2, 7
1516 784 : IF (sublist(i) /= sublist(j)) THEN
1517 112 : j = j + 1
1518 112 : sublist(j) = sublist(i)
1519 : END IF
1520 : END DO
1521 112 : maxcount = MAX(maxcount, j)
1522 672 : sublist(j + 1:ultimate_max) = HUGE(sublist)
1523 896 : list(1, :, icpu) = sublist
1524 952 : list(2, :, icpu) = 0
1525 : END DO
1526 :
1527 56 : CALL reallocate(list, 1, 2, 1, maxcount, 0, ncpu - 1)
1528 :
1529 56 : CALL timestop(handle)
1530 :
1531 112 : END SUBROUTINE create_destination_list
1532 :
1533 : ! **************************************************************************************************
1534 : !> \brief given a task list, compute the load of each process everywhere
1535 : !> giving this function the ability to loop over a (sub)set of rs_grids,
1536 : !> and do all the communication in one shot, would speed it up
1537 : !> \param loads ...
1538 : !> \param rs_descs ...
1539 : !> \param grid_level ...
1540 : !> \param ntasks ...
1541 : !> \param tasks ...
1542 : !> \param use_reordered_ranks ...
1543 : !> \par History
1544 : !> none
1545 : !> \author MattW 21/11/2007
1546 : ! **************************************************************************************************
1547 500 : SUBROUTINE get_current_loads(loads, rs_descs, grid_level, ntasks, tasks, use_reordered_ranks)
1548 : INTEGER(KIND=int_8), DIMENSION(:) :: loads
1549 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1550 : POINTER :: rs_descs
1551 : INTEGER :: grid_level, ntasks
1552 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1553 : LOGICAL, INTENT(IN) :: use_reordered_ranks
1554 :
1555 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_current_loads'
1556 :
1557 : INTEGER :: handle, i, iatom, ilevel, img, ipgf, &
1558 : iset, jatom, jpgf, jset
1559 : INTEGER(KIND=int_8) :: total_cost_local
1560 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: recv_buf_i, send_buf_i
1561 : TYPE(realspace_grid_desc_type), POINTER :: desc
1562 :
1563 500 : CALL timeset(routineN, handle)
1564 :
1565 500 : desc => rs_descs(grid_level)%rs_desc
1566 :
1567 : ! allocate local arrays
1568 1500 : ALLOCATE (send_buf_i(desc%group_size))
1569 1000 : ALLOCATE (recv_buf_i(desc%group_size))
1570 :
1571 : ! communication step 1 : compute the total local cost of the tasks
1572 : ! each proc needs to know the amount of work he will receive
1573 :
1574 : ! send buffer now contains for each target the cost of the tasks it will receive
1575 500 : send_buf_i = 0
1576 39092 : DO i = 1, ntasks
1577 38592 : ilevel = tasks(i)%grid_level
1578 38592 : img = tasks(i)%image
1579 38592 : iatom = tasks(i)%iatom
1580 38592 : jatom = tasks(i)%jatom
1581 38592 : iset = tasks(i)%iset
1582 38592 : jset = tasks(i)%jset
1583 38592 : ipgf = tasks(i)%ipgf
1584 38592 : jpgf = tasks(i)%jpgf
1585 38592 : IF (ilevel /= grid_level) CYCLE
1586 11144 : IF (use_reordered_ranks) THEN
1587 : send_buf_i(rs_descs(ilevel)%rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs))) + 1) = &
1588 : send_buf_i(rs_descs(ilevel)%rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs))) + 1) &
1589 3846 : + tasks(i)%cost
1590 : ELSE
1591 : send_buf_i(decode_rank(tasks(i)%destination, SIZE(rs_descs)) + 1) = &
1592 : send_buf_i(decode_rank(tasks(i)%destination, SIZE(rs_descs)) + 1) &
1593 6798 : + tasks(i)%cost
1594 : END IF
1595 : END DO
1596 500 : CALL desc%group%alltoall(send_buf_i, recv_buf_i, 1)
1597 :
1598 : ! communication step 2 : compute the global cost of the tasks
1599 1500 : total_cost_local = SUM(recv_buf_i)
1600 :
1601 : ! after this step, the recv buffer contains the local cost for each CPU
1602 1500 : CALL desc%group%allgather(total_cost_local, loads)
1603 :
1604 500 : CALL timestop(handle)
1605 :
1606 1000 : END SUBROUTINE get_current_loads
1607 : ! **************************************************************************************************
1608 : !> \brief performs load balancing shifting tasks on the replicated grids
1609 : !> this modifies the destination of some of the tasks on replicated
1610 : !> grids, and in this way balances the load
1611 : !> \param rs_descs ...
1612 : !> \param ntasks ...
1613 : !> \param tasks ...
1614 : !> \par History
1615 : !> none
1616 : !> \author MattW 21/11/2007
1617 : ! **************************************************************************************************
1618 50 : SUBROUTINE load_balance_replicated(rs_descs, ntasks, tasks)
1619 :
1620 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1621 : POINTER :: rs_descs
1622 : INTEGER :: ntasks
1623 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1624 :
1625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'load_balance_replicated'
1626 :
1627 : INTEGER :: handle, i, iatom, ilevel, img, ipgf, &
1628 : iset, j, jatom, jpgf, jset, &
1629 : no_overloaded, no_underloaded, &
1630 : proc_receiving
1631 : INTEGER(KIND=int_8) :: average_cost, cost_task_rep, count, &
1632 : offset, total_cost_global
1633 50 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: load_imbalance, loads, recv_buf_i
1634 50 : INTEGER, ALLOCATABLE, DIMENSION(:) :: index
1635 : TYPE(realspace_grid_desc_type), POINTER :: desc
1636 :
1637 50 : CALL timeset(routineN, handle)
1638 :
1639 50 : desc => rs_descs(1)%rs_desc
1640 :
1641 : ! allocate local arrays
1642 150 : ALLOCATE (recv_buf_i(desc%group_size))
1643 100 : ALLOCATE (loads(desc%group_size))
1644 :
1645 50 : recv_buf_i = 0
1646 244 : DO i = 1, SIZE(rs_descs)
1647 194 : CALL get_current_loads(loads, rs_descs, i, ntasks, tasks, use_reordered_ranks=.TRUE.)
1648 632 : recv_buf_i(:) = recv_buf_i + loads
1649 : END DO
1650 :
1651 150 : total_cost_global = SUM(recv_buf_i)
1652 50 : average_cost = total_cost_global/desc%group_size
1653 :
1654 : !
1655 : ! compute how to redistribute the replicated tasks so that the average cost is reached
1656 : !
1657 :
1658 : ! load imbalance measures the load of a given CPU relative
1659 : ! to the optimal load distribution (load=average)
1660 150 : ALLOCATE (load_imbalance(desc%group_size))
1661 150 : ALLOCATE (INDEX(desc%group_size))
1662 :
1663 150 : load_imbalance(:) = recv_buf_i - average_cost
1664 50 : no_overloaded = 0
1665 50 : no_underloaded = 0
1666 :
1667 150 : DO i = 1, desc%group_size
1668 100 : IF (load_imbalance(i) > 0) no_overloaded = no_overloaded + 1
1669 150 : IF (load_imbalance(i) < 0) no_underloaded = no_underloaded + 1
1670 : END DO
1671 :
1672 : ! sort the recv_buffer on number of tasks, gives us index which provides a
1673 : ! mapping between processor ranks and how overloaded the processor
1674 50 : CALL sort(recv_buf_i, SIZE(recv_buf_i), index)
1675 :
1676 : ! find out the number of replicated tasks each proc has
1677 : ! but only those tasks which have not yet been assigned
1678 50 : cost_task_rep = 0
1679 3896 : DO i = 1, ntasks
1680 : IF (tasks(i)%dist_type == 0 &
1681 3896 : .AND. decode_rank(tasks(i)%destination, SIZE(rs_descs)) == decode_rank(tasks(i)%source, SIZE(rs_descs))) THEN
1682 2370 : cost_task_rep = cost_task_rep + tasks(i)%cost
1683 : END IF
1684 : END DO
1685 :
1686 : ! now, correct the load imbalance for the overloaded CPUs
1687 : ! they will send away not more than the total load of replicated tasks
1688 50 : CALL desc%group%allgather(cost_task_rep, recv_buf_i)
1689 :
1690 150 : DO i = 1, desc%group_size
1691 : ! At the moment we can only offload replicated tasks
1692 150 : IF (load_imbalance(i) > 0) load_imbalance(i) = MIN(load_imbalance(i), recv_buf_i(i))
1693 : END DO
1694 :
1695 : ! simplest algorithm I can think of of is that the processor with the most
1696 : ! excess tasks fills up the process needing most, then moves on to next most.
1697 : ! At the moment if we've got less replicated tasks than we're overloaded then
1698 : ! task balancing will be incomplete
1699 :
1700 : ! only need to do anything if I've excess tasks
1701 50 : IF (load_imbalance(desc%my_pos + 1) > 0) THEN
1702 :
1703 23 : count = 0 ! weighted amount of tasks offloaded
1704 23 : offset = 0 ! no of underloaded processes already filled by other more overloaded procs
1705 :
1706 : ! calculate offset
1707 23 : DO i = desc%group_size, desc%group_size - no_overloaded + 1, -1
1708 23 : IF (INDEX(i) == desc%my_pos + 1) THEN
1709 : EXIT
1710 : ELSE
1711 0 : offset = offset + load_imbalance(INDEX(i))
1712 : END IF
1713 : END DO
1714 :
1715 : ! find my starting processor to send to
1716 23 : proc_receiving = HUGE(proc_receiving)
1717 23 : DO i = 1, no_underloaded
1718 23 : offset = offset + load_imbalance(INDEX(i))
1719 23 : IF (offset <= 0) THEN
1720 : proc_receiving = i
1721 : EXIT
1722 : END IF
1723 : END DO
1724 :
1725 : ! offset now contains minus the number of tasks proc_receiving requires
1726 : ! we fill this up by adjusting the destination of tasks on the replicated grid,
1727 : ! then move to next most underloaded proc
1728 1408 : DO j = 1, ntasks
1729 : IF (tasks(j)%dist_type == 0 &
1730 1408 : .AND. decode_rank(tasks(j)%destination, SIZE(rs_descs)) == decode_rank(tasks(j)%source, SIZE(rs_descs))) THEN
1731 : ! just avoid sending to non existing procs due to integer truncation
1732 : ! in the computation of the average
1733 819 : IF (proc_receiving > no_underloaded) EXIT
1734 : ! set new destination
1735 819 : ilevel = tasks(j)%grid_level
1736 819 : img = tasks(j)%image
1737 819 : iatom = tasks(j)%iatom
1738 819 : jatom = tasks(j)%jatom
1739 819 : iset = tasks(j)%iset
1740 819 : jset = tasks(j)%jset
1741 819 : ipgf = tasks(j)%ipgf
1742 819 : jpgf = tasks(j)%jpgf
1743 819 : tasks(j)%destination = encode_rank(INDEX(proc_receiving) - 1, ilevel, SIZE(rs_descs))
1744 819 : offset = offset + tasks(j)%cost
1745 819 : count = count + tasks(j)%cost
1746 819 : IF (count >= load_imbalance(desc%my_pos + 1)) EXIT
1747 796 : IF (offset > 0) THEN
1748 0 : proc_receiving = proc_receiving + 1
1749 : ! just avoid sending to non existing procs due to integer truncation
1750 : ! in the computation of the average
1751 0 : IF (proc_receiving > no_underloaded) EXIT
1752 0 : offset = load_imbalance(INDEX(proc_receiving))
1753 : END IF
1754 : END IF
1755 : END DO
1756 : END IF
1757 :
1758 50 : DEALLOCATE (index)
1759 50 : DEALLOCATE (load_imbalance)
1760 :
1761 50 : CALL timestop(handle)
1762 :
1763 100 : END SUBROUTINE load_balance_replicated
1764 :
1765 : ! **************************************************************************************************
1766 : !> \brief given an input task list, redistribute so that all tasks can be processed locally,
1767 : !> i.e. dest equals rank
1768 : !> \param rs_descs ...
1769 : !> \param ntasks ...
1770 : !> \param tasks ...
1771 : !> \param ntasks_recv ...
1772 : !> \param tasks_recv ...
1773 : !> \par History
1774 : !> none
1775 : !> \author MattW 21/11/2007
1776 : ! **************************************************************************************************
1777 50 : SUBROUTINE create_local_tasks(rs_descs, ntasks, tasks, ntasks_recv, tasks_recv)
1778 :
1779 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1780 : POINTER :: rs_descs
1781 : INTEGER :: ntasks
1782 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1783 : INTEGER :: ntasks_recv
1784 : TYPE(task_type), DIMENSION(:), POINTER :: tasks_recv
1785 :
1786 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_local_tasks'
1787 :
1788 : INTEGER :: handle, i, j, k, l, rank
1789 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: recv_buf, send_buf
1790 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_disps, recv_sizes, send_disps, &
1791 : send_sizes
1792 : TYPE(realspace_grid_desc_type), POINTER :: desc
1793 :
1794 50 : CALL timeset(routineN, handle)
1795 :
1796 50 : desc => rs_descs(1)%rs_desc
1797 :
1798 : ! allocate local arrays
1799 150 : ALLOCATE (send_sizes(desc%group_size))
1800 100 : ALLOCATE (recv_sizes(desc%group_size))
1801 100 : ALLOCATE (send_disps(desc%group_size))
1802 100 : ALLOCATE (recv_disps(desc%group_size))
1803 150 : ALLOCATE (send_buf(desc%group_size))
1804 100 : ALLOCATE (recv_buf(desc%group_size))
1805 :
1806 : ! fill send buffer, now counting how many tasks will be send (stored in an int8 array for convenience only).
1807 50 : send_buf = 0
1808 3896 : DO i = 1, ntasks
1809 : rank = rs_descs(decode_level(tasks(i)%destination, SIZE(rs_descs))) &
1810 3846 : %rs_desc%virtual2real(decode_rank(tasks(i)%destination, SIZE(rs_descs)))
1811 3896 : send_buf(rank + 1) = send_buf(rank + 1) + 1
1812 : END DO
1813 :
1814 50 : CALL desc%group%alltoall(send_buf, recv_buf, 1)
1815 :
1816 : ! pack the tasks, and send them around
1817 :
1818 50 : send_sizes = 0
1819 50 : send_disps = 0
1820 50 : recv_sizes = 0
1821 50 : recv_disps = 0
1822 :
1823 50 : send_sizes(1) = INT(send_buf(1)*task_size_in_int8)
1824 50 : recv_sizes(1) = INT(recv_buf(1)*task_size_in_int8)
1825 100 : DO i = 2, desc%group_size
1826 50 : send_sizes(i) = INT(send_buf(i)*task_size_in_int8)
1827 50 : recv_sizes(i) = INT(recv_buf(i)*task_size_in_int8)
1828 50 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1829 100 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1830 : END DO
1831 :
1832 : ! deallocate old send/recv buffers
1833 50 : DEALLOCATE (send_buf)
1834 50 : DEALLOCATE (recv_buf)
1835 :
1836 : ! allocate them with new sizes
1837 243 : ALLOCATE (send_buf(SUM(send_sizes)))
1838 248 : ALLOCATE (recv_buf(SUM(recv_sizes)))
1839 :
1840 : ! do packing
1841 50 : send_buf = 0
1842 50 : send_sizes = 0
1843 3896 : DO j = 1, ntasks
1844 : i = rs_descs(decode_level(tasks(j)%destination, SIZE(rs_descs))) &
1845 3846 : %rs_desc%virtual2real(decode_rank(tasks(j)%destination, SIZE(rs_descs))) + 1
1846 3846 : l = send_disps(i) + send_sizes(i)
1847 3846 : CALL serialize_task(tasks(j), send_buf(l + 1:l + task_size_in_int8))
1848 3896 : send_sizes(i) = send_sizes(i) + task_size_in_int8
1849 : END DO
1850 :
1851 : ! do communication
1852 50 : CALL desc%group%alltoall(send_buf, send_sizes, send_disps, recv_buf, recv_sizes, recv_disps)
1853 :
1854 50 : DEALLOCATE (send_buf)
1855 :
1856 150 : ntasks_recv = SUM(recv_sizes)/task_size_in_int8
1857 4144 : ALLOCATE (tasks_recv(ntasks_recv))
1858 :
1859 : ! do unpacking
1860 50 : l = 0
1861 150 : DO i = 1, desc%group_size
1862 3996 : DO j = 0, recv_sizes(i)/task_size_in_int8 - 1
1863 3846 : l = l + 1
1864 3846 : k = recv_disps(i) + j*task_size_in_int8
1865 3946 : CALL deserialize_task(tasks_recv(l), recv_buf(k + 1:k + task_size_in_int8))
1866 : END DO
1867 : END DO
1868 :
1869 50 : DEALLOCATE (recv_buf)
1870 50 : DEALLOCATE (send_sizes)
1871 50 : DEALLOCATE (recv_sizes)
1872 50 : DEALLOCATE (send_disps)
1873 50 : DEALLOCATE (recv_disps)
1874 :
1875 50 : CALL timestop(handle)
1876 :
1877 50 : END SUBROUTINE create_local_tasks
1878 :
1879 : ! **************************************************************************************************
1880 : !> \brief Assembles tasks to be performed on local grid
1881 : !> \param rs_descs the grids
1882 : !> \param ntasks Number of tasks for local processing
1883 : !> \param natoms ...
1884 : !> \param nimages ...
1885 : !> \param tasks the task set generated on this processor
1886 : !> \param rval ...
1887 : !> \param atom_pair_send ...
1888 : !> \param atom_pair_recv ...
1889 : !> \param symmetric ...
1890 : !> \param reorder_rs_grid_ranks ...
1891 : !> \param skip_load_balance_distributed ...
1892 : !> \par History
1893 : !> none
1894 : !> \author MattW 21/11/2007
1895 : ! **************************************************************************************************
1896 19941 : SUBROUTINE distribute_tasks(rs_descs, ntasks, natoms, &
1897 : tasks, atom_pair_send, atom_pair_recv, &
1898 : symmetric, reorder_rs_grid_ranks, skip_load_balance_distributed)
1899 :
1900 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1901 : POINTER :: rs_descs
1902 : INTEGER :: ntasks, natoms
1903 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1904 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_send, atom_pair_recv
1905 : LOGICAL, INTENT(IN) :: symmetric, reorder_rs_grid_ranks, &
1906 : skip_load_balance_distributed
1907 :
1908 : CHARACTER(LEN=*), PARAMETER :: routineN = 'distribute_tasks'
1909 :
1910 : INTEGER :: handle, igrid_level, irank, ntasks_recv
1911 : INTEGER(KIND=int_8) :: load_gap, max_load, replicated_load
1912 19941 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: total_loads, total_loads_tmp, trial_loads
1913 19941 : INTEGER(KIND=int_8), DIMENSION(:, :), POINTER :: loads
1914 19941 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indices, real2virtual, total_index
1915 : LOGICAL :: distributed_grids, fixed_first_grid
1916 : TYPE(realspace_grid_desc_type), POINTER :: desc
1917 19941 : TYPE(task_type), DIMENSION(:), POINTER :: tasks_recv
1918 :
1919 19941 : CALL timeset(routineN, handle)
1920 :
1921 19941 : CPASSERT(ASSOCIATED(tasks))
1922 :
1923 : ! *** figure out if we have distributed grids
1924 19941 : distributed_grids = .FALSE.
1925 98940 : DO igrid_level = 1, SIZE(rs_descs)
1926 98940 : IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
1927 56 : distributed_grids = .TRUE.
1928 : END IF
1929 : END DO
1930 19941 : desc => rs_descs(1)%rs_desc
1931 :
1932 19941 : IF (distributed_grids) THEN
1933 :
1934 200 : ALLOCATE (loads(0:desc%group_size - 1, SIZE(rs_descs)))
1935 150 : ALLOCATE (total_loads(0:desc%group_size - 1))
1936 :
1937 50 : total_loads = 0
1938 :
1939 : ! First round of balancing on the distributed grids
1940 : ! we just balance each of the distributed grids independently
1941 244 : DO igrid_level = 1, SIZE(rs_descs)
1942 244 : IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
1943 :
1944 56 : IF (.NOT. skip_load_balance_distributed) THEN
1945 56 : CALL load_balance_distributed(tasks, ntasks, rs_descs, igrid_level, natoms)
1946 : END IF
1947 :
1948 : CALL get_current_loads(loads(:, igrid_level), rs_descs, igrid_level, ntasks, &
1949 56 : tasks, use_reordered_ranks=.FALSE.)
1950 :
1951 168 : total_loads(:) = total_loads + loads(:, igrid_level)
1952 :
1953 : END IF
1954 : END DO
1955 :
1956 : ! calculate the total load of replicated tasks, so we can decide if it is worth
1957 : ! reordering the distributed grid levels
1958 :
1959 50 : replicated_load = 0
1960 244 : DO igrid_level = 1, SIZE(rs_descs)
1961 244 : IF (.NOT. rs_descs(igrid_level)%rs_desc%distributed) THEN
1962 : CALL get_current_loads(loads(:, igrid_level), rs_descs, igrid_level, ntasks, &
1963 138 : tasks, use_reordered_ranks=.FALSE.)
1964 414 : replicated_load = replicated_load + SUM(loads(:, igrid_level))
1965 : END IF
1966 : END DO
1967 :
1968 : ! Now we adjust the rank ordering based on the current loads
1969 : ! we leave the first distributed level and all the replicated levels in the default order
1970 50 : IF (reorder_rs_grid_ranks) THEN
1971 50 : fixed_first_grid = .FALSE.
1972 244 : DO igrid_level = 1, SIZE(rs_descs)
1973 244 : IF (rs_descs(igrid_level)%rs_desc%distributed) THEN
1974 56 : IF (fixed_first_grid .EQV. .FALSE.) THEN
1975 150 : total_loads(:) = loads(:, igrid_level)
1976 : fixed_first_grid = .TRUE.
1977 : ELSE
1978 18 : ALLOCATE (trial_loads(0:desc%group_size - 1))
1979 :
1980 18 : trial_loads(:) = total_loads + loads(:, igrid_level)
1981 18 : max_load = MAXVAL(trial_loads)
1982 : load_gap = 0
1983 18 : DO irank = 0, desc%group_size - 1
1984 18 : load_gap = load_gap + max_load - trial_loads(irank)
1985 : END DO
1986 :
1987 : ! If there is not enough replicated load to load balance well enough
1988 : ! then we will reorder this grid level
1989 6 : IF (load_gap > replicated_load*1.05_dp) THEN
1990 :
1991 18 : ALLOCATE (indices(0:desc%group_size - 1))
1992 12 : ALLOCATE (total_index(0:desc%group_size - 1))
1993 12 : ALLOCATE (total_loads_tmp(0:desc%group_size - 1))
1994 12 : ALLOCATE (real2virtual(0:desc%group_size - 1))
1995 :
1996 18 : total_loads_tmp(:) = total_loads
1997 6 : CALL sort(total_loads_tmp, desc%group_size, total_index)
1998 6 : CALL sort(loads(:, igrid_level), desc%group_size, indices)
1999 :
2000 : ! Reorder so that the rank with smallest load on this grid level is paired with
2001 : ! the highest load in total
2002 18 : DO irank = 0, desc%group_size - 1
2003 : total_loads(total_index(irank) - 1) = total_loads(total_index(irank) - 1) + &
2004 12 : loads(desc%group_size - irank - 1, igrid_level)
2005 18 : real2virtual(total_index(irank) - 1) = indices(desc%group_size - irank - 1) - 1
2006 : END DO
2007 :
2008 6 : CALL rs_grid_reorder_ranks(rs_descs(igrid_level)%rs_desc, real2virtual)
2009 :
2010 6 : DEALLOCATE (indices)
2011 6 : DEALLOCATE (total_index)
2012 6 : DEALLOCATE (total_loads_tmp)
2013 6 : DEALLOCATE (real2virtual)
2014 : ELSE
2015 0 : total_loads(:) = trial_loads
2016 : END IF
2017 :
2018 6 : DEALLOCATE (trial_loads)
2019 :
2020 : END IF
2021 : END IF
2022 : END DO
2023 : END IF
2024 :
2025 : ! Now we use the replicated tasks to balance out the rest of the load
2026 50 : CALL load_balance_replicated(rs_descs, ntasks, tasks)
2027 :
2028 : ! given a list of tasks, this will do the needed reshuffle so that all tasks will be local
2029 50 : CALL create_local_tasks(rs_descs, ntasks, tasks, ntasks_recv, tasks_recv)
2030 :
2031 : !
2032 : ! tasks list are complete, we can compute the list of atomic blocks (atom pairs)
2033 : ! we will be sending. These lists are needed for redistribute_matrix.
2034 : !
2035 50 : CALL get_atom_pair(atom_pair_send, tasks, ntasks=ntasks, send=.TRUE., symmetric=symmetric, rs_descs=rs_descs)
2036 :
2037 50 : CALL get_atom_pair(atom_pair_recv, tasks_recv, ntasks=ntasks_recv, send=.FALSE., symmetric=symmetric, rs_descs=rs_descs)
2038 :
2039 : ! cleanup, at this point we don't need the original tasks anymore
2040 50 : DEALLOCATE (tasks)
2041 50 : DEALLOCATE (loads)
2042 50 : DEALLOCATE (total_loads)
2043 :
2044 : ELSE
2045 19891 : tasks_recv => tasks
2046 19891 : ntasks_recv = ntasks
2047 19891 : CALL get_atom_pair(atom_pair_recv, tasks_recv, ntasks=ntasks_recv, send=.FALSE., symmetric=symmetric, rs_descs=rs_descs)
2048 : ! not distributed, hence atom_pair_send not needed
2049 : END IF
2050 :
2051 : ! here we sort the task list we will process locally.
2052 59433 : ALLOCATE (indices(ntasks_recv))
2053 19941 : CALL tasks_sort(tasks_recv, ntasks_recv, indices)
2054 19941 : DEALLOCATE (indices)
2055 :
2056 : !
2057 : ! final lists are ready
2058 : !
2059 :
2060 19941 : tasks => tasks_recv
2061 19941 : ntasks = ntasks_recv
2062 :
2063 19941 : CALL timestop(handle)
2064 :
2065 39882 : END SUBROUTINE distribute_tasks
2066 :
2067 : ! **************************************************************************************************
2068 : !> \brief ...
2069 : !> \param atom_pair ...
2070 : !> \param my_tasks ...
2071 : !> \param send ...
2072 : !> \param symmetric ...
2073 : !> \param natoms ...
2074 : !> \param nimages ...
2075 : !> \param rs_descs ...
2076 : ! **************************************************************************************************
2077 19991 : SUBROUTINE get_atom_pair(atom_pair, tasks, ntasks, send, symmetric, rs_descs)
2078 :
2079 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair
2080 : TYPE(task_type), DIMENSION(:), INTENT(INOUT) :: tasks
2081 : INTEGER, INTENT(IN) :: ntasks
2082 : LOGICAL, INTENT(IN) :: send, symmetric
2083 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), INTENT(IN) :: rs_descs
2084 :
2085 : INTEGER :: i, ilevel, iatom, jatom, npairs, virt_rank
2086 19991 : INTEGER, DIMENSION(:), ALLOCATABLE :: indices
2087 19991 : TYPE(atom_pair_type), DIMENSION(:), ALLOCATABLE :: atom_pair_tmp
2088 :
2089 19991 : CPASSERT(.NOT. ASSOCIATED(atom_pair))
2090 19991 : IF (ntasks == 0) THEN
2091 397 : ALLOCATE (atom_pair(0))
2092 : RETURN
2093 : END IF
2094 :
2095 : ! calculate list of atom pairs
2096 : ! fill pair list taking into account symmetry
2097 10886346 : ALLOCATE (atom_pair_tmp(ntasks))
2098 10847158 : DO i = 1, ntasks
2099 10827564 : atom_pair_tmp(i)%image = tasks(i)%image
2100 10827564 : iatom = tasks(i)%iatom
2101 10827564 : jatom = tasks(i)%jatom
2102 10827564 : IF (symmetric .AND. iatom > jatom) THEN
2103 : ! iatom / jatom swapped
2104 4144106 : atom_pair_tmp(i)%row = jatom
2105 4144106 : atom_pair_tmp(i)%col = iatom
2106 : ELSE
2107 6683458 : atom_pair_tmp(i)%row = iatom
2108 6683458 : atom_pair_tmp(i)%col = jatom
2109 : END IF
2110 :
2111 10847158 : IF (send) THEN
2112 : ! If sending, we need to use the 'real rank' as the pair has to be sent to the process which
2113 : ! actually has the correct part of the rs_grid to do the mapping
2114 3846 : ilevel = tasks(i)%grid_level
2115 3846 : virt_rank = decode_rank(tasks(i)%destination, SIZE(rs_descs))
2116 3846 : atom_pair_tmp(i)%rank = rs_descs(ilevel)%rs_desc%virtual2real(virt_rank)
2117 : ELSE
2118 : ! If we are receiving, then no conversion is needed as the rank is that of the process with the
2119 : ! required matrix block, and the ordering of the rs grid is irrelevant
2120 10823718 : atom_pair_tmp(i)%rank = decode_rank(tasks(i)%source, SIZE(rs_descs))
2121 : END IF
2122 : END DO
2123 :
2124 : ! find unique atom pairs that I'm sending/receiving
2125 58782 : ALLOCATE (indices(ntasks))
2126 19594 : CALL atom_pair_sort(atom_pair_tmp, ntasks, indices)
2127 19594 : npairs = 1
2128 19594 : tasks(indices(1))%pair_index = 1
2129 10827564 : DO i = 2, ntasks
2130 10807970 : IF (atom_pair_less_than(atom_pair_tmp(i - 1), atom_pair_tmp(i))) THEN
2131 698354 : npairs = npairs + 1
2132 698354 : atom_pair_tmp(npairs) = atom_pair_tmp(i)
2133 : END IF
2134 10827564 : tasks(indices(i))%pair_index = npairs
2135 : END DO
2136 19594 : DEALLOCATE (indices)
2137 :
2138 : ! Copy unique pairs to final location.
2139 776730 : ALLOCATE (atom_pair(npairs))
2140 737542 : atom_pair(:) = atom_pair_tmp(:npairs)
2141 19594 : DEALLOCATE (atom_pair_tmp)
2142 :
2143 : END SUBROUTINE get_atom_pair
2144 :
2145 : ! **************************************************************************************************
2146 : !> \brief redistributes the matrix so that it can be used in realspace operations
2147 : !> i.e. according to the task lists for collocate and integrate.
2148 : !> This routine can become a bottleneck in large calculations.
2149 : !> \param rs_descs ...
2150 : !> \param pmats ...
2151 : !> \param atom_pair_send ...
2152 : !> \param atom_pair_recv ...
2153 : !> \param natoms ...
2154 : !> \param nimages ...
2155 : !> \param scatter ...
2156 : !> \param hmats ...
2157 : ! **************************************************************************************************
2158 0 : SUBROUTINE rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, &
2159 : nimages, scatter, hmats)
2160 :
2161 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
2162 : POINTER :: rs_descs
2163 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmats
2164 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_send, atom_pair_recv
2165 : INTEGER :: nimages
2166 : LOGICAL :: scatter
2167 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
2168 : POINTER :: hmats
2169 :
2170 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rs_distribute_matrix'
2171 :
2172 : INTEGER :: acol, arow, handle, i, img, j, k, l, me, &
2173 : nblkcols_total, nblkrows_total, ncol, &
2174 : nrow, nthread, nthread_left
2175 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row, recv_disps, &
2176 0 : recv_pair_count, recv_pair_disps, recv_sizes, send_disps, send_pair_count, &
2177 0 : send_pair_disps, send_sizes
2178 0 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
2179 : LOGICAL :: found
2180 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: recv_buf_r, send_buf_r
2181 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block
2182 : TYPE(dbcsr_type), POINTER :: hmat, pmat
2183 : TYPE(realspace_grid_desc_type), POINTER :: desc
2184 0 : REAL(kind=dp), DIMENSION(:), POINTER :: vector
2185 :
2186 0 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
2187 :
2188 0 : CALL timeset(routineN, handle)
2189 :
2190 0 : IF (.NOT. scatter) THEN
2191 0 : CPASSERT(PRESENT(hmats))
2192 : END IF
2193 :
2194 0 : desc => rs_descs(1)%rs_desc
2195 0 : me = desc%my_pos + 1
2196 :
2197 : ! allocate local arrays
2198 0 : ALLOCATE (send_sizes(desc%group_size))
2199 0 : ALLOCATE (recv_sizes(desc%group_size))
2200 0 : ALLOCATE (send_disps(desc%group_size))
2201 0 : ALLOCATE (recv_disps(desc%group_size))
2202 0 : ALLOCATE (send_pair_count(desc%group_size))
2203 0 : ALLOCATE (recv_pair_count(desc%group_size))
2204 0 : ALLOCATE (send_pair_disps(desc%group_size))
2205 0 : ALLOCATE (recv_pair_disps(desc%group_size))
2206 :
2207 0 : pmat => pmats(1)%matrix
2208 : CALL dbcsr_get_info(pmat, &
2209 : row_blk_size=row_blk_size, &
2210 : col_blk_size=col_blk_size, &
2211 : nblkrows_total=nblkrows_total, &
2212 0 : nblkcols_total=nblkcols_total)
2213 0 : ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total), &
2214 0 : first_col(nblkcols_total), last_col(nblkcols_total))
2215 0 : CALL dbcsr_convert_sizes_to_offsets(row_blk_size, first_row, last_row)
2216 0 : CALL dbcsr_convert_sizes_to_offsets(col_blk_size, first_col, last_col)
2217 :
2218 : ! set up send buffer sizes
2219 0 : send_sizes = 0
2220 0 : send_pair_count = 0
2221 0 : DO i = 1, SIZE(atom_pair_send)
2222 0 : k = atom_pair_send(i)%rank + 1 ! proc we're sending this block to
2223 0 : arow = atom_pair_send(i)%row
2224 0 : acol = atom_pair_send(i)%col
2225 0 : nrow = last_row(arow) - first_row(arow) + 1
2226 0 : ncol = last_col(acol) - first_col(acol) + 1
2227 0 : send_sizes(k) = send_sizes(k) + nrow*ncol
2228 0 : send_pair_count(k) = send_pair_count(k) + 1
2229 : END DO
2230 :
2231 0 : send_disps = 0
2232 0 : send_pair_disps = 0
2233 0 : DO i = 2, desc%group_size
2234 0 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
2235 0 : send_pair_disps(i) = send_pair_disps(i - 1) + send_pair_count(i - 1)
2236 : END DO
2237 :
2238 0 : ALLOCATE (send_buf_r(SUM(send_sizes)))
2239 :
2240 : ! set up recv buffer
2241 :
2242 0 : recv_sizes = 0
2243 0 : recv_pair_count = 0
2244 0 : DO i = 1, SIZE(atom_pair_recv)
2245 0 : k = atom_pair_recv(i)%rank + 1 ! proc we're receiving this data from
2246 0 : arow = atom_pair_recv(i)%row
2247 0 : acol = atom_pair_recv(i)%col
2248 0 : nrow = last_row(arow) - first_row(arow) + 1
2249 0 : ncol = last_col(acol) - first_col(acol) + 1
2250 0 : recv_sizes(k) = recv_sizes(k) + nrow*ncol
2251 0 : recv_pair_count(k) = recv_pair_count(k) + 1
2252 : END DO
2253 :
2254 0 : recv_disps = 0
2255 0 : recv_pair_disps = 0
2256 0 : DO i = 2, desc%group_size
2257 0 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
2258 0 : recv_pair_disps(i) = recv_pair_disps(i - 1) + recv_pair_count(i - 1)
2259 : END DO
2260 0 : ALLOCATE (recv_buf_r(SUM(recv_sizes)))
2261 :
2262 : !$OMP PARALLEL DEFAULT(NONE) &
2263 : !$OMP SHARED(desc,send_pair_count,send_pair_disps,nimages),&
2264 : !$OMP SHARED(last_row,first_row,last_col,first_col),&
2265 : !$OMP SHARED(pmats,send_buf_r,send_disps,send_sizes),&
2266 : !$OMP SHARED(atom_pair_send,me,hmats,nblkrows_total),&
2267 : !$OMP SHARED(atom_pair_recv,recv_buf_r,scatter,recv_pair_disps), &
2268 : !$OMP SHARED(recv_sizes,recv_disps,recv_pair_count,locks), &
2269 : !$OMP PRIVATE(i,img,arow,acol,nrow,ncol,p_block,found,j,k,l),&
2270 0 : !$OMP PRIVATE(nthread,h_block,nthread_left,hmat,pmat,vector)
2271 :
2272 : nthread = 1
2273 : !$ nthread = omp_get_num_threads()
2274 : nthread_left = 1
2275 : !$ nthread_left = MAX(1, nthread - 1)
2276 :
2277 : ! do packing
2278 : !$OMP DO schedule(guided)
2279 : DO l = 1, desc%group_size
2280 : IF (l == me) CYCLE
2281 : send_sizes(l) = 0
2282 : DO i = 1, send_pair_count(l)
2283 : arow = atom_pair_send(send_pair_disps(l) + i)%row
2284 : acol = atom_pair_send(send_pair_disps(l) + i)%col
2285 : img = atom_pair_send(send_pair_disps(l) + i)%image
2286 : nrow = last_row(arow) - first_row(arow) + 1
2287 : ncol = last_col(acol) - first_col(acol) + 1
2288 : pmat => pmats(img)%matrix
2289 : CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, block=p_block, found=found)
2290 : CPASSERT(found)
2291 :
2292 : DO k = 1, ncol
2293 : DO j = 1, nrow
2294 : send_buf_r(send_disps(l) + send_sizes(l) + j + (k - 1)*nrow) = p_block(j, k)
2295 : END DO
2296 : END DO
2297 : send_sizes(l) = send_sizes(l) + nrow*ncol
2298 : END DO
2299 : END DO
2300 : !$OMP END DO
2301 :
2302 : IF (.NOT. scatter) THEN
2303 : ! We need locks to protect concurrent summation into H
2304 : !$OMP SINGLE
2305 : !$ ALLOCATE (locks(nthread*10))
2306 : !$OMP END SINGLE
2307 :
2308 : !$OMP DO
2309 : !$ do i = 1, nthread*10
2310 : !$ call omp_init_lock(locks(i))
2311 : !$ end do
2312 : !$OMP END DO
2313 : END IF
2314 :
2315 : !$OMP MASTER
2316 : ! do communication
2317 : CALL desc%group%alltoall(send_buf_r, send_sizes, send_disps, &
2318 : recv_buf_r, recv_sizes, recv_disps)
2319 : !$OMP END MASTER
2320 :
2321 : ! If this is a scatter, then no need to copy local blocks,
2322 : ! If not, we sum them directly into H (bypassing the alltoall)
2323 : IF (.NOT. scatter) THEN
2324 :
2325 : ! Distribute work over remaining threads assuming one is still in the alltoall
2326 : !$OMP DO schedule(dynamic,MAX(1,send_pair_count(me)/nthread_left))
2327 : DO i = 1, send_pair_count(me)
2328 : arow = atom_pair_send(send_pair_disps(me) + i)%row
2329 : acol = atom_pair_send(send_pair_disps(me) + i)%col
2330 : img = atom_pair_send(send_pair_disps(me) + i)%image
2331 : nrow = last_row(arow) - first_row(arow) + 1
2332 : ncol = last_col(acol) - first_col(acol) + 1
2333 : hmat => hmats(img)%matrix
2334 : pmat => pmats(img)%matrix
2335 : CALL dbcsr_get_block_p(matrix=hmat, row=arow, col=acol, BLOCK=h_block, found=found)
2336 : CPASSERT(found)
2337 : CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, BLOCK=p_block, found=found)
2338 : CPASSERT(found)
2339 :
2340 : !$ call omp_set_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2341 : DO k = 1, ncol
2342 : DO j = 1, nrow
2343 : h_block(j, k) = h_block(j, k) + p_block(j, k)
2344 : END DO
2345 : END DO
2346 : !$ call omp_unset_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2347 : END DO
2348 : !$OMP END DO
2349 : ELSE
2350 : ! We will insert new blocks into P, so create mutable work matrices
2351 : DO img = 1, nimages
2352 : pmat => pmats(img)%matrix
2353 : CALL dbcsr_work_create(pmat, work_mutable=.TRUE., &
2354 : nblks_guess=SIZE(atom_pair_recv)/nthread, sizedata_guess=SIZE(recv_buf_r)/nthread, &
2355 : n=nthread)
2356 : END DO
2357 : END IF
2358 :
2359 : ! wait for comm and setup to finish
2360 : !$OMP BARRIER
2361 :
2362 : !do unpacking
2363 : !$OMP DO schedule(guided)
2364 : DO l = 1, desc%group_size
2365 : IF (l == me) CYCLE
2366 : recv_sizes(l) = 0
2367 : DO i = 1, recv_pair_count(l)
2368 : arow = atom_pair_recv(recv_pair_disps(l) + i)%row
2369 : acol = atom_pair_recv(recv_pair_disps(l) + i)%col
2370 : img = atom_pair_recv(recv_pair_disps(l) + i)%image
2371 : nrow = last_row(arow) - first_row(arow) + 1
2372 : ncol = last_col(acol) - first_col(acol) + 1
2373 : pmat => pmats(img)%matrix
2374 : NULLIFY (p_block)
2375 : CALL dbcsr_get_block_p(matrix=pmat, row=arow, col=acol, BLOCK=p_block, found=found)
2376 :
2377 : IF (PRESENT(hmats)) THEN
2378 : hmat => hmats(img)%matrix
2379 : CALL dbcsr_get_block_p(matrix=hmat, row=arow, col=acol, BLOCK=h_block, found=found)
2380 : CPASSERT(found)
2381 : END IF
2382 :
2383 : IF (scatter .AND. .NOT. ASSOCIATED(p_block)) THEN
2384 : vector => recv_buf_r(recv_disps(l) + recv_sizes(l) + 1:recv_disps(l) + recv_sizes(l) + nrow*ncol)
2385 : CALL dbcsr_put_block(pmat, arow, acol, block=RESHAPE(vector, [nrow, ncol]))
2386 : END IF
2387 : IF (.NOT. scatter) THEN
2388 : !$ call omp_set_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2389 : DO k = 1, ncol
2390 : DO j = 1, nrow
2391 : h_block(j, k) = h_block(j, k) + recv_buf_r(recv_disps(l) + recv_sizes(l) + j + (k - 1)*nrow)
2392 : END DO
2393 : END DO
2394 : !$ call omp_unset_lock(locks((arow - 1)*nthread*10/nblkrows_total + 1))
2395 : END IF
2396 : recv_sizes(l) = recv_sizes(l) + nrow*ncol
2397 : END DO
2398 : END DO
2399 : !$OMP END DO
2400 :
2401 : !$ IF (.not. scatter) THEN
2402 : !$OMP DO
2403 : !$ do i = 1, nthread*10
2404 : !$ call omp_destroy_lock(locks(i))
2405 : !$ end do
2406 : !$OMP END DO
2407 : !$ END IF
2408 :
2409 : !$OMP SINGLE
2410 : !$ IF (.not. scatter) THEN
2411 : !$ DEALLOCATE (locks)
2412 : !$ END IF
2413 : !$OMP END SINGLE NOWAIT
2414 :
2415 : IF (scatter) THEN
2416 : ! Blocks were added to P
2417 : DO img = 1, nimages
2418 : pmat => pmats(img)%matrix
2419 : CALL dbcsr_finalize(pmat)
2420 : END DO
2421 : END IF
2422 : !$OMP END PARALLEL
2423 :
2424 0 : DEALLOCATE (send_buf_r)
2425 0 : DEALLOCATE (recv_buf_r)
2426 :
2427 0 : DEALLOCATE (send_sizes)
2428 0 : DEALLOCATE (recv_sizes)
2429 0 : DEALLOCATE (send_disps)
2430 0 : DEALLOCATE (recv_disps)
2431 0 : DEALLOCATE (send_pair_count)
2432 0 : DEALLOCATE (recv_pair_count)
2433 0 : DEALLOCATE (send_pair_disps)
2434 0 : DEALLOCATE (recv_pair_disps)
2435 :
2436 0 : DEALLOCATE (first_row, last_row, first_col, last_col)
2437 :
2438 0 : CALL timestop(handle)
2439 :
2440 0 : END SUBROUTINE rs_distribute_matrix
2441 :
2442 : ! **************************************************************************************************
2443 : !> \brief Calculates offsets and sizes for rs_scatter_matrix and rs_copy_matrix.
2444 : !> \author Ole Schuett
2445 : ! **************************************************************************************************
2446 17705 : SUBROUTINE rs_calc_offsets(pairs, nsgf, group_size, &
2447 : pair_offsets, rank_offsets, rank_sizes, buffer_size)
2448 : TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: pairs
2449 : INTEGER, DIMENSION(:), INTENT(IN) :: nsgf
2450 : INTEGER, INTENT(IN) :: group_size
2451 : INTEGER, DIMENSION(:), POINTER :: pair_offsets, rank_offsets, rank_sizes
2452 : INTEGER, INTENT(INOUT) :: buffer_size
2453 :
2454 : INTEGER :: acol, arow, i, block_size, total_size, k, prev_k
2455 :
2456 17705 : IF (ASSOCIATED(pair_offsets)) DEALLOCATE (pair_offsets)
2457 17705 : IF (ASSOCIATED(rank_offsets)) DEALLOCATE (rank_offsets)
2458 17705 : IF (ASSOCIATED(rank_sizes)) DEALLOCATE (rank_sizes)
2459 :
2460 : ! calculate buffer_size and pair_offsets
2461 52769 : ALLOCATE (pair_offsets(SIZE(pairs)))
2462 17705 : total_size = 0
2463 724631 : DO i = 1, SIZE(pairs)
2464 706926 : pair_offsets(i) = total_size
2465 706926 : arow = pairs(i)%row
2466 706926 : acol = pairs(i)%col
2467 706926 : block_size = nsgf(arow)*nsgf(acol)
2468 724631 : total_size = total_size + block_size
2469 : END DO
2470 17705 : buffer_size = total_size
2471 :
2472 : ! calculate rank_offsets and rank_sizes
2473 53115 : ALLOCATE (rank_offsets(group_size))
2474 35410 : ALLOCATE (rank_sizes(group_size))
2475 51889 : rank_offsets = 0
2476 51889 : rank_sizes = 0
2477 17705 : IF (SIZE(pairs) > 0) THEN
2478 17359 : prev_k = pairs(1)%rank + 1
2479 724285 : DO i = 1, SIZE(pairs)
2480 706926 : k = pairs(i)%rank + 1
2481 706926 : CPASSERT(k >= prev_k) ! expecting the pairs to be ordered by rank
2482 724285 : IF (k > prev_k) THEN
2483 65 : rank_offsets(k) = pair_offsets(i)
2484 65 : rank_sizes(prev_k) = rank_offsets(k) - rank_offsets(prev_k)
2485 65 : prev_k = k
2486 : END IF
2487 : END DO
2488 17359 : rank_sizes(k) = buffer_size - rank_offsets(k) ! complete last rank
2489 : END IF
2490 :
2491 17705 : END SUBROUTINE rs_calc_offsets
2492 :
2493 : ! **************************************************************************************************
2494 : !> \brief Group contributions to each matrix block without changing the communication order.
2495 : !> Rebuilt with the task list; contains indices only, independent of matrices and threads.
2496 : !> \param pairs atom pairs in communication order
2497 : !> \param block_starts first position of each block in pair_indices, plus an end sentinel
2498 : !> \param pair_indices original pair indices grouped by destination block
2499 : ! **************************************************************************************************
2500 17655 : SUBROUTINE rs_group_atom_pairs(pairs, block_starts, pair_indices)
2501 : TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: pairs
2502 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: block_starts, pair_indices
2503 :
2504 : INTEGER :: i, nblocks, npairs
2505 17655 : INTEGER, ALLOCATABLE, DIMENSION(:) :: starts
2506 17655 : TYPE(atom_pair_type), ALLOCATABLE, DIMENSION(:) :: sorted_pairs
2507 :
2508 17655 : npairs = SIZE(pairs)
2509 829656 : ALLOCATE (sorted_pairs(npairs), pair_indices(npairs), starts(npairs + 1))
2510 724414 : sorted_pairs(:) = pairs(:)
2511 : ! Ignore source rank for grouping. The stable sort preserves contribution order within each block.
2512 724414 : sorted_pairs(:)%rank = 0
2513 17655 : CALL atom_pair_sort(sorted_pairs, npairs, pair_indices)
2514 :
2515 17655 : nblocks = 0
2516 724414 : DO i = 1, npairs
2517 706759 : IF (i > 1) THEN
2518 689448 : IF (.NOT. atom_pair_less_than(sorted_pairs(i - 1), sorted_pairs(i))) CYCLE
2519 : END IF
2520 706686 : nblocks = nblocks + 1
2521 724414 : starts(nblocks) = i
2522 : END DO
2523 17655 : starts(nblocks + 1) = npairs + 1
2524 70620 : ALLOCATE (block_starts(nblocks + 1))
2525 741996 : block_starts(:) = starts(:nblocks + 1)
2526 :
2527 17655 : END SUBROUTINE rs_group_atom_pairs
2528 :
2529 : ! **************************************************************************************************
2530 : !> \brief Scatters dbcsr matrix blocks and receives them into a buffer as needed before collocation.
2531 : !> \author Ole Schuett
2532 : ! **************************************************************************************************
2533 298 : SUBROUTINE rs_scatter_matrices(src_matrices, dest_buffer, task_list, group)
2534 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2535 : TYPE(offload_buffer_type), INTENT(INOUT) :: dest_buffer
2536 : TYPE(task_list_type), INTENT(IN) :: task_list
2537 : TYPE(mp_comm_type), INTENT(IN) :: group
2538 :
2539 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rs_scatter_matrices'
2540 :
2541 : INTEGER :: handle
2542 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: buffer_send
2543 :
2544 298 : CALL timeset(routineN, handle)
2545 873 : ALLOCATE (buffer_send(task_list%buffer_size_send))
2546 :
2547 : ! pack dbcsr blocks into send buffer
2548 298 : CPASSERT(ASSOCIATED(task_list%atom_pair_send))
2549 : CALL rs_pack_buffer(src_matrices=src_matrices, &
2550 : dest_buffer=buffer_send, &
2551 : atom_pair=task_list%atom_pair_send, &
2552 298 : pair_offsets=task_list%pair_offsets_send)
2553 :
2554 : ! mpi all-to-all communication, receiving directly into blocks_recv%buffer.
2555 : CALL group%alltoall(buffer_send, task_list%rank_sizes_send, task_list%rank_offsets_send, &
2556 : dest_buffer%host_buffer, &
2557 2682 : task_list%rank_sizes_recv, task_list%rank_offsets_recv)
2558 :
2559 298 : DEALLOCATE (buffer_send)
2560 298 : CALL timestop(handle)
2561 :
2562 298 : END SUBROUTINE rs_scatter_matrices
2563 :
2564 : ! **************************************************************************************************
2565 : !> \brief Gather the dbcsr matrix blocks and receives them into a buffer as needed after integration.
2566 : !> \author Ole Schuett
2567 : ! **************************************************************************************************
2568 248 : SUBROUTINE rs_gather_matrices(src_buffer, dest_matrices, task_list, group)
2569 : TYPE(offload_buffer_type), INTENT(IN) :: src_buffer
2570 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2571 : TYPE(task_list_type), INTENT(IN) :: task_list
2572 : TYPE(mp_comm_type), INTENT(IN) :: group
2573 :
2574 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rs_gather_matrices'
2575 :
2576 : INTEGER :: handle
2577 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: buffer_send
2578 :
2579 248 : CALL timeset(routineN, handle)
2580 :
2581 : ! Caution: The meaning of send and recv are reversed in this routine.
2582 730 : ALLOCATE (buffer_send(task_list%buffer_size_send)) ! e.g. this is actually used for receiving
2583 :
2584 : ! mpi all-to-all communication
2585 : CALL group%alltoall(src_buffer%host_buffer, task_list%rank_sizes_recv, task_list%rank_offsets_recv, &
2586 2232 : buffer_send, task_list%rank_sizes_send, task_list%rank_offsets_send)
2587 :
2588 : ! unpack dbcsr blocks from send buffer
2589 248 : CPASSERT(ASSOCIATED(task_list%atom_pair_send))
2590 : CALL rs_unpack_buffer(src_buffer=buffer_send, &
2591 : dest_matrices=dest_matrices, &
2592 : atom_pair=task_list%atom_pair_send, &
2593 : pair_offsets=task_list%pair_offsets_send, &
2594 : block_starts=task_list%unpack_block_starts, &
2595 248 : pair_indices=task_list%unpack_pair_indices)
2596 :
2597 248 : DEALLOCATE (buffer_send)
2598 248 : CALL timestop(handle)
2599 :
2600 248 : END SUBROUTINE rs_gather_matrices
2601 :
2602 : ! **************************************************************************************************
2603 : !> \brief Copies the DBCSR blocks into buffer, replaces rs_scatter_matrix for non-distributed grids.
2604 : !> \author Ole Schuett
2605 : ! **************************************************************************************************
2606 289785 : SUBROUTINE rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
2607 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2608 : TYPE(offload_buffer_type), INTENT(INOUT) :: dest_buffer
2609 : TYPE(task_list_type), INTENT(IN) :: task_list
2610 :
2611 : CALL rs_pack_buffer(src_matrices=src_matrices, &
2612 : dest_buffer=dest_buffer%host_buffer, &
2613 : atom_pair=task_list%atom_pair_recv, &
2614 289785 : pair_offsets=task_list%pair_offsets_recv)
2615 :
2616 289785 : END SUBROUTINE rs_copy_to_buffer
2617 :
2618 : ! **************************************************************************************************
2619 : !> \brief Copies from buffer into DBCSR matrics, replaces rs_gather_matrix for non-distributed grids.
2620 : !> \author Ole Schuett
2621 : ! **************************************************************************************************
2622 242594 : SUBROUTINE rs_copy_to_matrices(src_buffer, dest_matrices, task_list)
2623 : TYPE(offload_buffer_type), INTENT(IN) :: src_buffer
2624 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2625 : TYPE(task_list_type), INTENT(IN) :: task_list
2626 :
2627 : CALL rs_unpack_buffer(src_buffer=src_buffer%host_buffer, &
2628 : dest_matrices=dest_matrices, &
2629 : atom_pair=task_list%atom_pair_recv, &
2630 : pair_offsets=task_list%pair_offsets_recv, &
2631 : block_starts=task_list%unpack_block_starts, &
2632 242594 : pair_indices=task_list%unpack_pair_indices)
2633 :
2634 242594 : END SUBROUTINE rs_copy_to_matrices
2635 :
2636 : ! **************************************************************************************************
2637 : !> \brief Helper routine for rs_scatter_matrix and rs_copy_to_buffer.
2638 : !> \author Ole Schuett
2639 : ! **************************************************************************************************
2640 290083 : SUBROUTINE rs_pack_buffer(src_matrices, dest_buffer, atom_pair, pair_offsets)
2641 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: src_matrices
2642 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: dest_buffer
2643 : TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: atom_pair
2644 : INTEGER, DIMENSION(:), INTENT(IN) :: pair_offsets
2645 :
2646 : INTEGER :: acol, arow, img, i, offset, block_size
2647 : LOGICAL :: found
2648 290083 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2649 :
2650 : !$OMP PARALLEL DEFAULT(NONE), &
2651 : !$OMP SHARED(src_matrices,atom_pair,pair_offsets,dest_buffer), &
2652 290083 : !$OMP PRIVATE(acol,arow,img,i,offset,block_size,found,block)
2653 : !$OMP DO schedule(guided)
2654 : DO i = 1, SIZE(atom_pair)
2655 : arow = atom_pair(i)%row
2656 : acol = atom_pair(i)%col
2657 : img = atom_pair(i)%image
2658 : CALL dbcsr_get_block_p(matrix=src_matrices(img)%matrix, row=arow, col=acol, &
2659 : block=block, found=found)
2660 : CPASSERT(found)
2661 : block_size = SIZE(block)
2662 : offset = pair_offsets(i)
2663 : dest_buffer(offset + 1:offset + block_size) = RESHAPE(block, shape=[block_size])
2664 : END DO
2665 : !$OMP END DO
2666 : !$OMP END PARALLEL
2667 :
2668 290083 : END SUBROUTINE rs_pack_buffer
2669 :
2670 : ! **************************************************************************************************
2671 : !> \brief Helper routine for rs_gather_matrix and rs_copy_to_matrices.
2672 : !> \author Ole Schuett
2673 : ! **************************************************************************************************
2674 242842 : SUBROUTINE rs_unpack_buffer(src_buffer, dest_matrices, atom_pair, pair_offsets, block_starts, pair_indices)
2675 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: src_buffer
2676 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dest_matrices
2677 : TYPE(atom_pair_type), DIMENSION(:), INTENT(IN) :: atom_pair
2678 : INTEGER, DIMENSION(:), INTENT(IN) :: pair_offsets
2679 : INTEGER, DIMENSION(:), INTENT(IN) :: block_starts, pair_indices
2680 :
2681 : INTEGER :: acol, arow, i, iblock, img, ipair, &
2682 : j, ncols, nrows, offset, pair
2683 : LOGICAL :: found
2684 242842 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
2685 :
2686 : ! One thread owns each destination block and accumulates all ranks' contributions.
2687 : !$OMP PARALLEL DEFAULT(NONE), &
2688 : !$OMP SHARED(src_buffer,atom_pair,pair_offsets,dest_matrices,block_starts,pair_indices), &
2689 242842 : !$OMP PRIVATE(acol,arow,i,iblock,img,ipair,j,ncols,nrows,offset,pair,found,block)
2690 : !$OMP DO schedule(guided)
2691 : DO iblock = 1, SIZE(block_starts) - 1
2692 : pair = pair_indices(block_starts(iblock))
2693 : arow = atom_pair(pair)%row
2694 : acol = atom_pair(pair)%col
2695 : img = atom_pair(pair)%image
2696 : CALL dbcsr_get_block_p(matrix=dest_matrices(img)%matrix, row=arow, col=acol, &
2697 : block=block, found=found)
2698 : CPASSERT(found)
2699 : nrows = SIZE(block, 1)
2700 : ncols = SIZE(block, 2)
2701 : DO ipair = block_starts(iblock), block_starts(iblock + 1) - 1
2702 : offset = pair_offsets(pair_indices(ipair))
2703 : DO j = 1, ncols
2704 : DO i = 1, nrows
2705 : block(i, j) = block(i, j) + src_buffer(offset + i + (j - 1)*nrows)
2706 : END DO
2707 : END DO
2708 : END DO
2709 : END DO
2710 : !$OMP END DO
2711 : !$OMP END PARALLEL
2712 :
2713 242842 : END SUBROUTINE rs_unpack_buffer
2714 :
2715 : ! **************************************************************************************************
2716 : !> \brief determines the rank of the processor who's real rs grid contains point
2717 : !> \param rs_desc ...
2718 : !> \param igrid_level ...
2719 : !> \param n_levels ...
2720 : !> \param cube_center ...
2721 : !> \param ntasks ...
2722 : !> \param tasks ...
2723 : !> \param lb_cube ...
2724 : !> \param ub_cube ...
2725 : !> \param added_tasks ...
2726 : !> \par History
2727 : !> 11.2007 created [MattW]
2728 : !> 10.2008 rewritten [Joost VandeVondele]
2729 : !> \author MattW
2730 : ! **************************************************************************************************
2731 1476 : SUBROUTINE rs_find_node(rs_desc, igrid_level, n_levels, cube_center, ntasks, tasks, &
2732 : lb_cube, ub_cube, added_tasks)
2733 :
2734 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
2735 : INTEGER, INTENT(IN) :: igrid_level, n_levels
2736 : INTEGER, DIMENSION(3), INTENT(IN) :: cube_center
2737 : INTEGER, INTENT(INOUT) :: ntasks
2738 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
2739 : INTEGER, DIMENSION(3), INTENT(IN) :: lb_cube, ub_cube
2740 : INTEGER, INTENT(OUT) :: added_tasks
2741 :
2742 : INTEGER, PARAMETER :: add_tasks = 1000
2743 : REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp
2744 :
2745 : INTEGER :: bit_index, coord(3), curr_tasks, dest, i, icoord(3), idest, itask, ix, iy, iz, &
2746 : lb_coord(3), lb_domain(3), lbc(3), ub_coord(3), ub_domain(3), ubc(3)
2747 : INTEGER :: bit_pattern
2748 : LOGICAL :: dir_periodic(3)
2749 :
2750 1476 : coord(1) = rs_desc%x2coord(cube_center(1))
2751 1476 : coord(2) = rs_desc%y2coord(cube_center(2))
2752 1476 : coord(3) = rs_desc%z2coord(cube_center(3))
2753 1476 : dest = rs_desc%coord2rank(coord(1), coord(2), coord(3))
2754 :
2755 : ! the real cube coordinates
2756 5904 : lbc = lb_cube + cube_center
2757 5904 : ubc = ub_cube + cube_center
2758 :
2759 10332 : IF (ALL((rs_desc%lb_global(:, dest) - rs_desc%border) <= lbc) .AND. &
2760 : ALL((rs_desc%ub_global(:, dest) + rs_desc%border) >= ubc)) THEN
2761 : !standard distributed collocation/integration
2762 1476 : tasks(ntasks)%destination = encode_rank(dest, igrid_level, n_levels)
2763 1476 : tasks(ntasks)%dist_type = 1
2764 1476 : tasks(ntasks)%subpatch_pattern = 0
2765 1476 : added_tasks = 1
2766 :
2767 : ! here we figure out if there is an alternate location for this task
2768 : ! i.e. even though the cube_center is not in the real local domain,
2769 : ! it might fully fit in the halo of the neighbor
2770 : ! if its radius is smaller than the maximum radius
2771 : ! the list of possible neighbors is stored here as a bit pattern
2772 : ! to reconstruct what the rank of the neigbor is,
2773 : ! we can use (note this requires the correct rs_grid)
2774 : ! IF (BTEST(bit_pattern,0)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[-1,0,0])
2775 : ! IF (BTEST(bit_pattern,1)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[+1,0,0])
2776 : ! IF (BTEST(bit_pattern,2)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,-1,0])
2777 : ! IF (BTEST(bit_pattern,3)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,+1,0])
2778 : ! IF (BTEST(bit_pattern,4)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,0,-1])
2779 : ! IF (BTEST(bit_pattern,5)) rank=rs_grid_locate_rank(rs_desc,tasks(ntasks)%destination,[0,0,+1])
2780 1476 : bit_index = 0
2781 1476 : bit_pattern = 0
2782 5904 : DO i = 1, 3
2783 5904 : IF (rs_desc%perd(i) == 1) THEN
2784 2952 : bit_pattern = IBCLR(bit_pattern, bit_index)
2785 2952 : bit_index = bit_index + 1
2786 2952 : bit_pattern = IBCLR(bit_pattern, bit_index)
2787 2952 : bit_index = bit_index + 1
2788 : ELSE
2789 : ! fits the left neighbor ?
2790 1476 : IF (ubc(i) <= rs_desc%lb_global(i, dest) - 1 + rs_desc%border) THEN
2791 695 : bit_pattern = IBSET(bit_pattern, bit_index)
2792 695 : bit_index = bit_index + 1
2793 : ELSE
2794 781 : bit_pattern = IBCLR(bit_pattern, bit_index)
2795 781 : bit_index = bit_index + 1
2796 : END IF
2797 : ! fits the right neighbor ?
2798 1476 : IF (lbc(i) >= rs_desc%ub_global(i, dest) + 1 - rs_desc%border) THEN
2799 250 : bit_pattern = IBSET(bit_pattern, bit_index)
2800 250 : bit_index = bit_index + 1
2801 : ELSE
2802 1226 : bit_pattern = IBCLR(bit_pattern, bit_index)
2803 1226 : bit_index = bit_index + 1
2804 : END IF
2805 : END IF
2806 : END DO
2807 1476 : tasks(ntasks)%subpatch_pattern = bit_pattern
2808 :
2809 : ELSE
2810 : ! generalised collocation/integration needed
2811 : ! first we figure out how many neighbors we have to add to include the lbc/ubc
2812 : ! in the available domains (inclusive of halo points)
2813 : ! first we 'ignore' periodic boundary conditions
2814 : ! i.e. ub_coord-lb_coord+1 might be larger than group_dim
2815 0 : lb_coord = coord
2816 0 : ub_coord = coord
2817 0 : lb_domain = rs_desc%lb_global(:, dest) - rs_desc%border
2818 0 : ub_domain = rs_desc%ub_global(:, dest) + rs_desc%border
2819 0 : DO i = 1, 3
2820 : ! only if the grid is not periodic in this direction we need to take care of adding neighbors
2821 0 : IF (rs_desc%perd(i) == 0) THEN
2822 : ! if the domain lower bound is greater than the lbc we need to add the size of the neighbor domain
2823 0 : DO
2824 0 : IF (lb_domain(i) > lbc(i)) THEN
2825 0 : lb_coord(i) = lb_coord(i) - 1
2826 0 : icoord = MODULO(lb_coord, rs_desc%group_dim)
2827 0 : idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2828 0 : lb_domain(i) = lb_domain(i) - (rs_desc%ub_global(i, idest) - rs_desc%lb_global(i, idest) + 1)
2829 : ELSE
2830 : EXIT
2831 : END IF
2832 : END DO
2833 : ! same for the upper bound
2834 0 : DO
2835 0 : IF (ub_domain(i) < ubc(i)) THEN
2836 0 : ub_coord(i) = ub_coord(i) + 1
2837 0 : icoord = MODULO(ub_coord, rs_desc%group_dim)
2838 0 : idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2839 0 : ub_domain(i) = ub_domain(i) + (rs_desc%ub_global(i, idest) - rs_desc%lb_global(i, idest) + 1)
2840 : ELSE
2841 : EXIT
2842 : END IF
2843 : END DO
2844 : END IF
2845 : END DO
2846 :
2847 : ! some care is needed for the periodic boundaries
2848 0 : DO i = 1, 3
2849 0 : IF (ub_domain(i) - lb_domain(i) + 1 >= rs_desc%npts(i)) THEN
2850 0 : dir_periodic(i) = .TRUE.
2851 0 : lb_coord(i) = 0
2852 0 : ub_coord(i) = rs_desc%group_dim(i) - 1
2853 : ELSE
2854 0 : dir_periodic(i) = .FALSE.
2855 : END IF
2856 : END DO
2857 :
2858 0 : added_tasks = PRODUCT(ub_coord - lb_coord + 1)
2859 0 : itask = ntasks
2860 0 : ntasks = ntasks + added_tasks - 1
2861 0 : IF (ntasks > SIZE(tasks)) THEN
2862 0 : curr_tasks = INT((SIZE(tasks) + add_tasks)*mult_tasks)
2863 0 : CALL reallocate_tasks(tasks, curr_tasks)
2864 : END IF
2865 0 : DO iz = lb_coord(3), ub_coord(3)
2866 0 : DO iy = lb_coord(2), ub_coord(2)
2867 0 : DO ix = lb_coord(1), ub_coord(1)
2868 0 : icoord = MODULO([ix, iy, iz], rs_desc%group_dim)
2869 0 : idest = rs_desc%coord2rank(icoord(1), icoord(2), icoord(3))
2870 0 : tasks(itask)%destination = encode_rank(idest, igrid_level, n_levels)
2871 0 : tasks(itask)%dist_type = 2
2872 0 : tasks(itask)%subpatch_pattern = 0
2873 : ! encode the domain size for this task
2874 : ! if the bit is set, we need to add the border in that direction
2875 0 : IF (ix == lb_coord(1) .AND. .NOT. dir_periodic(1)) THEN
2876 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 0)
2877 : END IF
2878 0 : IF (ix == ub_coord(1) .AND. .NOT. dir_periodic(1)) THEN
2879 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 1)
2880 : END IF
2881 0 : IF (iy == lb_coord(2) .AND. .NOT. dir_periodic(2)) THEN
2882 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 2)
2883 : END IF
2884 0 : IF (iy == ub_coord(2) .AND. .NOT. dir_periodic(2)) THEN
2885 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 3)
2886 : END IF
2887 0 : IF (iz == lb_coord(3) .AND. .NOT. dir_periodic(3)) THEN
2888 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 4)
2889 : END IF
2890 0 : IF (iz == ub_coord(3) .AND. .NOT. dir_periodic(3)) THEN
2891 0 : tasks(itask)%subpatch_pattern = IBSET(tasks(itask)%subpatch_pattern, 5)
2892 : END IF
2893 0 : itask = itask + 1
2894 : END DO
2895 : END DO
2896 : END DO
2897 : END IF
2898 :
2899 1476 : END SUBROUTINE rs_find_node
2900 :
2901 : ! **************************************************************************************************
2902 : !> \brief utility functions for encoding the grid level with a rank, allowing
2903 : !> different grid levels to maintain a different rank ordering without
2904 : !> losing information. These encoded_ints are stored in the tasks(1:2,:) array
2905 : !> \param rank ...
2906 : !> \param grid_level ...
2907 : !> \param n_levels ...
2908 : !> \return ...
2909 : !> \par History
2910 : !> 4.2009 created [Iain Bethune]
2911 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2912 : ! **************************************************************************************************
2913 21649731 : FUNCTION encode_rank(rank, grid_level, n_levels) RESULT(encoded_int)
2914 :
2915 : INTEGER, INTENT(IN) :: rank, grid_level, n_levels
2916 : INTEGER :: encoded_int
2917 :
2918 : ! ordered so can still sort by rank
2919 :
2920 21649731 : encoded_int = rank*n_levels + grid_level - 1
2921 :
2922 21649731 : END FUNCTION encode_rank
2923 :
2924 : ! **************************************************************************************************
2925 : !> \brief ...
2926 : !> \param encoded_int ...
2927 : !> \param n_levels ...
2928 : !> \return ...
2929 : ! **************************************************************************************************
2930 10873716 : FUNCTION decode_rank(encoded_int, n_levels) RESULT(rank)
2931 :
2932 : INTEGER, INTENT(IN) :: encoded_int
2933 : INTEGER, INTENT(IN) :: n_levels
2934 : INTEGER :: rank
2935 :
2936 10873716 : rank = INT(encoded_int/n_levels)
2937 :
2938 10873716 : END FUNCTION decode_rank
2939 :
2940 : ! **************************************************************************************************
2941 : !> \brief ...
2942 : !> \param encoded_int ...
2943 : !> \param n_levels ...
2944 : !> \return ...
2945 : ! **************************************************************************************************
2946 7692 : FUNCTION decode_level(encoded_int, n_levels) RESULT(grid_level)
2947 :
2948 : INTEGER, INTENT(IN) :: encoded_int
2949 : INTEGER, INTENT(IN) :: n_levels
2950 : INTEGER :: grid_level
2951 :
2952 7692 : grid_level = INT(MODULO(encoded_int, n_levels)) + 1
2953 :
2954 7692 : END FUNCTION decode_level
2955 :
2956 : ! **************************************************************************************************
2957 : !> \brief Sort pgf index pair (ipgf,iset,iatom),(jpgf,jset,jatom) for which all atom pairs are
2958 : !> grouped, and for each atom pair all set pairs are grouped and for each set pair,
2959 : !> all pgfs are grouped.
2960 : !> This yields the right order of the tasks for collocation after the sort
2961 : !> in distribute_tasks. E.g. for a atom pair, all sets and pgfs are computed in one go.
2962 : !> The exception is the gridlevel. Tasks are first ordered wrt to grid_level. This implies
2963 : !> that a given density matrix block will be decontracted several times,
2964 : !> but cache effects on the grid make up for this.
2965 : !> \param a ...
2966 : !> \param b ...
2967 : !> \return ...
2968 : !> \author Ole Schuett
2969 : ! **************************************************************************************************
2970 104134820 : PURE FUNCTION tasks_less_than(a, b) RESULT(res)
2971 : TYPE(task_type), INTENT(IN) :: a, b
2972 : LOGICAL :: res
2973 :
2974 104134820 : IF (a%grid_level /= b%grid_level) THEN
2975 25212610 : res = a%grid_level < b%grid_level
2976 :
2977 78922210 : ELSE IF (a%image /= b%image) THEN
2978 17755254 : res = a%image < b%image
2979 :
2980 61166956 : ELSE IF (a%iatom /= b%iatom) THEN
2981 8440495 : res = a%iatom < b%iatom
2982 :
2983 52726461 : ELSE IF (a%jatom /= b%jatom) THEN
2984 9274828 : res = a%jatom < b%jatom
2985 :
2986 43451633 : ELSE IF (a%iset /= b%iset) THEN
2987 6509869 : res = a%iset < b%iset
2988 :
2989 36941764 : ELSE IF (a%jset /= b%jset) THEN
2990 6854303 : res = a%jset < b%jset
2991 :
2992 30087461 : ELSE IF (a%ipgf /= b%ipgf) THEN
2993 8495656 : res = a%ipgf < b%ipgf
2994 :
2995 : ELSE
2996 21591805 : res = a%jpgf < b%jpgf
2997 :
2998 : END IF
2999 104134820 : END FUNCTION tasks_less_than
3000 :
3001 248287722 : #:call array_sort(prefix='tasks', type='TYPE(task_type)')
3002 : #:endcall
3003 :
3004 : ! **************************************************************************************************
3005 : !> \brief Order atom pairs to find duplicates.
3006 : !> \param a ...
3007 : !> \param b ...
3008 : !> \return ...
3009 : !> \author Ole Schuett
3010 : ! **************************************************************************************************
3011 54922940 : PURE FUNCTION atom_pair_less_than(a, b) RESULT(res)
3012 : TYPE(atom_pair_type), INTENT(IN) :: a, b
3013 : LOGICAL :: res
3014 :
3015 54922940 : IF (a%rank /= b%rank) THEN
3016 4524 : res = a%rank < b%rank
3017 :
3018 54918416 : ELSE IF (a%row /= b%row) THEN
3019 7675520 : res = a%row < b%row
3020 :
3021 47242896 : ELSE IF (a%col /= b%col) THEN
3022 7617174 : res = a%col < b%col
3023 :
3024 : ELSE
3025 39625722 : res = a%image < b%image
3026 :
3027 : END IF
3028 54922940 : END FUNCTION atom_pair_less_than
3029 :
3030 117832197 : #:call array_sort(prefix='atom_pair', type='TYPE(atom_pair_type)')
3031 : #:endcall
3032 :
3033 : END MODULE task_list_methods
|
