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 Represents a complex full matrix distributed on many processors.
10 : !> \author Joost VandeVondele, based on Fawzi's cp_fm_* routines
11 : ! **************************************************************************************************
12 : MODULE cp_cfm_types
13 : USE cp_blacs_env, ONLY: cp_blacs_env_type
14 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
15 : cp_fm_struct_equivalent,&
16 : cp_fm_struct_get,&
17 : cp_fm_struct_release,&
18 : cp_fm_struct_retain,&
19 : cp_fm_struct_type
20 : USE cp_fm_types, ONLY: cp_fm_type
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: z_one,&
23 : z_zero
24 : USE message_passing, ONLY: cp2k_is_parallel,&
25 : mp_any_source,&
26 : mp_para_env_type,&
27 : mp_proc_null,&
28 : mp_request_null,&
29 : mp_request_type,&
30 : mp_waitall
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 : PRIVATE
35 :
36 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_types'
38 : INTEGER, PARAMETER, PRIVATE :: src_tag = 3, dest_tag = 5, send_tag = 7, recv_tag = 11
39 :
40 : PUBLIC :: cp_cfm_type, cp_cfm_p_type, copy_cfm_info_type
41 : PUBLIC :: cp_cfm_cleanup_copy_general, &
42 : cp_cfm_create, &
43 : cp_cfm_finish_copy_general, &
44 : cp_cfm_get_element, &
45 : cp_cfm_get_info, &
46 : cp_cfm_get_submatrix, &
47 : cp_cfm_release, &
48 : cp_cfm_set_all, &
49 : cp_cfm_set_element, &
50 : cp_cfm_set_submatrix, &
51 : cp_cfm_start_copy_general, &
52 : cp_cfm_to_cfm, &
53 : cp_cfm_to_fm, &
54 : cp_fm_to_cfm
55 :
56 : INTERFACE cp_cfm_to_cfm
57 : MODULE PROCEDURE cp_cfm_to_cfm_matrix, & ! a full matrix
58 : cp_cfm_to_cfm_columns ! just a number of columns
59 : END INTERFACE
60 :
61 : ! **************************************************************************************************
62 : !> \brief Represent a complex full matrix.
63 : !> \param name the name of the matrix, used for printing
64 : !> \param matrix_struct structure of this matrix
65 : !> \param local_data array with the data of the matrix (its content depends on
66 : !> the matrix type used: in parallel run it will be in
67 : !> ScaLAPACK format, in sequential run it will simply contain the matrix)
68 : ! **************************************************************************************************
69 : TYPE cp_cfm_type
70 : CHARACTER(len=60) :: name = ""
71 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct => NULL()
72 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: local_data => NULL()
73 : END TYPE cp_cfm_type
74 :
75 : ! **************************************************************************************************
76 : !> \brief Just to build arrays of pointers to matrices.
77 : !> \param matrix the pointer to the matrix
78 : ! **************************************************************************************************
79 : TYPE cp_cfm_p_type
80 : TYPE(cp_cfm_type), POINTER :: matrix => NULL()
81 : END TYPE cp_cfm_p_type
82 :
83 : ! **************************************************************************************************
84 : !> \brief Stores the state of a copy between cp_cfm_start_copy_general
85 : !> and cp_cfm_finish_copy_general.
86 : !> \par History
87 : !> Jan 2017 derived type 'copy_info_type' has been created [Mark T]
88 : !> Jan 2018 the type 'copy_info_type' has been adapted for complex matrices [Sergey Chulkov]
89 : ! **************************************************************************************************
90 : TYPE copy_cfm_info_type
91 : !> number of MPI processes that send data
92 : INTEGER :: send_size = -1
93 : !> number of locally stored rows (1) and columns (2) of the destination matrix
94 : INTEGER, DIMENSION(2) :: nlocal_recv = -1
95 : !> number of rows (1) and columns (2) of the ScaLAPACK block of the source matrix
96 : INTEGER, DIMENSION(2) :: nblock_src = -1
97 : !> BLACS process grid shape of the source matrix: (1) nproc_row, (2) nproc_col
98 : INTEGER, DIMENSION(2) :: src_num_pe = -1
99 : !> displacements into recv_buf
100 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_disp
101 : !> MPI requests for non-blocking receive and send operations
102 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_request, send_request
103 : !> global column and row indices of locally stored elements of the destination matrix
104 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices => NULL(), recv_row_indices => NULL()
105 : !> rank of MPI process with BLACS coordinates (prow, pcol)
106 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: src_blacs2mpi
107 : !> receiving and sending buffers for non-blocking MPI communication
108 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buf, send_buf
109 : END TYPE copy_cfm_info_type
110 :
111 : CONTAINS
112 :
113 : ! **************************************************************************************************
114 : !> \brief Creates a new full matrix with the given structure.
115 : !> \param matrix matrix to be created
116 : !> \param matrix_struct structure of the matrix
117 : !> \param name name of the matrix
118 : !> \param nrow ...
119 : !> \param ncol ...
120 : !> \param set_zero ...
121 : !> \note
122 : !> preferred allocation routine
123 : ! **************************************************************************************************
124 560934 : SUBROUTINE cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
125 : TYPE(cp_cfm_type), INTENT(OUT) :: matrix
126 : TYPE(cp_fm_struct_type), INTENT(IN), TARGET :: matrix_struct
127 : CHARACTER(len=*), INTENT(in), OPTIONAL :: name
128 : INTEGER, INTENT(IN), OPTIONAL :: nrow, ncol
129 : LOGICAL, INTENT(IN), OPTIONAL :: set_zero
130 :
131 : INTEGER :: ncol_global, ncol_local, nrow_global, &
132 : nrow_local
133 : TYPE(cp_blacs_env_type), POINTER :: context
134 : TYPE(cp_fm_struct_type), POINTER :: cfm_struct
135 :
136 560946 : IF (PRESENT(nrow) .OR. PRESENT(ncol)) THEN
137 12 : CALL cp_fm_struct_get(matrix_struct, nrow_global=nrow_global, ncol_global=ncol_global)
138 12 : IF (PRESENT(nrow)) nrow_global = nrow
139 12 : IF (PRESENT(ncol)) ncol_global = ncol
140 : CALL cp_fm_struct_create(cfm_struct, template_fmstruct=matrix_struct, &
141 12 : nrow_global=nrow_global, ncol_global=ncol_global)
142 :
143 12 : context => cfm_struct%context
144 12 : matrix%matrix_struct => cfm_struct
145 12 : CALL cp_fm_struct_retain(matrix%matrix_struct)
146 :
147 12 : nrow_local = cfm_struct%local_leading_dimension
148 12 : ncol_local = MAX(1, cfm_struct%ncol_locals(context%mepos(2)))
149 :
150 12 : CALL cp_fm_struct_release(cfm_struct)
151 : ELSE
152 560922 : context => matrix_struct%context
153 560922 : matrix%matrix_struct => matrix_struct
154 560922 : CALL cp_fm_struct_retain(matrix%matrix_struct)
155 :
156 560922 : nrow_local = matrix_struct%local_leading_dimension
157 560922 : ncol_local = MAX(1, matrix_struct%ncol_locals(context%mepos(2)))
158 : END IF
159 :
160 560934 : NULLIFY (matrix%local_data)
161 2243736 : ALLOCATE (matrix%local_data(nrow_local, ncol_local))
162 :
163 560934 : IF (PRESENT(set_zero)) THEN
164 2170 : IF (set_zero) THEN
165 6044045 : matrix%local_data(1:nrow_local, 1:ncol_local) = z_zero
166 : END IF
167 : END IF
168 :
169 560934 : IF (PRESENT(name)) THEN
170 53338 : matrix%name = name
171 : ELSE
172 507596 : matrix%name = 'full complex matrix'
173 : END IF
174 :
175 560934 : END SUBROUTINE cp_cfm_create
176 :
177 : ! **************************************************************************************************
178 : !> \brief Releases a full matrix.
179 : !> \param matrix the matrix to release
180 : ! **************************************************************************************************
181 585982 : SUBROUTINE cp_cfm_release(matrix)
182 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix
183 :
184 585982 : IF (ASSOCIATED(matrix%local_data)) THEN
185 560934 : DEALLOCATE (matrix%local_data)
186 : END IF
187 585982 : matrix%name = ""
188 585982 : CALL cp_fm_struct_release(matrix%matrix_struct)
189 585982 : END SUBROUTINE cp_cfm_release
190 :
191 : ! **************************************************************************************************
192 : !> \brief Set all elements of the full matrix to alpha. Besides, set all
193 : !> diagonal matrix elements to beta (if given explicitly).
194 : !> \param matrix matrix to initialise
195 : !> \param alpha value of off-diagonal matrix elements
196 : !> \param beta value of diagonal matrix elements (equal to alpha if absent)
197 : !> \date 12.06.2001
198 : !> \author Matthias Krack
199 : !> \version 1.0
200 : ! **************************************************************************************************
201 146078 : SUBROUTINE cp_cfm_set_all(matrix, alpha, beta)
202 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
203 : COMPLEX(kind=dp), INTENT(in) :: alpha
204 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: beta
205 :
206 : INTEGER :: irow_local, nrow_local
207 : #if defined(__parallel)
208 : INTEGER :: icol_local, ncol_local
209 146078 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
210 : #endif
211 :
212 438234 : CALL zcopy(SIZE(matrix%local_data), alpha, 0, matrix%local_data(1, 1), 1)
213 :
214 146078 : IF (PRESENT(beta)) THEN
215 : #if defined(__parallel)
216 : CALL cp_cfm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
217 56336 : row_indices=row_indices, col_indices=col_indices)
218 :
219 56336 : icol_local = 1
220 56336 : irow_local = 1
221 :
222 228616 : DO WHILE (irow_local <= nrow_local .AND. icol_local <= ncol_local)
223 228616 : IF (row_indices(irow_local) < col_indices(icol_local)) THEN
224 0 : irow_local = irow_local + 1
225 172280 : ELSE IF (row_indices(irow_local) > col_indices(icol_local)) THEN
226 70057 : icol_local = icol_local + 1
227 : ELSE
228 102223 : matrix%local_data(irow_local, icol_local) = beta
229 102223 : irow_local = irow_local + 1
230 102223 : icol_local = icol_local + 1
231 : END IF
232 : END DO
233 : #else
234 : nrow_local = MIN(matrix%matrix_struct%nrow_global, matrix%matrix_struct%ncol_global)
235 :
236 : DO irow_local = 1, nrow_local
237 : matrix%local_data(irow_local, irow_local) = beta
238 : END DO
239 : #endif
240 : END IF
241 :
242 146078 : END SUBROUTINE cp_cfm_set_all
243 :
244 : ! **************************************************************************************************
245 : !> \brief Get the matrix element by its global index.
246 : !> \param matrix full matrix
247 : !> \param irow_global global row index
248 : !> \param icol_global global column index
249 : !> \param alpha matrix element
250 : !> \par History
251 : !> , TCH, created
252 : !> always return the answer
253 : ! **************************************************************************************************
254 1546674 : SUBROUTINE cp_cfm_get_element(matrix, irow_global, icol_global, alpha)
255 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
256 : INTEGER, INTENT(in) :: irow_global, icol_global
257 : COMPLEX(kind=dp), INTENT(out) :: alpha
258 :
259 : #if defined(__parallel)
260 : INTEGER :: icol_local, ipcol, iprow, irow_local, &
261 : mypcol, myprow, npcol, nprow
262 : INTEGER, DIMENSION(9) :: desca
263 : TYPE(cp_blacs_env_type), POINTER :: context
264 : #endif
265 :
266 : #if defined(__parallel)
267 1546674 : context => matrix%matrix_struct%context
268 1546674 : myprow = context%mepos(1)
269 1546674 : mypcol = context%mepos(2)
270 1546674 : nprow = context%num_pe(1)
271 1546674 : npcol = context%num_pe(2)
272 :
273 15466740 : desca(:) = matrix%matrix_struct%descriptor(:)
274 :
275 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
276 1546674 : irow_local, icol_local, iprow, ipcol)
277 :
278 1546674 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
279 773337 : alpha = matrix%local_data(irow_local, icol_local)
280 773337 : CALL context%ZGEBS2D('All', ' ', 1, 1, alpha, 1)
281 : ELSE
282 773337 : CALL context%ZGEBR2D('All', ' ', 1, 1, alpha, 1, iprow, ipcol)
283 : END IF
284 :
285 : #else
286 : alpha = matrix%local_data(irow_global, icol_global)
287 : #endif
288 :
289 1546674 : END SUBROUTINE cp_cfm_get_element
290 :
291 : ! **************************************************************************************************
292 : !> \brief Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
293 : !> \param matrix full matrix
294 : !> \param irow_global global row index
295 : !> \param icol_global global column index
296 : !> \param alpha value of the matrix element
297 : !> \date 12.06.2001
298 : !> \author Matthias Krack
299 : !> \version 1.0
300 : ! **************************************************************************************************
301 1205316 : SUBROUTINE cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
302 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
303 : INTEGER, INTENT(in) :: irow_global, icol_global
304 : COMPLEX(kind=dp), INTENT(in) :: alpha
305 :
306 : #if defined(__parallel)
307 : INTEGER :: icol_local, ipcol, iprow, irow_local, &
308 : mypcol, myprow, npcol, nprow
309 : INTEGER, DIMENSION(9) :: desca
310 : TYPE(cp_blacs_env_type), POINTER :: context
311 : #endif
312 :
313 : #if defined(__parallel)
314 1205316 : context => matrix%matrix_struct%context
315 1205316 : myprow = context%mepos(1)
316 1205316 : mypcol = context%mepos(2)
317 1205316 : nprow = context%num_pe(1)
318 1205316 : npcol = context%num_pe(2)
319 :
320 12053160 : desca(:) = matrix%matrix_struct%descriptor(:)
321 :
322 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
323 1205316 : irow_local, icol_local, iprow, ipcol)
324 :
325 1205316 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
326 602658 : matrix%local_data(irow_local, icol_local) = alpha
327 : END IF
328 :
329 : #else
330 : matrix%local_data(irow_global, icol_global) = alpha
331 : #endif
332 :
333 1205316 : END SUBROUTINE cp_cfm_set_element
334 :
335 : ! **************************************************************************************************
336 : !> \brief Extract a sub-matrix from the full matrix:
337 : !> op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n_rows,start_col:start_col+n_cols).
338 : !> Sub-matrix 'target_m' is replicated on each CPU. Using this call is expensive.
339 : !> \param fm full matrix you want to get the elements from
340 : !> \param target_m 2-D array to store the extracted sub-matrix
341 : !> \param start_row global row index of the matrix element target_m(1,1) (defaults to 1)
342 : !> \param start_col global column index of the matrix element target_m(1,1) (defaults to 1)
343 : !> \param n_rows number of rows to extract (defaults to size(op(target_m),1))
344 : !> \param n_cols number of columns to extract (defaults to size(op(target_m),2))
345 : !> \param transpose indicates that the extracted sub-matrix target_m should be transposed:
346 : !> op(target_m) = target_m^T if .TRUE.,
347 : !> op(target_m) = target_m if .FALSE. (defaults to false)
348 : !> \par History
349 : !> * 04.2016 created borrowing from Fawzi's cp_fm_get_submatrix [Lianheng Tong]
350 : !> * 01.2018 drop innermost conditional branching [Sergey Chulkov]
351 : !> \author Lianheng Tong
352 : !> \note
353 : !> Optimized for full column updates. The matrix target_m is replicated and valid on all CPUs.
354 : ! **************************************************************************************************
355 2606 : SUBROUTINE cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
356 : TYPE(cp_cfm_type), INTENT(IN) :: fm
357 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(out) :: target_m
358 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
359 : LOGICAL, INTENT(in), OPTIONAL :: transpose
360 :
361 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_get_submatrix'
362 :
363 2606 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
364 : INTEGER :: end_col_global, end_col_local, end_row_global, end_row_local, handle, i, j, &
365 : ncol_global, ncol_local, nrow_global, nrow_local, start_col_global, start_col_local, &
366 : start_row_global, start_row_local, this_col
367 2606 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
368 : LOGICAL :: do_zero, tr_a
369 : TYPE(mp_para_env_type), POINTER :: para_env
370 :
371 2606 : CALL timeset(routineN, handle)
372 :
373 7818 : IF (SIZE(target_m) /= 0) THEN
374 : #if defined(__parallel)
375 2606 : do_zero = .TRUE.
376 : #else
377 : do_zero = .FALSE.
378 : #endif
379 :
380 2606 : tr_a = .FALSE.
381 2606 : IF (PRESENT(transpose)) tr_a = transpose
382 :
383 : ! find out the first and last global row/column indices
384 2606 : start_row_global = 1
385 2606 : start_col_global = 1
386 2606 : IF (PRESENT(start_row)) start_row_global = start_row
387 2606 : IF (PRESENT(start_col)) start_col_global = start_col
388 :
389 2606 : IF (tr_a) THEN
390 0 : end_row_global = SIZE(target_m, 2)
391 0 : end_col_global = SIZE(target_m, 1)
392 : ELSE
393 2606 : end_row_global = SIZE(target_m, 1)
394 2606 : end_col_global = SIZE(target_m, 2)
395 : END IF
396 2606 : IF (PRESENT(n_rows)) end_row_global = n_rows
397 2606 : IF (PRESENT(n_cols)) end_col_global = n_cols
398 :
399 2606 : end_row_global = end_row_global + start_row_global - 1
400 2606 : end_col_global = end_col_global + start_col_global - 1
401 :
402 : CALL cp_cfm_get_info(matrix=fm, &
403 : nrow_global=nrow_global, ncol_global=ncol_global, &
404 : nrow_local=nrow_local, ncol_local=ncol_local, &
405 2606 : row_indices=row_indices, col_indices=col_indices)
406 2606 : IF (end_row_global > nrow_global) THEN
407 : end_row_global = nrow_global
408 : do_zero = .TRUE.
409 : END IF
410 2606 : IF (end_col_global > ncol_global) THEN
411 : end_col_global = ncol_global
412 : do_zero = .TRUE.
413 : END IF
414 :
415 : ! find out row/column indices of locally stored matrix elements that needs to be copied.
416 : ! Arrays row_indices and col_indices are assumed to be sorted in ascending order
417 13082 : DO start_row_local = 1, nrow_local
418 13082 : IF (row_indices(start_row_local) >= start_row_global) EXIT
419 : END DO
420 :
421 14468 : DO end_row_local = start_row_local, nrow_local
422 14468 : IF (row_indices(end_row_local) > end_row_global) EXIT
423 : END DO
424 2606 : end_row_local = end_row_local - 1
425 :
426 10928 : DO start_col_local = 1, ncol_local
427 10928 : IF (col_indices(start_col_local) >= start_col_global) EXIT
428 : END DO
429 :
430 43206 : DO end_col_local = start_col_local, ncol_local
431 43206 : IF (col_indices(end_col_local) > end_col_global) EXIT
432 : END DO
433 2606 : end_col_local = end_col_local - 1
434 :
435 2606 : para_env => fm%matrix_struct%para_env
436 2606 : local_data => fm%local_data
437 :
438 : ! wipe the content of the target matrix if:
439 : ! * the source matrix is distributed across a number of processes, or
440 : ! * not all elements of the target matrix will be assigned, e.g.
441 : ! when the target matrix is larger then the source matrix
442 : IF (do_zero) THEN
443 7818 : CALL zcopy(SIZE(target_m), z_zero, 0, target_m(1, 1), 1)
444 : END IF
445 :
446 2606 : IF (tr_a) THEN
447 0 : DO j = start_col_local, end_col_local
448 0 : this_col = col_indices(j) - start_col_global + 1
449 0 : DO i = start_row_local, end_row_local
450 0 : target_m(this_col, row_indices(i) - start_row_global + 1) = local_data(i, j)
451 : END DO
452 : END DO
453 : ELSE
454 43206 : DO j = start_col_local, end_col_local
455 40600 : this_col = col_indices(j) - start_col_global + 1
456 362037 : DO i = start_row_local, end_row_local
457 359431 : target_m(row_indices(i) - start_row_global + 1, this_col) = local_data(i, j)
458 : END DO
459 : END DO
460 : END IF
461 :
462 1039882 : CALL para_env%sum(target_m)
463 : END IF
464 :
465 2606 : CALL timestop(handle)
466 2606 : END SUBROUTINE cp_cfm_get_submatrix
467 :
468 : ! **************************************************************************************************
469 : !> \brief Set a sub-matrix of the full matrix:
470 : !> matrix(start_row:start_row+n_rows,start_col:start_col+n_cols)
471 : !> = alpha*op(new_values)(1:n_rows,1:n_cols) +
472 : !> beta*matrix(start_row:start_row+n_rows,start_col:start_col+n_cols)
473 : !> \param matrix full to update
474 : !> \param new_values replicated 2-D array that holds new elements of the updated sub-matrix
475 : !> \param start_row global row index of the matrix element new_values(1,1) (defaults to 1)
476 : !> \param start_col global column index of the matrix element new_values(1,1) (defaults to 1)
477 : !> \param n_rows number of rows to update (defaults to size(op(new_values),1))
478 : !> \param n_cols number of columns to update (defaults to size(op(new_values),2))
479 : !> \param alpha scale factor for the new values (defaults to (1.0,0.0))
480 : !> \param beta scale factor for the old values (defaults to (0.0,0.0))
481 : !> \param transpose indicates that the matrix new_values should be transposed:
482 : !> op(new_values) = new_values^T if .TRUE.,
483 : !> op(new_values) = new_values if .FALSE. (defaults to false)
484 : !> \par History
485 : !> * 04.2016 created borrowing from Fawzi's cp_fm_set_submatrix [Lianheng Tong]
486 : !> * 01.2018 drop innermost conditional branching [Sergey Chulkov]
487 : !> \author Lianheng Tong
488 : !> \note
489 : !> Optimized for alpha=(1.0,0.0), beta=(0.0,0.0)
490 : !> All matrix elements 'new_values' need to be valid on all CPUs
491 : ! **************************************************************************************************
492 1306 : SUBROUTINE cp_cfm_set_submatrix(matrix, new_values, start_row, &
493 : start_col, n_rows, n_cols, alpha, beta, transpose)
494 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
495 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(in) :: new_values
496 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
497 : COMPLEX(kind=dp), INTENT(in), OPTIONAL :: alpha, beta
498 : LOGICAL, INTENT(in), OPTIONAL :: transpose
499 :
500 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_set_submatrix'
501 :
502 : COMPLEX(kind=dp) :: al, be
503 1306 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
504 : INTEGER :: end_col_global, end_col_local, end_row_global, end_row_local, handle, i, j, &
505 : ncol_global, ncol_local, nrow_global, nrow_local, start_col_global, start_col_local, &
506 : start_row_global, start_row_local, this_col
507 1306 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
508 : LOGICAL :: tr_a
509 :
510 1306 : CALL timeset(routineN, handle)
511 :
512 1306 : al = z_one
513 1306 : be = z_zero
514 1306 : IF (PRESENT(alpha)) al = alpha
515 1306 : IF (PRESENT(beta)) be = beta
516 :
517 : ! find out the first and last global row/column indices
518 1306 : start_row_global = 1
519 1306 : start_col_global = 1
520 1306 : IF (PRESENT(start_row)) start_row_global = start_row
521 1306 : IF (PRESENT(start_col)) start_col_global = start_col
522 :
523 1306 : tr_a = .FALSE.
524 1306 : IF (PRESENT(transpose)) tr_a = transpose
525 :
526 0 : IF (tr_a) THEN
527 0 : end_row_global = SIZE(new_values, 2)
528 0 : end_col_global = SIZE(new_values, 1)
529 : ELSE
530 1306 : end_row_global = SIZE(new_values, 1)
531 1306 : end_col_global = SIZE(new_values, 2)
532 : END IF
533 1306 : IF (PRESENT(n_rows)) end_row_global = n_rows
534 1306 : IF (PRESENT(n_cols)) end_col_global = n_cols
535 :
536 1306 : end_row_global = end_row_global + start_row_global - 1
537 1306 : end_col_global = end_col_global + start_col_global - 1
538 :
539 : CALL cp_cfm_get_info(matrix=matrix, &
540 : nrow_global=nrow_global, ncol_global=ncol_global, &
541 : nrow_local=nrow_local, ncol_local=ncol_local, &
542 1306 : row_indices=row_indices, col_indices=col_indices)
543 1306 : IF (end_row_global > nrow_global) end_row_global = nrow_global
544 1306 : IF (end_col_global > ncol_global) end_col_global = ncol_global
545 :
546 : ! find out row/column indices of locally stored matrix elements that needs to be set.
547 : ! Arrays row_indices and col_indices are assumed to be sorted in ascending order
548 129036 : DO start_row_local = 1, nrow_local
549 129036 : IF (row_indices(start_row_local) >= start_row_global) EXIT
550 : END DO
551 :
552 9073 : DO end_row_local = start_row_local, nrow_local
553 9073 : IF (row_indices(end_row_local) > end_row_global) EXIT
554 : END DO
555 1306 : end_row_local = end_row_local - 1
556 :
557 242810 : DO start_col_local = 1, ncol_local
558 242810 : IF (col_indices(start_col_local) >= start_col_global) EXIT
559 : END DO
560 :
561 37300 : DO end_col_local = start_col_local, ncol_local
562 37300 : IF (col_indices(end_col_local) > end_col_global) EXIT
563 : END DO
564 1306 : end_col_local = end_col_local - 1
565 :
566 1306 : local_data => matrix%local_data
567 :
568 1306 : IF (al == z_one .AND. be == z_zero) THEN
569 1306 : IF (tr_a) THEN
570 0 : DO j = start_col_local, end_col_local
571 0 : this_col = col_indices(j) - start_col_global + 1
572 0 : DO i = start_row_local, end_row_local
573 0 : local_data(i, j) = new_values(this_col, row_indices(i) - start_row_global + 1)
574 : END DO
575 : END DO
576 : ELSE
577 37300 : DO j = start_col_local, end_col_local
578 35994 : this_col = col_indices(j) - start_col_global + 1
579 301908 : DO i = start_row_local, end_row_local
580 300602 : local_data(i, j) = new_values(row_indices(i) - start_row_global + 1, this_col)
581 : END DO
582 : END DO
583 : END IF
584 : ELSE
585 0 : IF (tr_a) THEN
586 0 : DO j = start_col_local, end_col_local
587 0 : this_col = col_indices(j) - start_col_global + 1
588 0 : DO i = start_row_local, end_row_local
589 : local_data(i, j) = al*new_values(this_col, row_indices(i) - start_row_global + 1) + &
590 0 : be*local_data(i, j)
591 : END DO
592 : END DO
593 : ELSE
594 0 : DO j = start_col_local, end_col_local
595 0 : this_col = col_indices(j) - start_col_global + 1
596 0 : DO i = start_row_local, end_row_local
597 : local_data(i, j) = al*new_values(row_indices(i) - start_row_global + 1, this_col) + &
598 0 : be*local_data(i, j)
599 : END DO
600 : END DO
601 : END IF
602 : END IF
603 :
604 1306 : CALL timestop(handle)
605 1306 : END SUBROUTINE cp_cfm_set_submatrix
606 :
607 : ! **************************************************************************************************
608 : !> \brief Returns information about a full matrix.
609 : !> \param matrix matrix
610 : !> \param name name of the matrix
611 : !> \param nrow_global total number of rows
612 : !> \param ncol_global total number of columns
613 : !> \param nrow_block number of rows per ScaLAPACK block
614 : !> \param ncol_block number of columns per ScaLAPACK block
615 : !> \param nrow_local number of locally stored rows
616 : !> \param ncol_local number of locally stored columns
617 : !> \param row_indices global indices of locally stored rows
618 : !> \param col_indices global indices of locally stored columns
619 : !> \param local_data locally stored matrix elements
620 : !> \param context BLACS context
621 : !> \param matrix_struct matrix structure
622 : !> \param para_env parallel environment
623 : !> \date 12.06.2001
624 : !> \author Matthias Krack
625 : !> \version 1.0
626 : ! **************************************************************************************************
627 1519230 : SUBROUTINE cp_cfm_get_info(matrix, name, nrow_global, ncol_global, &
628 : nrow_block, ncol_block, nrow_local, ncol_local, &
629 : row_indices, col_indices, local_data, context, &
630 : matrix_struct, para_env)
631 : TYPE(cp_cfm_type), INTENT(IN) :: matrix
632 : CHARACTER(len=*), INTENT(OUT), OPTIONAL :: name
633 : INTEGER, INTENT(OUT), OPTIONAL :: nrow_global, ncol_global, nrow_block, &
634 : ncol_block, nrow_local, ncol_local
635 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: row_indices, col_indices
636 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
637 : OPTIONAL, POINTER :: local_data
638 : TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: context
639 : TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: matrix_struct
640 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
641 :
642 0 : IF (PRESENT(name)) name = matrix%name
643 1519230 : IF (PRESENT(matrix_struct)) matrix_struct => matrix%matrix_struct
644 1519230 : IF (PRESENT(local_data)) local_data => matrix%local_data ! not hiding things anymore :-(
645 :
646 : CALL cp_fm_struct_get(matrix%matrix_struct, nrow_local=nrow_local, &
647 : ncol_local=ncol_local, nrow_global=nrow_global, &
648 : ncol_global=ncol_global, nrow_block=nrow_block, &
649 : ncol_block=ncol_block, context=context, &
650 1519230 : row_indices=row_indices, col_indices=col_indices, para_env=para_env)
651 :
652 1519230 : END SUBROUTINE cp_cfm_get_info
653 :
654 : ! **************************************************************************************************
655 : !> \brief Copy content of a full matrix into another full matrix of the same size.
656 : !> \param source source matrix
657 : !> \param destination destination matrix
658 : !> \author Joost VandeVondele
659 : ! **************************************************************************************************
660 556584 : SUBROUTINE cp_cfm_to_cfm_matrix(source, destination)
661 : TYPE(cp_cfm_type), INTENT(IN) :: source, destination
662 :
663 : INTEGER :: npcol, nprow
664 :
665 556584 : nprow = source%matrix_struct%context%num_pe(1)
666 556584 : npcol = source%matrix_struct%context%num_pe(2)
667 :
668 556584 : IF (.NOT. cp2k_is_parallel .OR. &
669 : cp_fm_struct_equivalent(source%matrix_struct, &
670 : destination%matrix_struct)) THEN
671 556584 : IF (SIZE(source%local_data, 1) /= SIZE(destination%local_data, 1) .OR. &
672 : SIZE(source%local_data, 2) /= SIZE(destination%local_data, 2)) THEN
673 0 : CPABORT("internal local_data has different sizes")
674 : END IF
675 1669752 : CALL zcopy(SIZE(source%local_data), source%local_data(1, 1), 1, destination%local_data(1, 1), 1)
676 : ELSE
677 0 : IF (source%matrix_struct%nrow_global /= destination%matrix_struct%nrow_global) THEN
678 0 : CPABORT("cannot copy between full matrixes of differen sizes")
679 : END IF
680 0 : IF (source%matrix_struct%ncol_global /= destination%matrix_struct%ncol_global) THEN
681 0 : CPABORT("cannot copy between full matrixes of differen sizes")
682 : END IF
683 : #if defined(__parallel)
684 : CALL pzcopy(source%matrix_struct%nrow_global* &
685 : source%matrix_struct%ncol_global, &
686 : source%local_data(1, 1), 1, 1, source%matrix_struct%descriptor, 1, &
687 0 : destination%local_data(1, 1), 1, 1, destination%matrix_struct%descriptor, 1)
688 : #else
689 : CPABORT("cp_cfm_to_cfm_matrix: pzcopy needs __parallel build flag")
690 : #endif
691 : END IF
692 556584 : END SUBROUTINE cp_cfm_to_cfm_matrix
693 :
694 : ! **************************************************************************************************
695 : !> \brief Copy a number of sequential columns of a full matrix into another full matrix.
696 : !> \param msource source matrix
697 : !> \param mtarget destination matrix
698 : !> \param ncol number of columns to copy
699 : !> \param source_start global index of the first column to copy within the source matrix
700 : !> \param target_start global index of the first column to copy within the destination matrix
701 : ! **************************************************************************************************
702 66616 : SUBROUTINE cp_cfm_to_cfm_columns(msource, mtarget, ncol, source_start, &
703 : target_start)
704 :
705 : TYPE(cp_cfm_type), INTENT(IN) :: msource, mtarget
706 : INTEGER, INTENT(IN) :: ncol
707 : INTEGER, INTENT(IN), OPTIONAL :: source_start, target_start
708 :
709 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_to_cfm_columns'
710 :
711 66616 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: a, b
712 : INTEGER :: handle, n, ss, ts
713 : #if defined(__parallel)
714 : INTEGER :: i
715 : INTEGER, DIMENSION(9) :: desca, descb
716 : #endif
717 :
718 66616 : CALL timeset(routineN, handle)
719 :
720 66616 : ss = 1
721 66616 : ts = 1
722 :
723 66616 : IF (PRESENT(source_start)) ss = source_start
724 66616 : IF (PRESENT(target_start)) ts = target_start
725 :
726 66616 : n = msource%matrix_struct%nrow_global
727 :
728 66616 : a => msource%local_data
729 66616 : b => mtarget%local_data
730 :
731 : #if defined(__parallel)
732 666160 : desca(:) = msource%matrix_struct%descriptor(:)
733 666160 : descb(:) = mtarget%matrix_struct%descriptor(:)
734 1312442 : DO i = 0, ncol - 1
735 1312442 : CALL pzcopy(n, a(1, 1), 1, ss + i, desca, 1, b(1, 1), 1, ts + i, descb, 1)
736 : END DO
737 : #else
738 : CALL zcopy(ncol*n, a(1, ss), 1, b(1, ts), 1)
739 : #endif
740 :
741 66616 : CALL timestop(handle)
742 :
743 66616 : END SUBROUTINE cp_cfm_to_cfm_columns
744 :
745 : ! **************************************************************************************************
746 : !> \brief Copy just a triangular matrix.
747 : !> \param msource source matrix
748 : !> \param mtarget target matrix
749 : !> \param uplo 'U' for upper triangular, 'L' for lower triangular
750 : ! **************************************************************************************************
751 0 : SUBROUTINE cp_cfm_to_cfm_triangular(msource, mtarget, uplo)
752 : TYPE(cp_cfm_type), INTENT(IN) :: msource, mtarget
753 : CHARACTER(len=*), INTENT(IN) :: uplo
754 :
755 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_to_cfm_triangular'
756 :
757 0 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa, bb
758 : INTEGER :: handle, ncol, nrow
759 : #if defined(__parallel)
760 : INTEGER, DIMENSION(9) :: desca, descb
761 : #endif
762 :
763 0 : CALL timeset(routineN, handle)
764 :
765 0 : nrow = msource%matrix_struct%nrow_global
766 0 : ncol = msource%matrix_struct%ncol_global
767 :
768 0 : aa => msource%local_data
769 0 : bb => mtarget%local_data
770 :
771 : #if defined(__parallel)
772 0 : desca(:) = msource%matrix_struct%descriptor(:)
773 0 : descb(:) = mtarget%matrix_struct%descriptor(:)
774 0 : CALL pzlacpy(uplo, nrow, ncol, aa(1, 1), 1, 1, desca, bb(1, 1), 1, 1, descb)
775 : #else
776 : CALL zlacpy(uplo, nrow, ncol, aa(1, 1), nrow, bb(1, 1), nrow)
777 : #endif
778 :
779 0 : CALL timestop(handle)
780 0 : END SUBROUTINE cp_cfm_to_cfm_triangular
781 :
782 : ! **************************************************************************************************
783 : !> \brief Copy real and imaginary parts of a complex full matrix into
784 : !> separate real-value full matrices.
785 : !> \param msource complex matrix
786 : !> \param mtargetr (optional) real part of the source matrix
787 : !> \param mtargeti (optional) imaginary part of the source matrix
788 : !> \note
789 : !> Matrix structures are assumed to be equivalent.
790 : ! **************************************************************************************************
791 350948 : SUBROUTINE cp_cfm_to_fm(msource, mtargetr, mtargeti)
792 :
793 : TYPE(cp_cfm_type), INTENT(IN) :: msource
794 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: mtargetr, mtargeti
795 :
796 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_to_fm'
797 :
798 350948 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: zmat
799 : INTEGER :: handle
800 350948 : REAL(kind=dp), DIMENSION(:, :), POINTER :: imat, rmat
801 :
802 350948 : CALL timeset(routineN, handle)
803 :
804 350948 : zmat => msource%local_data
805 350948 : IF (PRESENT(mtargetr)) THEN
806 340652 : rmat => mtargetr%local_data
807 : IF ((.NOT. cp_fm_struct_equivalent(mtargetr%matrix_struct, msource%matrix_struct)) .OR. &
808 340652 : (SIZE(rmat, 1) /= SIZE(zmat, 1)) .OR. &
809 : (SIZE(rmat, 2) /= SIZE(zmat, 2))) THEN
810 0 : CPABORT("size of local_data of mtargetr differ to msource")
811 : END IF
812 : ! copy local data
813 466684770 : rmat = REAL(zmat, kind=dp)
814 : ELSE
815 350948 : NULLIFY (rmat)
816 : END IF
817 350948 : IF (PRESENT(mtargeti)) THEN
818 283444 : imat => mtargeti%local_data
819 : IF ((.NOT. cp_fm_struct_equivalent(mtargeti%matrix_struct, msource%matrix_struct)) .OR. &
820 283444 : (SIZE(imat, 1) /= SIZE(zmat, 1)) .OR. &
821 : (SIZE(imat, 2) /= SIZE(zmat, 2))) THEN
822 0 : CPABORT("size of local_data of mtargeti differ to msource")
823 : END IF
824 : ! copy local data
825 464586266 : imat = REAL(AIMAG(zmat), kind=dp)
826 : ELSE
827 350948 : NULLIFY (imat)
828 : END IF
829 :
830 350948 : CALL timestop(handle)
831 :
832 350948 : END SUBROUTINE cp_cfm_to_fm
833 :
834 : ! **************************************************************************************************
835 : !> \brief Construct a complex full matrix by taking its real and imaginary parts from
836 : !> two separate real-value full matrices.
837 : !> \param msourcer (optional) real part of the complex matrix (defaults to 0.0)
838 : !> \param msourcei (optional) imaginary part of the complex matrix (defaults to 0.0)
839 : !> \param mtarget resulting complex matrix
840 : !> \note
841 : !> Matrix structures are assumed to be equivalent.
842 : ! **************************************************************************************************
843 366620 : SUBROUTINE cp_fm_to_cfm(msourcer, msourcei, mtarget)
844 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: msourcer, msourcei
845 : TYPE(cp_cfm_type), INTENT(IN) :: mtarget
846 :
847 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_to_cfm'
848 :
849 366620 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: zmat
850 : INTEGER :: handle, mode
851 366620 : REAL(kind=dp), DIMENSION(:, :), POINTER :: imat, rmat
852 :
853 366620 : CALL timeset(routineN, handle)
854 :
855 366620 : mode = 0
856 366620 : zmat => mtarget%local_data
857 366620 : IF (PRESENT(msourcer)) THEN
858 360996 : rmat => msourcer%local_data
859 : IF ((.NOT. cp_fm_struct_equivalent(msourcer%matrix_struct, mtarget%matrix_struct)) .OR. &
860 360996 : (SIZE(rmat, 1) /= SIZE(zmat, 1)) .OR. &
861 : (SIZE(rmat, 2) /= SIZE(zmat, 2))) THEN
862 0 : CPABORT("size of local_data of msourcer differ to mtarget")
863 : END IF
864 : mode = mode + 1
865 : ELSE
866 : NULLIFY (rmat)
867 : END IF
868 366620 : IF (PRESENT(msourcei)) THEN
869 286116 : imat => msourcei%local_data
870 : IF ((.NOT. cp_fm_struct_equivalent(msourcei%matrix_struct, mtarget%matrix_struct)) .OR. &
871 286116 : (SIZE(imat, 1) /= SIZE(zmat, 1)) .OR. &
872 : (SIZE(imat, 2) /= SIZE(zmat, 2))) THEN
873 0 : CPABORT("size of local_data of msourcei differ to mtarget")
874 : END IF
875 286116 : mode = mode + 2
876 : ELSE
877 : NULLIFY (imat)
878 : END IF
879 : ! copy local data
880 : SELECT CASE (mode)
881 : CASE (0)
882 0 : zmat(:, :) = z_zero
883 : CASE (1)
884 11661060 : zmat(:, :) = CMPLX(rmat(:, :), 0.0_dp, kind=dp)
885 : CASE (2)
886 138736 : zmat(:, :) = CMPLX(0.0_dp, imat(:, :), kind=dp)
887 : CASE (3)
888 452617508 : zmat(:, :) = CMPLX(rmat(:, :), imat(:, :), kind=dp)
889 : END SELECT
890 :
891 366620 : CALL timestop(handle)
892 :
893 366620 : END SUBROUTINE cp_fm_to_cfm
894 :
895 : ! **************************************************************************************************
896 : !> \brief Initiate the copy operation: get distribution data, post MPI isend and irecvs.
897 : !> \param source input complex-valued fm matrix
898 : !> \param destination output complex-valued fm matrix
899 : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
900 : !> of the input and output matrices
901 : !> \param info all of the data that will be needed to complete the copy operation
902 : !> \note a slightly modified version of the subroutine cp_fm_start_copy_general() that uses
903 : !> allocatable arrays instead of pointers wherever possible.
904 : ! **************************************************************************************************
905 123920 : SUBROUTINE cp_cfm_start_copy_general(source, destination, para_env, info)
906 : TYPE(cp_cfm_type), INTENT(IN) :: source, destination
907 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
908 : TYPE(copy_cfm_info_type), INTENT(out) :: info
909 :
910 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_start_copy_general'
911 :
912 : INTEGER :: dest_p_i, dest_q_j, global_rank, global_size, handle, i, j, k, mpi_rank, &
913 : ncol_block_dest, ncol_block_src, ncol_local_recv, ncol_local_send, ncols, &
914 : nrow_block_dest, nrow_block_src, nrow_local_recv, nrow_local_send, nrows, p, q, &
915 : recv_rank, recv_size, send_rank, send_size
916 12392 : INTEGER, ALLOCATABLE, DIMENSION(:) :: all_ranks, dest2global, dest_p, dest_q, &
917 24784 : recv_count, send_count, send_disp, &
918 12392 : source2global, src_p, src_q
919 12392 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dest_blacs2mpi
920 : INTEGER, DIMENSION(2) :: dest_block, dest_block_tmp, dest_num_pe, &
921 : src_block, src_block_tmp, src_num_pe
922 24784 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices, &
923 24784 : send_col_indices, send_row_indices
924 : TYPE(cp_fm_struct_type), POINTER :: recv_dist, send_dist
925 173488 : TYPE(mp_request_type), DIMENSION(6) :: recv_req, send_req
926 :
927 12392 : CALL timeset(routineN, handle)
928 :
929 : IF (.NOT. cp2k_is_parallel) THEN
930 : ! Just copy all of the matrix data into a 'send buffer', to be unpacked later
931 : nrow_local_send = SIZE(source%local_data, 1)
932 : ncol_local_send = SIZE(source%local_data, 2)
933 : ALLOCATE (info%send_buf(nrow_local_send*ncol_local_send))
934 : k = 0
935 : DO j = 1, ncol_local_send
936 : DO i = 1, nrow_local_send
937 : k = k + 1
938 : info%send_buf(k) = source%local_data(i, j)
939 : END DO
940 : END DO
941 : ELSE
942 12392 : NULLIFY (recv_dist, send_dist)
943 12392 : NULLIFY (recv_col_indices, recv_row_indices, send_col_indices, send_row_indices)
944 :
945 : ! The 'global' communicator contains both the source and destination decompositions
946 12392 : global_size = para_env%num_pe
947 12392 : global_rank = para_env%mepos
948 :
949 : ! The source/send decomposition and destination/recv decompositions may only exist on
950 : ! on a subset of the processes involved in the communication
951 : ! Check if the source and/or destination arguments are .not. ASSOCIATED():
952 : ! if so, skip the send / recv parts (since these processes do not participate in the sending/receiving distribution)
953 12392 : IF (ASSOCIATED(destination%matrix_struct)) THEN
954 12392 : recv_dist => destination%matrix_struct
955 12392 : recv_rank = recv_dist%para_env%mepos
956 : ELSE
957 0 : recv_rank = mp_proc_null
958 : END IF
959 :
960 12392 : IF (ASSOCIATED(source%matrix_struct)) THEN
961 6196 : send_dist => source%matrix_struct
962 6196 : send_rank = send_dist%para_env%mepos
963 : ELSE
964 6196 : send_rank = mp_proc_null
965 : END IF
966 :
967 : ! Map the rank in the source/dest communicator to the global rank
968 37176 : ALLOCATE (all_ranks(0:global_size - 1))
969 :
970 12392 : CALL para_env%allgather(send_rank, all_ranks)
971 12392 : IF (ASSOCIATED(destination%matrix_struct)) THEN
972 61960 : ALLOCATE (source2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
973 37176 : DO i = 0, global_size - 1
974 37176 : IF (all_ranks(i) /= mp_proc_null) THEN
975 12392 : source2global(all_ranks(i)) = i
976 : END IF
977 : END DO
978 : END IF
979 :
980 12392 : CALL para_env%allgather(recv_rank, all_ranks)
981 12392 : IF (ASSOCIATED(source%matrix_struct)) THEN
982 30980 : ALLOCATE (dest2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
983 18588 : DO i = 0, global_size - 1
984 18588 : IF (all_ranks(i) /= mp_proc_null) THEN
985 12392 : dest2global(all_ranks(i)) = i
986 : END IF
987 : END DO
988 : END IF
989 12392 : DEALLOCATE (all_ranks)
990 :
991 : ! Some data from the two decompositions will be needed by all processes in the global group :
992 : ! process grid shape, block size, and the BLACS-to-MPI mapping
993 :
994 : ! The global root process will receive the data (from the root process in each decomposition)
995 86744 : send_req(:) = mp_request_null
996 12392 : IF (global_rank == 0) THEN
997 43372 : recv_req(:) = mp_request_null
998 6196 : CALL para_env%irecv(src_block, mp_any_source, recv_req(1), tag=src_tag)
999 6196 : CALL para_env%irecv(dest_block, mp_any_source, recv_req(2), tag=dest_tag)
1000 6196 : CALL para_env%irecv(src_num_pe, mp_any_source, recv_req(3), tag=src_tag)
1001 6196 : CALL para_env%irecv(dest_num_pe, mp_any_source, recv_req(4), tag=dest_tag)
1002 : END IF
1003 :
1004 12392 : IF (ASSOCIATED(source%matrix_struct)) THEN
1005 6196 : IF ((send_rank == 0)) THEN
1006 : ! need to use separate buffers here in case this is actually global rank 0
1007 18588 : src_block_tmp = [send_dist%nrow_block, send_dist%ncol_block]
1008 6196 : CALL para_env%isend(src_block_tmp, 0, send_req(1), tag=src_tag)
1009 6196 : CALL para_env%isend(send_dist%context%num_pe, 0, send_req(2), tag=src_tag)
1010 : END IF
1011 : END IF
1012 :
1013 12392 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1014 12392 : IF ((recv_rank == 0)) THEN
1015 18588 : dest_block_tmp = [recv_dist%nrow_block, recv_dist%ncol_block]
1016 6196 : CALL para_env%isend(dest_block_tmp, 0, send_req(3), tag=dest_tag)
1017 6196 : CALL para_env%isend(recv_dist%context%num_pe, 0, send_req(4), tag=dest_tag)
1018 : END IF
1019 : END IF
1020 :
1021 12392 : IF (global_rank == 0) THEN
1022 6196 : CALL mp_waitall(recv_req(1:4))
1023 : ! Now we know the process decomposition, we can allocate the arrays to hold the blacs2mpi mapping
1024 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1025 43372 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1))
1026 6196 : CALL para_env%irecv(info%src_blacs2mpi, mp_any_source, recv_req(5), tag=src_tag)
1027 6196 : CALL para_env%irecv(dest_blacs2mpi, mp_any_source, recv_req(6), tag=dest_tag)
1028 : END IF
1029 :
1030 12392 : IF (ASSOCIATED(source%matrix_struct)) THEN
1031 6196 : IF ((send_rank == 0)) THEN
1032 6196 : CALL para_env%isend(send_dist%context%blacs2mpi(:, :), 0, send_req(5), tag=src_tag)
1033 : END IF
1034 : END IF
1035 :
1036 12392 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1037 12392 : IF ((recv_rank == 0)) THEN
1038 6196 : CALL para_env%isend(recv_dist%context%blacs2mpi(:, :), 0, send_req(6), tag=dest_tag)
1039 : END IF
1040 : END IF
1041 :
1042 12392 : IF (global_rank == 0) THEN
1043 6196 : CALL mp_waitall(recv_req(5:6))
1044 : END IF
1045 :
1046 : ! Finally, broadcast the data to all processes in the global communicator
1047 12392 : CALL para_env%bcast(src_block, 0)
1048 12392 : CALL para_env%bcast(dest_block, 0)
1049 12392 : CALL para_env%bcast(src_num_pe, 0)
1050 12392 : CALL para_env%bcast(dest_num_pe, 0)
1051 37176 : info%src_num_pe(1:2) = src_num_pe(1:2)
1052 37176 : info%nblock_src(1:2) = src_block(1:2)
1053 12392 : IF (global_rank /= 0) THEN
1054 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1055 43372 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1))
1056 : END IF
1057 12392 : CALL para_env%bcast(info%src_blacs2mpi, 0)
1058 12392 : CALL para_env%bcast(dest_blacs2mpi, 0)
1059 :
1060 12392 : recv_size = dest_num_pe(1)*dest_num_pe(2)
1061 12392 : send_size = src_num_pe(1)*src_num_pe(2)
1062 12392 : info%send_size = send_size
1063 12392 : CALL mp_waitall(send_req(:))
1064 :
1065 : ! Setup is now complete, we can start the actual communication here.
1066 : ! The order implemented here is:
1067 : ! DEST_1
1068 : ! compute recv sizes
1069 : ! call irecv
1070 : ! SRC_1
1071 : ! compute send sizes
1072 : ! pack send buffers
1073 : ! call isend
1074 : ! DEST_2
1075 : ! wait for the recvs and unpack buffers (this part eventually will go into another
1076 : ! routine to allow comms to run concurrently)
1077 : ! SRC_2
1078 : ! wait for the sends
1079 :
1080 : ! DEST_1
1081 12392 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1082 : CALL cp_fm_struct_get(recv_dist, row_indices=recv_row_indices, &
1083 12392 : col_indices=recv_col_indices)
1084 12392 : info%recv_col_indices => recv_col_indices
1085 12392 : info%recv_row_indices => recv_row_indices
1086 12392 : nrow_block_src = src_block(1)
1087 12392 : ncol_block_src = src_block(2)
1088 74352 : ALLOCATE (recv_count(0:send_size - 1), info%recv_disp(0:send_size - 1), info%recv_request(0:send_size - 1))
1089 :
1090 : ! Determine the recv counts, allocate the receive buffers, call mpi_irecv for all the non-zero sized receives
1091 12392 : nrow_local_recv = recv_dist%nrow_locals(recv_dist%context%mepos(1))
1092 12392 : ncol_local_recv = recv_dist%ncol_locals(recv_dist%context%mepos(2))
1093 12392 : info%nlocal_recv(1) = nrow_local_recv
1094 12392 : info%nlocal_recv(2) = ncol_local_recv
1095 : ! Initialise src_p, src_q arrays (sized using number of rows/cols in the receiving distribution)
1096 61960 : ALLOCATE (src_p(nrow_local_recv), src_q(ncol_local_recv))
1097 76456 : DO i = 1, nrow_local_recv
1098 : ! For each local row we will receive, we look up its global row (in recv_row_indices),
1099 : ! then work out which row block it comes from, and which process row that row block comes from.
1100 76456 : src_p(i) = MOD(((recv_row_indices(i) - 1)/nrow_block_src), src_num_pe(1))
1101 : END DO
1102 140520 : DO j = 1, ncol_local_recv
1103 : ! Similarly for the columns
1104 140520 : src_q(j) = MOD(((recv_col_indices(j) - 1)/ncol_block_src), src_num_pe(2))
1105 : END DO
1106 : ! src_p/q now contains the process row/column ID that will send data to that row/column
1107 :
1108 24784 : DO q = 0, src_num_pe(2) - 1
1109 140520 : ncols = COUNT(src_q == q)
1110 37176 : DO p = 0, src_num_pe(1) - 1
1111 76456 : nrows = COUNT(src_p == p)
1112 : ! Use the send_dist here as we are looking up the processes where the data comes from
1113 24784 : recv_count(info%src_blacs2mpi(p, q)) = nrows*ncols
1114 : END DO
1115 : END DO
1116 12392 : DEALLOCATE (src_p, src_q)
1117 :
1118 : ! Use one long buffer (and displacements into that buffer)
1119 : ! this prevents the need for a rectangular array where not all elements will be populated
1120 49568 : ALLOCATE (info%recv_buf(SUM(recv_count(:))))
1121 12392 : info%recv_disp(0) = 0
1122 12392 : DO i = 1, send_size - 1
1123 12392 : info%recv_disp(i) = info%recv_disp(i - 1) + recv_count(i - 1)
1124 : END DO
1125 :
1126 : ! Issue receive calls on ranks which expect data
1127 24784 : DO k = 0, send_size - 1
1128 24784 : IF (recv_count(k) > 0) THEN
1129 : CALL para_env%irecv(info%recv_buf(info%recv_disp(k) + 1:info%recv_disp(k) + recv_count(k)), &
1130 12392 : source2global(k), info%recv_request(k))
1131 : ELSE
1132 0 : info%recv_request(k) = mp_request_null
1133 : END IF
1134 : END DO
1135 12392 : DEALLOCATE (source2global)
1136 : END IF ! ASSOCIATED(destination)
1137 :
1138 : ! SRC_1
1139 12392 : IF (ASSOCIATED(source%matrix_struct)) THEN
1140 : CALL cp_fm_struct_get(send_dist, row_indices=send_row_indices, &
1141 6196 : col_indices=send_col_indices)
1142 6196 : nrow_block_dest = dest_block(1)
1143 6196 : ncol_block_dest = dest_block(2)
1144 43372 : ALLOCATE (send_count(0:recv_size - 1), send_disp(0:recv_size - 1), info%send_request(0:recv_size - 1))
1145 :
1146 : ! Determine the send counts, allocate the send buffers
1147 6196 : nrow_local_send = send_dist%nrow_locals(send_dist%context%mepos(1))
1148 6196 : ncol_local_send = send_dist%ncol_locals(send_dist%context%mepos(2))
1149 :
1150 : ! Initialise dest_p, dest_q arrays (sized nrow_local, ncol_local)
1151 : ! i.e. number of rows,cols in the sending distribution
1152 30980 : ALLOCATE (dest_p(nrow_local_send), dest_q(ncol_local_send))
1153 :
1154 70260 : DO i = 1, nrow_local_send
1155 : ! Use the send_dist%row_indices() here (we are looping over the local rows we will send)
1156 70260 : dest_p(i) = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1157 : END DO
1158 70260 : DO j = 1, ncol_local_send
1159 70260 : dest_q(j) = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1160 : END DO
1161 : ! dest_p/q now contain the process row/column ID that will receive data from that row/column
1162 :
1163 12392 : DO q = 0, dest_num_pe(2) - 1
1164 70260 : ncols = COUNT(dest_q == q)
1165 24784 : DO p = 0, dest_num_pe(1) - 1
1166 140520 : nrows = COUNT(dest_p == p)
1167 18588 : send_count(dest_blacs2mpi(p, q)) = nrows*ncols
1168 : END DO
1169 : END DO
1170 6196 : DEALLOCATE (dest_p, dest_q)
1171 :
1172 : ! Allocate the send buffer using send_count -- and calculate the offset into the buffer for each process
1173 30980 : ALLOCATE (info%send_buf(SUM(send_count(:))))
1174 6196 : send_disp(0) = 0
1175 12392 : DO k = 1, recv_size - 1
1176 12392 : send_disp(k) = send_disp(k - 1) + send_count(k - 1)
1177 : END DO
1178 :
1179 : ! Loop over the smat, pack the send buffers
1180 6196 : send_count(:) = 0
1181 70260 : DO j = 1, ncol_local_send
1182 : ! Use send_col_indices and row_indices here, as we are looking up the global row/column number of local rows.
1183 64064 : dest_q_j = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1184 797876 : DO i = 1, nrow_local_send
1185 727616 : dest_p_i = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1186 727616 : mpi_rank = dest_blacs2mpi(dest_p_i, dest_q_j)
1187 727616 : send_count(mpi_rank) = send_count(mpi_rank) + 1
1188 791680 : info%send_buf(send_disp(mpi_rank) + send_count(mpi_rank)) = source%local_data(i, j)
1189 : END DO
1190 : END DO
1191 :
1192 : ! For each non-zero send_count, call mpi_isend
1193 18588 : DO k = 0, recv_size - 1
1194 18588 : IF (send_count(k) > 0) THEN
1195 : CALL para_env%isend(info%send_buf(send_disp(k) + 1:send_disp(k) + send_count(k)), &
1196 12392 : dest2global(k), info%send_request(k))
1197 : ELSE
1198 0 : info%send_request(k) = mp_request_null
1199 : END IF
1200 : END DO
1201 6196 : DEALLOCATE (send_count, send_disp, dest2global)
1202 : END IF ! ASSOCIATED(source)
1203 12392 : DEALLOCATE (dest_blacs2mpi)
1204 :
1205 : END IF !IF (.NOT. cp2k_is_parallel)
1206 :
1207 12392 : CALL timestop(handle)
1208 :
1209 49568 : END SUBROUTINE cp_cfm_start_copy_general
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief Complete the copy operation: wait for comms, unpack, clean up MPI state.
1213 : !> \param destination output cfm matrix
1214 : !> \param info all of the data that will be needed to complete the copy operation
1215 : !> \note a slightly modified version of the subroutine cp_fm_finish_copy_general() that uses
1216 : !> allocatable arrays instead of pointers wherever possible.
1217 : ! **************************************************************************************************
1218 12392 : SUBROUTINE cp_cfm_finish_copy_general(destination, info)
1219 : TYPE(cp_cfm_type), INTENT(IN) :: destination
1220 : TYPE(copy_cfm_info_type), INTENT(inout) :: info
1221 :
1222 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_finish_copy_general'
1223 :
1224 : INTEGER :: handle, i, j, k, mpi_rank, ni, nj, &
1225 : src_q_j
1226 12392 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_count, src_p_i
1227 12392 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices
1228 :
1229 12392 : CALL timeset(routineN, handle)
1230 :
1231 : IF (.NOT. cp2k_is_parallel) THEN
1232 : ! Now unpack the data from the 'send buffer'
1233 : k = 0
1234 : DO j = 1, SIZE(destination%local_data, 2)
1235 : DO i = 1, SIZE(destination%local_data, 1)
1236 : k = k + 1
1237 : destination%local_data(i, j) = info%send_buf(k)
1238 : END DO
1239 : END DO
1240 : DEALLOCATE (info%send_buf)
1241 : ELSE
1242 : ! Set up local variables ...
1243 12392 : recv_col_indices => info%recv_col_indices
1244 12392 : recv_row_indices => info%recv_row_indices
1245 :
1246 : ! ... use the local variables to do the work
1247 : ! DEST_2
1248 12392 : CALL mp_waitall(info%recv_request(:))
1249 :
1250 12392 : nj = info%nlocal_recv(2)
1251 12392 : ni = info%nlocal_recv(1)
1252 61960 : ALLOCATE (recv_count(0:info%send_size - 1), src_p_i(ni))
1253 : ! Loop over the rmat, filling it in with data from the recv buffers
1254 : ! (here the block sizes, num_pes refer to the distribution of the source matrix)
1255 12392 : recv_count(:) = 0
1256 76456 : DO i = 1, ni
1257 76456 : src_p_i(i) = MOD(((recv_row_indices(i) - 1)/info%nblock_src(1)), info%src_num_pe(1))
1258 : END DO
1259 :
1260 140520 : DO j = 1, nj
1261 128128 : src_q_j = MOD(((recv_col_indices(j) - 1)/info%nblock_src(2)), info%src_num_pe(2))
1262 868136 : DO i = 1, ni
1263 727616 : mpi_rank = info%src_blacs2mpi(src_p_i(i), src_q_j)
1264 727616 : recv_count(mpi_rank) = recv_count(mpi_rank) + 1
1265 855744 : destination%local_data(i, j) = info%recv_buf(info%recv_disp(mpi_rank) + recv_count(mpi_rank))
1266 : END DO
1267 : END DO
1268 :
1269 12392 : DEALLOCATE (recv_count, src_p_i)
1270 : ! Invalidate the stored state
1271 12392 : NULLIFY (info%recv_col_indices, info%recv_row_indices)
1272 12392 : DEALLOCATE (info%recv_disp, info%recv_request, info%recv_buf, info%src_blacs2mpi)
1273 : END IF
1274 :
1275 12392 : CALL timestop(handle)
1276 :
1277 12392 : END SUBROUTINE cp_cfm_finish_copy_general
1278 :
1279 : ! **************************************************************************************************
1280 : !> \brief Complete the copy operation: wait for comms clean up MPI state.
1281 : !> \param info all of the data that will be needed to complete the copy operation
1282 : !> \note a slightly modified version of the subroutine cp_fm_cleanup_copy_general() that uses
1283 : !> allocatable arrays instead of pointers wherever possible.
1284 : ! **************************************************************************************************
1285 6196 : SUBROUTINE cp_cfm_cleanup_copy_general(info)
1286 : TYPE(copy_cfm_info_type), INTENT(inout) :: info
1287 :
1288 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cleanup_copy_general'
1289 :
1290 : INTEGER :: handle
1291 :
1292 6196 : CALL timeset(routineN, handle)
1293 :
1294 : IF (cp2k_is_parallel) THEN
1295 : ! SRC_2
1296 : ! If this process is also in the destination decomposition, this deallocate
1297 : ! Was already done in cp_fm_finish_copy_general
1298 6196 : IF (ALLOCATED(info%src_blacs2mpi)) DEALLOCATE (info%src_blacs2mpi)
1299 6196 : CALL mp_waitall(info%send_request(:))
1300 6196 : DEALLOCATE (info%send_request, info%send_buf)
1301 : END IF
1302 :
1303 6196 : CALL timestop(handle)
1304 6196 : END SUBROUTINE cp_cfm_cleanup_copy_general
1305 0 : END MODULE cp_cfm_types
|
