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 Distribution methods for atoms, particles, or molecules
10 : !> \par History
11 : !> - 1d-distribution of molecules and particles (Sep. 2003, MK)
12 : !> - 2d-distribution for Quickstep updated with molecules (Oct. 2003, MK)
13 : !> \author MK (22.08.2003)
14 : ! **************************************************************************************************
15 : MODULE distribution_methods
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind,&
18 : get_atomic_kind_set
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_type
21 : USE cell_types, ONLY: cell_type,&
22 : pbc,&
23 : real_to_scaled,&
24 : scaled_to_real
25 : USE cp_array_utils, ONLY: cp_1d_i_p_type
26 : USE cp_blacs_env, ONLY: cp_blacs_env_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_distribution_get_num_images
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_get_default_unit_nr,&
31 : cp_logger_type
32 : USE cp_min_heap, ONLY: cp_heap_fill,&
33 : cp_heap_get_first,&
34 : cp_heap_new,&
35 : cp_heap_release,&
36 : cp_heap_reset_first,&
37 : cp_heap_type
38 : USE cp_output_handling, ONLY: cp_p_file,&
39 : cp_print_key_finished_output,&
40 : cp_print_key_should_output,&
41 : cp_print_key_unit_nr
42 : USE distribution_1d_types, ONLY: distribution_1d_create,&
43 : distribution_1d_type
44 : USE distribution_2d_types, ONLY: distribution_2d_create,&
45 : distribution_2d_type,&
46 : distribution_2d_write
47 : USE input_constants, ONLY: model_block_count,&
48 : model_block_lmax
49 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
50 : section_vals_type,&
51 : section_vals_val_get
52 : USE kinds, ONLY: dp,&
53 : int_8
54 : USE machine, ONLY: m_flush
55 : USE mathconstants, ONLY: pi
56 : USE mathlib, ONLY: gcd,&
57 : lcm
58 : USE molecule_kind_types, ONLY: get_molecule_kind,&
59 : get_molecule_kind_set,&
60 : molecule_kind_type
61 : USE molecule_types, ONLY: molecule_type
62 : USE parallel_rng_types, ONLY: UNIFORM,&
63 : rng_stream_type
64 : USE particle_types, ONLY: particle_type
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : qs_kind_type
67 : USE util, ONLY: sort
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : ! *** Global parameters (in this module) ***
75 :
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'distribution_methods'
77 :
78 : ! *** Public subroutines ***
79 :
80 : PUBLIC :: distribute_molecules_1d, &
81 : distribute_molecules_2d
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief Distribute molecules and particles
87 : !> \param atomic_kind_set particle (atomic) kind information
88 : !> \param particle_set particle information
89 : !> \param local_particles distribution of particles created by this routine
90 : !> \param molecule_kind_set molecule kind information
91 : !> \param molecule_set molecule information
92 : !> \param local_molecules distribution of molecules created by this routine
93 : !> \param force_env_section ...
94 : !> \param prev_molecule_kind_set previous molecule kind information, used with
95 : !> prev_local_molecules
96 : !> \param prev_local_molecules previous distribution of molecules, new one will
97 : !> be identical if all the prev_* arguments are present and associated
98 : !> \par History
99 : !> none
100 : !> \author MK (Jun. 2003)
101 : ! **************************************************************************************************
102 12551 : SUBROUTINE distribute_molecules_1d(atomic_kind_set, particle_set, &
103 : local_particles, &
104 : molecule_kind_set, molecule_set, &
105 : local_molecules, force_env_section, &
106 : prev_molecule_kind_set, &
107 : prev_local_molecules)
108 :
109 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
110 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
111 : TYPE(distribution_1d_type), POINTER :: local_particles
112 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
113 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
114 : TYPE(distribution_1d_type), POINTER :: local_molecules
115 : TYPE(section_vals_type), POINTER :: force_env_section
116 : TYPE(molecule_kind_type), DIMENSION(:), OPTIONAL, &
117 : POINTER :: prev_molecule_kind_set
118 : TYPE(distribution_1d_type), OPTIONAL, POINTER :: prev_local_molecules
119 :
120 : CHARACTER(len=*), PARAMETER :: routineN = 'distribute_molecules_1d'
121 :
122 : INTEGER :: atom_a, bin, handle, iatom, imolecule, imolecule_kind, imolecule_local, &
123 : imolecule_prev_kind, iparticle_kind, ipe, iw, kind_a, molecule_a, n, natom, nbins, nload, &
124 : nmolecule, nmolecule_kind, nparticle_kind, nsgf, output_unit
125 : INTEGER(int_8) :: bin_price
126 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: workload_count, workload_fill
127 12551 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nmolecule_local, nparticle_local, work
128 12551 : INTEGER, DIMENSION(:), POINTER :: molecule_list
129 : LOGICAL :: found, has_prev_subsys_info, is_local
130 12551 : TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: local_molecule
131 : TYPE(cp_heap_type) :: bin_heap_count, bin_heap_fill
132 : TYPE(cp_logger_type), POINTER :: logger
133 : TYPE(molecule_kind_type), POINTER :: molecule_kind
134 :
135 12551 : CALL timeset(routineN, handle)
136 :
137 12551 : has_prev_subsys_info = .FALSE.
138 12551 : IF (PRESENT(prev_local_molecules) .AND. &
139 : PRESENT(prev_molecule_kind_set)) THEN
140 2637 : IF (ASSOCIATED(prev_local_molecules) .AND. &
141 : ASSOCIATED(prev_molecule_kind_set)) THEN
142 44 : has_prev_subsys_info = .TRUE.
143 : END IF
144 : END IF
145 :
146 12551 : logger => cp_get_default_logger()
147 :
148 : ASSOCIATE (group => logger%para_env, mype => logger%para_env%mepos + 1, &
149 : npe => logger%para_env%num_pe)
150 :
151 37653 : ALLOCATE (workload_count(npe))
152 12551 : workload_count(:) = 0
153 :
154 25102 : ALLOCATE (workload_fill(npe))
155 12551 : workload_fill(:) = 0
156 :
157 12551 : nmolecule_kind = SIZE(molecule_kind_set)
158 :
159 37653 : ALLOCATE (nmolecule_local(nmolecule_kind))
160 12551 : nmolecule_local(:) = 0
161 :
162 200135 : ALLOCATE (local_molecule(nmolecule_kind))
163 :
164 12551 : nparticle_kind = SIZE(atomic_kind_set)
165 :
166 37653 : ALLOCATE (nparticle_local(nparticle_kind))
167 12551 : nparticle_local(:) = 0
168 :
169 12551 : nbins = npe
170 :
171 12551 : CALL cp_heap_new(bin_heap_count, nbins)
172 12551 : CALL cp_heap_fill(bin_heap_count, workload_count)
173 :
174 12551 : CALL cp_heap_new(bin_heap_fill, nbins)
175 12551 : CALL cp_heap_fill(bin_heap_fill, workload_fill)
176 :
177 162482 : DO imolecule_kind = 1, nmolecule_kind
178 :
179 149931 : molecule_kind => molecule_kind_set(imolecule_kind)
180 :
181 149931 : NULLIFY (molecule_list)
182 :
183 : ! *** Get the number of molecules and the number of ***
184 : ! *** atoms in each molecule of that molecular kind ***
185 :
186 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
187 : molecule_list=molecule_list, &
188 : natom=natom, &
189 149931 : nsgf=nsgf)
190 :
191 : ! *** Consider the number of atoms or basis ***
192 : ! *** functions which depends on the method ***
193 :
194 149931 : nload = MAX(natom, nsgf)
195 149931 : nmolecule = SIZE(molecule_list)
196 :
197 : ! *** Get the number of local molecules of the current molecule kind ***
198 :
199 468705 : DO imolecule = 1, nmolecule
200 468705 : IF (has_prev_subsys_info) THEN
201 297660 : DO imolecule_prev_kind = 1, SIZE(prev_molecule_kind_set)
202 3188004 : IF (ANY(prev_local_molecules%list(imolecule_prev_kind)%array( &
203 5060 : 1:prev_local_molecules%n_el(imolecule_prev_kind)) == molecule_list(imolecule))) THEN
204 : ! molecule used to be local
205 2530 : nmolecule_local(imolecule_kind) = nmolecule_local(imolecule_kind) + 1
206 : END IF
207 : END DO
208 : ELSE
209 313714 : CALL cp_heap_get_first(bin_heap_count, bin, bin_price, found)
210 313714 : IF (.NOT. found) THEN
211 0 : CPABORT("No topmost heap element found.")
212 : END IF
213 :
214 313714 : ipe = bin
215 313714 : IF (bin_price /= workload_count(ipe)) THEN
216 0 : CPABORT("inconsistent heap")
217 : END IF
218 :
219 313714 : workload_count(ipe) = workload_count(ipe) + nload
220 313714 : IF (ipe == mype) THEN
221 172041 : nmolecule_local(imolecule_kind) = nmolecule_local(imolecule_kind) + 1
222 : END IF
223 :
224 313714 : bin_price = workload_count(ipe)
225 313714 : CALL cp_heap_reset_first(bin_heap_count, bin_price)
226 : END IF
227 : END DO
228 :
229 : ! *** Distribute the molecules ***
230 149931 : n = nmolecule_local(imolecule_kind)
231 :
232 149931 : IF (n > 0) THEN
233 240306 : ALLOCATE (local_molecule(imolecule_kind)%array(n))
234 : ELSE
235 69829 : NULLIFY (local_molecule(imolecule_kind)%array)
236 : END IF
237 :
238 : imolecule_local = 0
239 631187 : DO imolecule = 1, nmolecule
240 318774 : is_local = .FALSE.
241 318774 : IF (has_prev_subsys_info) THEN
242 297660 : DO imolecule_prev_kind = 1, SIZE(prev_molecule_kind_set)
243 3188004 : IF (ANY(prev_local_molecules%list(imolecule_prev_kind)%array( &
244 5060 : 1:prev_local_molecules%n_el(imolecule_prev_kind)) == molecule_list(imolecule))) THEN
245 2530 : is_local = .TRUE.
246 : END IF
247 : END DO
248 : ELSE
249 313714 : CALL cp_heap_get_first(bin_heap_fill, bin, bin_price, found)
250 313714 : IF (.NOT. found) THEN
251 0 : CPABORT("No topmost heap element found.")
252 : END IF
253 :
254 313714 : ipe = bin
255 313714 : IF (bin_price /= workload_fill(ipe)) THEN
256 0 : CPABORT("inconsistent heap")
257 : END IF
258 :
259 313714 : workload_fill(ipe) = workload_fill(ipe) + nload
260 313714 : is_local = (ipe == mype)
261 : END IF
262 318774 : IF (is_local) THEN
263 174571 : imolecule_local = imolecule_local + 1
264 174571 : molecule_a = molecule_list(imolecule)
265 174571 : local_molecule(imolecule_kind)%array(imolecule_local) = molecule_a
266 595747 : DO iatom = 1, natom
267 421176 : atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
268 :
269 : CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
270 421176 : kind_number=kind_a)
271 595747 : nparticle_local(kind_a) = nparticle_local(kind_a) + 1
272 : END DO
273 : END IF
274 468705 : IF (.NOT. has_prev_subsys_info) THEN
275 313714 : bin_price = workload_fill(ipe)
276 313714 : CALL cp_heap_reset_first(bin_heap_fill, bin_price)
277 : END IF
278 : END DO
279 :
280 : END DO
281 :
282 36594 : IF (ANY(workload_fill /= workload_count)) THEN
283 0 : CPABORT("Inconsistent heaps encountered")
284 : END IF
285 :
286 12551 : CALL cp_heap_release(bin_heap_count)
287 12551 : CALL cp_heap_release(bin_heap_fill)
288 :
289 : ! *** Create the local molecule structure ***
290 :
291 : CALL distribution_1d_create(local_molecules, &
292 : n_el=nmolecule_local, &
293 12551 : para_env=logger%para_env)
294 :
295 : ! *** Create the local particle structure ***
296 :
297 : CALL distribution_1d_create(local_particles, &
298 : n_el=nparticle_local, &
299 12551 : para_env=logger%para_env)
300 :
301 : ! *** Store the generated local molecule and particle distributions ***
302 :
303 12551 : nparticle_local(:) = 0
304 :
305 162482 : DO imolecule_kind = 1, nmolecule_kind
306 :
307 149931 : IF (nmolecule_local(imolecule_kind) == 0) CYCLE
308 :
309 : local_molecules%list(imolecule_kind)%array(:) = &
310 254673 : local_molecule(imolecule_kind)%array(:)
311 :
312 80102 : molecule_kind => molecule_kind_set(imolecule_kind)
313 :
314 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
315 80102 : natom=natom)
316 :
317 267224 : DO imolecule = 1, nmolecule_local(imolecule_kind)
318 174571 : molecule_a = local_molecule(imolecule_kind)%array(imolecule)
319 745678 : DO iatom = 1, natom
320 421176 : atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
321 : CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
322 421176 : kind_number=kind_a)
323 421176 : nparticle_local(kind_a) = nparticle_local(kind_a) + 1
324 595747 : local_particles%list(kind_a)%array(nparticle_local(kind_a)) = atom_a
325 : END DO
326 : END DO
327 :
328 : END DO
329 :
330 : ! *** Print distribution, if requested ***
331 :
332 12551 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
333 37653 : force_env_section, "PRINT%DISTRIBUTION1D"), cp_p_file)) THEN
334 :
335 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION1D", &
336 158 : extension=".Log")
337 :
338 158 : iw = output_unit
339 158 : IF (output_unit < 0) iw = cp_logger_get_default_unit_nr(logger, LOCAL=.TRUE.)
340 :
341 : ! *** Print molecule distribution ***
342 :
343 474 : ALLOCATE (work(npe))
344 158 : work(:) = 0
345 :
346 586 : work(mype) = SUM(nmolecule_local)
347 158 : CALL group%sum(work)
348 :
349 158 : IF (output_unit > 0) THEN
350 : WRITE (UNIT=output_unit, &
351 : FMT="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
352 79 : "DISTRIBUTION OF THE MOLECULES", &
353 79 : "Process Number of molecules", &
354 316 : (ipe - 1, work(ipe), ipe=1, npe)
355 : WRITE (UNIT=output_unit, FMT="(T55, A3, T73, I8)") &
356 237 : "Sum", SUM(work)
357 79 : CALL m_flush(output_unit)
358 : END IF
359 :
360 158 : CALL group%sync()
361 :
362 474 : DO ipe = 1, npe
363 316 : IF (ipe == mype) THEN
364 : WRITE (UNIT=iw, FMT="(/, T3, A)") &
365 158 : "Process Kind Local molecules (global indices)"
366 586 : DO imolecule_kind = 1, nmolecule_kind
367 586 : IF (imolecule_kind == 1) THEN
368 : WRITE (UNIT=iw, FMT="(T4, I6, 2X, I5, (T21, 10I6))") &
369 158 : ipe - 1, imolecule_kind, &
370 372 : (local_molecules%list(imolecule_kind)%array(imolecule), &
371 688 : imolecule=1, nmolecule_local(imolecule_kind))
372 : ELSE
373 : WRITE (UNIT=iw, FMT="(T12, I5, (T21, 10I6))") &
374 270 : imolecule_kind, &
375 405 : (local_molecules%list(imolecule_kind)%array(imolecule), &
376 945 : imolecule=1, nmolecule_local(imolecule_kind))
377 : END IF
378 : END DO
379 : END IF
380 316 : CALL m_flush(iw)
381 474 : CALL group%sync()
382 : END DO
383 :
384 : ! *** Print particle distribution ***
385 :
386 158 : work(:) = 0
387 :
388 530 : work(mype) = SUM(nparticle_local)
389 158 : CALL group%sum(work)
390 :
391 158 : IF (output_unit > 0) THEN
392 : WRITE (UNIT=output_unit, &
393 : FMT="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
394 79 : "DISTRIBUTION OF THE PARTICLES", &
395 79 : "Process Number of particles", &
396 316 : (ipe - 1, work(ipe), ipe=1, npe)
397 : WRITE (UNIT=output_unit, FMT="(T55, A3, T73, I8)") &
398 237 : "Sum", SUM(work)
399 79 : CALL m_flush(output_unit)
400 : END IF
401 :
402 158 : CALL group%sync()
403 :
404 474 : DO ipe = 1, npe
405 316 : IF (ipe == mype) THEN
406 : WRITE (UNIT=iw, FMT="(/, T3, A)") &
407 158 : "Process Kind Local particles (global indices)"
408 530 : DO iparticle_kind = 1, nparticle_kind
409 530 : IF (iparticle_kind == 1) THEN
410 : WRITE (UNIT=iw, FMT="(T4, I6, 2X, I5, (T20, 10I6))") &
411 158 : ipe - 1, iparticle_kind, &
412 658 : (local_particles%list(iparticle_kind)%array(iatom), &
413 974 : iatom=1, nparticle_local(iparticle_kind))
414 : ELSE
415 : WRITE (UNIT=iw, FMT="(T12, I5, (T20, 10I6))") &
416 214 : iparticle_kind, &
417 1149 : (local_particles%list(iparticle_kind)%array(iatom), &
418 1577 : iatom=1, nparticle_local(iparticle_kind))
419 : END IF
420 : END DO
421 : END IF
422 316 : CALL m_flush(iw)
423 474 : CALL group%sync()
424 : END DO
425 158 : DEALLOCATE (work)
426 :
427 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
428 158 : "PRINT%DISTRIBUTION1D")
429 : END IF
430 : END ASSOCIATE
431 : ! *** Release work storage ***
432 :
433 12551 : DEALLOCATE (workload_count)
434 :
435 12551 : DEALLOCATE (workload_fill)
436 :
437 12551 : DEALLOCATE (nmolecule_local)
438 :
439 12551 : DEALLOCATE (nparticle_local)
440 :
441 162482 : DO imolecule_kind = 1, nmolecule_kind
442 162482 : IF (ASSOCIATED(local_molecule(imolecule_kind)%array)) THEN
443 80102 : DEALLOCATE (local_molecule(imolecule_kind)%array)
444 : END IF
445 : END DO
446 12551 : DEALLOCATE (local_molecule)
447 :
448 12551 : CALL timestop(handle)
449 :
450 25102 : END SUBROUTINE distribute_molecules_1d
451 :
452 : ! **************************************************************************************************
453 : !> \brief Distributes the particle pairs creating a 2d distribution optimally
454 : !> suited for quickstep
455 : !> \param cell ...
456 : !> \param atomic_kind_set ...
457 : !> \param particle_set ...
458 : !> \param qs_kind_set ...
459 : !> \param molecule_kind_set ...
460 : !> \param molecule_set ...
461 : !> \param distribution_2d the distribution that will be created by this
462 : !> method
463 : !> \param blacs_env the parallel environment at the basis of the
464 : !> distribution
465 : !> \param force_env_section ...
466 : !> \par History
467 : !> - local_rows & cols blocksize optimizations (Aug. 2003, MK)
468 : !> - cleanup of distribution_2d (Sep. 2003, fawzi)
469 : !> - update for molecules (Oct. 2003, MK)
470 : !> \author fawzi (Feb. 2003)
471 : !> \note
472 : !> Intermediate generation of a 2d distribution of the molecules, but
473 : !> only the corresponding particle (atomic) distribution is currently
474 : !> used. The 2d distribution of the molecules is deleted, but may easily
475 : !> be recovered (MK).
476 : ! **************************************************************************************************
477 9724 : SUBROUTINE distribute_molecules_2d(cell, atomic_kind_set, particle_set, &
478 : qs_kind_set, molecule_kind_set, molecule_set, &
479 : distribution_2d, blacs_env, force_env_section)
480 : TYPE(cell_type), POINTER :: cell
481 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
482 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
483 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
484 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
485 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
486 : TYPE(distribution_2d_type), POINTER :: distribution_2d
487 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
488 : TYPE(section_vals_type), POINTER :: force_env_section
489 :
490 : CHARACTER(len=*), PARAMETER :: routineN = 'distribute_molecules_2d'
491 :
492 : INTEGER :: cluster_price, cost_model, handle, iatom, iatom_mol, iatom_one, ikind, imol, &
493 : imolecule, imolecule_kind, iparticle_kind, ipcol, iprow, iw, kind_a, n, natom, natom_mol, &
494 : nclusters, nmolecule, nmolecule_kind, nparticle_kind, nsgf, output_unit
495 9724 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_list, cluster_prices, &
496 9724 : nparticle_local_col, &
497 9724 : nparticle_local_row, work
498 9724 : INTEGER, DIMENSION(:), POINTER :: lmax_basis, molecule_list
499 9724 : INTEGER, DIMENSION(:, :), POINTER :: cluster_col_distribution, &
500 9724 : cluster_row_distribution, &
501 9724 : col_distribution, row_distribution
502 : LOGICAL :: basic_cluster_optimization, basic_optimization, basic_spatial_optimization, &
503 : molecular_distribution, skip_optimization
504 9724 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coords, pbc_scaled_coords
505 : REAL(KIND=dp), DIMENSION(3) :: center
506 9724 : TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: local_particle_col, local_particle_row
507 : TYPE(cp_logger_type), POINTER :: logger
508 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
509 : TYPE(molecule_kind_type), POINTER :: molecule_kind
510 : TYPE(section_vals_type), POINTER :: distribution_section
511 :
512 : !...
513 :
514 9724 : CALL timeset(routineN, handle)
515 :
516 9724 : logger => cp_get_default_logger()
517 :
518 9724 : distribution_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DISTRIBUTION")
519 :
520 9724 : CALL section_vals_val_get(distribution_section, "2D_MOLECULAR_DISTRIBUTION", l_val=molecular_distribution)
521 9724 : CALL section_vals_val_get(distribution_section, "SKIP_OPTIMIZATION", l_val=skip_optimization)
522 9724 : CALL section_vals_val_get(distribution_section, "BASIC_OPTIMIZATION", l_val=basic_optimization)
523 9724 : CALL section_vals_val_get(distribution_section, "BASIC_SPATIAL_OPTIMIZATION", l_val=basic_spatial_optimization)
524 9724 : CALL section_vals_val_get(distribution_section, "BASIC_CLUSTER_OPTIMIZATION", l_val=basic_cluster_optimization)
525 :
526 9724 : CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
527 : !
528 :
529 : ASSOCIATE (group => blacs_env%para_env, myprow => blacs_env%mepos(1) + 1, mypcol => blacs_env%mepos(2) + 1, &
530 : nprow => blacs_env%num_pe(1), npcol => blacs_env%num_pe(2))
531 :
532 9724 : nmolecule_kind = SIZE(molecule_kind_set)
533 9724 : CALL get_molecule_kind_set(molecule_kind_set, nmolecule=nmolecule)
534 :
535 9724 : nparticle_kind = SIZE(atomic_kind_set)
536 9724 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
537 :
538 : !
539 : ! we need to generate two representations of the distribution, one as a straight array with global particles
540 : ! one ordered wrt to kinds and only listing the local particles
541 : !
542 29172 : ALLOCATE (row_distribution(natom, 2))
543 19448 : ALLOCATE (col_distribution(natom, 2))
544 : ! Initialize the distributions to -1, as the second dimension only gets set with cluster optimization
545 : ! but the information is needed by dbcsr
546 270008 : row_distribution = -1; col_distribution = -1
547 :
548 47433 : ALLOCATE (local_particle_col(nparticle_kind))
549 37709 : ALLOCATE (local_particle_row(nparticle_kind))
550 29172 : ALLOCATE (nparticle_local_row(nparticle_kind))
551 19448 : ALLOCATE (nparticle_local_col(nparticle_kind))
552 :
553 9724 : IF (basic_optimization .OR. basic_spatial_optimization .OR. basic_cluster_optimization) THEN
554 :
555 9724 : IF (molecular_distribution) THEN
556 2 : nclusters = nmolecule
557 : ELSE
558 : nclusters = natom
559 : END IF
560 :
561 29172 : ALLOCATE (cluster_list(nclusters))
562 19448 : ALLOCATE (cluster_prices(nclusters))
563 29172 : ALLOCATE (cluster_row_distribution(nclusters, 2))
564 19448 : ALLOCATE (cluster_col_distribution(nclusters, 2))
565 279700 : cluster_row_distribution = -1; cluster_col_distribution = -1
566 :
567 : ! Fill in the clusters and their prices
568 9724 : CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
569 9724 : IF (.NOT. molecular_distribution) THEN
570 65057 : DO iatom = 1, natom
571 55335 : IF (iatom > nclusters) THEN
572 0 : CPABORT("Bounds error")
573 : END IF
574 55335 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
575 55335 : cluster_list(iatom) = iatom
576 55335 : SELECT CASE (cost_model)
577 : CASE (model_block_count)
578 55335 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
579 55335 : cluster_price = nsgf
580 : CASE (model_block_lmax)
581 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
582 0 : CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
583 0 : cluster_price = MAXVAL(lmax_basis)
584 : CASE default
585 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
586 0 : CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
587 55335 : cluster_price = 8 + (MAXVAL(lmax_basis)**2)
588 : END SELECT
589 120392 : cluster_prices(iatom) = cluster_price
590 : END DO
591 : ELSE
592 : imol = 0
593 4 : DO imolecule_kind = 1, nmolecule_kind
594 2 : molecule_kind => molecule_kind_set(imolecule_kind)
595 2 : CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
596 8 : DO imolecule = 1, SIZE(molecule_list)
597 4 : imol = imol + 1
598 4 : cluster_list(imol) = imol
599 4 : cluster_price = 0
600 16 : DO iatom_mol = 1, natom_mol
601 12 : iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
602 12 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
603 16 : SELECT CASE (cost_model)
604 : CASE (model_block_count)
605 12 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
606 12 : cluster_price = cluster_price + nsgf
607 : CASE (model_block_lmax)
608 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
609 0 : CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
610 0 : cluster_price = cluster_price + MAXVAL(lmax_basis)
611 : CASE default
612 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
613 0 : CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
614 12 : cluster_price = cluster_price + 8 + (MAXVAL(lmax_basis)**2)
615 : END SELECT
616 : END DO
617 6 : cluster_prices(imol) = cluster_price
618 : END DO
619 : END DO
620 : END IF
621 :
622 : ! And distribute
623 9724 : IF (basic_optimization) THEN
624 : CALL make_basic_distribution(cluster_list, cluster_prices, &
625 9556 : nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
626 : ELSE
627 168 : IF (basic_cluster_optimization) THEN
628 4 : IF (molecular_distribution) THEN
629 0 : CPABORT("clustering and molecular blocking NYI")
630 : END IF
631 16 : ALLOCATE (pbc_scaled_coords(3, natom), coords(3, natom))
632 118 : DO iatom = 1, natom
633 114 : CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
634 118 : coords(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
635 : END DO
636 : CALL make_cluster_distribution(coords, pbc_scaled_coords, cell, cluster_prices, &
637 4 : nprow, cluster_row_distribution, npcol, cluster_col_distribution)
638 : ELSE ! basic_spatial_optimization
639 492 : ALLOCATE (pbc_scaled_coords(3, nclusters))
640 164 : IF (.NOT. molecular_distribution) THEN
641 : ! just scaled coords
642 2658 : DO iatom = 1, natom
643 2658 : CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
644 : END DO
645 : ELSE
646 : ! use scaled coords of geometric center, folding when appropriate
647 : imol = 0
648 0 : DO imolecule_kind = 1, nmolecule_kind
649 0 : molecule_kind => molecule_kind_set(imolecule_kind)
650 0 : CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
651 0 : DO imolecule = 1, SIZE(molecule_list)
652 0 : imol = imol + 1
653 0 : iatom_one = molecule_set(molecule_list(imolecule))%first_atom
654 0 : center = 0.0_dp
655 0 : DO iatom_mol = 1, natom_mol
656 0 : iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
657 : center = center + &
658 0 : pbc(particle_set(iatom)%r(:) - particle_set(iatom_one)%r(:), cell) + particle_set(iatom_one)%r(:)
659 : END DO
660 0 : center = center/natom_mol
661 0 : CALL real_to_scaled(pbc_scaled_coords(:, imol), pbc(center, cell), cell)
662 : END DO
663 : END DO
664 : END IF
665 :
666 : CALL make_basic_spatial_distribution(pbc_scaled_coords, cluster_prices, &
667 164 : nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
668 :
669 164 : DEALLOCATE (pbc_scaled_coords)
670 : END IF
671 : END IF
672 :
673 : ! And assign back
674 9724 : IF (.NOT. molecular_distribution) THEN
675 279672 : row_distribution = cluster_row_distribution
676 279672 : col_distribution = cluster_col_distribution
677 : ELSE
678 : imol = 0
679 4 : DO imolecule_kind = 1, nmolecule_kind
680 2 : molecule_kind => molecule_kind_set(imolecule_kind)
681 2 : CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
682 8 : DO imolecule = 1, SIZE(molecule_list)
683 4 : imol = imol + 1
684 18 : DO iatom_mol = 1, natom_mol
685 12 : iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
686 72 : row_distribution(iatom, :) = cluster_row_distribution(imol, :)
687 76 : col_distribution(iatom, :) = cluster_col_distribution(imol, :)
688 : END DO
689 : END DO
690 : END DO
691 : END IF
692 :
693 : ! cleanup
694 9724 : DEALLOCATE (cluster_list)
695 9724 : DEALLOCATE (cluster_prices)
696 9724 : DEALLOCATE (cluster_row_distribution)
697 19448 : DEALLOCATE (cluster_col_distribution)
698 :
699 : ELSE
700 : ! expects nothing else
701 0 : CPABORT("Invalid optimization for DFT%QS%DISTRIBUTION")
702 : END IF
703 :
704 : ! prepare the lists of local particles
705 :
706 : ! count local particles of a given kind
707 9724 : nparticle_local_col = 0
708 9724 : nparticle_local_row = 0
709 65071 : DO iatom = 1, natom
710 55347 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
711 55347 : IF (row_distribution(iatom, 1) == myprow) nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
712 120418 : IF (col_distribution(iatom, 1) == mypcol) nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
713 : END DO
714 :
715 : ! allocate space
716 27985 : DO iparticle_kind = 1, nparticle_kind
717 18261 : n = nparticle_local_row(iparticle_kind)
718 49435 : ALLOCATE (local_particle_row(iparticle_kind)%array(n))
719 :
720 18261 : n = nparticle_local_col(iparticle_kind)
721 64507 : ALLOCATE (local_particle_col(iparticle_kind)%array(n))
722 : END DO
723 :
724 : ! store
725 9724 : nparticle_local_col = 0
726 9724 : nparticle_local_row = 0
727 65071 : DO iatom = 1, natom
728 55347 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
729 55347 : IF (row_distribution(iatom, 1) == myprow) THEN
730 28749 : nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
731 28749 : local_particle_row(kind_a)%array(nparticle_local_row(kind_a)) = iatom
732 : END IF
733 120418 : IF (col_distribution(iatom, 1) == mypcol) THEN
734 55347 : nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
735 55347 : local_particle_col(kind_a)%array(nparticle_local_col(kind_a)) = iatom
736 : END IF
737 : END DO
738 :
739 : ! *** Generate the 2d distribution structure but take care of the zero offsets required
740 65071 : row_distribution(:, 1) = row_distribution(:, 1) - 1
741 65071 : col_distribution(:, 1) = col_distribution(:, 1) - 1
742 : CALL distribution_2d_create(distribution_2d, &
743 : row_distribution_ptr=row_distribution, &
744 : col_distribution_ptr=col_distribution, &
745 : local_rows_ptr=local_particle_row, &
746 : local_cols_ptr=local_particle_col, &
747 9724 : blacs_env=blacs_env)
748 :
749 9724 : NULLIFY (local_particle_row)
750 9724 : NULLIFY (local_particle_col)
751 9724 : NULLIFY (row_distribution)
752 9724 : NULLIFY (col_distribution)
753 :
754 : ! *** Print distribution, if requested ***
755 9724 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
756 9724 : force_env_section, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
757 :
758 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION", &
759 108 : extension=".Log")
760 :
761 : ! *** Print row distribution ***
762 :
763 324 : ALLOCATE (work(nprow))
764 108 : work(:) = 0
765 :
766 322 : IF (mypcol == 1) work(myprow) = SUM(distribution_2d%n_local_rows)
767 :
768 108 : CALL group%sum(work)
769 :
770 108 : IF (output_unit > 0) THEN
771 : WRITE (UNIT=output_unit, &
772 : FMT="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
773 54 : "DISTRIBUTION OF THE PARTICLES (ROWS)", &
774 54 : "Process row Number of particles Number of matrix rows", &
775 214 : (iprow - 1, work(iprow), -1, iprow=1, nprow)
776 : WRITE (UNIT=output_unit, FMT="(T23, A3, T41, I10, T71, I10)") &
777 160 : "Sum", SUM(work), -1
778 54 : CALL m_flush(output_unit)
779 : END IF
780 :
781 108 : DEALLOCATE (work)
782 :
783 : ! *** Print column distribution ***
784 :
785 324 : ALLOCATE (work(npcol))
786 108 : work(:) = 0
787 :
788 217 : IF (myprow == 1) work(mypcol) = SUM(distribution_2d%n_local_cols)
789 :
790 108 : CALL group%sum(work)
791 :
792 108 : IF (output_unit > 0) THEN
793 : WRITE (UNIT=output_unit, &
794 : FMT="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
795 54 : "DISTRIBUTION OF THE PARTICLES (COLUMNS)", &
796 54 : "Process col Number of particles Number of matrix columns", &
797 162 : (ipcol - 1, work(ipcol), -1, ipcol=1, npcol)
798 : WRITE (UNIT=output_unit, FMT="(T23, A3, T41, I10, T71, I10)") &
799 108 : "Sum", SUM(work), -1
800 54 : CALL m_flush(output_unit)
801 : END IF
802 :
803 108 : DEALLOCATE (work)
804 :
805 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
806 108 : "PRINT%DISTRIBUTION")
807 : END IF
808 : END ASSOCIATE
809 :
810 9724 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
811 : force_env_section, "PRINT%DISTRIBUTION2D"), cp_p_file)) THEN
812 :
813 98 : iw = cp_logger_get_default_unit_nr(logger, LOCAL=.TRUE.)
814 : CALL distribution_2d_write(distribution_2d, &
815 : unit_nr=iw, &
816 : local=.TRUE., &
817 98 : long_description=.TRUE.)
818 :
819 : END IF
820 :
821 : ! *** Release work storage ***
822 :
823 9724 : DEALLOCATE (nparticle_local_row)
824 :
825 9724 : DEALLOCATE (nparticle_local_col)
826 :
827 9724 : CALL timestop(handle)
828 :
829 29172 : END SUBROUTINE distribute_molecules_2d
830 :
831 : ! **************************************************************************************************
832 : !> \brief Creates a basic distribution
833 : !> \param cluster_list ...
834 : !> \param cluster_prices ...
835 : !> \param nprows ...
836 : !> \param row_distribution ...
837 : !> \param npcols ...
838 : !> \param col_distribution ...
839 : !> \par History
840 : !> - Created 2010-08-06 UB
841 : ! **************************************************************************************************
842 9556 : SUBROUTINE make_basic_distribution(cluster_list, cluster_prices, &
843 9556 : nprows, row_distribution, npcols, col_distribution)
844 : INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_list, cluster_prices
845 : INTEGER, INTENT(IN) :: nprows
846 : INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
847 : INTEGER, INTENT(IN) :: npcols
848 : INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
849 :
850 : CHARACTER(len=*), PARAMETER :: routineN = 'make_basic_distribution'
851 :
852 : INTEGER :: bin, cluster, cluster_index, &
853 : cluster_price, nbins, nclusters, pcol, &
854 : pgrid_gcd, prow, timing_handle
855 : INTEGER(int_8) :: bin_price
856 : LOGICAL :: found
857 : TYPE(cp_heap_type) :: bin_heap
858 :
859 : ! ---------------------------------------------------------------------------
860 :
861 9556 : CALL timeset(routineN, timing_handle)
862 9556 : nbins = lcm(nprows, npcols)
863 9556 : pgrid_gcd = gcd(nprows, npcols)
864 9556 : CALL sort(cluster_prices, SIZE(cluster_list), cluster_list)
865 9556 : CALL cp_heap_new(bin_heap, nbins)
866 46288 : CALL cp_heap_fill(bin_heap, [(0_int_8, bin=1, nbins)])
867 : !
868 9556 : nclusters = SIZE(cluster_list)
869 : ! Put the most expensive cluster in the bin with the smallest
870 : ! price and repeat.
871 62287 : DO cluster_index = nclusters, 1, -1
872 52731 : cluster = cluster_list(cluster_index)
873 52731 : CALL cp_heap_get_first(bin_heap, bin, bin_price, found)
874 52731 : IF (.NOT. found) THEN
875 0 : CPABORT("No topmost heap element found.")
876 : END IF
877 : !
878 52731 : prow = INT((bin - 1)*pgrid_gcd/npcols)
879 52731 : IF (prow >= nprows) THEN
880 0 : CPABORT("Invalid process row.")
881 : END IF
882 52731 : pcol = INT((bin - 1)*pgrid_gcd/nprows)
883 52731 : IF (pcol >= npcols) THEN
884 0 : CPABORT("Invalid process column.")
885 : END IF
886 52731 : row_distribution(cluster) = prow + 1
887 52731 : col_distribution(cluster) = pcol + 1
888 : !
889 52731 : cluster_price = cluster_prices(cluster_index)
890 52731 : bin_price = bin_price + cluster_price
891 115018 : CALL cp_heap_reset_first(bin_heap, bin_price)
892 : END DO
893 9556 : CALL cp_heap_release(bin_heap)
894 9556 : CALL timestop(timing_handle)
895 9556 : END SUBROUTINE make_basic_distribution
896 :
897 : ! **************************************************************************************************
898 : !> \brief Creates a basic spatial distribution
899 : !> that tries to make the corresponding blocks as homogeneous as possible
900 : !> \param pbc_scaled_coords ...
901 : !> \param costs ...
902 : !> \param nprows ...
903 : !> \param row_distribution ...
904 : !> \param npcols ...
905 : !> \param col_distribution ...
906 : !> \par History
907 : !> - Created 2010-11-11 Joost VandeVondele
908 : ! **************************************************************************************************
909 164 : SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, &
910 164 : nprows, row_distribution, npcols, col_distribution)
911 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
912 : INTEGER, DIMENSION(:), INTENT(IN) :: costs
913 : INTEGER, INTENT(IN) :: nprows
914 : INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
915 : INTEGER, INTENT(IN) :: npcols
916 : INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
917 :
918 : CHARACTER(len=*), PARAMETER :: routineN = 'make_basic_spatial_distribution'
919 :
920 : INTEGER :: handle, iatom, natoms, nbins, pgrid_gcd
921 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_costs, distribution
922 :
923 164 : CALL timeset(routineN, handle)
924 :
925 164 : natoms = SIZE(costs)
926 164 : nbins = lcm(nprows, npcols)
927 164 : pgrid_gcd = gcd(nprows, npcols)
928 820 : ALLOCATE (bin_costs(nbins), distribution(natoms))
929 164 : bin_costs = 0
930 :
931 5152 : CALL spatial_recurse(pbc_scaled_coords, costs, [(iatom, iatom=1, natoms)], bin_costs, distribution, 0)
932 :
933 : ! WRITE(*, *) "Final bin costs: ", bin_costs
934 :
935 : ! final row_distribution / col_distribution
936 2658 : DO iatom = 1, natoms
937 2494 : row_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/npcols + 1
938 2658 : col_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/nprows + 1
939 : END DO
940 :
941 164 : DEALLOCATE (bin_costs, distribution)
942 :
943 164 : CALL timestop(handle)
944 :
945 164 : END SUBROUTINE make_basic_spatial_distribution
946 :
947 : ! **************************************************************************************************
948 : !> \brief ...
949 : !> \param pbc_scaled_coords ...
950 : !> \param costs ...
951 : !> \param indices ...
952 : !> \param bin_costs ...
953 : !> \param distribution ...
954 : !> \param level ...
955 : ! **************************************************************************************************
956 3168 : RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_costs, distribution, level)
957 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
958 : INTEGER, DIMENSION(:), INTENT(IN) :: costs, indices
959 : INTEGER, DIMENSION(:), INTENT(INOUT) :: bin_costs, distribution
960 : INTEGER, INTENT(IN) :: level
961 :
962 : INTEGER :: iatom, ibin, natoms, nbins, nhalf
963 3168 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_costs_sorted, atom_permutation, &
964 3168 : bin_costs_sorted, permutation
965 3168 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coord
966 :
967 3168 : natoms = SIZE(costs)
968 3168 : nbins = SIZE(bin_costs)
969 3168 : nhalf = (natoms + 1)/2
970 :
971 3168 : IF (natoms <= nbins) THEN
972 : ! assign the most expensive atom to the least costly bin
973 6664 : ALLOCATE (bin_costs_sorted(nbins), permutation(nbins))
974 4998 : bin_costs_sorted(:) = bin_costs
975 1666 : CALL sort(bin_costs_sorted, nbins, permutation)
976 6664 : ALLOCATE (atom_costs_sorted(natoms), atom_permutation(natoms))
977 4160 : atom_costs_sorted(:) = costs
978 1666 : CALL sort(atom_costs_sorted, natoms, atom_permutation)
979 1666 : ibin = 0
980 : ! WRITE(*, *) "Dealing with a new bunch of atoms "
981 4160 : DO iatom = natoms, 1, -1
982 2494 : ibin = ibin + 1
983 : ! WRITE(*, *) "atom", indices(atom_permutation(iatom)), "cost", atom_costs_sorted(iatom), &
984 : ! "bin", permutation(ibin), "its cost", bin_costs(permutation(ibin))
985 : ! WRITE(100, '(A, I0, 3F12.6)') "A", permutation(ibin), pbc_scaled_coords(:, atom_permutation(iatom))
986 2494 : bin_costs(permutation(ibin)) = bin_costs(permutation(ibin)) + atom_costs_sorted(iatom)
987 4160 : distribution(indices(atom_permutation(iatom))) = permutation(ibin)
988 : END DO
989 1666 : DEALLOCATE (bin_costs_sorted, permutation, atom_costs_sorted, atom_permutation)
990 : ELSE
991 : ! divide atoms in two subsets, sorting according to their coordinates, alternatively x, y, z
992 : ! recursively do this for both subsets
993 7510 : ALLOCATE (coord(natoms), permutation(natoms))
994 12216 : coord(:) = pbc_scaled_coords(MOD(level, 3) + 1, :)
995 1502 : CALL sort(coord, natoms, permutation)
996 : CALL spatial_recurse(pbc_scaled_coords(:, permutation(1:nhalf)), costs(permutation(1:nhalf)), &
997 36098 : indices(permutation(1:nhalf)), bin_costs, distribution, level + 1)
998 : CALL spatial_recurse(pbc_scaled_coords(:, permutation(nhalf + 1:)), costs(permutation(nhalf + 1:)), &
999 31190 : indices(permutation(nhalf + 1:)), bin_costs, distribution, level + 1)
1000 1502 : DEALLOCATE (coord, permutation)
1001 : END IF
1002 :
1003 3168 : END SUBROUTINE spatial_recurse
1004 :
1005 : ! **************************************************************************************************
1006 : !> \brief creates a distribution placing close by atoms into clusters and
1007 : !> putting them on the same processors. Load balancing is
1008 : !> performed by balancing sum of the cluster costs per processor
1009 : !> \param coords coordinates of the system
1010 : !> \param scaled_coords scaled coordinates
1011 : !> \param cell the cell_type
1012 : !> \param costs costs per atomic block
1013 : !> \param nprows number of precessors per row on the 2d grid
1014 : !> \param row_distribution the resulting distribution over proc_rows of atomic blocks
1015 : !> \param npcols number of precessors per col on the 2d grid
1016 : !> \param col_distribution the resulting distribution over proc_cols of atomic blocks
1017 : ! **************************************************************************************************
1018 4 : SUBROUTINE make_cluster_distribution(coords, scaled_coords, cell, costs, &
1019 4 : nprows, row_distribution, npcols, col_distribution)
1020 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coords, scaled_coords
1021 : TYPE(cell_type), POINTER :: cell
1022 : INTEGER, DIMENSION(:), INTENT(IN) :: costs
1023 : INTEGER, INTENT(IN) :: nprows
1024 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: row_distribution
1025 : INTEGER, INTENT(IN) :: npcols
1026 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: col_distribution
1027 :
1028 : CHARACTER(len=*), PARAMETER :: routineN = 'make_cluster_distribution'
1029 :
1030 : INTEGER :: handle, i, icluster, level, natom, &
1031 : output_unit
1032 : INTEGER(KIND=int_8) :: ncluster
1033 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_to_cluster, cluster_cost, &
1034 4 : cluster_count, cluster_to_col, &
1035 : cluster_to_row, piv_cost, proc_cost, &
1036 4 : sorted_cost
1037 : REAL(KIND=dp) :: fold(3)
1038 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cluster_center, cluster_high, cluster_low
1039 :
1040 4 : CALL timeset(routineN, handle)
1041 :
1042 4 : output_unit = cp_logger_get_default_io_unit()
1043 :
1044 4 : natom = SIZE(costs)
1045 118 : ncluster = dbcsr_distribution_get_num_images(SUM(costs), natom, nprows, npcols)
1046 12 : ALLOCATE (atom_to_cluster(natom))
1047 12 : ALLOCATE (cluster_cost(ncluster))
1048 8 : ALLOCATE (cluster_to_row(ncluster))
1049 8 : ALLOCATE (cluster_to_col(ncluster))
1050 8 : ALLOCATE (sorted_cost(ncluster))
1051 8 : ALLOCATE (piv_cost(ncluster))
1052 4 : cluster_cost(:) = 0
1053 :
1054 4 : icluster = 0
1055 4 : CALL cluster_recurse(coords, scaled_coords, cell, costs, atom_to_cluster, ncluster, icluster, cluster_cost)
1056 :
1057 16 : sorted_cost(:) = cluster_cost(:)
1058 4 : CALL sort(sorted_cost, INT(ncluster), piv_cost)
1059 :
1060 12 : ALLOCATE (proc_cost(nprows))
1061 4 : proc_cost = 0; level = 1
1062 4 : CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_row, nprows)
1063 :
1064 12 : DEALLOCATE (proc_cost); ALLOCATE (proc_cost(npcols))
1065 4 : proc_cost = 0; level = 1
1066 4 : CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_col, npcols)
1067 :
1068 118 : DO i = 1, natom
1069 114 : row_distribution(i, 1) = cluster_to_row(atom_to_cluster(i))
1070 114 : row_distribution(i, 2) = atom_to_cluster(i)
1071 114 : col_distribution(i, 1) = cluster_to_col(atom_to_cluster(i))
1072 118 : col_distribution(i, 2) = atom_to_cluster(i)
1073 : END DO
1074 :
1075 : ! generate some statistics on clusters
1076 12 : ALLOCATE (cluster_center(3, ncluster))
1077 8 : ALLOCATE (cluster_low(3, ncluster))
1078 8 : ALLOCATE (cluster_high(3, ncluster))
1079 8 : ALLOCATE (cluster_count(ncluster))
1080 4 : cluster_count = 0
1081 118 : DO i = 1, natom
1082 114 : cluster_count(atom_to_cluster(i)) = cluster_count(atom_to_cluster(i)) + 1
1083 460 : cluster_center(:, atom_to_cluster(i)) = coords(:, i)
1084 : END DO
1085 52 : cluster_low = HUGE(0.0_dp)/2
1086 52 : cluster_high = -HUGE(0.0_dp)/2
1087 118 : DO i = 1, natom
1088 798 : fold = pbc(coords(:, i) - cluster_center(:, atom_to_cluster(i)), cell) + cluster_center(:, atom_to_cluster(i))
1089 456 : cluster_low(:, atom_to_cluster(i)) = MIN(cluster_low(:, atom_to_cluster(i)), fold(:))
1090 460 : cluster_high(:, atom_to_cluster(i)) = MAX(cluster_high(:, atom_to_cluster(i)), fold(:))
1091 : END DO
1092 4 : IF (output_unit > 0) THEN
1093 2 : WRITE (output_unit, *)
1094 2 : WRITE (output_unit, '(T2,A)') "Cluster distribution information"
1095 2 : WRITE (output_unit, '(T2,A,T48,I8)') "Number of atoms", natom
1096 2 : WRITE (output_unit, '(T2,A,T48,I8)') "Number of clusters", ncluster
1097 8 : WRITE (output_unit, '(T2,A,T48,I8)') "Largest cluster in atoms", MAXVAL(cluster_count)
1098 8 : WRITE (output_unit, '(T2,A,T48,I8)') "Smallest cluster in atoms", MINVAL(cluster_count)
1099 2 : WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster x=", &
1100 8 : MAXVAL(cluster_high(1, :) - cluster_low(1, :), MASK=(cluster_count > 0)), &
1101 14 : MAXLOC(cluster_high(1, :) - cluster_low(1, :), MASK=(cluster_count > 0))
1102 2 : WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster y=", &
1103 8 : MAXVAL(cluster_high(2, :) - cluster_low(2, :), MASK=(cluster_count > 0)), &
1104 14 : MAXLOC(cluster_high(2, :) - cluster_low(2, :), MASK=(cluster_count > 0))
1105 2 : WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster z=", &
1106 8 : MAXVAL(cluster_high(3, :) - cluster_low(3, :), MASK=(cluster_count > 0)), &
1107 14 : MAXLOC(cluster_high(3, :) - cluster_low(3, :), MASK=(cluster_count > 0))
1108 : END IF
1109 :
1110 4 : DEALLOCATE (atom_to_cluster, cluster_cost, cluster_to_row, cluster_to_col, sorted_cost, piv_cost, proc_cost)
1111 4 : CALL timestop(handle)
1112 :
1113 8 : END SUBROUTINE make_cluster_distribution
1114 :
1115 : ! **************************************************************************************************
1116 : !> \brief assigns the clusters to processors, tryimg to balance the cost on the nodes
1117 : !> \param cluster_cost vector with the cost of each cluster
1118 : !> \param piv_cost pivoting vector sorting the cluster_cost
1119 : !> \param proc_cost cost per processor, on input 0 everywhere
1120 : !> \param cluster_assign assgnment of clusters on proc
1121 : !> \param nproc number of processor over which clusters are distributed
1122 : ! **************************************************************************************************
1123 8 : SUBROUTINE assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_assign, nproc)
1124 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, piv_cost, proc_cost, &
1125 : cluster_assign
1126 : INTEGER :: nproc
1127 :
1128 : CHARACTER(len=*), PARAMETER :: routineN = 'assign_clusters'
1129 :
1130 : INTEGER :: handle, i, ilevel, offset, &
1131 16 : piv_pcost(nproc), sort_proc_cost(nproc)
1132 :
1133 8 : CALL timeset(routineN, handle)
1134 :
1135 26 : DO ilevel = 1, SIZE(cluster_cost)/nproc
1136 42 : sort_proc_cost(:) = proc_cost(:)
1137 18 : CALL sort(sort_proc_cost, nproc, piv_pcost)
1138 :
1139 18 : offset = (SIZE(cluster_cost)/nproc - ilevel + 1)*nproc + 1
1140 50 : DO i = 1, nproc
1141 24 : cluster_assign(piv_cost(offset - i)) = piv_pcost(i)
1142 42 : proc_cost(piv_pcost(i)) = proc_cost(piv_pcost(i)) + cluster_cost(piv_cost(offset - i))
1143 : END DO
1144 : END DO
1145 :
1146 8 : CALL timestop(handle)
1147 :
1148 8 : END SUBROUTINE assign_clusters
1149 :
1150 : ! **************************************************************************************************
1151 : !> \brief recursive routine to cluster atoms.
1152 : !> Low level uses a modified KMEANS algorithm
1153 : !> recursion is used to reduce cost.
1154 : !> each level will subdivide a cluster into smaller clusters
1155 : !> If only a single split is necessary atoms are assigned to the current cluster
1156 : !> \param coord coordinates of the system
1157 : !> \param scaled_coord scaled coordinates
1158 : !> \param cell the cell_type
1159 : !> \param costs costs per atomic block
1160 : !> \param cluster_inds the atom_to cluster mapping
1161 : !> \param ncluster number of clusters still to be created on a given recursion level
1162 : !> \param icluster the index of the current cluster to be created
1163 : !> \param fin_cluster_cost total cost of the final clusters
1164 : ! **************************************************************************************************
1165 16 : RECURSIVE SUBROUTINE cluster_recurse(coord, scaled_coord, cell, costs, cluster_inds, ncluster, icluster, fin_cluster_cost)
1166 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coord, scaled_coord
1167 : TYPE(cell_type), POINTER :: cell
1168 : INTEGER, DIMENSION(:), INTENT(IN) :: costs
1169 : INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_inds
1170 : INTEGER(KIND=int_8), INTENT(INOUT) :: ncluster
1171 : INTEGER, INTENT(INOUT) :: icluster
1172 : INTEGER, DIMENSION(:), INTENT(INOUT) :: fin_cluster_cost
1173 :
1174 : INTEGER :: i, ibeg, iend, maxv(1), min_seed, &
1175 : natoms, nleft, nsplits, seed, tot_cost
1176 16 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: ncluster_new
1177 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, inds_tmp, nat_cluster, piv
1178 : LOGICAL :: found
1179 : REAL(KIND=dp) :: balance, balance_new, conv
1180 :
1181 16 : natoms = SIZE(coord, 2)
1182 : ! This is a bit of an arbitrary choice, simply a try to avoid too many clusters on large systems and too few for balancing on
1183 : ! small systems or subclusters
1184 16 : IF (natoms <= 1) THEN
1185 2 : nsplits = 1
1186 : ELSE
1187 14 : nsplits = MIN(INT(MIN(INT(MAX(6, INT(60.00/LOG(REAL(natoms, KIND=dp)))), KIND=int_8), ncluster)), natoms)
1188 : END IF
1189 16 : IF (nsplits == 1) THEN
1190 12 : icluster = icluster + 1
1191 126 : cluster_inds = icluster
1192 126 : fin_cluster_cost(icluster) = SUM(costs)
1193 : ELSE
1194 36 : ALLOCATE (cluster_cost(nsplits), ncluster_new(nsplits), inds_tmp(natoms), piv(natoms), nat_cluster(nsplits))
1195 : ! initialise some values
1196 4 : cluster_cost = 0; seed = 300; found = .TRUE.; min_seed = seed
1197 4 : CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed, conv)
1198 28 : balance = MAXVAL(REAL(nat_cluster, KIND=dp))/MINVAL(REAL(nat_cluster, KIND=dp))
1199 :
1200 : ! If the system is small enough try to do better in terms of balancing number of atoms per cluster
1201 : ! by changing the seed for the initial guess
1202 4 : IF (natoms < 1000 .AND. balance > 1.1) THEN
1203 24 : found = .FALSE.
1204 24 : DO i = 1, 5
1205 24 : IF (balance > 1.1) THEN
1206 20 : CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed + i*40, conv)
1207 140 : balance_new = MAXVAL(REAL(nat_cluster, KIND=dp))/MINVAL(REAL(nat_cluster, KIND=dp))
1208 20 : IF (balance_new < balance) THEN
1209 0 : balance = balance_new
1210 0 : min_seed = seed + i*40
1211 : END IF
1212 : ELSE
1213 : found = .TRUE.
1214 : EXIT
1215 : END IF
1216 : END DO
1217 : END IF
1218 : !If we do not match the convergence than recompute at least the best assignment
1219 4 : IF (.NOT. found) CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, min_seed, conv)
1220 :
1221 : ! compute the cost of each cluster to decide how many splits have to be performed on the next lower level
1222 118 : DO i = 1, natoms
1223 118 : cluster_cost(cluster_inds(i)) = cluster_cost(cluster_inds(i)) + costs(i)
1224 : END DO
1225 16 : tot_cost = SUM(cluster_cost)
1226 : ! compute new splitting, can be done more elegant
1227 16 : ncluster_new(:) = ncluster*cluster_cost(:)/tot_cost
1228 16 : nleft = INT(ncluster - SUM(ncluster_new))
1229 : ! As we won't have empty clusters, we can not have 0 as new size, so we correct for this at first
1230 16 : DO i = 1, nsplits
1231 16 : IF (ncluster_new(i) == 0) THEN
1232 6 : ncluster_new(i) = 1
1233 6 : nleft = nleft - 1
1234 : END IF
1235 : END DO
1236 : ! now comes the next part that the number of clusters will not match anymore, so try to correct in a meaningful way without
1237 : ! introducing 0 sized blocks again
1238 4 : IF (nleft /= 0) THEN
1239 0 : DO i = 1, ABS(nleft)
1240 0 : IF (nleft < 0) THEN
1241 0 : maxv = MINLOC(cluster_cost/ncluster_new)
1242 0 : IF (ncluster_new(maxv(1)) /= 1) THEN
1243 0 : ncluster_new(maxv) = ncluster_new(maxv) - 1
1244 : ELSE
1245 0 : maxv = MAXLOC(ncluster_new)
1246 0 : ncluster_new(maxv) = ncluster_new(maxv) - 1
1247 : END IF
1248 : ELSE
1249 0 : maxv = MAXLOC(cluster_cost/ncluster_new)
1250 0 : ncluster_new(maxv) = ncluster_new(maxv) + 1
1251 : END IF
1252 : END DO
1253 : END IF
1254 :
1255 : !Now get the permutations to sort the atoms in the nsplits clusters for the next level of iteration
1256 118 : inds_tmp(:) = cluster_inds(:)
1257 4 : CALL sort(inds_tmp, natoms, piv)
1258 :
1259 4 : ibeg = 1; iend = 0
1260 16 : DO i = 1, nsplits
1261 12 : IF (nat_cluster(i) == 0) CYCLE
1262 12 : iend = iend + nat_cluster(i)
1263 : CALL cluster_recurse(coord(:, piv(ibeg:iend)), scaled_coord(:, piv(ibeg:iend)), cell, costs(piv(ibeg:iend)), &
1264 1038 : inds_tmp(ibeg:iend), ncluster_new(i), icluster, fin_cluster_cost)
1265 16 : ibeg = ibeg + nat_cluster(i)
1266 : END DO
1267 : ! copy the sorted cluster IDs on the old layout, inds_tmp gets set at the lowest level of recursion
1268 118 : cluster_inds(piv(:)) = inds_tmp
1269 4 : DEALLOCATE (cluster_cost, ncluster_new, inds_tmp, piv, nat_cluster)
1270 :
1271 : END IF
1272 :
1273 16 : END SUBROUTINE cluster_recurse
1274 :
1275 : ! **************************************************************************************************
1276 : !> \brief A modified version of the kmeans algorithm.
1277 : !> The assignment has a penalty function in case clusters become
1278 : !> larger than average. Like this more even sized clusters are created
1279 : !> trading it for locality
1280 : !> \param ncent number of centers to be created
1281 : !> \param coord coordinates
1282 : !> \param scaled_coord scaled coord
1283 : !> \param cell the cell_type
1284 : !> \param cluster atom to cluster assignment
1285 : !> \param nat_cl atoms per cluster
1286 : !> \param seed seed for the RNG. Algorithm might need multiple tries to deliver best results
1287 : !> \param tot_var the total variance of the clusters around the centers
1288 : ! **************************************************************************************************
1289 28 : SUBROUTINE kmeans(ncent, coord, scaled_coord, cell, cluster, nat_cl, seed, tot_var)
1290 : INTEGER :: ncent
1291 : REAL(KIND=dp), DIMENSION(:, :) :: coord, scaled_coord
1292 : TYPE(cell_type), POINTER :: cell
1293 : INTEGER, DIMENSION(:) :: cluster, nat_cl
1294 : INTEGER :: seed
1295 : REAL(KIND=dp) :: tot_var
1296 :
1297 : CHARACTER(len=*), PARAMETER :: routineN = 'kmeans'
1298 :
1299 : INTEGER :: handle, i, ind, itn, j, nat, oldc
1300 : LOGICAL :: changed
1301 56 : REAL(KIND=dp) :: average(3, ncent, 2), cent_coord(3, ncent), devi, deviat(ncent), dist, &
1302 56 : dvec(3), old_var, rn, scaled_cent(3, ncent), var_cl(ncent)
1303 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat
1304 : REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed
1305 : TYPE(rng_stream_type) :: rng_stream
1306 :
1307 28 : CALL timeset(routineN, handle)
1308 :
1309 252 : initial_seed = REAL(seed, dp); nat = SIZE(coord, 2)
1310 112 : ALLOCATE (dmat(ncent, nat))
1311 :
1312 : rng_stream = rng_stream_type(name="kmeans uniform distribution [0,1]", &
1313 28 : distribution_type=UNIFORM, seed=initial_seed)
1314 :
1315 : ! try to find a clever initial guess with centers being somewhat distributed
1316 28 : rn = rng_stream%next()
1317 28 : ind = CEILING(rn*nat)
1318 112 : cent_coord(:, 1) = coord(:, ind)
1319 84 : DO i = 2, ncent
1320 28 : DO
1321 928 : rn = rng_stream%next()
1322 928 : ind = CEILING(rn*nat)
1323 3712 : cent_coord(:, i) = coord(:, ind)
1324 : devi = HUGE(1.0_dp)
1325 2004 : DO j = 1, i - 1
1326 1076 : dvec = pbc(cent_coord(:, j), cent_coord(:, i), cell)
1327 4304 : dist = NORM2(dvec)
1328 928 : IF (dist < devi) devi = dist
1329 : END DO
1330 928 : rn = rng_stream%next()
1331 928 : IF (rn < devi**2/169.0) EXIT
1332 : END DO
1333 : END DO
1334 :
1335 : ! Now start the KMEANS but penalise it in case it starts packing too many atoms into a single set
1336 : ! Unfoirtunatelz as this is dependent on what happened before it cant be parallel
1337 826 : cluster = 0; old_var = HUGE(1.0_dp)
1338 106 : DO itn = 1, 1000
1339 1210 : changed = .FALSE.; var_cl = 0.0_dp; tot_var = 0.0_dp; nat_cl = 0; deviat = 0.0_dp
1340 : ! !$OMP PARALLEL DO PRIVATE(i,j,dvec)
1341 4402 : DO i = 1, nat
1342 21418 : DO j = 1, ncent
1343 17016 : dvec = pbc(cent_coord(:, j), coord(:, i), cell)
1344 72360 : dmat(j, i) = DOT_PRODUCT(dvec, dvec)
1345 : END DO
1346 : END DO
1347 4402 : DO i = 1, nat
1348 4296 : devi = HUGE(1.0_dp); oldc = cluster(i)
1349 21312 : DO j = 1, ncent
1350 17016 : dist = dmat(j, i) + MAX(nat_cl(j)**2/nat*ncent, nat/ncent)
1351 21312 : IF (dist < devi) THEN
1352 8760 : devi = dist; cluster(i) = j
1353 : END IF
1354 : END DO
1355 4296 : deviat(cluster(i)) = deviat(cluster(i)) + SQRT(devi)
1356 4296 : nat_cl(cluster(i)) = nat_cl(cluster(i)) + 1
1357 4296 : tot_var = tot_var + devi
1358 4402 : IF (oldc /= cluster(i)) changed = .TRUE.
1359 : END DO
1360 : ! get the update of the centers done, add a new one in case one center lost all its atoms
1361 : ! the algorithm would survive, but its nice to really create what you demand
1362 106 : IF (tot_var >= old_var) EXIT
1363 106 : IF (changed) THEN
1364 : ! Here misery of computing the center of geometry of the clusters in PBC.
1365 : ! The mapping on the unit circle allows to circumvent all problems
1366 2506 : average = 0.0_dp
1367 3576 : DO i = 1, SIZE(coord, 2)
1368 13992 : average(:, cluster(i), 1) = average(:, cluster(i), 1) + COS(scaled_coord(:, i)*2.0_dp*pi)
1369 14070 : average(:, cluster(i), 2) = average(:, cluster(i), 2) + SIN(scaled_coord(:, i)*2.0_dp*pi)
1370 : END DO
1371 :
1372 362 : DO i = 1, ncent
1373 362 : IF (nat_cl(i) == 0) THEN
1374 0 : rn = rng_stream%next()
1375 0 : scaled_cent(:, i) = scaled_coord(:, CEILING(rn*nat))
1376 : ELSE
1377 1136 : average(:, i, 1) = average(:, i, 1)/REAL(nat_cl(i), dp)
1378 1136 : average(:, i, 2) = average(:, i, 2)/REAL(nat_cl(i), dp)
1379 1136 : scaled_cent(:, i) = (ATAN2(-average(:, i, 2), -average(:, i, 1)) + pi)/(2.0_dp*pi)
1380 284 : CALL scaled_to_real(cent_coord(:, i), scaled_cent(:, i), cell)
1381 : END IF
1382 : END DO
1383 : ELSE
1384 : EXIT
1385 : END IF
1386 : END DO
1387 :
1388 28 : CALL timestop(handle)
1389 :
1390 728 : END SUBROUTINE kmeans
1391 :
1392 : END MODULE distribution_methods
|
