(git:4cc4584)
Loading...
Searching...
No Matches
cp_fm_diag_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Auxiliary tools to redistribute cp_fm_type and cp_cfm_type matrices before and after
10!> diagonalization. Heuristics are used to determine the optimal number of CPUs for
11!> diagonalization and the input matrices are redistributed if necessary
12!> \par History
13!> - [01.2018] moved redistribution related code from cp_fm_syevd here
14!> - [08.2026] added complex matrix (cp_cfm_type) counterparts
15!> \author Nico Holmberg [01.2018]
16! **************************************************************************************************
17MODULE cp_fm_diag_utils
18 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
19 cp_blacs_env_release,&
20 cp_blacs_env_type
21 USE cp_cfm_types, ONLY: cp_cfm_create,&
22 cp_cfm_get_info,&
23 cp_cfm_release,&
24 cp_cfm_type
25 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 cp_fm_struct_release,&
27 cp_fm_struct_type
28 USE cp_fm_types, ONLY: cp_fm_create,&
29 cp_fm_get_info,&
30 cp_fm_release,&
31 cp_fm_type
32 USE cp_log_handling, ONLY: cp_get_default_logger,&
33 cp_logger_get_default_io_unit,&
34 cp_logger_type
35 USE kinds, ONLY: dp
36 USE mathlib, ONLY: gcd
37 USE message_passing, ONLY: mp_para_env_type
38#include "../base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_diag_utils'
45
46 ! Information on redistribution
47 TYPE, PUBLIC :: cp_fm_redistribute_info
48 INTEGER :: matrix_order = -1
49 INTEGER :: num_pe_old = -1 ! number of processes before a potential redistribute
50 INTEGER :: num_pe_new = -1 ! number of processes after a potential redistribute
51 INTEGER :: num_pe_opt = -1 ! optimal number of processes for the given matrix
52 INTEGER :: num_pe_max_nz_col = -1 ! the maximal number of processes s.t. no column has zero width, may be < 0 if ignored
53 LOGICAL :: redistribute = .false. ! whether or not the matrix was actually redistributed
54 CONTAINS
55 PROCEDURE, pass(self) :: write => cp_fm_redistribute_info_write
56 END TYPE cp_fm_redistribute_info
57
58 ! Container for redistribution settings and temporary work structs
59 TYPE cp_fm_redistribute_type
60 ! Settings
61 INTEGER :: a = -1, x = -1
62 LOGICAL :: should_print = .false.
63 LOGICAL :: elpa_force_redistribute = .false.
64 ! Temporaries
65 INTEGER, DIMENSION(:), POINTER :: group_distribution => null(), &
66 group_partition => null()
67 TYPE(cp_blacs_env_type), POINTER :: blacs_env_new => null()
68 TYPE(mp_para_env_type), POINTER :: para_env_new => null()
69 END TYPE cp_fm_redistribute_type
70
71 ! Permanent instance of the redistribute type
72 TYPE(cp_fm_redistribute_type), PRIVATE, &
73 SAVE :: work_redistribute
74
75 ! Public subroutines
76
77 PUBLIC :: cp_fm_redistribute_start, &
78 cp_fm_redistribute_end, &
79 cp_fm_redistribute_init, &
80 cp_cfm_redistribute_start, &
81 cp_cfm_redistribute_end
82#if defined(__parallel)
83 PUBLIC :: cp_fm_max_ncpu_non_zero_column
84#endif
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief Write the redistribute info nicely formatted to the given I/O unit
90!> \param self reference to the cp_fm_redistribute_info instance
91!> \param io_unit I/O unit to use for writing
92! **************************************************************************************************
93 SUBROUTINE cp_fm_redistribute_info_write(self, io_unit)
94 CLASS(cp_fm_redistribute_info), INTENT(IN) :: self
95 INTEGER, INTENT(IN) :: io_unit
96
97 WRITE (unit=io_unit, fmt="(A)") ""
98 WRITE (unit=io_unit, fmt="(T2,A,T71,I10)") &
99 "CP_FM_DIAG| Number of processes over which the matrix is distributed ", self%num_pe_old, &
100 "CP_FM_DIAG| Matrix order ", self%matrix_order
101 WRITE (unit=io_unit, fmt="(T2,A,T71,I10)") &
102 "CP_FM_DIAG| Optimal number of CPUs ", self%num_pe_opt
103 IF (self%num_pe_max_nz_col < 0) THEN
104 WRITE (unit=io_unit, fmt="(T2,A,T71,A10)") &
105 "CP_FM_DIAG| Maximum number of CPUs (with non-zero columns) ", "<N/A>"
106 ELSE
107 WRITE (unit=io_unit, fmt="(T2,A,T71,I10)") &
108 "CP_FM_DIAG| Maximum number of CPUs (with non-zero columns): ", self%num_pe_max_nz_col
109 END IF
110 IF (self%redistribute) THEN
111 WRITE (unit=io_unit, fmt="(T2,A,T71,I10)") &
112 "CP_FM_DIAG| Number of processes for the redistribution ", self%num_pe_new
113 ELSE
114 WRITE (unit=io_unit, fmt="(T2,A)") &
115 "CP_FM_DIAG| The matrix will NOT be redistributed"
116 END IF
117 WRITE (unit=io_unit, fmt="(A)") ""
118
119 END SUBROUTINE cp_fm_redistribute_info_write
120
121! **************************************************************************************************
122!> \brief Releases the temporary storage needed when redistributing arrays
123!> \param has_redistributed flag that determines if the processors holds a part of the
124!> redistributed array
125!> \author Nico Holmberg [01.2018]
126! **************************************************************************************************
127 SUBROUTINE cp_fm_redistribute_work_finalize(has_redistributed)
128 LOGICAL, INTENT(IN) :: has_redistributed
129
130 IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
131 IF (has_redistributed) THEN
132 CALL cp_blacs_env_release(work_redistribute%blacs_env_new)
133 END IF
134 CALL work_redistribute%para_env_new%free()
135 DEALLOCATE (work_redistribute%para_env_new)
136 DEALLOCATE (work_redistribute%group_distribution)
137 DEALLOCATE (work_redistribute%group_partition)
138 END IF
139 ! Return the temporaries to their initial state; the settings written by
140 ! cp_fm_redistribute_init have to survive every redistribution
141 NULLIFY (work_redistribute%group_distribution)
142 NULLIFY (work_redistribute%group_partition)
143 NULLIFY (work_redistribute%blacs_env_new)
144 NULLIFY (work_redistribute%para_env_new)
145
146 END SUBROUTINE cp_fm_redistribute_work_finalize
147
148! **************************************************************************************************
149!> \brief Initializes the parameters that determine how to calculate the optimal number of CPUs
150!> for diagonalizing a matrix. The parameters are read from the GLOBAL input section.
151!> \param a integer parameter used to define the rule for determining the optimal
152!> number of CPUs for diagonalization
153!> \param x integer parameter used to define the rule for determining the optimal
154!> number of CPUs for diagonalization
155!> \param should_print flag that determines if information about the redistribution process
156!> should be printed
157!> \param elpa_force_redistribute flag that if redistribution should always be performed when
158!> the ELPA diagonalization library is in use
159!> \author Nico Holmberg [01.2018]
160! **************************************************************************************************
161 SUBROUTINE cp_fm_redistribute_init(a, x, should_print, elpa_force_redistribute)
162 INTEGER, INTENT(IN) :: a, x
163 LOGICAL, INTENT(IN) :: should_print, elpa_force_redistribute
164
165 cpassert(a > 0 .AND. x > 0)
166
167 ! Init work first: the structure constructor resets every component, so the
168 ! settings have to be written after it, not before
169 work_redistribute = cp_fm_redistribute_type()
170 work_redistribute%a = a
171 work_redistribute%x = x
172 work_redistribute%should_print = should_print
173 work_redistribute%elpa_force_redistribute = elpa_force_redistribute
174
175 END SUBROUTINE cp_fm_redistribute_init
176
177! **************************************************************************************************
178!> \brief Calculates the optimal number of CPUs for diagonalizing a matrix.
179!> \param size the size of the diagonalized matrix
180!> \return the optimal number of CPUs
181!> \author Nico Holmberg [01.2018]
182! **************************************************************************************************
183 PURE FUNCTION cp_fm_diag_get_optimal_ncpu(size) RESULT(ncpu)
184 INTEGER, INTENT(IN) :: size
185 INTEGER :: ncpu
186
187 ncpu = ((size + work_redistribute%a*work_redistribute%x - 1)/ &
188 (work_redistribute%a*work_redistribute%x))*work_redistribute%a
189
190 END FUNCTION cp_fm_diag_get_optimal_ncpu
191
192#if defined(__parallel)
193! **************************************************************************************************
194!> \brief Finds a safe ELPA layout at or below the requested number of CPUs.
195!> \param matrix_struct the original matrix layout
196!> \param ncpu_max the upper limit on the number of CPUs
197!> \param block_size the block size of the redistributed matrix
198!> \return a number of CPUs whose layout has no empty process rows or columns
199!> \author Nico Holmberg [01.2018]
200! **************************************************************************************************
201 FUNCTION cp_fm_max_ncpu_non_zero_column(matrix_struct, ncpu_max, block_size) RESULT(ncpu)
202 TYPE(cp_fm_struct_type), INTENT(IN) :: matrix_struct
203 INTEGER, INTENT(IN) :: ncpu_max, block_size
204 INTEGER :: ncpu
205
206 INTEGER :: gcd_max, ipe, jpe, npcol, nprow, nzero
207 INTEGER, EXTERNAL :: numroc
208
209 ncpu = max(1, min(ncpu_max, matrix_struct%para_env%num_pe))
210 ! Avoid layouts with odd number of CPUs (blacs grid layout will be square).
211 IF (ncpu > 2) ncpu = ncpu - modulo(ncpu, 2)
212
213 ! Keep a valid original layout only when no redistribution is needed.
214 IF (ncpu == matrix_struct%para_env%num_pe .AND. &
215 block_size == matrix_struct%nrow_block .AND. block_size == matrix_struct%ncol_block) THEN
216 nzero = count(matrix_struct%nrow_locals == 0) + count(matrix_struct%ncol_locals == 0)
217 IF (nzero == 0) RETURN
218 ! A smaller rank count ensures that an invalid original grid is replaced.
219 ncpu = ncpu - nzero
220 IF (ncpu > 2) ncpu = ncpu - modulo(ncpu, 2)
221 END IF
222
223 ! Fewer CPUs can mean more process columns (e.g. 32: 4x8, 20: 2x10).
224 ! Check the actual candidate grid with the final block size.
225 DO WHILE (ncpu > 1)
226 ! Determine layout of new blacs grid with ncpu CPUs
227 ! (snippet copied from cp_blacs_env.F:cp_blacs_env_create)
228 gcd_max = -1
229 DO ipe = 1, ceiling(sqrt(real(ncpu, dp)))
230 jpe = ncpu/ipe
231 IF (ipe*jpe /= ncpu) THEN
232 cycle
233 END IF
234 IF (gcd(ipe, jpe) >= gcd_max) THEN
235 nprow = ipe
236 npcol = jpe
237 gcd_max = gcd(ipe, jpe)
238 END IF
239 END DO
240
241 ! Count empty process rows and columns.
242 ! (snippet copied from cp_fm_struct.F:cp_fm_struct_create)
243 nzero = 0
244 DO ipe = 0, npcol - 1
245 IF (numroc(matrix_struct%ncol_global, block_size, ipe, 0, npcol) == 0) THEN
246 nzero = nzero + 1
247 END IF
248 END DO
249 DO ipe = 0, nprow - 1
250 IF (numroc(matrix_struct%nrow_global, block_size, ipe, 0, nprow) == 0) THEN
251 nzero = nzero + 1
252 END IF
253 END DO
254
255 IF (nzero == 0) THEN
256 EXIT
257 END IF
258
259 ncpu = ncpu - nzero
260
261 IF (ncpu > 2) THEN
262 ncpu = ncpu - modulo(ncpu, 2)
263 END IF
264 END DO
265
266 END FUNCTION cp_fm_max_ncpu_non_zero_column
267#endif
268
269! **************************************************************************************************
270!> \brief Determines the optimal number of CPUs for matrix diagonalization and redistributes
271!> the input matrices if necessary
272!> \param matrix the input cp_fm_type matrix to be diagonalized
273!> \param eigenvectors the cp_fm_type matrix that will hold the eigenvectors of the input matrix
274!> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
275!> or a pointer to the original matrix if no redistribution is required
276!> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
277!> matrix if no redistribution is required
278!> \param caller_is_elpa flag that determines if ELPA is used for diagonalization
279!> \param redist_info get info about the redistribution
280!> \par History
281!> - [01.2018] created by moving redistribution related code from cp_fm_syevd here
282!> \author Nico Holmberg [01.2018]
283! **************************************************************************************************
284 SUBROUTINE cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
285 caller_is_elpa, redist_info)
286
287 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
288 TYPE(cp_fm_type), INTENT(OUT) :: matrix_new, eigenvectors_new
289 LOGICAL, OPTIONAL, INTENT(IN) :: caller_is_elpa
290
291 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_redistribute_start'
292
293 INTEGER :: handle
294 LOGICAL :: is_elpa
295 TYPE(cp_fm_redistribute_info), OPTIONAL, INTENT(OUT) :: redist_info
296
297#if defined(__parallel)
298 REAL(kind=dp) :: fake_local_data(1, 1)
299 INTEGER :: fake_descriptor(9), mepos_old, &
300 io_unit, ngroups, ncol_block, blksize, nrow_block
301 TYPE(cp_fm_struct_type), POINTER :: fm_struct_new
302 TYPE(mp_para_env_type), POINTER :: para_env
303 TYPE(cp_logger_type), POINTER :: logger
304 TYPE(cp_fm_redistribute_info) :: rdinfo
305#endif
306
307 CALL timeset(routinen, handle)
308 is_elpa = .false.
309 IF (PRESENT(caller_is_elpa)) THEN
310#if defined(__ELPA)
311 is_elpa = caller_is_elpa
312#else
313 cpabort("CP2K compiled without the ELPA library.")
314#endif
315 END IF
316
317#if defined(__parallel)
318
319 logger => cp_get_default_logger()
320 io_unit = cp_logger_get_default_io_unit(logger)
321
322 ! first figure out the optimal number of cpus
323 ! this is pure heuristics, the defaults are based on rosa timings
324 ! that demonstrate that timings go up sharply if too many tasks are used
325 ! we take a multiple of 4, and approximately n/60
326 para_env => matrix%matrix_struct%para_env
327 mepos_old = para_env%mepos
328 ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
329 nrow_block = -1
330 blksize = -1
331
332 rdinfo%matrix_order = matrix%matrix_struct%nrow_global
333 rdinfo%num_pe_old = para_env%num_pe
334 rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
335 rdinfo%num_pe_new = rdinfo%num_pe_opt
336 rdinfo%num_pe_max_nz_col = -1
337 rdinfo%redistribute = .false.
338
339 IF (is_elpa) THEN
340 CALL cp_fm_get_info(matrix, ncol_block=ncol_block, nrow_block=nrow_block)
341
342 ! On GPUs, ELPA requires the block size to be a power of 2.
343 blksize = 1
344 DO WHILE (2*blksize <= min(nrow_block, ncol_block))
345 blksize = blksize*2
346 END DO
347 nrow_block = blksize
348 ncol_block = blksize
349
350 ! Only redistribute when necessary unless the user requests the heuristic.
351 rdinfo%num_pe_max_nz_col = cp_fm_max_ncpu_non_zero_column( &
352 matrix%matrix_struct, rdinfo%num_pe_old, blksize)
353 rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
354 IF (work_redistribute%elpa_force_redistribute .AND. rdinfo%num_pe_opt < rdinfo%num_pe_new) THEN
355 rdinfo%num_pe_new = cp_fm_max_ncpu_non_zero_column( &
356 matrix%matrix_struct, rdinfo%num_pe_opt, blksize)
357 END IF
358 END IF
359
360 ! finally, only redistribute if we're going to use less CPUs than before or changed the block size
361 rdinfo%redistribute = (rdinfo%num_pe_old > rdinfo%num_pe_new) .OR. (blksize >= 0 .AND. &
362 ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
363
364 IF (work_redistribute%should_print .AND. io_unit > 0) THEN
365 IF (is_elpa) THEN
366 IF (work_redistribute%elpa_force_redistribute) THEN
367 WRITE (unit=io_unit, fmt="(T2,A,T78,A3)") &
368 "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
369 ELSE
370 WRITE (unit=io_unit, fmt="(T2,A,T79,A2)") &
371 "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
372 END IF
373 END IF
374 CALL rdinfo%write(io_unit)
375 END IF
376 CALL para_env%sync()
377
378 ! if the optimal is smaller than num_pe, we will redistribute the input matrix
379 IF (rdinfo%redistribute) THEN
380 ! split comm, the first num_pe_new tasks will do the work
381 ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
382 ALLOCATE (work_redistribute%group_partition(0:1))
383 work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
384 ALLOCATE (work_redistribute%para_env_new)
385 CALL work_redistribute%para_env_new%from_split( &
386 comm=para_env, ngroups=ngroups, group_distribution=work_redistribute%group_distribution, &
387 n_subgroups=2, group_partition=work_redistribute%group_partition)
388
389 IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
390
391 ! create blacs, should inherit the preferences for the layout and so on, from the higher level
392 NULLIFY (work_redistribute%blacs_env_new)
393 CALL cp_blacs_env_create(blacs_env=work_redistribute%blacs_env_new, para_env=work_redistribute%para_env_new)
394
395 ! create new matrix
396 NULLIFY (fm_struct_new)
397 IF (nrow_block == -1 .OR. ncol_block == -1) THEN
398 CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
399 para_env=work_redistribute%para_env_new, &
400 context=work_redistribute%blacs_env_new, &
401 nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
402 ncol_block=ncol_block, nrow_block=nrow_block)
403 ELSE
404 CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
405 para_env=work_redistribute%para_env_new, &
406 context=work_redistribute%blacs_env_new, &
407 nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
408 ncol_block=ncol_block, nrow_block=nrow_block, force_block=.true.)
409 END IF
410 CALL cp_fm_create(matrix_new, matrix_struct=fm_struct_new, name="yevd_new_mat")
411 CALL cp_fm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="yevd_new_vec")
412 CALL cp_fm_struct_release(fm_struct_new)
413
414 ! redistribute old
415 CALL pdgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
416 matrix%matrix_struct%descriptor, &
417 matrix_new%local_data(1, 1), 1, 1, matrix_new%matrix_struct%descriptor, &
418 matrix%matrix_struct%context)
419 ELSE
420 ! these tasks must help redistribute (they own part of the data),
421 ! but need fake 'new' data, and their descriptor must indicate this with -1
422 ! see also scalapack comments on pdgemr2d
423 fake_descriptor = -1
424 CALL pdgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
425 matrix%matrix_struct%descriptor, &
426 fake_local_data(1, 1), 1, 1, fake_descriptor, &
427 matrix%matrix_struct%context)
428 END IF
429 ELSE
430 ! No need to redistribute, just return pointers to the original arrays
431 matrix_new = matrix
432 eigenvectors_new = eigenvectors
433 END IF
434
435 IF (PRESENT(redist_info)) THEN
436 redist_info = rdinfo
437 END IF
438#else
439
440 mark_used(matrix)
441 mark_used(eigenvectors)
442 mark_used(matrix_new)
443 mark_used(eigenvectors_new)
444 mark_used(redist_info)
445 cpabort("Routine called in non-parallel case.")
446#endif
447
448 CALL timestop(handle)
449
450 END SUBROUTINE cp_fm_redistribute_start
451
452! **************************************************************************************************
453!> \brief Redistributes eigenvectors and eigenvalues back to the original communicator group
454!> \param matrix the input cp_fm_type matrix to be diagonalized
455!> \param eigenvectors the cp_fm_type matrix that will hold the eigenvectors of the input matrix
456!> \param eig global array holding the eigenvalues of the input matrixmatrix
457!> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
458!> or a pointer to the original matrix if no redistribution is required
459!> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
460!> matrix if no redistribution is required
461!> \par History
462!> - [01.2018] created by moving redistribution related code from cp_fm_syevd here
463!> \author Nico Holmberg [01.2018]
464! **************************************************************************************************
465 SUBROUTINE cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
466
467 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
468 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eig
469 TYPE(cp_fm_type), INTENT(INOUT) :: matrix_new, eigenvectors_new
470
471 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_redistribute_end'
472
473 INTEGER :: handle
474#if defined(__parallel)
475 REAL(kind=dp) :: fake_local_data(1, 1)
476 INTEGER :: fake_descriptor(9), mepos_old, n
477 TYPE(mp_para_env_type), POINTER :: para_env
478#endif
479
480 CALL timeset(routinen, handle)
481
482#if defined(__parallel)
483
484 ! Check if matrix was redistributed
485 IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
486 n = matrix%matrix_struct%nrow_global
487 para_env => matrix%matrix_struct%para_env
488 mepos_old = para_env%mepos
489
490 IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
491 ! redistribute results on CPUs that hold the redistributed matrix
492 CALL pdgemr2d(n, n, eigenvectors_new%local_data(1, 1), 1, 1, eigenvectors_new%matrix_struct%descriptor, &
493 eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
494 eigenvectors%matrix_struct%context)
495 CALL cp_fm_release(matrix_new)
496 CALL cp_fm_release(eigenvectors_new)
497 ELSE
498 ! these tasks must help redistribute (they own part of the data),
499 ! but need fake 'new' data, and their descriptor must indicate this with -1
500 ! see also scalapack comments on pdgemr2d
501 fake_descriptor = -1
502 CALL pdgemr2d(n, n, fake_local_data(1, 1), 1, 1, fake_descriptor, &
503 eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
504 eigenvectors%matrix_struct%context)
505 END IF
506 ! free work
507 CALL cp_fm_redistribute_work_finalize(work_redistribute%group_distribution(mepos_old) == 0)
508
509 ! finally, also the eigenvalues need to end up on the non-group member tasks
510 CALL para_env%bcast(eig, 0)
511 END IF
512
513#else
514
515 mark_used(matrix)
516 mark_used(eigenvectors)
517 mark_used(eig)
518 mark_used(matrix_new)
519 mark_used(eigenvectors_new)
520 cpabort("Routine called in non-parallel case.")
521#endif
522
523 CALL timestop(handle)
524
525 END SUBROUTINE cp_fm_redistribute_end
526
527! **************************************************************************************************
528!> \brief Determines the optimal number of CPUs for matrix diagonalization and redistributes
529!> the input matrices if necessary
530!> \param matrix the input cp_cfm_type matrix to be diagonalized
531!> \param eigenvectors the cp_cfm_type matrix that will hold the eigenvectors of the input matrix
532!> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
533!> or a pointer to the original matrix if no redistribution is required
534!> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
535!> matrix if no redistribution is required
536!> \param caller_is_elpa flag that determines if ELPA is used for diagonalization
537!> \param redist_info get info about the redistribution
538!> \par History
539!> - [08.2026] created by mirroring cp_fm_redistribute_start for complex matrices
540! **************************************************************************************************
541 SUBROUTINE cp_cfm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
542 caller_is_elpa, redist_info)
543
544 TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
545 TYPE(cp_cfm_type), INTENT(OUT) :: matrix_new, eigenvectors_new
546 LOGICAL, OPTIONAL, INTENT(IN) :: caller_is_elpa
547
548 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_redistribute_start'
549
550 INTEGER :: handle
551 LOGICAL :: is_elpa
552 TYPE(cp_fm_redistribute_info), OPTIONAL, INTENT(OUT) :: redist_info
553
554#if defined(__parallel)
555 COMPLEX(KIND=dp) :: fake_local_data(1, 1)
556 INTEGER :: fake_descriptor(9), mepos_old, &
557 io_unit, ngroups, ncol_block, blksize, nrow_block
558 TYPE(cp_fm_struct_type), POINTER :: fm_struct_new
559 TYPE(mp_para_env_type), POINTER :: para_env
560 TYPE(cp_logger_type), POINTER :: logger
561 TYPE(cp_fm_redistribute_info) :: rdinfo
562#endif
563
564 CALL timeset(routinen, handle)
565 is_elpa = .false.
566 IF (PRESENT(caller_is_elpa)) THEN
567#if defined(__ELPA)
568 is_elpa = caller_is_elpa
569#else
570 cpabort("CP2K compiled without the ELPA library.")
571#endif
572 END IF
573
574#if defined(__parallel)
575
576 logger => cp_get_default_logger()
577 io_unit = cp_logger_get_default_io_unit(logger)
578
579 ! first figure out the optimal number of cpus
580 ! this is pure heuristics, the defaults are based on rosa timings
581 ! that demonstrate that timings go up sharply if too many tasks are used
582 ! we take a multiple of 4, and approximately n/60
583 para_env => matrix%matrix_struct%para_env
584 mepos_old = para_env%mepos
585 ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
586 nrow_block = -1
587 blksize = -1
588
589 rdinfo%matrix_order = matrix%matrix_struct%nrow_global
590 rdinfo%num_pe_old = para_env%num_pe
591 rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
592 rdinfo%num_pe_new = rdinfo%num_pe_opt
593 rdinfo%num_pe_max_nz_col = -1
594 rdinfo%redistribute = .false.
595
596 IF (is_elpa) THEN
597 CALL cp_cfm_get_info(matrix, ncol_block=ncol_block, nrow_block=nrow_block)
598
599 ! On GPUs, ELPA requires the block size to be a power of 2.
600 blksize = 1
601 DO WHILE (2*blksize <= min(nrow_block, ncol_block))
602 blksize = blksize*2
603 END DO
604 nrow_block = blksize
605 ncol_block = blksize
606
607 ! Only redistribute when necessary unless the user requests the heuristic.
608 rdinfo%num_pe_max_nz_col = cp_fm_max_ncpu_non_zero_column( &
609 matrix%matrix_struct, rdinfo%num_pe_old, blksize)
610 rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
611 IF (work_redistribute%elpa_force_redistribute .AND. rdinfo%num_pe_opt < rdinfo%num_pe_new) THEN
612 rdinfo%num_pe_new = cp_fm_max_ncpu_non_zero_column( &
613 matrix%matrix_struct, rdinfo%num_pe_opt, blksize)
614 END IF
615 END IF
616
617 ! finally, only redistribute if we're going to use less CPUs than before or changed the block size
618 rdinfo%redistribute = (rdinfo%num_pe_old > rdinfo%num_pe_new) .OR. (blksize >= 0 .AND. &
619 ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
620
621 IF (work_redistribute%should_print .AND. io_unit > 0) THEN
622 IF (is_elpa) THEN
623 IF (work_redistribute%elpa_force_redistribute) THEN
624 WRITE (unit=io_unit, fmt="(T2,A,T78,A3)") &
625 "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
626 ELSE
627 WRITE (unit=io_unit, fmt="(T2,A,T79,A2)") &
628 "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
629 END IF
630 END IF
631 CALL rdinfo%write(io_unit)
632 END IF
633 CALL para_env%sync()
634
635 ! if the optimal is smaller than num_pe, we will redistribute the input matrix
636 IF (rdinfo%redistribute) THEN
637 ! split comm, the first num_pe_new tasks will do the work
638 ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
639 ALLOCATE (work_redistribute%group_partition(0:1))
640 work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
641 ALLOCATE (work_redistribute%para_env_new)
642 CALL work_redistribute%para_env_new%from_split( &
643 comm=para_env, ngroups=ngroups, group_distribution=work_redistribute%group_distribution, &
644 n_subgroups=2, group_partition=work_redistribute%group_partition)
645
646 IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
647
648 ! create blacs, should inherit the preferences for the layout and so on, from the higher level
649 NULLIFY (work_redistribute%blacs_env_new)
650 CALL cp_blacs_env_create(blacs_env=work_redistribute%blacs_env_new, para_env=work_redistribute%para_env_new)
651
652 ! create new matrix
653 NULLIFY (fm_struct_new)
654 IF (nrow_block == -1 .OR. ncol_block == -1) THEN
655 CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
656 para_env=work_redistribute%para_env_new, &
657 context=work_redistribute%blacs_env_new, &
658 nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
659 ncol_block=ncol_block, nrow_block=nrow_block)
660 ELSE
661 CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
662 para_env=work_redistribute%para_env_new, &
663 context=work_redistribute%blacs_env_new, &
664 nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
665 ncol_block=ncol_block, nrow_block=nrow_block, force_block=.true.)
666 END IF
667 CALL cp_cfm_create(matrix_new, matrix_struct=fm_struct_new, name="zheevd_new_mat")
668 CALL cp_cfm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="zheevd_new_vec")
669 CALL cp_fm_struct_release(fm_struct_new)
670
671 ! redistribute old
672 CALL pzgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
673 matrix%matrix_struct%descriptor, &
674 matrix_new%local_data(1, 1), 1, 1, matrix_new%matrix_struct%descriptor, &
675 matrix%matrix_struct%context)
676 ELSE
677 ! these tasks must help redistribute (they own part of the data),
678 ! but need fake 'new' data, and their descriptor must indicate this with -1
679 ! see also scalapack comments on pzgemr2d
680 fake_descriptor = -1
681 CALL pzgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
682 matrix%matrix_struct%descriptor, &
683 fake_local_data(1, 1), 1, 1, fake_descriptor, &
684 matrix%matrix_struct%context)
685 END IF
686 ELSE
687 ! No need to redistribute, just return pointers to the original arrays
688 matrix_new = matrix
689 eigenvectors_new = eigenvectors
690 END IF
691
692 IF (PRESENT(redist_info)) THEN
693 redist_info = rdinfo
694 END IF
695#else
696
697 mark_used(matrix)
698 mark_used(eigenvectors)
699 mark_used(matrix_new)
700 mark_used(eigenvectors_new)
701 mark_used(redist_info)
702 cpabort("Routine called in non-parallel case.")
703#endif
704
705 CALL timestop(handle)
706
707 END SUBROUTINE cp_cfm_redistribute_start
708
709! **************************************************************************************************
710!> \brief Redistributes eigenvectors and eigenvalues back to the original communicator group
711!> \param matrix the input cp_cfm_type matrix to be diagonalized
712!> \param eigenvectors the cp_cfm_type matrix that will hold the eigenvectors of the input matrix
713!> \param eig global array holding the eigenvalues of the input matrixmatrix
714!> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
715!> or a pointer to the original matrix if no redistribution is required
716!> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
717!> matrix if no redistribution is required
718!> \par History
719!> - [08.2026] created by mirroring cp_fm_redistribute_end for complex matrices
720! **************************************************************************************************
721 SUBROUTINE cp_cfm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
722
723 TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
724 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eig
725 TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_new, eigenvectors_new
726
727 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_redistribute_end'
728
729 INTEGER :: handle
730#if defined(__parallel)
731 COMPLEX(KIND=dp) :: fake_local_data(1, 1)
732 INTEGER :: fake_descriptor(9), mepos_old, n
733 TYPE(mp_para_env_type), POINTER :: para_env
734#endif
735
736 CALL timeset(routinen, handle)
737
738#if defined(__parallel)
739
740 ! Check if matrix was redistributed
741 IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
742 n = matrix%matrix_struct%nrow_global
743 para_env => matrix%matrix_struct%para_env
744 mepos_old = para_env%mepos
745
746 IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
747 ! redistribute results on CPUs that hold the redistributed matrix
748 CALL pzgemr2d(n, n, eigenvectors_new%local_data(1, 1), 1, 1, eigenvectors_new%matrix_struct%descriptor, &
749 eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
750 eigenvectors%matrix_struct%context)
751 CALL cp_cfm_release(matrix_new)
752 CALL cp_cfm_release(eigenvectors_new)
753 ELSE
754 ! these tasks must help redistribute (they own part of the data),
755 ! but need fake 'new' data, and their descriptor must indicate this with -1
756 ! see also scalapack comments on pzgemr2d
757 fake_descriptor = -1
758 CALL pzgemr2d(n, n, fake_local_data(1, 1), 1, 1, fake_descriptor, &
759 eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
760 eigenvectors%matrix_struct%context)
761 END IF
762 ! free work
763 CALL cp_fm_redistribute_work_finalize(work_redistribute%group_distribution(mepos_old) == 0)
764
765 ! finally, also the eigenvalues need to end up on the non-group member tasks
766 CALL para_env%bcast(eig, 0)
767 END IF
768
769#else
770
771 mark_used(matrix)
772 mark_used(eigenvectors)
773 mark_used(eig)
774 mark_used(matrix_new)
775 mark_used(eigenvectors_new)
776 cpabort("Routine called in non-parallel case.")
777#endif
778
779 CALL timestop(handle)
780
781 END SUBROUTINE cp_cfm_redistribute_end
782
783END MODULE cp_fm_diag_utils
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_fm_diag_utils
Auxiliary tools to redistribute cp_fm_type and cp_cfm_type matrices before and after diagonalization....
Definition cp_fm_diag_utils.F:17
cp_fm_diag_utils::cp_fm_redistribute_end
subroutine, public cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
Redistributes eigenvectors and eigenvalues back to the original communicator group.
Definition cp_fm_diag_utils.F:466
cp_fm_diag_utils::cp_cfm_redistribute_start
subroutine, public cp_cfm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, caller_is_elpa, redist_info)
Determines the optimal number of CPUs for matrix diagonalization and redistributes the input matrices...
Definition cp_fm_diag_utils.F:543
cp_fm_diag_utils::cp_cfm_redistribute_end
subroutine, public cp_cfm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
Redistributes eigenvectors and eigenvalues back to the original communicator group.
Definition cp_fm_diag_utils.F:722
cp_fm_diag_utils::cp_fm_redistribute_start
subroutine, public cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, caller_is_elpa, redist_info)
Determines the optimal number of CPUs for matrix diagonalization and redistributes the input matrices...
Definition cp_fm_diag_utils.F:286
cp_fm_diag_utils::cp_fm_redistribute_init
subroutine, public cp_fm_redistribute_init(a, x, should_print, elpa_force_redistribute)
Initializes the parameters that determine how to calculate the optimal number of CPUs for diagonalizi...
Definition cp_fm_diag_utils.F:162
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1046
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1289
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_fm_diag_utils::cp_fm_redistribute_info
Definition cp_fm_diag_utils.F:47
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768