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_get,&
27 : cp_fm_struct_release,&
28 : cp_fm_struct_type
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_get_info,&
31 : cp_fm_release,&
32 : cp_fm_type
33 : USE cp_log_handling, ONLY: cp_get_default_logger,&
34 : cp_logger_get_default_io_unit,&
35 : cp_logger_type
36 : USE kinds, ONLY: dp
37 : USE mathlib, ONLY: gcd
38 : USE message_passing, ONLY: mp_para_env_type
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_diag_utils'
46 :
47 : ! Information on redistribution
48 : TYPE, PUBLIC :: cp_fm_redistribute_info
49 : INTEGER :: matrix_order = -1
50 : INTEGER :: num_pe_old = -1 ! number of processes before a potential redistribute
51 : INTEGER :: num_pe_new = -1 ! number of processes after a potential redistribute
52 : INTEGER :: num_pe_opt = -1 ! optimal number of processes for the given matrix
53 : INTEGER :: num_pe_max_nz_col = -1 ! the maximal number of processes s.t. no column has zero width, may be < 0 if ignored
54 : LOGICAL :: redistribute = .FALSE. ! whether or not the matrix was actually redistributed
55 : CONTAINS
56 : PROCEDURE, PASS(self) :: write => cp_fm_redistribute_info_write
57 : END TYPE cp_fm_redistribute_info
58 :
59 : ! Container for redistribution settings and temporary work structs
60 : TYPE cp_fm_redistribute_type
61 : ! Settings
62 : INTEGER :: a = -1, x = -1
63 : LOGICAL :: should_print = .FALSE.
64 : LOGICAL :: elpa_force_redistribute = .FALSE.
65 : ! Temporaries
66 : INTEGER, DIMENSION(:), POINTER :: group_distribution => NULL(), &
67 : group_partition => NULL()
68 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_new => NULL()
69 : TYPE(mp_para_env_type), POINTER :: para_env_new => NULL()
70 : END TYPE cp_fm_redistribute_type
71 :
72 : ! Permanent instance of the redistribute type
73 : TYPE(cp_fm_redistribute_type), PRIVATE, &
74 : SAVE :: work_redistribute
75 :
76 : ! Public subroutines
77 :
78 : PUBLIC :: cp_fm_redistribute_start, &
79 : cp_fm_redistribute_end, &
80 : cp_fm_redistribute_init, &
81 : cp_cfm_redistribute_start, &
82 : cp_cfm_redistribute_end
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Write the redistribute info nicely formatted to the given I/O unit
88 : !> \param self reference to the cp_fm_redistribute_info instance
89 : !> \param io_unit I/O unit to use for writing
90 : ! **************************************************************************************************
91 0 : SUBROUTINE cp_fm_redistribute_info_write(self, io_unit)
92 : CLASS(cp_fm_redistribute_info), INTENT(IN) :: self
93 : INTEGER, INTENT(IN) :: io_unit
94 :
95 0 : WRITE (UNIT=io_unit, FMT="(A)") ""
96 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
97 0 : "CP_FM_DIAG| Number of processes over which the matrix is distributed ", self%num_pe_old, &
98 0 : "CP_FM_DIAG| Matrix order ", self%matrix_order
99 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
100 0 : "CP_FM_DIAG| Optimal number of CPUs ", self%num_pe_opt
101 0 : IF (self%num_pe_max_nz_col < 0) THEN
102 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,A10)") &
103 0 : "CP_FM_DIAG| Maximum number of CPUs (with non-zero columns) ", "<N/A>"
104 : ELSE
105 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
106 0 : "CP_FM_DIAG| Maximum number of CPUs (with non-zero columns): ", self%num_pe_max_nz_col
107 : END IF
108 0 : IF (self%redistribute) THEN
109 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
110 0 : "CP_FM_DIAG| Number of processes for the redistribution ", self%num_pe_new
111 : ELSE
112 : WRITE (UNIT=io_unit, FMT="(T2,A)") &
113 0 : "CP_FM_DIAG| The matrix will NOT be redistributed"
114 : END IF
115 0 : WRITE (UNIT=io_unit, FMT="(A)") ""
116 :
117 0 : END SUBROUTINE cp_fm_redistribute_info_write
118 :
119 : ! **************************************************************************************************
120 : !> \brief Releases the temporary storage needed when redistributing arrays
121 : !> \param has_redistributed flag that determines if the processors holds a part of the
122 : !> redistributed array
123 : !> \author Nico Holmberg [01.2018]
124 : ! **************************************************************************************************
125 275518 : SUBROUTINE cp_fm_redistribute_work_finalize(has_redistributed)
126 : LOGICAL, INTENT(IN) :: has_redistributed
127 :
128 275518 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
129 275518 : IF (has_redistributed) THEN
130 139645 : CALL cp_blacs_env_release(work_redistribute%blacs_env_new)
131 : END IF
132 275518 : CALL work_redistribute%para_env_new%free()
133 275518 : DEALLOCATE (work_redistribute%para_env_new)
134 275518 : DEALLOCATE (work_redistribute%group_distribution)
135 275518 : DEALLOCATE (work_redistribute%group_partition)
136 : END IF
137 : ! Return work to its initial state
138 275518 : work_redistribute = cp_fm_redistribute_type()
139 :
140 275518 : END SUBROUTINE cp_fm_redistribute_work_finalize
141 :
142 : ! **************************************************************************************************
143 : !> \brief Initializes the parameters that determine how to calculate the optimal number of CPUs
144 : !> for diagonalizing a matrix. The parameters are read from the GLOBAL input section.
145 : !> \param a integer parameter used to define the rule for determining the optimal
146 : !> number of CPUs for diagonalization
147 : !> \param x integer parameter used to define the rule for determining the optimal
148 : !> number of CPUs for diagonalization
149 : !> \param should_print flag that determines if information about the redistribution process
150 : !> should be printed
151 : !> \param elpa_force_redistribute flag that if redistribution should always be performed when
152 : !> the ELPA diagonalization library is in use
153 : !> \author Nico Holmberg [01.2018]
154 : ! **************************************************************************************************
155 11451 : SUBROUTINE cp_fm_redistribute_init(a, x, should_print, elpa_force_redistribute)
156 : INTEGER, INTENT(IN) :: a, x
157 : LOGICAL, INTENT(IN) :: should_print, elpa_force_redistribute
158 :
159 : work_redistribute%a = a
160 : work_redistribute%x = x
161 : work_redistribute%should_print = should_print
162 : work_redistribute%elpa_force_redistribute = elpa_force_redistribute
163 : ! Init work
164 11451 : work_redistribute = cp_fm_redistribute_type()
165 :
166 11451 : END SUBROUTINE cp_fm_redistribute_init
167 :
168 : ! **************************************************************************************************
169 : !> \brief Calculates the optimal number of CPUs for diagonalizing a matrix.
170 : !> \param size the size of the diagonalized matrix
171 : !> \return the optimal number of CPUs
172 : !> \author Nico Holmberg [01.2018]
173 : ! **************************************************************************************************
174 279668 : PURE FUNCTION cp_fm_diag_get_optimal_ncpu(size) RESULT(ncpu)
175 : INTEGER, INTENT(IN) :: size
176 : INTEGER :: ncpu
177 :
178 : ncpu = ((size + work_redistribute%a*work_redistribute%x - 1)/ &
179 279668 : (work_redistribute%a*work_redistribute%x))*work_redistribute%a
180 :
181 279668 : END FUNCTION cp_fm_diag_get_optimal_ncpu
182 :
183 : #if defined(__parallel)
184 : ! **************************************************************************************************
185 : !> \brief Determines the largest number of CPUs a matrix can be distributed on without any of the
186 : !> processors getting a zero-width column (currently only needed for ELPA).
187 : !> \param matrix the matrix that will be diagonalized
188 : !> \return the maximum number of CPUs for ELPA
189 : !> \author Nico Holmberg [01.2018]
190 : ! **************************************************************************************************
191 14422 : FUNCTION cp_fm_max_ncpu_non_zero_column(matrix) RESULT(ncpu)
192 : TYPE(cp_fm_type), INTENT(IN) :: matrix
193 : INTEGER :: ncpu
194 :
195 : INTEGER :: gcd_max, ipe, jpe, ncol_block, &
196 : ncol_global, npcol, nrow_block, &
197 : nrow_global, num_pe_old, nzero
198 14422 : INTEGER, DIMENSION(:), POINTER :: ncol_locals
199 : INTEGER, EXTERNAL :: numroc
200 :
201 14422 : NULLIFY (ncol_locals)
202 : ! First check if there are any zero width columns in current layout
203 : CALL cp_fm_get_info(matrix, ncol_locals=ncol_locals, &
204 : nrow_global=nrow_global, ncol_global=ncol_global, &
205 14422 : nrow_block=nrow_block, ncol_block=ncol_block)
206 28850 : nzero = COUNT(ncol_locals == 0)
207 14422 : num_pe_old = matrix%matrix_struct%para_env%num_pe
208 14422 : ncpu = num_pe_old - nzero
209 :
210 : ! Avoid layouts with odd number of CPUs (blacs grid layout will be square)
211 14422 : IF (ncpu > 2) THEN
212 0 : ncpu = ncpu - MODULO(ncpu, 2)
213 : END IF
214 :
215 : ! if there are no zero-width columns and the number of processors was even, leave it at that
216 14422 : IF (ncpu == num_pe_old) THEN
217 : RETURN
218 : END IF
219 :
220 : ! Iteratively search for the maximum number of CPUs for ELPA
221 : ! On each step, we test whether the blacs grid created with ncpu processes
222 : ! contains any columns with zero width
223 0 : DO WHILE (ncpu > 1)
224 : ! Determine layout of new blacs grid with ncpu CPUs
225 : ! (snippet copied from cp_blacs_env.F:cp_blacs_env_create)
226 0 : gcd_max = -1
227 0 : DO ipe = 1, CEILING(SQRT(REAL(ncpu, dp)))
228 0 : jpe = ncpu/ipe
229 0 : IF (ipe*jpe /= ncpu) THEN
230 : CYCLE
231 : END IF
232 0 : IF (gcd(ipe, jpe) >= gcd_max) THEN
233 0 : npcol = jpe
234 0 : gcd_max = gcd(ipe, jpe)
235 : END IF
236 : END DO
237 :
238 : ! Count the number of processors without any columns
239 : ! (snippet copied from cp_fm_struct.F:cp_fm_struct_create)
240 0 : nzero = 0
241 0 : DO ipe = 0, npcol - 1
242 0 : IF (numroc(ncol_global, ncol_block, ipe, 0, npcol) == 0) THEN
243 0 : nzero = nzero + 1
244 : END IF
245 : END DO
246 :
247 0 : IF (nzero == 0) THEN
248 : EXIT
249 : END IF
250 :
251 0 : ncpu = ncpu - nzero
252 :
253 0 : IF (ncpu > 2) THEN
254 0 : ncpu = ncpu - MODULO(ncpu, 2)
255 : END IF
256 : END DO
257 :
258 14422 : END FUNCTION cp_fm_max_ncpu_non_zero_column
259 : #endif
260 :
261 : ! **************************************************************************************************
262 : !> \brief Determines the optimal number of CPUs for matrix diagonalization and redistributes
263 : !> the input matrices if necessary
264 : !> \param matrix the input cp_fm_type matrix to be diagonalized
265 : !> \param eigenvectors the cp_fm_type matrix that will hold the eigenvectors of the input matrix
266 : !> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
267 : !> or a pointer to the original matrix if no redistribution is required
268 : !> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
269 : !> matrix if no redistribution is required
270 : !> \param caller_is_elpa flag that determines if ELPA is used for diagonalization
271 : !> \param redist_info get info about the redistribution
272 : !> \par History
273 : !> - [01.2018] created by moving redistribution related code from cp_fm_syevd here
274 : !> \author Nico Holmberg [01.2018]
275 : ! **************************************************************************************************
276 14422 : SUBROUTINE cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
277 : caller_is_elpa, redist_info)
278 :
279 : TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
280 : TYPE(cp_fm_type), INTENT(OUT) :: matrix_new, eigenvectors_new
281 : LOGICAL, OPTIONAL, INTENT(IN) :: caller_is_elpa
282 :
283 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_redistribute_start'
284 :
285 : INTEGER :: handle
286 : LOGICAL :: is_elpa
287 : TYPE(cp_fm_redistribute_info), OPTIONAL, INTENT(OUT) :: redist_info
288 :
289 : #if defined(__parallel)
290 : REAL(KIND=dp) :: fake_local_data(1, 1)
291 : INTEGER :: fake_descriptor(9), mepos_old, &
292 : io_unit, ngroups, ncol_block, blksize, nrow_block
293 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_new
294 : TYPE(mp_para_env_type), POINTER :: para_env
295 : TYPE(cp_logger_type), POINTER :: logger
296 : TYPE(cp_fm_redistribute_info) :: rdinfo
297 : #endif
298 :
299 279596 : CALL timeset(routineN, handle)
300 279596 : is_elpa = .FALSE.
301 279596 : IF (PRESENT(caller_is_elpa)) THEN
302 : #if defined(__ELPA)
303 14422 : is_elpa = caller_is_elpa
304 : #else
305 : CPABORT("CP2K compiled without the ELPA library.")
306 : #endif
307 : END IF
308 :
309 : #if defined(__parallel)
310 :
311 279596 : logger => cp_get_default_logger()
312 279596 : io_unit = cp_logger_get_default_io_unit(logger)
313 :
314 : ! first figure out the optimal number of cpus
315 : ! this is pure heuristics, the defaults are based on rosa timings
316 : ! that demonstrate that timings go up sharply if too many tasks are used
317 : ! we take a multiple of 4, and approximately n/60
318 279596 : para_env => matrix%matrix_struct%para_env
319 279596 : mepos_old = para_env%mepos
320 279596 : ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
321 279596 : nrow_block = -1
322 279596 : blksize = -1
323 :
324 279596 : rdinfo%matrix_order = matrix%matrix_struct%nrow_global
325 279596 : rdinfo%num_pe_old = para_env%num_pe
326 279596 : rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
327 279596 : rdinfo%num_pe_new = rdinfo%num_pe_opt
328 : rdinfo%num_pe_max_nz_col = -1
329 : rdinfo%redistribute = .FALSE.
330 :
331 279596 : IF (is_elpa) THEN
332 : ! with ELPA we don't have to redistribute if not necessary (scales, unlike ScaLAPACK)
333 14422 : rdinfo%num_pe_new = rdinfo%num_pe_old
334 :
335 : ! BUT: Diagonalization with ELPA fails when a processor column has zero width
336 : ! Determine the maximum number of CPUs the matrix can be distributed without zero-width columns
337 : ! for the current block size.
338 14422 : rdinfo%num_pe_max_nz_col = cp_fm_max_ncpu_non_zero_column(matrix)
339 :
340 : ! if the user wants to redistribute to the ScaLAPACK optimal number of CPUs anyway, let him if it's safe.
341 14422 : IF (work_redistribute%elpa_force_redistribute .AND. rdinfo%num_pe_opt < rdinfo%num_pe_max_nz_col) THEN
342 : ! Use heuristics to determine the need for redistribution (when num_pe_opt is smaller than the safe maximum)
343 : ! in this case we can also take the block size used for ScaLAPACK
344 0 : rdinfo%num_pe_new = rdinfo%num_pe_opt
345 14422 : ELSE IF (rdinfo%num_pe_old > rdinfo%num_pe_max_nz_col) THEN
346 : ! Otherwise, only redistribute if we have to
347 0 : rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
348 : ! do NOT let cp_fm_struct_create automatically adjust the block size because the
349 : ! calculated number of processors such that no block has 0 columns wouldn't match (see #578):
350 : ! if the automatically chosen block size is larger than the present one we would still end
351 : ! up with empty processors
352 : END IF
353 :
354 14422 : CALL cp_fm_get_info(matrix, ncol_block=ncol_block, nrow_block=nrow_block)
355 :
356 : ! On GPUs, ELPA requires the block size to be a power of 2
357 14422 : blksize = 1
358 88704 : DO WHILE (2*blksize <= MIN(nrow_block, ncol_block))
359 14422 : blksize = blksize*2
360 : END DO
361 14422 : nrow_block = blksize
362 14422 : ncol_block = blksize
363 : END IF
364 :
365 : ! finally, only redistribute if we're going to use less CPUs than before or changed the block size
366 : rdinfo%redistribute = (rdinfo%num_pe_old > rdinfo%num_pe_new) .OR. (blksize >= 0 .AND. &
367 279596 : ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
368 :
369 279596 : IF (work_redistribute%should_print .AND. io_unit > 0) THEN
370 0 : IF (is_elpa) THEN
371 0 : IF (work_redistribute%elpa_force_redistribute) THEN
372 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
373 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
374 : ELSE
375 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
376 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
377 : END IF
378 : END IF
379 0 : CALL rdinfo%write(io_unit)
380 : END IF
381 279596 : CALL para_env%sync()
382 :
383 : ! if the optimal is smaller than num_pe, we will redistribute the input matrix
384 279596 : IF (rdinfo%redistribute) THEN
385 : ! split comm, the first num_pe_new tasks will do the work
386 826338 : ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
387 275446 : ALLOCATE (work_redistribute%group_partition(0:1))
388 826338 : work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
389 275446 : ALLOCATE (work_redistribute%para_env_new)
390 : CALL work_redistribute%para_env_new%from_split( &
391 : comm=para_env, ngroups=ngroups, group_distribution=work_redistribute%group_distribution, &
392 275446 : n_subgroups=2, group_partition=work_redistribute%group_partition)
393 :
394 275446 : IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
395 :
396 : ! create blacs, should inherit the preferences for the layout and so on, from the higher level
397 139609 : NULLIFY (work_redistribute%blacs_env_new)
398 139609 : CALL cp_blacs_env_create(blacs_env=work_redistribute%blacs_env_new, para_env=work_redistribute%para_env_new)
399 :
400 : ! create new matrix
401 139609 : NULLIFY (fm_struct_new)
402 139609 : IF (nrow_block == -1 .OR. ncol_block == -1) THEN
403 : CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
404 : para_env=work_redistribute%para_env_new, &
405 : context=work_redistribute%blacs_env_new, &
406 : nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
407 134473 : ncol_block=ncol_block, nrow_block=nrow_block)
408 : ELSE
409 : CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
410 : para_env=work_redistribute%para_env_new, &
411 : context=work_redistribute%blacs_env_new, &
412 : nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
413 5136 : ncol_block=ncol_block, nrow_block=nrow_block, force_block=.TRUE.)
414 : END IF
415 139609 : CALL cp_fm_create(matrix_new, matrix_struct=fm_struct_new, name="yevd_new_mat")
416 139609 : CALL cp_fm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="yevd_new_vec")
417 139609 : CALL cp_fm_struct_release(fm_struct_new)
418 :
419 : ! redistribute old
420 : CALL pdgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
421 : matrix%matrix_struct%descriptor, &
422 : matrix_new%local_data(1, 1), 1, 1, matrix_new%matrix_struct%descriptor, &
423 139609 : matrix%matrix_struct%context)
424 : ELSE
425 : ! these tasks must help redistribute (they own part of the data),
426 : ! but need fake 'new' data, and their descriptor must indicate this with -1
427 : ! see also scalapack comments on pdgemr2d
428 1358370 : fake_descriptor = -1
429 : CALL pdgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
430 : matrix%matrix_struct%descriptor, &
431 : fake_local_data(1, 1), 1, 1, fake_descriptor, &
432 135837 : matrix%matrix_struct%context)
433 : END IF
434 : ELSE
435 : ! No need to redistribute, just return pointers to the original arrays
436 4150 : matrix_new = matrix
437 4150 : eigenvectors_new = eigenvectors
438 : END IF
439 :
440 279596 : IF (PRESENT(redist_info)) THEN
441 14422 : redist_info = rdinfo
442 : END IF
443 : #else
444 :
445 : MARK_USED(matrix)
446 : MARK_USED(eigenvectors)
447 : MARK_USED(matrix_new)
448 : MARK_USED(eigenvectors_new)
449 : MARK_USED(redist_info)
450 : CPABORT("Routine called in non-parallel case.")
451 : #endif
452 :
453 279596 : CALL timestop(handle)
454 :
455 279596 : END SUBROUTINE cp_fm_redistribute_start
456 :
457 : ! **************************************************************************************************
458 : !> \brief Redistributes eigenvectors and eigenvalues back to the original communicator group
459 : !> \param matrix the input cp_fm_type matrix to be diagonalized
460 : !> \param eigenvectors the cp_fm_type matrix that will hold the eigenvectors of the input matrix
461 : !> \param eig global array holding the eigenvalues of the input matrixmatrix
462 : !> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
463 : !> or a pointer to the original matrix if no redistribution is required
464 : !> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
465 : !> matrix if no redistribution is required
466 : !> \par History
467 : !> - [01.2018] created by moving redistribution related code from cp_fm_syevd here
468 : !> \author Nico Holmberg [01.2018]
469 : ! **************************************************************************************************
470 279596 : SUBROUTINE cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
471 :
472 : TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
473 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eig
474 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_new, eigenvectors_new
475 :
476 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_redistribute_end'
477 :
478 : INTEGER :: handle
479 : #if defined(__parallel)
480 : REAL(KIND=dp) :: fake_local_data(1, 1)
481 : INTEGER :: fake_descriptor(9), mepos_old, n
482 : TYPE(mp_para_env_type), POINTER :: para_env
483 : #endif
484 :
485 279596 : CALL timeset(routineN, handle)
486 :
487 : #if defined(__parallel)
488 :
489 : ! Check if matrix was redistributed
490 279596 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
491 275446 : n = matrix%matrix_struct%nrow_global
492 275446 : para_env => matrix%matrix_struct%para_env
493 275446 : mepos_old = para_env%mepos
494 :
495 275446 : IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
496 : ! redistribute results on CPUs that hold the redistributed matrix
497 : CALL pdgemr2d(n, n, eigenvectors_new%local_data(1, 1), 1, 1, eigenvectors_new%matrix_struct%descriptor, &
498 : eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
499 139609 : eigenvectors%matrix_struct%context)
500 139609 : CALL cp_fm_release(matrix_new)
501 139609 : CALL cp_fm_release(eigenvectors_new)
502 : ELSE
503 : ! these tasks must help redistribute (they own part of the data),
504 : ! but need fake 'new' data, and their descriptor must indicate this with -1
505 : ! see also scalapack comments on pdgemr2d
506 1358370 : fake_descriptor = -1
507 : CALL pdgemr2d(n, n, fake_local_data(1, 1), 1, 1, fake_descriptor, &
508 : eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
509 135837 : eigenvectors%matrix_struct%context)
510 : END IF
511 : ! free work
512 275446 : CALL cp_fm_redistribute_work_finalize(work_redistribute%group_distribution(mepos_old) == 0)
513 :
514 : ! finally, also the eigenvalues need to end up on the non-group member tasks
515 7662682 : CALL para_env%bcast(eig, 0)
516 : END IF
517 :
518 : #else
519 :
520 : MARK_USED(matrix)
521 : MARK_USED(eigenvectors)
522 : MARK_USED(eig)
523 : MARK_USED(matrix_new)
524 : MARK_USED(eigenvectors_new)
525 : CPABORT("Routine called in non-parallel case.")
526 : #endif
527 :
528 279596 : CALL timestop(handle)
529 :
530 279596 : END SUBROUTINE cp_fm_redistribute_end
531 :
532 : #if defined(__parallel)
533 : ! **************************************************************************************************
534 : !> \brief Determines the largest number of CPUs a matrix can be distributed on without any of the
535 : !> processors getting a zero-width column (currently only needed for ELPA).
536 : !> \param matrix the matrix that will be diagonalized
537 : !> \return the maximum number of CPUs for ELPA
538 : ! **************************************************************************************************
539 72 : FUNCTION cp_cfm_max_ncpu_non_zero_column(matrix) RESULT(ncpu)
540 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
541 : INTEGER :: ncpu
542 :
543 : INTEGER :: gcd_max, ipe, jpe, ncol_block, &
544 : ncol_global, npcol, nrow_block, &
545 : nrow_global, num_pe_old, nzero
546 72 : INTEGER, DIMENSION(:), POINTER :: ncol_locals
547 : INTEGER, EXTERNAL :: numroc
548 :
549 72 : NULLIFY (ncol_locals)
550 : ! First check if there are any zero width columns in current layout
551 : CALL cp_fm_struct_get(matrix%matrix_struct, ncol_locals=ncol_locals, &
552 : nrow_global=nrow_global, ncol_global=ncol_global, &
553 72 : nrow_block=nrow_block, ncol_block=ncol_block)
554 144 : nzero = COUNT(ncol_locals == 0)
555 72 : num_pe_old = matrix%matrix_struct%para_env%num_pe
556 72 : ncpu = num_pe_old - nzero
557 :
558 : ! Avoid layouts with odd number of CPUs (blacs grid layout will be square)
559 72 : IF (ncpu > 2) THEN
560 0 : ncpu = ncpu - MODULO(ncpu, 2)
561 : END IF
562 :
563 : ! if there are no zero-width columns and the number of processors was even, leave it at that
564 72 : IF (ncpu == num_pe_old) THEN
565 : RETURN
566 : END IF
567 :
568 : ! Iteratively search for the maximum number of CPUs for ELPA
569 : ! On each step, we test whether the blacs grid created with ncpu processes
570 : ! contains any columns with zero width
571 0 : DO WHILE (ncpu > 1)
572 : ! Determine layout of new blacs grid with ncpu CPUs
573 : ! (snippet copied from cp_blacs_env.F:cp_blacs_env_create)
574 0 : gcd_max = -1
575 0 : DO ipe = 1, CEILING(SQRT(REAL(ncpu, dp)))
576 0 : jpe = ncpu/ipe
577 0 : IF (ipe*jpe /= ncpu) THEN
578 : CYCLE
579 : END IF
580 0 : IF (gcd(ipe, jpe) >= gcd_max) THEN
581 0 : npcol = jpe
582 0 : gcd_max = gcd(ipe, jpe)
583 : END IF
584 : END DO
585 :
586 : ! Count the number of processors without any columns
587 : ! (snippet copied from cp_fm_struct.F:cp_fm_struct_create)
588 0 : nzero = 0
589 0 : DO ipe = 0, npcol - 1
590 0 : IF (numroc(ncol_global, ncol_block, ipe, 0, npcol) == 0) THEN
591 0 : nzero = nzero + 1
592 : END IF
593 : END DO
594 :
595 0 : IF (nzero == 0) THEN
596 : EXIT
597 : END IF
598 :
599 0 : ncpu = ncpu - nzero
600 :
601 0 : IF (ncpu > 2) THEN
602 0 : ncpu = ncpu - MODULO(ncpu, 2)
603 : END IF
604 : END DO
605 :
606 72 : END FUNCTION cp_cfm_max_ncpu_non_zero_column
607 : #endif
608 :
609 : ! **************************************************************************************************
610 : !> \brief Determines the optimal number of CPUs for matrix diagonalization and redistributes
611 : !> the input matrices if necessary
612 : !> \param matrix the input cp_cfm_type matrix to be diagonalized
613 : !> \param eigenvectors the cp_cfm_type matrix that will hold the eigenvectors of the input matrix
614 : !> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
615 : !> or a pointer to the original matrix if no redistribution is required
616 : !> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
617 : !> matrix if no redistribution is required
618 : !> \param caller_is_elpa flag that determines if ELPA is used for diagonalization
619 : !> \param redist_info get info about the redistribution
620 : !> \par History
621 : !> - [08.2026] created by mirroring cp_fm_redistribute_start for complex matrices
622 : ! **************************************************************************************************
623 72 : SUBROUTINE cp_cfm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
624 : caller_is_elpa, redist_info)
625 :
626 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
627 : TYPE(cp_cfm_type), INTENT(OUT) :: matrix_new, eigenvectors_new
628 : LOGICAL, OPTIONAL, INTENT(IN) :: caller_is_elpa
629 :
630 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_redistribute_start'
631 :
632 : INTEGER :: handle
633 : LOGICAL :: is_elpa
634 : TYPE(cp_fm_redistribute_info), OPTIONAL, INTENT(OUT) :: redist_info
635 :
636 : #if defined(__parallel)
637 : COMPLEX(KIND=dp) :: fake_local_data(1, 1)
638 : INTEGER :: fake_descriptor(9), mepos_old, &
639 : io_unit, ngroups, ncol_block, blksize, nrow_block
640 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_new
641 : TYPE(mp_para_env_type), POINTER :: para_env
642 : TYPE(cp_logger_type), POINTER :: logger
643 : TYPE(cp_fm_redistribute_info) :: rdinfo
644 : #endif
645 :
646 72 : CALL timeset(routineN, handle)
647 72 : is_elpa = .FALSE.
648 72 : IF (PRESENT(caller_is_elpa)) THEN
649 : #if defined(__ELPA)
650 72 : is_elpa = caller_is_elpa
651 : #else
652 : CPABORT("CP2K compiled without the ELPA library.")
653 : #endif
654 : END IF
655 :
656 : #if defined(__parallel)
657 :
658 72 : logger => cp_get_default_logger()
659 72 : io_unit = cp_logger_get_default_io_unit(logger)
660 :
661 : ! first figure out the optimal number of cpus
662 : ! this is pure heuristics, the defaults are based on rosa timings
663 : ! that demonstrate that timings go up sharply if too many tasks are used
664 : ! we take a multiple of 4, and approximately n/60
665 72 : para_env => matrix%matrix_struct%para_env
666 72 : mepos_old = para_env%mepos
667 72 : ncol_block = -1 ! normally we also want to adjust the block size according to the optimal # of CPUs
668 72 : nrow_block = -1
669 72 : blksize = -1
670 :
671 72 : rdinfo%matrix_order = matrix%matrix_struct%nrow_global
672 72 : rdinfo%num_pe_old = para_env%num_pe
673 72 : rdinfo%num_pe_opt = cp_fm_diag_get_optimal_ncpu(rdinfo%matrix_order)
674 72 : rdinfo%num_pe_new = rdinfo%num_pe_opt
675 : rdinfo%num_pe_max_nz_col = -1
676 : rdinfo%redistribute = .FALSE.
677 :
678 72 : IF (is_elpa) THEN
679 : ! with ELPA we don't have to redistribute if not necessary (scales, unlike ScaLAPACK)
680 72 : rdinfo%num_pe_new = rdinfo%num_pe_old
681 :
682 : ! BUT: Diagonalization with ELPA fails when a processor column has zero width
683 : ! Determine the maximum number of CPUs the matrix can be distributed without zero-width columns
684 : ! for the current block size.
685 72 : rdinfo%num_pe_max_nz_col = cp_cfm_max_ncpu_non_zero_column(matrix)
686 :
687 : ! if the user wants to redistribute to the ScaLAPACK optimal number of CPUs anyway, let him if it's safe.
688 72 : IF (work_redistribute%elpa_force_redistribute .AND. rdinfo%num_pe_opt < rdinfo%num_pe_max_nz_col) THEN
689 : ! Use heuristics to determine the need for redistribution (when num_pe_opt is smaller than the safe maximum)
690 : ! in this case we can also take the block size used for ScaLAPACK
691 0 : rdinfo%num_pe_new = rdinfo%num_pe_opt
692 72 : ELSE IF (rdinfo%num_pe_old > rdinfo%num_pe_max_nz_col) THEN
693 : ! Otherwise, only redistribute if we have to
694 0 : rdinfo%num_pe_new = rdinfo%num_pe_max_nz_col
695 : ! do NOT let cp_fm_struct_create automatically adjust the block size because the
696 : ! calculated number of processors such that no block has 0 columns wouldn't match (see #578):
697 : ! if the automatically chosen block size is larger than the present one we would still end
698 : ! up with empty processors
699 : END IF
700 :
701 72 : CALL cp_cfm_get_info(matrix, ncol_block=ncol_block, nrow_block=nrow_block)
702 :
703 : ! On GPUs, ELPA requires the block size to be a power of 2
704 72 : blksize = 1
705 360 : DO WHILE (2*blksize <= MIN(nrow_block, ncol_block))
706 72 : blksize = blksize*2
707 : END DO
708 72 : nrow_block = blksize
709 72 : ncol_block = blksize
710 : END IF
711 :
712 : ! finally, only redistribute if we're going to use less CPUs than before or changed the block size
713 : rdinfo%redistribute = (rdinfo%num_pe_old > rdinfo%num_pe_new) .OR. (blksize >= 0 .AND. &
714 72 : ((blksize /= matrix%matrix_struct%ncol_block) .OR. (blksize /= matrix%matrix_struct%nrow_block)))
715 :
716 72 : IF (work_redistribute%should_print .AND. io_unit > 0) THEN
717 0 : IF (is_elpa) THEN
718 0 : IF (work_redistribute%elpa_force_redistribute) THEN
719 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
720 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "YES"
721 : ELSE
722 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
723 0 : "CP_FM_DIAG| Force redistribute (ELPA):", "NO"
724 : END IF
725 : END IF
726 0 : CALL rdinfo%write(io_unit)
727 : END IF
728 72 : CALL para_env%sync()
729 :
730 : ! if the optimal is smaller than num_pe, we will redistribute the input matrix
731 72 : IF (rdinfo%redistribute) THEN
732 : ! split comm, the first num_pe_new tasks will do the work
733 216 : ALLOCATE (work_redistribute%group_distribution(0:rdinfo%num_pe_old - 1))
734 72 : ALLOCATE (work_redistribute%group_partition(0:1))
735 216 : work_redistribute%group_partition = [rdinfo%num_pe_new, rdinfo%num_pe_old - rdinfo%num_pe_new]
736 72 : ALLOCATE (work_redistribute%para_env_new)
737 : CALL work_redistribute%para_env_new%from_split( &
738 : comm=para_env, ngroups=ngroups, group_distribution=work_redistribute%group_distribution, &
739 72 : n_subgroups=2, group_partition=work_redistribute%group_partition)
740 :
741 72 : IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
742 :
743 : ! create blacs, should inherit the preferences for the layout and so on, from the higher level
744 36 : NULLIFY (work_redistribute%blacs_env_new)
745 36 : CALL cp_blacs_env_create(blacs_env=work_redistribute%blacs_env_new, para_env=work_redistribute%para_env_new)
746 :
747 : ! create new matrix
748 36 : NULLIFY (fm_struct_new)
749 36 : IF (nrow_block == -1 .OR. ncol_block == -1) THEN
750 : CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
751 : para_env=work_redistribute%para_env_new, &
752 : context=work_redistribute%blacs_env_new, &
753 : nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
754 0 : ncol_block=ncol_block, nrow_block=nrow_block)
755 : ELSE
756 : CALL cp_fm_struct_create(fmstruct=fm_struct_new, &
757 : para_env=work_redistribute%para_env_new, &
758 : context=work_redistribute%blacs_env_new, &
759 : nrow_global=rdinfo%matrix_order, ncol_global=rdinfo%matrix_order, &
760 36 : ncol_block=ncol_block, nrow_block=nrow_block, force_block=.TRUE.)
761 : END IF
762 36 : CALL cp_cfm_create(matrix_new, matrix_struct=fm_struct_new, name="zheevd_new_mat")
763 36 : CALL cp_cfm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="zheevd_new_vec")
764 36 : CALL cp_fm_struct_release(fm_struct_new)
765 :
766 : ! redistribute old
767 : CALL pzgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
768 : matrix%matrix_struct%descriptor, &
769 : matrix_new%local_data(1, 1), 1, 1, matrix_new%matrix_struct%descriptor, &
770 36 : matrix%matrix_struct%context)
771 : ELSE
772 : ! these tasks must help redistribute (they own part of the data),
773 : ! but need fake 'new' data, and their descriptor must indicate this with -1
774 : ! see also scalapack comments on pzgemr2d
775 360 : fake_descriptor = -1
776 : CALL pzgemr2d(rdinfo%matrix_order, rdinfo%matrix_order, matrix%local_data(1, 1), 1, 1, &
777 : matrix%matrix_struct%descriptor, &
778 : fake_local_data(1, 1), 1, 1, fake_descriptor, &
779 36 : matrix%matrix_struct%context)
780 : END IF
781 : ELSE
782 : ! No need to redistribute, just return pointers to the original arrays
783 0 : matrix_new = matrix
784 0 : eigenvectors_new = eigenvectors
785 : END IF
786 :
787 72 : IF (PRESENT(redist_info)) THEN
788 72 : redist_info = rdinfo
789 : END IF
790 : #else
791 :
792 : MARK_USED(matrix)
793 : MARK_USED(eigenvectors)
794 : MARK_USED(matrix_new)
795 : MARK_USED(eigenvectors_new)
796 : MARK_USED(redist_info)
797 : CPABORT("Routine called in non-parallel case.")
798 : #endif
799 :
800 72 : CALL timestop(handle)
801 :
802 72 : END SUBROUTINE cp_cfm_redistribute_start
803 :
804 : ! **************************************************************************************************
805 : !> \brief Redistributes eigenvectors and eigenvalues back to the original communicator group
806 : !> \param matrix the input cp_cfm_type matrix to be diagonalized
807 : !> \param eigenvectors the cp_cfm_type matrix that will hold the eigenvectors of the input matrix
808 : !> \param eig global array holding the eigenvalues of the input matrixmatrix
809 : !> \param matrix_new the redistributed input matrix which will subsequently be diagonalized,
810 : !> or a pointer to the original matrix if no redistribution is required
811 : !> \param eigenvectors_new the redistributed eigenvectors matrix, or a pointer to the original
812 : !> matrix if no redistribution is required
813 : !> \par History
814 : !> - [08.2026] created by mirroring cp_fm_redistribute_end for complex matrices
815 : ! **************************************************************************************************
816 72 : SUBROUTINE cp_cfm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
817 :
818 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
819 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: eig
820 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_new, eigenvectors_new
821 :
822 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_redistribute_end'
823 :
824 : INTEGER :: handle
825 : #if defined(__parallel)
826 : COMPLEX(KIND=dp) :: fake_local_data(1, 1)
827 : INTEGER :: fake_descriptor(9), mepos_old, n
828 : TYPE(mp_para_env_type), POINTER :: para_env
829 : #endif
830 :
831 72 : CALL timeset(routineN, handle)
832 :
833 : #if defined(__parallel)
834 :
835 : ! Check if matrix was redistributed
836 72 : IF (ASSOCIATED(work_redistribute%group_distribution)) THEN
837 72 : n = matrix%matrix_struct%nrow_global
838 72 : para_env => matrix%matrix_struct%para_env
839 72 : mepos_old = para_env%mepos
840 :
841 72 : IF (work_redistribute%group_distribution(mepos_old) == 0) THEN
842 : ! redistribute results on CPUs that hold the redistributed matrix
843 : CALL pzgemr2d(n, n, eigenvectors_new%local_data(1, 1), 1, 1, eigenvectors_new%matrix_struct%descriptor, &
844 : eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
845 36 : eigenvectors%matrix_struct%context)
846 36 : CALL cp_cfm_release(matrix_new)
847 36 : CALL cp_cfm_release(eigenvectors_new)
848 : ELSE
849 : ! these tasks must help redistribute (they own part of the data),
850 : ! but need fake 'new' data, and their descriptor must indicate this with -1
851 : ! see also scalapack comments on pzgemr2d
852 360 : fake_descriptor = -1
853 : CALL pzgemr2d(n, n, fake_local_data(1, 1), 1, 1, fake_descriptor, &
854 : eigenvectors%local_data(1, 1), 1, 1, eigenvectors%matrix_struct%descriptor, &
855 36 : eigenvectors%matrix_struct%context)
856 : END IF
857 : ! free work
858 72 : CALL cp_fm_redistribute_work_finalize(work_redistribute%group_distribution(mepos_old) == 0)
859 :
860 : ! finally, also the eigenvalues need to end up on the non-group member tasks
861 5688 : CALL para_env%bcast(eig, 0)
862 : END IF
863 :
864 : #else
865 :
866 : MARK_USED(matrix)
867 : MARK_USED(eigenvectors)
868 : MARK_USED(eig)
869 : MARK_USED(matrix_new)
870 : MARK_USED(eigenvectors_new)
871 : CPABORT("Routine called in non-parallel case.")
872 : #endif
873 :
874 72 : CALL timestop(handle)
875 :
876 72 : END SUBROUTINE cp_cfm_redistribute_end
877 :
878 0 : END MODULE cp_fm_diag_utils
|
