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 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 : ! **************************************************************************************************
17 : MODULE 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 :
86 : CONTAINS
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 16 : 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 16 : WRITE (UNIT=io_unit, FMT="(A)") ""
98 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
99 16 : "CP_FM_DIAG| Number of processes over which the matrix is distributed ", self%num_pe_old, &
100 32 : "CP_FM_DIAG| Matrix order ", self%matrix_order
101 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
102 16 : "CP_FM_DIAG| Optimal number of CPUs ", self%num_pe_opt
103 16 : IF (self%num_pe_max_nz_col < 0) THEN
104 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,A10)") &
105 16 : "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 0 : "CP_FM_DIAG| Maximum number of CPUs (with non-zero columns): ", self%num_pe_max_nz_col
109 : END IF
110 16 : IF (self%redistribute) THEN
111 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
112 0 : "CP_FM_DIAG| Number of processes for the redistribution ", self%num_pe_new
113 : ELSE
114 : WRITE (UNIT=io_unit, FMT="(T2,A)") &
115 16 : "CP_FM_DIAG| The matrix will NOT be redistributed"
116 : END IF
117 16 : WRITE (UNIT=io_unit, FMT="(A)") ""
118 :
119 16 : 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 10408 : SUBROUTINE cp_fm_redistribute_work_finalize(has_redistributed)
128 : LOGICAL, INTENT(IN) :: has_redistributed
129 :
130 10408 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
131 10408 : IF (has_redistributed) THEN
132 5204 : CALL cp_blacs_env_release(work_redistribute%blacs_env_new)
133 : END IF
134 10408 : CALL work_redistribute%para_env_new%free()
135 10408 : DEALLOCATE (work_redistribute%para_env_new)
136 10408 : DEALLOCATE (work_redistribute%group_distribution)
137 10408 : 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 10408 : NULLIFY (work_redistribute%group_distribution)
142 10408 : NULLIFY (work_redistribute%group_partition)
143 10408 : NULLIFY (work_redistribute%blacs_env_new)
144 10408 : NULLIFY (work_redistribute%para_env_new)
145 :
146 10408 : 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 11722 : 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 11722 : 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 11722 : work_redistribute = cp_fm_redistribute_type()
170 11722 : work_redistribute%a = a
171 11722 : work_redistribute%x = x
172 11722 : work_redistribute%should_print = should_print
173 11722 : work_redistribute%elpa_force_redistribute = elpa_force_redistribute
174 :
175 11722 : 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 286096 : 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 286096 : (work_redistribute%a*work_redistribute%x))*work_redistribute%a
189 :
190 286096 : 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 14614 : 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 14614 : 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 14614 : 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 14614 : block_size == matrix_struct%nrow_block .AND. block_size == matrix_struct%ncol_block) THEN
216 16788 : nzero = COUNT(matrix_struct%nrow_locals == 0) + COUNT(matrix_struct%ncol_locals == 0)
217 4188 : IF (nzero == 0) RETURN
218 : ! A smaller rank count ensures that an invalid original grid is replaced.
219 4 : ncpu = ncpu - nzero
220 4 : 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 10444 : 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 10442 : gcd_max = -1
229 31402 : DO ipe = 1, CEILING(SQRT(REAL(ncpu, dp)))
230 20960 : jpe = ncpu/ipe
231 20960 : IF (ipe*jpe /= ncpu) THEN
232 : CYCLE
233 : END IF
234 31364 : IF (gcd(ipe, jpe) >= gcd_max) THEN
235 20892 : nprow = ipe
236 20892 : npcol = jpe
237 20892 : 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 10442 : nzero = 0
244 21076 : DO ipe = 0, npcol - 1
245 21076 : IF (numroc(matrix_struct%ncol_global, block_size, ipe, 0, npcol) == 0) THEN
246 40 : nzero = nzero + 1
247 : END IF
248 : END DO
249 31336 : DO ipe = 0, nprow - 1
250 31336 : IF (numroc(matrix_struct%nrow_global, block_size, ipe, 0, nprow) == 0) THEN
251 8 : nzero = nzero + 1
252 : END IF
253 : END DO
254 :
255 10442 : IF (nzero == 0) THEN
256 : EXIT
257 : END IF
258 :
259 14 : ncpu = ncpu - nzero
260 :
261 10444 : IF (ncpu > 2) THEN
262 10 : 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 14516 : 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 286024 : CALL timeset(routineN, handle)
308 286024 : is_elpa = .FALSE.
309 286024 : IF (PRESENT(caller_is_elpa)) THEN
310 : #if defined(__ELPA)
311 14516 : 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 286024 : logger => cp_get_default_logger()
320 286024 : 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 286024 : para_env => matrix%matrix_struct%para_env
327 286024 : mepos_old = para_env%mepos
328 286024 : ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
329 286024 : nrow_block = -1
330 286024 : blksize = -1
331 :
332 286024 : rdinfo%matrix_order = matrix%matrix_struct%nrow_global
333 286024 : rdinfo%num_pe_old = para_env%num_pe
334 286024 : rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
335 286024 : rdinfo%num_pe_new = rdinfo%num_pe_opt
336 286024 : rdinfo%num_pe_max_nz_col = -1
337 286024 : rdinfo%redistribute = .FALSE.
338 :
339 286024 : IF (is_elpa) THEN
340 14516 : 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 14516 : blksize = 1
344 89296 : DO WHILE (2*blksize <= MIN(nrow_block, ncol_block))
345 74780 : blksize = blksize*2
346 : END DO
347 14516 : nrow_block = blksize
348 14516 : 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 14516 : matrix%matrix_struct, rdinfo%num_pe_old, blksize)
353 14516 : rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
354 14516 : 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 0 : 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 286024 : ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
363 :
364 286024 : IF (work_redistribute%should_print .AND. io_unit > 0) THEN
365 16 : IF (is_elpa) THEN
366 0 : IF (work_redistribute%elpa_force_redistribute) THEN
367 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
368 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
369 : ELSE
370 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
371 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
372 : END IF
373 : END IF
374 16 : CALL rdinfo%write(io_unit)
375 : END IF
376 286024 : CALL para_env%sync()
377 :
378 : ! if the optimal is smaller than num_pe, we will redistribute the input matrix
379 286024 : IF (rdinfo%redistribute) THEN
380 : ! split comm, the first num_pe_new tasks will do the work
381 31008 : ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
382 10336 : ALLOCATE (work_redistribute%group_partition(0:1))
383 31008 : work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
384 10336 : 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 10336 : n_subgroups=2, group_partition=work_redistribute%group_partition)
388 :
389 10336 : 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 5168 : NULLIFY (work_redistribute%blacs_env_new)
393 5168 : 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 5168 : NULLIFY (fm_struct_new)
397 5168 : 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 0 : 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 5168 : ncol_block=ncol_block, nrow_block=nrow_block, force_block=.TRUE.)
409 : END IF
410 5168 : CALL cp_fm_create(matrix_new, matrix_struct=fm_struct_new, name="yevd_new_mat")
411 5168 : CALL cp_fm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="yevd_new_vec")
412 5168 : 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 5168 : 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 51680 : 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 5168 : matrix%matrix_struct%context)
428 : END IF
429 : ELSE
430 : ! No need to redistribute, just return pointers to the original arrays
431 275688 : matrix_new = matrix
432 275688 : eigenvectors_new = eigenvectors
433 : END IF
434 :
435 286024 : IF (PRESENT(redist_info)) THEN
436 14516 : 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 286024 : CALL timestop(handle)
449 :
450 286024 : 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 286024 : 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 286024 : CALL timeset(routineN, handle)
481 :
482 : #if defined(__parallel)
483 :
484 : ! Check if matrix was redistributed
485 286024 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
486 10336 : n = matrix%matrix_struct%nrow_global
487 10336 : para_env => matrix%matrix_struct%para_env
488 10336 : mepos_old = para_env%mepos
489 :
490 10336 : 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 5168 : eigenvectors%matrix_struct%context)
495 5168 : CALL cp_fm_release(matrix_new)
496 5168 : 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 51680 : 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 5168 : eigenvectors%matrix_struct%context)
505 : END IF
506 : ! free work
507 10336 : 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 1319296 : 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 286024 : CALL timestop(handle)
524 :
525 286024 : 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 72 : 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 72 : CALL timeset(routineN, handle)
565 72 : is_elpa = .FALSE.
566 72 : IF (PRESENT(caller_is_elpa)) THEN
567 : #if defined(__ELPA)
568 72 : 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 72 : logger => cp_get_default_logger()
577 72 : 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 72 : para_env => matrix%matrix_struct%para_env
584 72 : mepos_old = para_env%mepos
585 72 : ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
586 72 : nrow_block = -1
587 72 : blksize = -1
588 :
589 72 : rdinfo%matrix_order = matrix%matrix_struct%nrow_global
590 72 : rdinfo%num_pe_old = para_env%num_pe
591 72 : rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
592 72 : rdinfo%num_pe_new = rdinfo%num_pe_opt
593 72 : rdinfo%num_pe_max_nz_col = -1
594 72 : rdinfo%redistribute = .FALSE.
595 :
596 72 : IF (is_elpa) THEN
597 72 : 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 72 : blksize = 1
601 360 : DO WHILE (2*blksize <= MIN(nrow_block, ncol_block))
602 288 : blksize = blksize*2
603 : END DO
604 72 : nrow_block = blksize
605 72 : 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 72 : matrix%matrix_struct, rdinfo%num_pe_old, blksize)
610 72 : rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
611 72 : 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 0 : 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 72 : ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
620 :
621 72 : IF (work_redistribute%should_print .AND. io_unit > 0) THEN
622 0 : IF (is_elpa) THEN
623 0 : IF (work_redistribute%elpa_force_redistribute) THEN
624 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
625 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
626 : ELSE
627 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
628 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
629 : END IF
630 : END IF
631 0 : CALL rdinfo%write(io_unit)
632 : END IF
633 72 : CALL para_env%sync()
634 :
635 : ! if the optimal is smaller than num_pe, we will redistribute the input matrix
636 72 : IF (rdinfo%redistribute) THEN
637 : ! split comm, the first num_pe_new tasks will do the work
638 216 : ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
639 72 : ALLOCATE (work_redistribute%group_partition(0:1))
640 216 : work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
641 72 : 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 72 : n_subgroups=2, group_partition=work_redistribute%group_partition)
645 :
646 72 : 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 36 : NULLIFY (work_redistribute%blacs_env_new)
650 36 : 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 36 : NULLIFY (fm_struct_new)
654 36 : 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 0 : 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 36 : ncol_block=ncol_block, nrow_block=nrow_block, force_block=.TRUE.)
666 : END IF
667 36 : CALL cp_cfm_create(matrix_new, matrix_struct=fm_struct_new, name="zheevd_new_mat")
668 36 : CALL cp_cfm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="zheevd_new_vec")
669 36 : 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 36 : 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 360 : 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 36 : matrix%matrix_struct%context)
685 : END IF
686 : ELSE
687 : ! No need to redistribute, just return pointers to the original arrays
688 0 : matrix_new = matrix
689 0 : eigenvectors_new = eigenvectors
690 : END IF
691 :
692 72 : IF (PRESENT(redist_info)) THEN
693 72 : 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 72 : CALL timestop(handle)
706 :
707 72 : 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 72 : 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 72 : CALL timeset(routineN, handle)
737 :
738 : #if defined(__parallel)
739 :
740 : ! Check if matrix was redistributed
741 72 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
742 72 : n = matrix%matrix_struct%nrow_global
743 72 : para_env => matrix%matrix_struct%para_env
744 72 : mepos_old = para_env%mepos
745 :
746 72 : 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 36 : eigenvectors%matrix_struct%context)
751 36 : CALL cp_cfm_release(matrix_new)
752 36 : 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 360 : 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 36 : eigenvectors%matrix_struct%context)
761 : END IF
762 : ! free work
763 72 : 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 5688 : 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 72 : CALL timestop(handle)
780 :
781 72 : END SUBROUTINE cp_cfm_redistribute_end
782 :
783 0 : END MODULE cp_fm_diag_utils
|
