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