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 Routines to convert sparse matrices between DBCSR (distributed-blocks compressed sparse rows)
10 : !> and SIESTA (distributed compressed sparse columns) formats.
11 : !> \author Sergey Chulkov
12 : !> \author Christian Ahart
13 : !> \author Clotilde Cucinotta
14 : ! **************************************************************************************************
15 : MODULE smeagol_matrix_utils
16 : USE cell_types, ONLY: cell_type, &
17 : real_to_scaled, &
18 : scaled_to_real
19 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p, &
20 : dbcsr_get_info, &
21 : dbcsr_p_type, &
22 : dbcsr_set
23 : USE kinds, ONLY: dp, &
24 : dp_size, &
25 : int_8
26 : USE message_passing, ONLY: mp_para_env_type, &
27 : mp_request_type, &
28 : mp_waitall
29 : USE negf_matrix_utils, ONLY: get_index_by_cell
30 : #if defined(__SMEAGOL)
31 : USE parallel, ONLY: GetNodeOrbs, &
32 : GlobalToLocalOrb, &
33 : LocalToGlobalOrb, &
34 : WhichNodeOrb
35 : #endif
36 : USE particle_types, ONLY: particle_type
37 : USE qs_neighbor_list_types, ONLY: get_iterator_info, &
38 : get_neighbor_list_set_p, &
39 : neighbor_list_iterate, &
40 : neighbor_list_iterator_create, &
41 : neighbor_list_iterator_p_type, &
42 : neighbor_list_iterator_release, &
43 : neighbor_list_set_p_type
44 : USE qs_subsys_types, ONLY: qs_subsys_get, &
45 : qs_subsys_type
46 : USE util, ONLY: sort
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 : PRIVATE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'smeagol_matrix_utils'
53 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
54 :
55 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_iatom_index = 1
56 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_jatom_index = 2
57 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_dbcsr_image_index = 3
58 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_siesta_image_index = 4
59 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_siesta_transp_image_index = 5
60 : INTEGER, PARAMETER, PRIVATE :: neighbor_list_dim1 = neighbor_list_siesta_transp_image_index
61 :
62 : PUBLIC :: siesta_distrib_csc_struct_type
63 : PUBLIC :: siesta_struct_create, siesta_struct_release
64 : PUBLIC :: convert_dbcsr_to_distributed_siesta, convert_distributed_siesta_to_dbcsr
65 :
66 : PRIVATE :: get_negf_cell_ijk, index_in_canonical_enumeration, number_from_canonical_enumeration, pbc_0_1
67 : PRIVATE :: get_number_of_mpi_sendrecv_requests, assign_nonzero_elements_to_requests
68 :
69 : !> number of DBCSR matrix elements to receive from a given rank
70 : INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_recv = 1
71 : !> number of DBCSR matrix elements to send to a given rank
72 : INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_send = 2
73 : INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_dim2 = nelements_dbcsr_send
74 :
75 : ! 128 MiB (to limit memory usage for matrix redistribution)
76 : INTEGER(kind=int_8), PARAMETER, PRIVATE :: max_mpi_packet_size_bytes = 134217728
77 : INTEGER(kind=int_8), PARAMETER, PRIVATE :: max_mpi_packet_size_dp = max_mpi_packet_size_bytes/INT(dp_size, kind=int_8)
78 :
79 : ! a portable way to determine the upper bound for tag value is to call
80 : ! MPI_COMM_GET_ATTR(comm, MPI_TAG_UB, max_mpi_rank, flag, ierror).
81 : ! The MPI specification guarantees a value of 32767
82 : INTEGER, PARAMETER, PRIVATE :: max_mpi_rank = 32767
83 :
84 : ! **************************************************************************************************
85 : !> \brief Sparsity pattern of replicated SIESTA compressed sparse column (CSC) matrices
86 : ! **************************************************************************************************
87 : TYPE siesta_distrib_csc_struct_type
88 : !> gather all non-zero matrix elements on the given MPI process.
89 : !> Distribute the elements across MPI processes if gather_root < 0.
90 : !> The matrix elements should be located on I/O process in case of bulk transport,
91 : !> and should be distributed in case of SMEAGOL calculation.
92 : INTEGER :: gather_root = 0
93 : !> Based of full (.FALSE.) or upper-triangular (.TRUE.) DBCSR matrix.
94 : !> It is used in CPASSERTs to allow access to lower-triangular matrix elements of non-symmetric DBCSR matrices
95 : !> In case there is no bugs in this module, these CPASSERTs should newer trigger.
96 : !> Therefore the 'symmetric' variable alongside with relevant CPASSERT calls are excessive and can be removed.
97 : LOGICAL :: symmetric = .TRUE.
98 : !> number of neighbour list nodes for each MPI process (0:num_pe-1).
99 : !> If do_merge == .TRUE., nodes for different cell images along k cell vector are merged into one node
100 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nnodes_per_proc
101 :
102 : !> replicated neighbour list (1:neighbor_list_dim1, 1:SUM(nnodes_per_proc)).
103 : !> Neighbour list nodes are ordered according to their MPI ranks.
104 : !> Thus, the first nnodes_per_proc(0) nodes are stored on MPI rank 0,
105 : !> the next nnodes_per_proc(1) nodes reside on MPI rank 1, etc
106 : !> Nodes for cell images along transport direction are merged into one node.
107 : !> The number of non-zero DBCSR matrix blocks and their DBCSR cell image indices
108 : !> are stored into 'n_dbcsr_cell_images_to_merge' and 'dbcsr_cell_image_to_merge' arrays
109 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl_repl
110 :
111 : !> number of DBCSR images for each local merged neighbour list node (1:nnodes_per_proc(para_env%mepos))
112 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_dbcsr_cell_images_to_merge
113 : !> list of DBCSR image indices to merge; (1:SUM(n_dbcsr_cell_images_to_merge))
114 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dbcsr_cell_image_to_merge
115 :
116 : !> number of DBCSR non-zero matrix elements that should be received/sent from each MPI rank
117 : !> (0:num_pe-1, 1:nelements_dbcsr_dim2)
118 : INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:, :) :: nelements_per_proc
119 :
120 : !> number of non-zero matrix elements local to this MPI rank.
121 : INTEGER(kind=int_8) :: n_nonzero_elements = 0_int_8
122 : INTEGER :: nrows = 0, ncols = 0
123 : !> Number of non-zero matrix elements (columns) on each row
124 : !> n_nonzero_cols(1:nrows); same as 'numh' in SMEAGOL code.
125 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_nonzero_cols
126 : !> offset of the first non-zero matrix elements on each row.
127 : !> column_offset(1:nrows); same as 'listhptr' in SMEAGOL code.
128 : !> It should be declared as INTEGER(kind=int_8), but SMEAGOL expects it to be INTEGER
129 : INTEGER, ALLOCATABLE, DIMENSION(:) :: row_offset
130 : !> column index of each non-zero matrix element.
131 : !> col_index(1:n_nonzero_elements); same as 'listh' in SMEAGOL code
132 : !> col index of the first non-zero matrix elements of irow row is col_index(row_offset(irow)+1)
133 : INTEGER, ALLOCATABLE, DIMENSION(:) :: col_index
134 : !> index of the non-zero matrix element in a communication buffer for each rank
135 : INTEGER, ALLOCATABLE, DIMENSION(:) :: packed_index
136 : !> R_atom_row - R_atom_col
137 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xij
138 : !> equivalent atomic orbitals
139 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indxuo
140 : !> atomic index on which the orbital is centred
141 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iaorb
142 : !> coordinates of all atoms in the supercell
143 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xa
144 : END TYPE siesta_distrib_csc_struct_type
145 :
146 : CONTAINS
147 :
148 : ! **************************************************************************************************
149 : !> \brief Map non-zero matrix blocks between sparse matrices in DBCSR and SIESTA formats.
150 : !> \param siesta_struct structure that stores metadata (sparsity pattern) of sparse SIESTA matrices
151 : !> \param matrix_dbcsr_kp DBCSR matrices for each cell image
152 : !> \param subsys QuickStep molecular system
153 : !> \param cell_to_index array to convert 3-D cell indices to 1-D DBCSR image indices
154 : !> \param sab_nl pair-wise neighbour list
155 : !> \param para_env MPI parallel environment
156 : !> \param max_ij_cell_image largest index of cell images along i and j cell vectors (e.g. (2,0) in
157 : !> case of 5 cell images (0,0), (1,0), (-1,0), (2,0), and (-2,0))
158 : !> \param do_merge merge DBCSR images along transport direction (k cell vector)
159 : !> \param gather_root distribute non-zero matrix elements of SIESTA matrices across all
160 : !> parallel processes (-1), or gather them on the given MPI rank (>= 0).
161 : ! **************************************************************************************************
162 4 : SUBROUTINE siesta_struct_create(siesta_struct, matrix_dbcsr_kp, subsys, cell_to_index, &
163 : sab_nl, para_env, max_ij_cell_image, do_merge, gather_root)
164 : TYPE(siesta_distrib_csc_struct_type), &
165 : INTENT(inout) :: siesta_struct
166 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
167 : TYPE(qs_subsys_type), POINTER :: subsys
168 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
169 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
170 : POINTER :: sab_nl
171 : TYPE(mp_para_env_type), POINTER :: para_env
172 : INTEGER, DIMENSION(2), INTENT(inout) :: max_ij_cell_image
173 : LOGICAL, INTENT(in) :: do_merge
174 : INTEGER, INTENT(in) :: gather_root
175 :
176 : CHARACTER(len=*), PARAMETER :: routineN = 'siesta_struct_create'
177 :
178 : CHARACTER(len=20) :: str_nelem, str_nelem_max
179 : INTEGER :: handle, iatom, icol, icol_blk, icol_local, image, image_j, image_k, irow, &
180 : irow_local, natoms, ncells_siesta_total, ncols_blk, ncols_total, nrows_local, &
181 : nrows_total, offset
182 : INTEGER(kind=int_8) :: n_nonzero_elements_local
183 : INTEGER, DIMENSION(3) :: max_ijk_cell_image, ncells_siesta
184 4 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size
185 : LOGICAL :: do_distribute, is_root_rank
186 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: particle_coords
187 : REAL(kind=dp), DIMENSION(3) :: real_cell_shift, scaled_cell_shift
188 : TYPE(cell_type), POINTER :: cell
189 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
190 :
191 4 : CALL timeset(routineN, handle)
192 4 : do_distribute = gather_root < 0
193 4 : is_root_rank = gather_root == para_env%mepos
194 :
195 : ! here row_blk_offset / col_blk_offset are global indices of the first row / column of a given non-zero block.
196 : ! They are not offsets (index-1) but the actual indices.
197 : CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
198 : nfullrows_total=nrows_total, nfullcols_total=ncols_total, &
199 4 : nblkcols_total=ncols_blk, col_blk_size=col_blk_size, col_blk_offset=col_blk_offset)
200 : IF (debug_this_module) THEN
201 : CPASSERT(nrows_total == ncols_total)
202 : CPASSERT(gather_root < para_env%num_pe)
203 : END IF
204 :
205 4 : siesta_struct%gather_root = gather_root
206 4 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=siesta_struct%symmetric)
207 :
208 : ! apply periodic boundary conditions to atomic coordinates
209 4 : CALL qs_subsys_get(subsys, cell=cell, particle_set=particle_set, nparticle=natoms)
210 12 : ALLOCATE (particle_coords(3, natoms))
211 140 : DO iatom = 1, natoms
212 140 : CALL pbc_0_1(particle_coords(1:3, iatom), particle_set(iatom)%r(1:3), cell)
213 : END DO
214 :
215 : ! Note: in case we would like to limit the number of cell images along transport direction (k cell vector)
216 : ! by enabling 'BulkTransvCellSizeZ' keyword, we need to pass max_ijk_cell_image(1:3) vector
217 : ! to the subroutine instead of the reduced vector max_ij_cell_image(1:2)
218 12 : max_ijk_cell_image(1:2) = max_ij_cell_image(1:2)
219 :
220 : ! determine the actual number of cell images along k cell vector. Unless the third element is also passed
221 : ! via subroutine arguments, an extra MPI_Allreduce operation is needed each time we call
222 : ! replicate_neighbour_list() / get_nnodes_local().
223 4 : max_ijk_cell_image(3) = -1
224 : ! bulk-transport calculation expects exactly 3 cell images along transport direction
225 4 : IF (.NOT. do_merge) max_ijk_cell_image(3) = 1
226 :
227 : ! replicate pair-wise neighbour list. Identical non-zero matrix blocks from cell image along transport direction
228 : ! are grouped together if do_merge == .TRUE.
229 12 : ALLOCATE (siesta_struct%nnodes_per_proc(0:para_env%num_pe - 1))
230 : CALL replicate_neighbour_list(siesta_struct%nl_repl, &
231 : siesta_struct%n_dbcsr_cell_images_to_merge, &
232 : siesta_struct%dbcsr_cell_image_to_merge, &
233 : siesta_struct%nnodes_per_proc, &
234 4 : max_ijk_cell_image, sab_nl, para_env, particle_coords, cell, cell_to_index, do_merge)
235 12 : max_ij_cell_image(1:2) = max_ijk_cell_image(1:2)
236 :
237 : ! count number of non-zero matrix elements that need to be send to and received from other parallel processes
238 16 : ALLOCATE (siesta_struct%nelements_per_proc(0:para_env%num_pe - 1, nelements_dbcsr_dim2))
239 : CALL count_remote_dbcsr_elements(siesta_struct%nelements_per_proc, siesta_struct%nnodes_per_proc, &
240 4 : siesta_struct%nl_repl, matrix_dbcsr_kp, siesta_struct%symmetric, para_env, gather_root)
241 :
242 : ! number of SIESTA non-zero matrix elements that are going to be stored on this parallel process
243 12 : n_nonzero_elements_local = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
244 4 : siesta_struct%n_nonzero_elements = n_nonzero_elements_local
245 :
246 : ! as SMEAGOL uses 32-bits integers, the number of non-zero matrix elements is limited by 2^31 per MPI rank.
247 : ! Abort CP2K if we are about to exceed this limit.
248 4 : IF (n_nonzero_elements_local > INT(HUGE(0), kind=int_8)) THEN
249 0 : WRITE (str_nelem, '(I0)') n_nonzero_elements_local
250 0 : WRITE (str_nelem_max, '(I0)') HUGE(0)
251 : CALL cp_abort(__LOCATION__, &
252 : "The number of non-zero matrix elements per MPI process "//TRIM(str_nelem)// &
253 : " cannot exceed "//TRIM(str_nelem_max)// &
254 0 : ". Please increase the number of MPI processes to satisfy this SMEAGOL limitation.")
255 : END IF
256 :
257 : ! in case there is no nonzero matrix element stored on this process, allocate arrays with one element to avoid SEGFAULT
258 4 : IF (n_nonzero_elements_local == 0) n_nonzero_elements_local = 1
259 :
260 : ! number of SIESTA-matrix rows local to the given parallel process
261 4 : IF (do_distribute) THEN
262 : #if defined(__SMEAGOL)
263 0 : CALL GetNodeOrbs(nrows_total, para_env%mepos, para_env%num_pe, nrows_local)
264 : #else
265 : CALL cp_abort(__LOCATION__, &
266 : "CP2K was compiled with no SMEAGOL support.")
267 : #endif
268 : ELSE
269 4 : IF (is_root_rank) THEN
270 2 : nrows_local = nrows_total
271 : ELSE
272 2 : nrows_local = 0
273 : END IF
274 : END IF
275 :
276 : ! number of cell images along each cell vector. It is 2*m+1, as SIESTA images are ordered as 0, 1, -1, ..., m, -m
277 16 : ncells_siesta(1:3) = 2*max_ijk_cell_image(1:3) + 1
278 : ! in case of merged cell images along the transport direction, there will be just 1 'merged' cell image along it
279 4 : IF (do_merge) ncells_siesta(3) = 1
280 :
281 4 : ncells_siesta_total = ncells_siesta(1)*ncells_siesta(2)*ncells_siesta(3)
282 :
283 : ! number of rows local to the given parallel process. Rows are distributed in a block-cyclic manner
284 4 : siesta_struct%nrows = nrows_local
285 :
286 : ! number of columns of the matrix in its dense form. SIESTA uses 1-D (rows) block-cyclic distribution.
287 : ! All non-zero matrix elements on a given row are stored on the same parallel process
288 4 : siesta_struct%ncols = nrows_total*ncells_siesta_total
289 :
290 : ! allocate at least one array element to avoid SIGFAULT when passing unallocated arrays to subroutines
291 4 : IF (nrows_local == 0) nrows_local = 1
292 :
293 12 : ALLOCATE (siesta_struct%n_nonzero_cols(nrows_local))
294 8 : ALLOCATE (siesta_struct%row_offset(nrows_local))
295 12 : ALLOCATE (siesta_struct%col_index(n_nonzero_elements_local))
296 8 : ALLOCATE (siesta_struct%packed_index(n_nonzero_elements_local))
297 :
298 : ! restore the actual number of local rows
299 : nrows_local = siesta_struct%nrows
300 :
301 : ! get number of non-zero matrix element on each local row (n_nonzero_cols),
302 : ! offset of the first non-zero matrix element for each local row (row_offset),
303 : ! global column indices of all local non-zero matrix elements (col_index), and
304 : ! the indices of all local non-zero matrix elements in the communication buffer (packed_index)
305 : CALL get_nonzero_element_indices(siesta_struct%n_nonzero_cols, siesta_struct%row_offset, &
306 : siesta_struct%col_index, siesta_struct%packed_index, &
307 : siesta_struct%nl_repl, matrix_dbcsr_kp, &
308 4 : siesta_struct%symmetric, para_env, gather_root)
309 :
310 : ! indices of equivalent atomic orbitals
311 12 : ALLOCATE (siesta_struct%indxuo(siesta_struct%ncols))
312 1228 : DO icol = 1, ncols_total
313 1228 : siesta_struct%indxuo(icol) = icol
314 : END DO
315 300 : DO image = 2, ncells_siesta_total
316 181452 : siesta_struct%indxuo((image - 1)*ncols_total + 1:image*ncols_total) = siesta_struct%indxuo(1:ncols_total)
317 : END DO
318 :
319 : ! particle index on which the orbital is centred
320 8 : ALLOCATE (siesta_struct%iaorb(siesta_struct%ncols))
321 140 : DO icol_blk = 1, ncols_blk
322 : ! col_blk_offset() is not an offset but the index of the first atomic orbital in the column block
323 1364 : siesta_struct%iaorb(col_blk_offset(icol_blk):col_blk_offset(icol_blk) + col_blk_size(icol_blk) - 1) = icol_blk
324 : END DO
325 300 : DO image = 2, ncells_siesta_total
326 181452 : siesta_struct%iaorb((image - 1)*ncols_total + 1:image*ncols_total) = siesta_struct%iaorb(1:ncols_total) + (image - 1)*natoms
327 : END DO
328 :
329 : ! coordinates of all particles in each cell images
330 12 : ALLOCATE (siesta_struct%xa(3, natoms*ncells_siesta_total))
331 16 : DO image_k = 1, ncells_siesta(3)
332 : !icell_siesta(3) = image_k
333 12 : scaled_cell_shift(3) = REAL(number_from_canonical_enumeration(image_k), kind=dp) ! SIESTA -> actual cell index
334 76 : DO image_j = 1, ncells_siesta(2)
335 60 : scaled_cell_shift(2) = REAL(number_from_canonical_enumeration(image_j), kind=dp)
336 372 : DO image = 1, ncells_siesta(1)
337 300 : scaled_cell_shift(1) = REAL(number_from_canonical_enumeration(image), kind=dp)
338 300 : CALL scaled_to_real(real_cell_shift, scaled_cell_shift, cell)
339 300 : offset = (((image_k - 1)*ncells_siesta(2) + image_j - 1)*ncells_siesta(1) + image - 1)*natoms
340 10560 : DO iatom = 1, natoms
341 41100 : siesta_struct%xa(1:3, offset + iatom) = particle_set(iatom)%r(1:3) + real_cell_shift(1:3)
342 : END DO
343 : END DO
344 : END DO
345 : END DO
346 :
347 : ! inter-atomic distance
348 12 : ALLOCATE (siesta_struct%xij(3, n_nonzero_elements_local))
349 616 : DO irow_local = 1, nrows_local
350 612 : IF (do_distribute) THEN
351 : #if defined(__SMEAGOL)
352 0 : CALL LocalToGlobalOrb(irow_local, para_env%mepos, para_env%num_pe, irow)
353 : #else
354 : CALL cp_abort(__LOCATION__, &
355 : "CP2K was compiled with no SMEAGOL support.")
356 : #endif
357 : ELSE
358 612 : irow = irow_local
359 : IF (debug_this_module) THEN
360 : CPASSERT(is_root_rank)
361 : END IF
362 : END IF
363 612 : offset = siesta_struct%row_offset(irow_local)
364 2319484 : DO icol_local = offset + 1, offset + siesta_struct%n_nonzero_cols(irow_local)
365 2318868 : icol = siesta_struct%col_index(icol_local)
366 : siesta_struct%xij(1:3, icol_local) = siesta_struct%xa(1:3, siesta_struct%iaorb(icol)) - &
367 9276084 : siesta_struct%xa(1:3, siesta_struct%iaorb(irow))
368 : END DO
369 : END DO
370 :
371 4 : DEALLOCATE (particle_coords)
372 :
373 4 : CALL timestop(handle)
374 8 : END SUBROUTINE siesta_struct_create
375 :
376 : ! **************************************************************************************************
377 : !> \brief Release a SIESTA matrix structure
378 : !> \param siesta_struct structure to release
379 : ! **************************************************************************************************
380 4 : SUBROUTINE siesta_struct_release(siesta_struct)
381 : TYPE(siesta_distrib_csc_struct_type), &
382 : INTENT(inout) :: siesta_struct
383 :
384 : CHARACTER(len=*), PARAMETER :: routineN = 'siesta_struct_release'
385 :
386 : INTEGER :: handle
387 :
388 4 : CALL timeset(routineN, handle)
389 :
390 4 : siesta_struct%gather_root = -1
391 :
392 4 : IF (ALLOCATED(siesta_struct%nnodes_per_proc)) DEALLOCATE (siesta_struct%nnodes_per_proc)
393 4 : IF (ALLOCATED(siesta_struct%nl_repl)) DEALLOCATE (siesta_struct%nl_repl)
394 4 : IF (ALLOCATED(siesta_struct%n_dbcsr_cell_images_to_merge)) DEALLOCATE (siesta_struct%n_dbcsr_cell_images_to_merge)
395 4 : IF (ALLOCATED(siesta_struct%dbcsr_cell_image_to_merge)) DEALLOCATE (siesta_struct%dbcsr_cell_image_to_merge)
396 4 : IF (ALLOCATED(siesta_struct%nelements_per_proc)) DEALLOCATE (siesta_struct%nelements_per_proc)
397 :
398 4 : siesta_struct%n_nonzero_elements = 0
399 4 : siesta_struct%nrows = 0
400 4 : siesta_struct%ncols = 0
401 :
402 4 : IF (ALLOCATED(siesta_struct%n_nonzero_cols)) DEALLOCATE (siesta_struct%n_nonzero_cols)
403 4 : IF (ALLOCATED(siesta_struct%row_offset)) DEALLOCATE (siesta_struct%row_offset)
404 4 : IF (ALLOCATED(siesta_struct%col_index)) DEALLOCATE (siesta_struct%col_index)
405 4 : IF (ALLOCATED(siesta_struct%packed_index)) DEALLOCATE (siesta_struct%packed_index)
406 :
407 4 : IF (ALLOCATED(siesta_struct%xij)) DEALLOCATE (siesta_struct%xij)
408 4 : IF (ALLOCATED(siesta_struct%indxuo)) DEALLOCATE (siesta_struct%indxuo)
409 4 : IF (ALLOCATED(siesta_struct%iaorb)) DEALLOCATE (siesta_struct%iaorb)
410 4 : IF (ALLOCATED(siesta_struct%xa)) DEALLOCATE (siesta_struct%xa)
411 :
412 4 : CALL timestop(handle)
413 4 : END SUBROUTINE siesta_struct_release
414 :
415 : ! **************************************************************************************************
416 : !> \brief Convert matrix from DBCSR to sparse SIESTA format.
417 : !> \param matrix_siesta matrix in SIESTA format [out]
418 : !> \param matrix_dbcsr_kp DBCSR matrix [in]
419 : !> \param siesta_struct structure to map matrix blocks between formats
420 : !> \param para_env MPI parallel environment
421 : ! **************************************************************************************************
422 16 : SUBROUTINE convert_dbcsr_to_distributed_siesta(matrix_siesta, matrix_dbcsr_kp, siesta_struct, para_env)
423 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: matrix_siesta
424 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
425 : TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
426 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
427 :
428 : CHARACTER(len=*), PARAMETER :: routineN = 'convert_dbcsr_to_distributed_siesta'
429 :
430 : INTEGER :: first_col_minus_one, first_row_minus_one, handle, icol_blk, icol_local, &
431 : image_dbcsr, image_ind, image_ind_offset, image_siesta, image_siesta_transp, inode, &
432 : inode_proc, iproc, irequest, irow_blk, irow_local, irow_proc, mepos, n_image_ind, &
433 : ncols_blk, ncols_local, nnodes_proc, node_offset, nprocs, nrequests_recv, &
434 : nrequests_total, nrows_blk, nrows_local
435 : INTEGER(kind=int_8) :: n_nonzero_elements_dbcsr, &
436 : n_nonzero_elements_siesta, &
437 : offset_recv_mepos, offset_send_mepos
438 16 : INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: n_packed_elements_per_proc, &
439 16 : nelements_per_request, &
440 16 : offset_per_proc, offset_per_request
441 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: next_nonzero_element_offset, peer_rank, &
442 16 : request_tag
443 16 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
444 16 : row_blk_offset, row_blk_size
445 : LOGICAL :: do_distribute, found, is_root_rank, &
446 : symmetric
447 16 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buffer, reorder_recv_buffer, &
448 16 : send_buffer
449 16 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block, sm_block_merged
450 16 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: requests
451 :
452 16 : CALL timeset(routineN, handle)
453 9275496 : matrix_siesta(:) = 0.0_dp
454 :
455 16 : mepos = para_env%mepos
456 16 : nprocs = para_env%num_pe
457 16 : do_distribute = siesta_struct%gather_root < 0
458 16 : is_root_rank = siesta_struct%gather_root == mepos
459 :
460 : CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
461 : nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
462 : row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
463 16 : row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
464 16 : symmetric = siesta_struct%symmetric
465 :
466 : ! number of locally stored SIESTA non-zero matrix elements
467 48 : n_nonzero_elements_siesta = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
468 : ! number of locally stored DBCSR non-zero matrix elements
469 48 : n_nonzero_elements_dbcsr = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_send))
470 :
471 : ! number of concurrent MPI isend / irecv operations
472 : nrequests_recv = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
473 16 : max_mpi_packet_size_dp)
474 : nrequests_total = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
475 16 : max_mpi_packet_size_dp) + nrequests_recv
476 :
477 16 : IF (nrequests_total > 0) THEN
478 : ! allocate MPI-related arrays. request_tag is not actually needed, as MPI standard guarantees the order of
479 : ! peer-to-peer messages with the same tag between same processes
480 64 : ALLOCATE (requests(nrequests_total))
481 48 : ALLOCATE (peer_rank(nrequests_total), request_tag(nrequests_total))
482 64 : ALLOCATE (offset_per_request(nrequests_total), nelements_per_request(nrequests_total))
483 : !requests(:) = mp_request_null
484 :
485 : ! split large messages into a number of smaller messages. It is not really needed
486 : ! unless we are going to send > 2^31 matrix elements per MPI request
487 16 : IF (nrequests_recv > 0) THEN
488 : CALL assign_nonzero_elements_to_requests(offset_per_request(1:nrequests_recv), &
489 : nelements_per_request(1:nrequests_recv), &
490 : peer_rank(1:nrequests_recv), &
491 : request_tag(1:nrequests_recv), &
492 : mepos, &
493 : siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
494 8 : max_mpi_packet_size_dp)
495 : END IF
496 16 : IF (nrequests_total > nrequests_recv) THEN
497 : CALL assign_nonzero_elements_to_requests(offset_per_request(nrequests_recv + 1:nrequests_total), &
498 : nelements_per_request(nrequests_recv + 1:nrequests_total), &
499 : peer_rank(nrequests_recv + 1:nrequests_total), &
500 : request_tag(nrequests_recv + 1:nrequests_total), &
501 : mepos, &
502 : siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
503 8 : max_mpi_packet_size_dp)
504 : END IF
505 : END IF
506 :
507 : ! point-to-point recv/send can be replaced with alltoallv, if data to distribute per rank is < 2^31 elements
508 : ! (should be OK due to SMEAGOL limitation)
509 : ! in principle, it is possible to overcome this limit by using derived datatypes
510 : ! (which will require additional wrapper functions, indeed).
511 : !
512 : ! pre-post non-blocking receive operations
513 16 : IF (n_nonzero_elements_siesta > 0) THEN
514 24 : ALLOCATE (recv_buffer(n_nonzero_elements_siesta))
515 : END IF
516 24 : DO irequest = 1, nrequests_recv
517 : CALL para_env%irecv(recv_buffer(offset_per_request(irequest) + 1: &
518 : offset_per_request(irequest) + nelements_per_request(irequest)), &
519 24 : peer_rank(irequest), requests(irequest), request_tag(irequest))
520 : END DO
521 :
522 : ! pack local DBCSR non-zero matrix elements ordering by their target parallel process
523 64 : ALLOCATE (offset_per_proc(0:nprocs - 1), n_packed_elements_per_proc(0:nprocs - 1))
524 16 : offset_per_proc(0) = 0
525 32 : DO iproc = 1, nprocs - 1
526 32 : offset_per_proc(iproc) = offset_per_proc(iproc - 1) + siesta_struct%nelements_per_proc(iproc - 1, nelements_dbcsr_send)
527 : END DO
528 16 : n_packed_elements_per_proc(:) = 0
529 :
530 : ! number of local neighbour-list nodes and offset of the first local neighbour-list node
531 16 : nnodes_proc = siesta_struct%nnodes_per_proc(mepos)
532 : !node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos)) - siesta_struct%nnodes_per_proc(mepos)
533 40 : node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos)) - nnodes_proc
534 :
535 : ! if do_distribute == .FALSE., send all matrix elements to MPI process with rank gather_root
536 : ! in case of do_distribute == .TRUE., iproc is determined by calling WhichNodeOrb()
537 16 : iproc = siesta_struct%gather_root
538 :
539 16 : IF (n_nonzero_elements_dbcsr > 0) THEN
540 48 : ALLOCATE (send_buffer(n_nonzero_elements_dbcsr))
541 16 : send_buffer(:) = 0.0_dp
542 :
543 : ! iterate over locally-stored DBCSR matrix blocks.
544 : ! inode_proc is the target parallel process (where data are going to be sent)
545 16 : image_ind_offset = 0
546 59568 : DO inode_proc = 1, nnodes_proc
547 59552 : n_image_ind = siesta_struct%n_dbcsr_cell_images_to_merge(inode_proc)
548 59552 : IF (n_image_ind > 0) THEN
549 59552 : inode = node_offset + inode_proc
550 :
551 59552 : irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
552 59552 : icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
553 59552 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
554 59552 : image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
555 59552 : image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
556 :
557 59552 : nrows_local = row_blk_size(irow_blk)
558 59552 : ncols_local = col_blk_size(icol_blk)
559 59552 : first_row_minus_one = row_blk_offset(irow_blk) - 1
560 59552 : first_col_minus_one = col_blk_offset(icol_blk) - 1
561 :
562 : ! merging cell images along transport direction
563 59552 : IF (n_image_ind == 1) THEN
564 : ! the most common case. Nothing to merge, so there is no need to allocate memory for a merged block
565 59552 : image_dbcsr = siesta_struct%dbcsr_cell_image_to_merge(image_ind_offset + 1)
566 : CALL dbcsr_get_block_p(matrix=matrix_dbcsr_kp(image_dbcsr)%matrix, &
567 59552 : row=irow_blk, col=icol_blk, block=sm_block_merged, found=found)
568 59552 : CPASSERT(found)
569 : ELSE ! n_image_ind > 1
570 0 : ALLOCATE (sm_block_merged(nrows_local, ncols_local))
571 :
572 0 : DO image_ind = 1, n_image_ind
573 0 : image_dbcsr = siesta_struct%dbcsr_cell_image_to_merge(image_ind + image_ind_offset)
574 :
575 : CALL dbcsr_get_block_p(matrix=matrix_dbcsr_kp(image_dbcsr)%matrix, &
576 0 : row=irow_blk, col=icol_blk, block=sm_block, found=found)
577 0 : CPASSERT(found)
578 : sm_block_merged(1:nrows_local, 1:ncols_local) = sm_block_merged(1:nrows_local, 1:ncols_local) + &
579 0 : sm_block(1:nrows_local, 1:ncols_local)
580 : END DO
581 : END IF
582 :
583 : ! pack matrix elements for the 'normal' SIESTA matrix block
584 59552 : IF (image_siesta > 0) THEN
585 595520 : DO irow_local = 1, nrows_local
586 535968 : IF (do_distribute) THEN
587 : #if defined(__SMEAGOL)
588 0 : CALL WhichNodeOrb(irow_local + first_row_minus_one, nprocs, iproc)
589 : #else
590 : CALL cp_abort(__LOCATION__, &
591 : "CP2K was compiled with no SMEAGOL support.")
592 : #endif
593 : END IF
594 :
595 : ! CPASSERT
596 : IF (debug_this_module) THEN
597 : CPASSERT(iproc >= 0 .AND. iproc < nprocs)
598 : IF (n_packed_elements_per_proc(iproc) + ncols_local > &
599 : siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
600 : CALL cp__a(__SHORT_FILE__, __LINE__)
601 : END IF
602 : END IF
603 :
604 : offset_send_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
605 535968 : send_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = sm_block_merged(irow_local, 1:ncols_local)
606 5359680 :
607 : n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + ncols_local
608 595520 : END DO
609 : END IF
610 :
611 : ! pack matrix elements of the transposed SIESTA matrix block
612 : IF (image_siesta_transp > 0) THEN
613 59552 : DO icol_local = 1, ncols_local
614 549600 : IF (do_distribute) THEN
615 494640 : #if defined(__SMEAGOL)
616 : CALL WhichNodeOrb(icol_local + first_col_minus_one, nprocs, iproc) ! iproc_orb
617 0 : #else
618 : CALL cp_abort(__LOCATION__, &
619 : "CP2K was compiled with no SMEAGOL support.")
620 : #endif
621 : END IF
622 :
623 : ! CPASSERT
624 : IF (debug_this_module) THEN
625 : CPASSERT(iproc >= 0 .AND. iproc < nprocs)
626 : IF (n_packed_elements_per_proc(iproc) + nrows_local > &
627 : siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
628 : CALL cp__a(__SHORT_FILE__, __LINE__)
629 : END IF
630 : END IF
631 :
632 : offset_send_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
633 : send_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = sm_block_merged(1:nrows_local, icol_local)
634 494640 :
635 4946400 : n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + nrows_local
636 : END DO
637 549600 : END IF
638 :
639 : IF (n_image_ind > 1) THEN
640 : DEALLOCATE (sm_block_merged)
641 59552 : END IF
642 0 : END IF
643 :
644 : image_ind_offset = image_ind_offset + siesta_struct%n_dbcsr_cell_images_to_merge(inode_proc)
645 : END DO
646 59568 :
647 : IF (debug_this_module) THEN
648 : DO iproc = 0, nprocs - 1
649 : IF (n_packed_elements_per_proc(iproc) /= siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
650 : CALL cp__a(__SHORT_FILE__, __LINE__)
651 : END IF
652 : END DO
653 : END IF
654 :
655 : ! send packed data to other parallel processes
656 : DO irequest = nrequests_recv + 1, nrequests_total
657 : CALL para_env%isend(send_buffer(offset_per_request(irequest) + 1: &
658 24 : offset_per_request(irequest) + nelements_per_request(irequest)), &
659 : peer_rank(irequest), requests(irequest), request_tag(irequest))
660 : END DO
661 24 :
662 : ! copy data locally that stay on the same process.
663 : IF (mepos > 0) THEN
664 : offset_recv_mepos = SUM(siesta_struct%nelements_per_proc(0:mepos - 1, nelements_dbcsr_recv))
665 16 : ELSE
666 16 : offset_recv_mepos = 0
667 : END IF
668 : offset_send_mepos = offset_per_proc(mepos)
669 :
670 16 : IF (debug_this_module) THEN
671 : IF (n_packed_elements_per_proc(mepos) /= siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_recv)) THEN
672 : CALL cp__a(__SHORT_FILE__, __LINE__)
673 : END IF
674 : END IF
675 :
676 : IF (n_packed_elements_per_proc(mepos) > 0) THEN
677 : recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + n_packed_elements_per_proc(mepos)) = &
678 16 : send_buffer(offset_send_mepos + 1:offset_send_mepos + n_packed_elements_per_proc(mepos))
679 : END IF
680 4591088 : END IF
681 :
682 : IF (nrequests_total > 0) THEN
683 : ! wait for pending isend/irecv requests
684 16 : CALL mp_waitall(requests)
685 : DEALLOCATE (nelements_per_request, offset_per_request, peer_rank, requests, request_tag)
686 16 : END IF
687 16 :
688 : ! release send buffers
689 : IF (ALLOCATED(send_buffer)) DEALLOCATE (send_buffer)
690 : DEALLOCATE (offset_per_proc, n_packed_elements_per_proc)
691 16 :
692 16 : ! non-zero matrix elements in 'recv_buffer' array are grouped by their source MPI rank,
693 : ! local row index, and column index (in this order).
694 : ! Reorder the matrix elements ('reorder_recv_buffer') so they are grouped by their local row index,
695 : ! source MPI rank, and column index.
696 : ! (column indices are in the ascending order within each (row index, source MPI rank) block).
697 : ! The array 'packed_index' allows mapping matrix element between these intermediate order and SIESTA order :
698 : ! local row index, column index
699 : IF (n_nonzero_elements_siesta > 0) THEN
700 : ALLOCATE (reorder_recv_buffer(n_nonzero_elements_siesta))
701 16 : ALLOCATE (next_nonzero_element_offset(siesta_struct%nrows))
702 24 : next_nonzero_element_offset(:) = 0
703 24 : offset_recv_mepos = 0
704 8 :
705 8 : DO inode = 1, SIZE(siesta_struct%nl_repl, 2)
706 : irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
707 59560 : icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
708 59552 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
709 59552 : image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
710 59552 : image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
711 59552 :
712 59552 : nrows_local = row_blk_size(irow_blk)
713 : ncols_local = col_blk_size(icol_blk)
714 59552 : first_row_minus_one = row_blk_offset(irow_blk) - 1
715 59552 : first_col_minus_one = col_blk_offset(icol_blk) - 1
716 59552 :
717 59552 : ! normal block
718 : IF (image_siesta > 0) THEN
719 : DO irow_local = 1, nrows_local
720 59552 : IF (do_distribute) THEN
721 595520 : #if defined(__SMEAGOL)
722 535968 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
723 : #else
724 0 : CALL cp_abort(__LOCATION__, &
725 : "CP2K was compiled with no SMEAGOL support.")
726 : #endif
727 : ELSE
728 : IF (is_root_rank) THEN
729 : irow_proc = irow_local + first_row_minus_one
730 535968 : ELSE
731 535968 : irow_proc = 0
732 : END IF
733 0 : END IF
734 : IF (irow_proc > 0) THEN
735 : offset_send_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
736 595520 : reorder_recv_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = &
737 535968 : recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + ncols_local)
738 : offset_recv_mepos = offset_recv_mepos + ncols_local
739 5359680 : next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + ncols_local
740 535968 : END IF
741 535968 : END DO
742 : END IF
743 :
744 : ! transposed block
745 : IF (image_siesta_transp > 0) THEN
746 : DO icol_local = 1, ncols_local
747 59560 : IF (do_distribute) THEN
748 549600 : #if defined(__SMEAGOL)
749 494640 : CALL GlobalToLocalOrb(icol_local + first_col_minus_one, mepos, nprocs, irow_proc)
750 : #else
751 0 : CALL cp_abort(__LOCATION__, &
752 : "CP2K was compiled with no SMEAGOL support.")
753 : #endif
754 : ELSE
755 : IF (is_root_rank) THEN
756 : irow_proc = icol_local + first_col_minus_one
757 494640 : ELSE
758 494640 : irow_proc = 0
759 : END IF
760 0 : END IF
761 : IF (irow_proc > 0) THEN
762 : offset_send_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
763 549600 : reorder_recv_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = &
764 494640 : recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + nrows_local)
765 : offset_recv_mepos = offset_recv_mepos + nrows_local
766 4946400 : next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + nrows_local
767 494640 : END IF
768 494640 : END DO
769 : END IF
770 : END DO
771 :
772 : IF (debug_this_module) THEN
773 : DO irow_local = 1, siesta_struct%nrows
774 : IF (siesta_struct%n_nonzero_cols(irow_local) /= next_nonzero_element_offset(irow_local)) THEN
775 : CALL cp__a(__SHORT_FILE__, __LINE__)
776 : END IF
777 : END DO
778 : END IF
779 :
780 : DEALLOCATE (next_nonzero_element_offset)
781 : DEALLOCATE (recv_buffer)
782 8 :
783 8 : ! Map non-zero matrix element between the intermediate order and SIESTA order
784 : DO irow_local = 1, siesta_struct%nrows
785 : offset_recv_mepos = siesta_struct%row_offset(irow_local)
786 2456 : DO icol_local = 1, siesta_struct%n_nonzero_cols(irow_local)
787 2448 : matrix_siesta(offset_recv_mepos + icol_local) = &
788 9277928 : reorder_recv_buffer(offset_recv_mepos + siesta_struct%packed_index(offset_recv_mepos + icol_local))
789 : END DO
790 9277920 : END DO
791 : DEALLOCATE (reorder_recv_buffer)
792 : END IF
793 8 :
794 : CALL timestop(handle)
795 : END SUBROUTINE convert_dbcsr_to_distributed_siesta
796 16 :
797 48 : ! **************************************************************************************************
798 : !> \brief Convert matrix from DBCSR to sparse SIESTA format.
799 : !> \param matrix_dbcsr_kp DBCSR matrix [out]. The matrix is declared as INTENT(in) as pointers to
800 : !> dbcsr matrices remain intact. However we have intention to update matrix elements
801 : !> \param matrix_siesta matrix in SIESTA format [in]
802 : !> \param siesta_struct structure to map matrix blocks between formats
803 : !> \param para_env MPI parallel environment
804 : ! **************************************************************************************************
805 : SUBROUTINE convert_distributed_siesta_to_dbcsr(matrix_dbcsr_kp, matrix_siesta, siesta_struct, para_env)
806 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
807 0 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: matrix_siesta
808 : TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
809 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
810 :
811 : CHARACTER(len=*), PARAMETER :: routineN = 'convert_distributed_siesta_to_dbcsr'
812 :
813 : INTEGER :: first_col_minus_one, first_row_minus_one, handle, icol_blk, icol_local, &
814 : image_dbcsr, image_siesta, image_siesta_transp, inode, inode_proc, iproc, irequest, &
815 : irow_blk, irow_local, irow_proc, mepos, n_image_ind, ncols_blk, ncols_local, nnodes_proc, &
816 : node_offset, nprocs, nrequests_recv, nrequests_total, nrows_blk, nrows_local
817 : INTEGER(kind=int_8) :: n_nonzero_elements_dbcsr, &
818 : n_nonzero_elements_siesta, &
819 : offset_recv_mepos, offset_send_mepos
820 : INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: n_packed_elements_per_proc, &
821 : nelements_per_request, &
822 0 : offset_per_proc, offset_per_request
823 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: next_nonzero_element_offset, peer_rank, &
824 0 : request_tag
825 0 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
826 0 : row_blk_offset, row_blk_size
827 0 : LOGICAL :: do_distribute, found, is_root_rank, &
828 0 : symmetric
829 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buffer, reorder_send_buffer, &
830 : send_buffer
831 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
832 0 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: requests
833 0 :
834 0 : CALL timeset(routineN, handle)
835 : DO image_dbcsr = 1, SIZE(matrix_dbcsr_kp)
836 0 : CALL dbcsr_set(matrix_dbcsr_kp(image_dbcsr)%matrix, 0.0_dp)
837 0 : END DO
838 0 :
839 : mepos = para_env%mepos
840 : nprocs = para_env%num_pe
841 0 : do_distribute = siesta_struct%gather_root < 0
842 0 : is_root_rank = siesta_struct%gather_root == mepos
843 0 :
844 0 : CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
845 : nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
846 : row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
847 : row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
848 : symmetric = siesta_struct%symmetric
849 0 :
850 0 : n_nonzero_elements_siesta = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
851 : n_nonzero_elements_dbcsr = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_send))
852 0 :
853 0 : nrequests_recv = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
854 : max_mpi_packet_size_dp)
855 : nrequests_total = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
856 0 : max_mpi_packet_size_dp) + nrequests_recv
857 : IF (nrequests_total > 0) THEN
858 0 : ALLOCATE (requests(nrequests_total))
859 0 : ALLOCATE (peer_rank(nrequests_total), request_tag(nrequests_total))
860 0 : ALLOCATE (offset_per_request(nrequests_total), nelements_per_request(nrequests_total))
861 0 : !requests(:) = mp_request_null
862 0 : IF (nrequests_recv > 0) THEN
863 : CALL assign_nonzero_elements_to_requests(offset_per_request(1:nrequests_recv), &
864 0 : nelements_per_request(1:nrequests_recv), &
865 : peer_rank(1:nrequests_recv), &
866 : request_tag(1:nrequests_recv), &
867 : mepos, &
868 : siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
869 : max_mpi_packet_size_dp)
870 : END IF
871 0 : IF (nrequests_total > nrequests_recv) THEN
872 : CALL assign_nonzero_elements_to_requests(offset_per_request(nrequests_recv + 1:nrequests_total), &
873 0 : nelements_per_request(nrequests_recv + 1:nrequests_total), &
874 : peer_rank(nrequests_recv + 1:nrequests_total), &
875 : request_tag(nrequests_recv + 1:nrequests_total), &
876 : mepos, &
877 : siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
878 : max_mpi_packet_size_dp)
879 : END IF
880 0 : END IF
881 :
882 : IF (n_nonzero_elements_dbcsr > 0) THEN
883 : ALLOCATE (recv_buffer(n_nonzero_elements_dbcsr))
884 0 : END IF
885 0 : DO irequest = 1, nrequests_recv
886 : CALL para_env%irecv(recv_buffer(offset_per_request(irequest) + 1: &
887 0 : offset_per_request(irequest) + nelements_per_request(irequest)), &
888 : peer_rank(irequest), requests(irequest), request_tag(irequest))
889 : END DO
890 0 :
891 : ALLOCATE (offset_per_proc(0:nprocs - 1), n_packed_elements_per_proc(0:nprocs - 1))
892 : offset_per_proc(0) = 0
893 0 : DO iproc = 1, nprocs - 1
894 0 : offset_per_proc(iproc) = offset_per_proc(iproc - 1) + siesta_struct%nelements_per_proc(iproc - 1, nelements_dbcsr_send)
895 0 : END DO
896 0 : n_packed_elements_per_proc(:) = 0
897 :
898 0 : IF (mepos > 0) THEN
899 : node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos - 1))
900 0 : ELSE
901 0 : node_offset = 0
902 : END IF
903 : nnodes_proc = siesta_struct%nnodes_per_proc(mepos)
904 :
905 0 : IF (n_nonzero_elements_siesta > 0) THEN
906 : ALLOCATE (send_buffer(n_nonzero_elements_siesta))
907 0 :
908 0 : ALLOCATE (reorder_send_buffer(n_nonzero_elements_siesta))
909 : DO irow_local = 1, siesta_struct%nrows
910 0 : offset_send_mepos = siesta_struct%row_offset(irow_local)
911 0 : DO icol_local = 1, siesta_struct%n_nonzero_cols(irow_local)
912 0 : reorder_send_buffer(offset_send_mepos + siesta_struct%packed_index(offset_send_mepos + icol_local)) = &
913 0 : matrix_siesta(offset_send_mepos + icol_local)
914 : END DO
915 0 : END DO
916 :
917 : ALLOCATE (next_nonzero_element_offset(siesta_struct%nrows))
918 : next_nonzero_element_offset(:) = 0
919 0 : offset_send_mepos = 0
920 0 :
921 0 : DO inode = 1, SIZE(siesta_struct%nl_repl, 2)
922 : irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
923 0 : icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
924 0 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
925 0 : image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
926 0 : image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
927 0 :
928 0 : nrows_local = row_blk_size(irow_blk)
929 : ncols_local = col_blk_size(icol_blk)
930 0 : first_row_minus_one = row_blk_offset(irow_blk) - 1
931 0 : first_col_minus_one = col_blk_offset(icol_blk) - 1
932 0 :
933 0 : IF (image_siesta > 0) THEN
934 : DO irow_local = 1, nrows_local
935 0 : IF (do_distribute) THEN
936 0 : #if defined(__SMEAGOL)
937 0 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
938 : #else
939 0 : CALL cp_abort(__LOCATION__, &
940 : "CP2K was compiled with no SMEAGOL support.")
941 : #endif
942 : ELSE
943 : IF (is_root_rank) THEN
944 : irow_proc = irow_local + first_row_minus_one
945 0 : ELSE
946 0 : irow_proc = 0
947 : END IF
948 0 : END IF
949 : IF (irow_proc > 0) THEN
950 : offset_recv_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
951 0 : send_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = &
952 0 : reorder_send_buffer(offset_recv_mepos + 1:offset_recv_mepos + ncols_local)
953 : offset_send_mepos = offset_send_mepos + ncols_local
954 0 : next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + ncols_local
955 0 : END IF
956 0 : END DO
957 : END IF
958 :
959 : ! transposed block
960 : IF (image_siesta_transp > 0) THEN
961 : DO icol_local = 1, ncols_local
962 0 : IF (do_distribute) THEN
963 0 : #if defined(__SMEAGOL)
964 0 : CALL GlobalToLocalOrb(icol_local + first_col_minus_one, mepos, nprocs, irow_proc)
965 : #else
966 0 : CALL cp_abort(__LOCATION__, &
967 : "CP2K was compiled with no SMEAGOL support.")
968 : #endif
969 : ELSE
970 : IF (is_root_rank) THEN
971 : irow_proc = icol_local + first_col_minus_one
972 0 : ELSE
973 0 : irow_proc = 0
974 : END IF
975 0 : END IF
976 : IF (irow_proc > 0) THEN
977 : offset_recv_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
978 0 : send_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = &
979 0 : reorder_send_buffer(offset_recv_mepos + 1:offset_recv_mepos + nrows_local)
980 : offset_send_mepos = offset_send_mepos + nrows_local
981 0 : next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + nrows_local
982 0 : END IF
983 0 : END DO
984 : END IF
985 : END DO
986 :
987 : IF (debug_this_module) THEN
988 : DO irow_local = 1, siesta_struct%nrows
989 : IF (siesta_struct%n_nonzero_cols(irow_local) /= next_nonzero_element_offset(irow_local)) THEN
990 : CALL cp__a(__SHORT_FILE__, __LINE__)
991 : END IF
992 : END DO
993 : END IF
994 :
995 : DEALLOCATE (next_nonzero_element_offset)
996 : DEALLOCATE (reorder_send_buffer)
997 0 :
998 0 : DO irequest = nrequests_recv + 1, nrequests_total
999 : CALL para_env%isend(send_buffer(offset_per_request(irequest) + 1: &
1000 0 : offset_per_request(irequest) + nelements_per_request(irequest)), &
1001 : peer_rank(irequest), requests(irequest), request_tag(irequest))
1002 : END DO
1003 0 :
1004 : ! copy data locally that stay on the same process.
1005 : IF (mepos > 0) THEN
1006 : offset_send_mepos = SUM(siesta_struct%nelements_per_proc(0:mepos - 1, nelements_dbcsr_recv))
1007 0 : ELSE
1008 0 : offset_send_mepos = 0
1009 : END IF
1010 : offset_recv_mepos = offset_per_proc(mepos)
1011 :
1012 0 : IF (debug_this_module) THEN
1013 : IF (siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_recv) /= &
1014 : siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_send)) THEN
1015 : CALL cp__a(__SHORT_FILE__, __LINE__)
1016 : END IF
1017 : END IF
1018 :
1019 : IF (siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_send) > 0) THEN
1020 : recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_send)) = &
1021 0 : send_buffer(offset_send_mepos + 1:offset_send_mepos + siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_send))
1022 : END IF
1023 0 : END IF
1024 :
1025 : IF (nrequests_total > 0) THEN
1026 : ! wait for pending isend/irecv requests
1027 0 : CALL mp_waitall(requests)
1028 : DEALLOCATE (nelements_per_request, offset_per_request, peer_rank, requests, request_tag)
1029 0 : END IF
1030 0 :
1031 : IF (ALLOCATED(send_buffer)) DEALLOCATE (send_buffer)
1032 :
1033 0 : ! if do_distribute == .FALSE., collect matrix elements from MPI process with rank gather_root
1034 : iproc = siesta_struct%gather_root
1035 : IF (n_nonzero_elements_dbcsr > 0) THEN
1036 0 : DO inode_proc = 1, nnodes_proc
1037 0 : n_image_ind = siesta_struct%n_dbcsr_cell_images_to_merge(inode_proc)
1038 0 : IF (n_image_ind > 0) THEN
1039 0 : inode = node_offset + inode_proc
1040 0 :
1041 0 : irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
1042 : icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
1043 0 : image_dbcsr = siesta_struct%nl_repl(neighbor_list_dbcsr_image_index, inode)
1044 0 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
1045 0 : image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
1046 0 : image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
1047 0 :
1048 0 : nrows_local = row_blk_size(irow_blk)
1049 : ncols_local = col_blk_size(icol_blk)
1050 0 : first_row_minus_one = row_blk_offset(irow_blk) - 1
1051 0 : first_col_minus_one = col_blk_offset(icol_blk) - 1
1052 0 :
1053 0 : CALL dbcsr_get_block_p(matrix=matrix_dbcsr_kp(image_dbcsr)%matrix, &
1054 : row=irow_blk, col=icol_blk, block=sm_block, found=found)
1055 : CPASSERT(found)
1056 0 :
1057 0 : IF (image_siesta > 0) THEN
1058 : DO irow_local = 1, nrows_local
1059 0 : IF (do_distribute) THEN
1060 0 : #if defined(__SMEAGOL)
1061 0 : CALL WhichNodeOrb(irow_local + first_row_minus_one, nprocs, iproc) ! iproc_orb
1062 : #else
1063 0 : CALL cp_abort(__LOCATION__, &
1064 : "CP2K was compiled with no SMEAGOL support.")
1065 : #endif
1066 : END IF
1067 : ! CPASSERT
1068 : IF (debug_this_module) THEN
1069 : CPASSERT(iproc >= 0 .AND. iproc < nprocs)
1070 : IF (n_packed_elements_per_proc(iproc) + ncols_local > &
1071 : siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
1072 : CALL cp__a(__SHORT_FILE__, __LINE__)
1073 : END IF
1074 : END IF
1075 :
1076 : offset_recv_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
1077 : sm_block(irow_local, 1:ncols_local) = recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + ncols_local)
1078 :
1079 0 : n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + ncols_local
1080 0 : END DO
1081 : END IF
1082 0 :
1083 : ! transposed block
1084 : IF (image_siesta_transp > 0) THEN
1085 : DO icol_local = 1, ncols_local
1086 : IF (do_distribute) THEN
1087 0 : #if defined(__SMEAGOL)
1088 0 : CALL WhichNodeOrb(icol_local + first_col_minus_one, nprocs, iproc) ! iproc_orb
1089 0 : #else
1090 : CALL cp_abort(__LOCATION__, &
1091 0 : "CP2K was compiled with no SMEAGOL support.")
1092 : #endif
1093 : END IF
1094 : ! CPASSERT
1095 : IF (debug_this_module) THEN
1096 : CPASSERT(iproc >= 0 .AND. iproc < nprocs)
1097 : IF (n_packed_elements_per_proc(iproc) + nrows_local > &
1098 : siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
1099 : CALL cp__a(__SHORT_FILE__, __LINE__)
1100 : END IF
1101 : END IF
1102 :
1103 : offset_recv_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
1104 : sm_block(1:nrows_local, icol_local) = recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + nrows_local)
1105 :
1106 : n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + nrows_local
1107 0 : END DO
1108 0 : END IF
1109 : END IF
1110 0 : END DO
1111 :
1112 : IF (debug_this_module) THEN
1113 : DO iproc = 0, nprocs - 1
1114 : IF (n_packed_elements_per_proc(iproc) /= siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
1115 : CALL cp__a(__SHORT_FILE__, __LINE__)
1116 : END IF
1117 : END DO
1118 : END IF
1119 :
1120 : DEALLOCATE (recv_buffer)
1121 : END IF
1122 :
1123 : DEALLOCATE (offset_per_proc, n_packed_elements_per_proc)
1124 0 :
1125 : CALL timestop(handle)
1126 : END SUBROUTINE convert_distributed_siesta_to_dbcsr
1127 0 :
1128 : ! *** PRIVATE SUBROUTINES ***
1129 0 :
1130 0 : ! **************************************************************************************************
1131 : !> \brief Computes number of neighbour-list nodes on the current parallel process.
1132 : !> \param nnodes_local number of nodes [out]
1133 : !> \param max_ijk_cell_image_local largest index of cell images along i, j and k cell vectors
1134 : !> on this parallel process [out]
1135 : !> \param max_ijk_cell_image largest index of cell images along i, j and k cell vectors [inout]
1136 : !> \param sab_nl pair-wise neighbour list [in]
1137 : !> \param para_env MPI parallel environment [in]
1138 : !> \param particle_coords list of atomic coordinates subject to periodic boundary conditions [in]
1139 : !> \param cell simulation unit cell [in]
1140 : ! **************************************************************************************************
1141 : SUBROUTINE get_nnodes_local(nnodes_local, max_ijk_cell_image_local, max_ijk_cell_image, sab_nl, para_env, particle_coords, cell)
1142 : INTEGER, INTENT(out) :: nnodes_local
1143 : INTEGER, DIMENSION(3), INTENT(out) :: max_ijk_cell_image_local
1144 : INTEGER, DIMENSION(3), INTENT(inout) :: max_ijk_cell_image
1145 4 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1146 : INTENT(in), POINTER :: sab_nl
1147 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
1148 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1149 : INTENT(in) :: particle_coords
1150 : TYPE(cell_type), INTENT(in), POINTER :: cell
1151 :
1152 : CHARACTER(len=*), PARAMETER :: routineN = 'get_nnodes_local'
1153 :
1154 : INTEGER :: handle, iatom, icoord, jatom
1155 : INTEGER, DIMENSION(3) :: cell_ijk, max_ijk_cell_image_tmp
1156 : LOGICAL :: update_ncells
1157 : REAL(kind=dp), DIMENSION(3) :: r_ij
1158 : TYPE(neighbor_list_iterator_p_type), &
1159 : DIMENSION(:), POINTER :: nl_iterator
1160 :
1161 : CALL timeset(routineN, handle)
1162 :
1163 4 : update_ncells = .FALSE.
1164 : DO icoord = 1, 3 ! x, y, z
1165 4 : IF (max_ijk_cell_image(icoord) >= 0) THEN
1166 : max_ijk_cell_image_tmp(icoord) = max_ijk_cell_image(icoord)
1167 4 : ELSE
1168 16 : max_ijk_cell_image_tmp(icoord) = HUGE(max_ijk_cell_image_tmp(icoord))
1169 16 : update_ncells = .TRUE.
1170 4 : END IF
1171 : END DO
1172 8 :
1173 8 : nnodes_local = 0
1174 : max_ijk_cell_image_local(:) = 0
1175 :
1176 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1177 4 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1178 4 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=r_ij)
1179 : CALL get_negf_cell_ijk(cell_ijk, r_ij, r_i=particle_coords(1:3, iatom), r_j=particle_coords(1:3, jatom), cell=cell)
1180 4 : cell_ijk(1:3) = ABS(cell_ijk(1:3))
1181 15164 :
1182 15160 : IF (cell_ijk(1) <= max_ijk_cell_image_tmp(1) .AND. cell_ijk(2) <= max_ijk_cell_image_tmp(2) .AND. &
1183 15160 : cell_ijk(3) <= max_ijk_cell_image_tmp(3)) THEN
1184 60640 : nnodes_local = nnodes_local + 1
1185 : max_ijk_cell_image_local(1:3) = MAX(max_ijk_cell_image_local(1:3), cell_ijk(1:3))
1186 15160 : END IF
1187 4 : END DO
1188 14888 : CALL neighbor_list_iterator_release(nl_iterator)
1189 59552 :
1190 : IF (update_ncells) THEN
1191 : max_ijk_cell_image_tmp(1:3) = max_ijk_cell_image_local(1:3)
1192 4 : CALL para_env%max(max_ijk_cell_image_tmp)
1193 : DO icoord = 1, 3
1194 4 : IF (max_ijk_cell_image(icoord) < 0) THEN
1195 4 : max_ijk_cell_image(icoord) = max_ijk_cell_image_tmp(icoord)
1196 4 : END IF
1197 16 : END DO
1198 16 : END IF
1199 8 :
1200 : CALL timestop(handle)
1201 : END SUBROUTINE get_nnodes_local
1202 :
1203 : ! **************************************************************************************************
1204 4 : !> \brief Construct list of neighbour-list's nodes on the current parallel process.
1205 4 : !> \param nl_local non-merged local neighbour-list's nodes [out]
1206 : !> \param max_ijk_cell_image_local largest index of cell images along i, j and k cell vectors
1207 : !> on this parallel process [in]
1208 : !> \param max_ijk_cell_image largest index of cell images along i, j and k cell vectors [in]
1209 : !> \param sab_nl pair-wise neighbour list [in]
1210 : !> \param particle_coords list of atomic coordinates subject to periodic boundary conditions [in]
1211 : !> \param cell simulation unit cell [in]
1212 : !> \param cell_to_index array to convert 3-D cell indices to 1-D DBCSR image indices
1213 : !> \param do_merge merge cell images along transport direction [in]
1214 : !> \param node_merged_indices nodes-related indices. Nodes with identical indices will be merged [out]
1215 : !> \param k_cells list of cell image indices along transport direction. Nodes to
1216 : !> be merged have the same 'node_merged_indices' but different 'k_cells' indices [out]
1217 : ! **************************************************************************************************
1218 : SUBROUTINE get_nl_nodes_local(nl_local, max_ijk_cell_image_local, max_ijk_cell_image, sab_nl, particle_coords, &
1219 : cell, cell_to_index, do_merge, node_merged_indices, k_cells)
1220 : INTEGER, DIMENSION(:, :), INTENT(out) :: nl_local
1221 : INTEGER, DIMENSION(3), INTENT(in) :: max_ijk_cell_image_local, &
1222 4 : max_ijk_cell_image
1223 4 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1224 : INTENT(in), POINTER :: sab_nl
1225 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1226 : INTENT(in) :: particle_coords
1227 : TYPE(cell_type), INTENT(in), POINTER :: cell
1228 : INTEGER, DIMENSION(:, :, :), INTENT(in), POINTER :: cell_to_index
1229 : LOGICAL, INTENT(in) :: do_merge
1230 : INTEGER(kind=int_8), DIMENSION(:), INTENT(out) :: node_merged_indices
1231 : INTEGER, DIMENSION(:), INTENT(out) :: k_cells
1232 :
1233 : CHARACTER(len=*), PARAMETER :: routineN = 'get_nl_nodes_local'
1234 :
1235 : INTEGER :: handle, iatom, icol_blk, image, inode, &
1236 : irow_blk, jatom, natoms
1237 : INTEGER(kind=8), DIMENSION(2) :: ncells_siesta_local
1238 : INTEGER, DIMENSION(2) :: ncells_siesta
1239 : INTEGER, DIMENSION(3) :: cell_ijk_abs, cell_ijk_dbcsr, &
1240 : cell_ijk_siesta
1241 : LOGICAL :: do_symmetric
1242 : REAL(kind=dp), DIMENSION(3) :: r_ij
1243 : TYPE(neighbor_list_iterator_p_type), &
1244 : DIMENSION(:), POINTER :: nl_iterator
1245 :
1246 : CALL timeset(routineN, handle)
1247 : ! natoms only used to compute a merged 1D index of each DBCSR block
1248 4 : natoms = SIZE(particle_coords, 2)
1249 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1250 4 : ncells_siesta(1:2) = 2*max_ijk_cell_image(1:2) + 1
1251 :
1252 4 : ncells_siesta_local(1:2) = INT(2*max_ijk_cell_image_local(1:2) + 1, kind=int_8)
1253 4 :
1254 12 : inode = 0
1255 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
1256 12 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1257 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell_ijk_dbcsr, r=r_ij)
1258 4 : CALL get_negf_cell_ijk(cell_ijk_abs, r_ij, r_i=particle_coords(1:3, iatom), r_j=particle_coords(1:3, jatom), cell=cell)
1259 4 :
1260 15164 : IF (ABS(cell_ijk_abs(1)) <= max_ijk_cell_image(1) .AND. ABS(cell_ijk_abs(2)) <= max_ijk_cell_image(2) .AND. &
1261 15160 : ABS(cell_ijk_abs(3)) <= max_ijk_cell_image(3)) THEN
1262 15160 :
1263 : inode = inode + 1
1264 15160 :
1265 4 : image = get_index_by_cell(cell_ijk_dbcsr, cell_to_index)
1266 : CPASSERT(image > 0)
1267 14888 : nl_local(neighbor_list_dbcsr_image_index, inode) = image
1268 :
1269 14888 : IF (do_symmetric .AND. iatom > jatom) THEN
1270 14888 : irow_blk = jatom
1271 14888 : icol_blk = iatom
1272 : cell_ijk_abs(1:3) = -cell_ijk_abs(1:3)
1273 14888 : ELSE
1274 28224 : irow_blk = iatom
1275 28224 : icol_blk = jatom
1276 28224 : END IF
1277 :
1278 : nl_local(neighbor_list_iatom_index, inode) = irow_blk
1279 : nl_local(neighbor_list_jatom_index, inode) = icol_blk
1280 :
1281 : cell_ijk_siesta(1:3) = index_in_canonical_enumeration(cell_ijk_abs(1:3)) ! absolute -> SIESTA
1282 14888 :
1283 14888 : IF (do_merge) THEN
1284 : node_merged_indices(inode) = (((cell_ijk_siesta(2) - 1)*ncells_siesta_local(1) + cell_ijk_siesta(1) - 1)* &
1285 59552 : INT(natoms, kind=int_8) + icol_blk - 1)*INT(natoms, kind=int_8) + INT(irow_blk - 1, kind=int_8)
1286 : image = cell_ijk_siesta(1) + ncells_siesta(1)*(cell_ijk_siesta(2) - 1)
1287 14888 : ELSE
1288 : node_merged_indices(inode) = ((((cell_ijk_siesta(3) - 1)*ncells_siesta_local(2) + &
1289 0 : cell_ijk_siesta(2) - 1)*ncells_siesta_local(1) + cell_ijk_siesta(1) - 1)* &
1290 0 : INT(natoms, kind=int_8) + icol_blk - 1)*INT(natoms, kind=int_8) + INT(irow_blk - 1, kind=int_8)
1291 : image = cell_ijk_siesta(1) + ncells_siesta(1)*(cell_ijk_siesta(2) - 1 + ncells_siesta(2)*(cell_ijk_siesta(3) - 1))
1292 : END IF
1293 : k_cells(inode) = cell_ijk_siesta(3)
1294 14888 : nl_local(neighbor_list_siesta_image_index, inode) = image
1295 14888 :
1296 : IF (do_symmetric .AND. irow_blk /= icol_blk) THEN
1297 14888 : cell_ijk_abs(1:3) = -cell_ijk_abs(1:3)
1298 14888 : cell_ijk_siesta(1:3) = index_in_canonical_enumeration(cell_ijk_abs(1:3)) ! absolute -> SIESTA
1299 : IF (do_merge) cell_ijk_siesta(3) = 1
1300 14888 : nl_local(neighbor_list_siesta_transp_image_index, inode) = &
1301 54960 : cell_ijk_siesta(1) + ncells_siesta(1)*(cell_ijk_siesta(2) - 1 + ncells_siesta(2)*(cell_ijk_siesta(3) - 1))
1302 54960 : ELSE
1303 13740 : nl_local(neighbor_list_siesta_transp_image_index, inode) = 0
1304 : END IF
1305 13740 : END IF
1306 : END DO
1307 1148 : CALL neighbor_list_iterator_release(nl_iterator)
1308 :
1309 : IF (debug_this_module) THEN
1310 : CPASSERT(SIZE(nl_local, 2) == inode)
1311 4 : END IF
1312 :
1313 : CALL timestop(handle)
1314 : END SUBROUTINE get_nl_nodes_local
1315 :
1316 : ! **************************************************************************************************
1317 4 : !> \brief Replicate (and optionally merge) pair-wise neighbour list.
1318 4 : !> \param repl_nl replicated neighbour list. It needs to be deallocated elsewhere [allocated]
1319 : !> \param n_dbcsr_cell_images_to_merge number of merged blocks per neighbour-list node [allocated]
1320 : !> \param dbcsr_cell_image_to_merge list of DBCSR image indices to merge [allocated]
1321 : !> \param nnodes_per_proc number of merged nodes on each parallel processes [out]
1322 : !> \param max_ijk_cell_image largest index of cell images along i, j and k cell vectors [inout]
1323 : !> \param sab_nl pair-wise neighbour list [in]
1324 : !> \param para_env MPI parallel environment [in]
1325 : !> \param particle_coords list of atomic coordinates subject to periodic boundary conditions [in]
1326 : !> \param cell simulation unit cell [in]
1327 : !> \param cell_to_index array to convert 3-D cell indices to 1-D DBCSR image indices [in]
1328 : !> \param do_merge merge cell images along transport direction [in]
1329 : ! **************************************************************************************************
1330 : SUBROUTINE replicate_neighbour_list(repl_nl, n_dbcsr_cell_images_to_merge, dbcsr_cell_image_to_merge, &
1331 : nnodes_per_proc, max_ijk_cell_image, sab_nl, para_env, particle_coords, &
1332 : cell, cell_to_index, do_merge)
1333 : INTEGER, ALLOCATABLE, DIMENSION(:, :), &
1334 4 : INTENT(inout) :: repl_nl
1335 4 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: n_dbcsr_cell_images_to_merge, &
1336 : dbcsr_cell_image_to_merge
1337 : INTEGER, DIMENSION(0:), INTENT(out) :: nnodes_per_proc
1338 : INTEGER, DIMENSION(3), INTENT(inout) :: max_ijk_cell_image
1339 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1340 : INTENT(in), POINTER :: sab_nl
1341 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
1342 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1343 : INTENT(in) :: particle_coords
1344 : TYPE(cell_type), INTENT(in), POINTER :: cell
1345 : INTEGER, DIMENSION(:, :, :), INTENT(in), POINTER :: cell_to_index
1346 : LOGICAL, INTENT(in) :: do_merge
1347 :
1348 : CHARACTER(len=*), PARAMETER :: routineN = 'replicate_neighbour_list'
1349 :
1350 : INTEGER :: handle, inode, iproc, kcell_closest, &
1351 : nnodes_local, nnodes_merged, &
1352 : nnodes_repl, offset_inode
1353 : INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: node_merged_indices
1354 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inodes_orig, k_cells
1355 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl_local
1356 : INTEGER, DIMENSION(3) :: max_ijk_cell_image_local
1357 4 :
1358 4 : CALL timeset(routineN, handle)
1359 4 : CPASSERT(.NOT. ALLOCATED(repl_nl))
1360 : CPASSERT(.NOT. ALLOCATED(n_dbcsr_cell_images_to_merge))
1361 : CPASSERT(.NOT. ALLOCATED(dbcsr_cell_image_to_merge))
1362 4 :
1363 4 : CALL get_nnodes_local(nnodes_local, max_ijk_cell_image_local, max_ijk_cell_image, sab_nl, para_env, particle_coords, cell)
1364 4 :
1365 4 : nnodes_per_proc(:) = 0
1366 :
1367 4 : IF (nnodes_local > 0) THEN
1368 : ALLOCATE (nl_local(neighbor_list_dim1, nnodes_local))
1369 12 : ALLOCATE (node_merged_indices(nnodes_local))
1370 : ALLOCATE (k_cells(nnodes_local))
1371 4 : CALL get_nl_nodes_local(nl_local, max_ijk_cell_image_local, max_ijk_cell_image, sab_nl, particle_coords, cell, &
1372 12 : cell_to_index, do_merge, node_merged_indices, k_cells)
1373 12 :
1374 12 : ALLOCATE (inodes_orig(nnodes_local))
1375 : CALL sort(node_merged_indices, nnodes_local, inodes_orig)
1376 4 :
1377 : nnodes_merged = 1
1378 8 : DO inode = 2, nnodes_local
1379 4 : IF (node_merged_indices(inode) > node_merged_indices(inode - 1)) nnodes_merged = nnodes_merged + 1
1380 : END DO
1381 4 : ELSE
1382 14888 : nnodes_merged = 0
1383 14888 : END IF
1384 :
1385 : nnodes_per_proc(para_env%mepos) = nnodes_merged
1386 : CALL para_env%sum(nnodes_per_proc)
1387 :
1388 : nnodes_repl = SUM(nnodes_per_proc(:))
1389 4 : ALLOCATE (repl_nl(neighbor_list_dim1, nnodes_repl))
1390 20 :
1391 : IF (nnodes_local > 0) THEN
1392 12 : IF (para_env%mepos > 0) THEN
1393 12 : offset_inode = SUM(nnodes_per_proc(0:para_env%mepos - 1))
1394 : ELSE
1395 4 : offset_inode = 0
1396 4 : END IF
1397 4 :
1398 : ALLOCATE (n_dbcsr_cell_images_to_merge(nnodes_merged))
1399 : ALLOCATE (dbcsr_cell_image_to_merge(nnodes_local))
1400 : n_dbcsr_cell_images_to_merge(:) = 0
1401 :
1402 12 : nnodes_merged = 1 !offset_inode + 1
1403 12 : repl_nl(:, offset_inode + 1) = nl_local(:, inodes_orig(1))
1404 4 : n_dbcsr_cell_images_to_merge(1) = 1
1405 : dbcsr_cell_image_to_merge(1) = nl_local(neighbor_list_dbcsr_image_index, inodes_orig(1))
1406 4 : kcell_closest = k_cells(inodes_orig(1))
1407 24 : DO inode = 2, nnodes_local
1408 4 : IF (node_merged_indices(inode) > node_merged_indices(inode - 1)) THEN
1409 4 : nnodes_merged = nnodes_merged + 1
1410 4 : repl_nl(:, offset_inode + nnodes_merged) = nl_local(:, inodes_orig(inode))
1411 14888 : !n_dbcsr_cell_images_to_merge(nnodes_merged) = 1
1412 14884 : kcell_closest = k_cells(inodes_orig(inode))
1413 14884 : ELSE
1414 89304 : IF (ABS(k_cells(inodes_orig(inode))) < ABS(kcell_closest) .OR. &
1415 : (ABS(k_cells(inodes_orig(inode))) == ABS(kcell_closest) .AND. kcell_closest < 0)) THEN
1416 14884 : repl_nl(:, offset_inode + nnodes_merged) = nl_local(:, inodes_orig(inode))
1417 : kcell_closest = k_cells(inodes_orig(inode))
1418 0 : END IF
1419 0 : END IF
1420 0 : dbcsr_cell_image_to_merge(inode) = nl_local(neighbor_list_dbcsr_image_index, inodes_orig(inode))
1421 0 : n_dbcsr_cell_images_to_merge(nnodes_merged) = n_dbcsr_cell_images_to_merge(nnodes_merged) + 1
1422 : END DO
1423 :
1424 14884 : IF (debug_this_module) THEN
1425 14888 : CPASSERT(SUM(n_dbcsr_cell_images_to_merge) == nnodes_local)
1426 : END IF
1427 :
1428 : DEALLOCATE (inodes_orig)
1429 : DEALLOCATE (node_merged_indices, k_cells)
1430 : DEALLOCATE (nl_local)
1431 : END IF
1432 4 :
1433 4 : IF (para_env%num_pe > 1) THEN
1434 4 : offset_inode = 0
1435 : DO iproc = 0, para_env%num_pe - 1
1436 : IF (nnodes_per_proc(iproc) > 0) THEN
1437 4 : CALL para_env%bcast(repl_nl(:, offset_inode + 1:offset_inode + nnodes_per_proc(iproc)), iproc)
1438 4 : offset_inode = offset_inode + nnodes_per_proc(iproc)
1439 12 : END IF
1440 12 : END DO
1441 8 : END IF
1442 8 :
1443 : CALL timestop(handle)
1444 : END SUBROUTINE replicate_neighbour_list
1445 :
1446 : ! **************************************************************************************************
1447 4 : !> \brief Count number of DBCSR matrix elements that should be received from (*,nelements_dbcsr_recv)
1448 8 : !> and send to (*,nelements_dbcsr_send) each parallel process.
1449 : !> \param nelements_per_proc number of non-zero matrix elements for each MPI process
1450 : !> \param nnodes_per_proc number of non-zero DBCSR matrix blocks (neighbour-list nodes)
1451 : !> \param nl_repl replicated neighbour list
1452 : !> \param matrix_dbcsr_kp DBCSR matrix
1453 : !> \param symmetric whether the DBCSR matrix is a symmetric one
1454 : !> \param para_env parallel environment
1455 : !> \param gather_root if >=0, gather all non-zero matrix element on the MPI process with
1456 : !> gather_root rank (useful for bulk transport calculation).
1457 : !> If <0, distribute non-zero matrix element across all MPI processes
1458 : ! **************************************************************************************************
1459 : SUBROUTINE count_remote_dbcsr_elements(nelements_per_proc, nnodes_per_proc, nl_repl, matrix_dbcsr_kp, &
1460 : symmetric, para_env, gather_root)
1461 : INTEGER(kind=int_8), DIMENSION(0:, :), INTENT(out) :: nelements_per_proc
1462 : INTEGER, DIMENSION(0:), INTENT(in) :: nnodes_per_proc
1463 4 : INTEGER, DIMENSION(:, :), INTENT(in) :: nl_repl
1464 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
1465 : LOGICAL, INTENT(in) :: symmetric
1466 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
1467 : INTEGER, INTENT(in) :: gather_root
1468 :
1469 : CHARACTER(len=*), PARAMETER :: routineN = 'count_remote_dbcsr_elements'
1470 :
1471 : INTEGER :: first_row_minus_one, handle, icol_blk, image, image_transp, inode, inode_proc, &
1472 : iproc, iproc_orb, irow_blk, irow_local, mepos, ncols_blk, ncols_local, nnodes_proc, &
1473 : nprocs, nrows_blk, nrows_local, offset_inode
1474 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1475 : row_blk_offset, row_blk_size
1476 : LOGICAL :: do_distribute
1477 :
1478 4 : CALL timeset(routineN, handle)
1479 4 : nelements_per_proc(:, :) = 0
1480 : mepos = para_env%mepos
1481 : nprocs = para_env%num_pe
1482 4 : do_distribute = gather_root < 0
1483 28 : IF (debug_this_module) THEN
1484 4 : CPASSERT(SIZE(nnodes_per_proc) == nprocs)
1485 4 : END IF
1486 4 :
1487 : CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
1488 : nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
1489 : row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
1490 : row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
1491 :
1492 : offset_inode = 0
1493 : iproc_orb = gather_root ! if do_distribute == .FALSE., send all matrix elements to MPI process with rank gather_root
1494 4 : DO iproc = LBOUND(nnodes_per_proc, 1), UBOUND(nnodes_per_proc, 1)
1495 : nnodes_proc = nnodes_per_proc(iproc)
1496 4 : DO inode_proc = 1, nnodes_proc
1497 4 : inode = inode_proc + offset_inode
1498 16 :
1499 8 : irow_blk = nl_repl(neighbor_list_iatom_index, inode)
1500 29784 : icol_blk = nl_repl(neighbor_list_jatom_index, inode)
1501 29776 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
1502 : image = nl_repl(neighbor_list_siesta_image_index, inode)
1503 29776 : image_transp = nl_repl(neighbor_list_siesta_transp_image_index, inode)
1504 29776 :
1505 29776 : IF (image > 0) THEN
1506 29776 : nrows_local = row_blk_size(irow_blk)
1507 29776 : first_row_minus_one = row_blk_offset(irow_blk) - 1
1508 : ncols_local = col_blk_size(icol_blk)
1509 29776 : DO irow_local = 1, nrows_local
1510 29776 : IF (do_distribute) THEN
1511 29776 : #if defined(__SMEAGOL)
1512 29776 : CALL WhichNodeOrb(irow_local + first_row_minus_one, nprocs, iproc_orb)
1513 297760 : #else
1514 267984 : CALL cp_abort(__LOCATION__, &
1515 : "CP2K was compiled with no SMEAGOL support.")
1516 0 : #endif
1517 : END IF
1518 : IF (iproc_orb == mepos) THEN
1519 : nelements_per_proc(iproc, nelements_dbcsr_recv) = nelements_per_proc(iproc, nelements_dbcsr_recv) + &
1520 : ncols_local
1521 : END IF
1522 267984 :
1523 : IF (iproc == mepos) THEN
1524 133992 : nelements_per_proc(iproc_orb, nelements_dbcsr_send) = nelements_per_proc(iproc_orb, nelements_dbcsr_send) + &
1525 : ncols_local
1526 : END IF
1527 297760 : END DO
1528 : END IF
1529 133992 :
1530 : ! transposed block
1531 : IF (image_transp > 0) THEN
1532 : nrows_local = col_blk_size(icol_blk)
1533 : first_row_minus_one = col_blk_offset(icol_blk) - 1
1534 : ncols_local = row_blk_size(irow_blk)
1535 29784 : DO irow_local = 1, nrows_local
1536 27480 : IF (do_distribute) THEN
1537 27480 : #if defined(__SMEAGOL)
1538 27480 : CALL WhichNodeOrb(irow_local + first_row_minus_one, nprocs, iproc_orb)
1539 274800 : #else
1540 247320 : CALL cp_abort(__LOCATION__, &
1541 : "CP2K was compiled with no SMEAGOL support.")
1542 0 : #endif
1543 : END IF
1544 : IF (iproc_orb == mepos) THEN
1545 : nelements_per_proc(iproc, nelements_dbcsr_recv) = nelements_per_proc(iproc, nelements_dbcsr_recv) + &
1546 : ncols_local
1547 : END IF
1548 247320 :
1549 : IF (iproc == mepos) THEN
1550 123660 : nelements_per_proc(iproc_orb, nelements_dbcsr_send) = nelements_per_proc(iproc_orb, nelements_dbcsr_send) + &
1551 : ncols_local
1552 : END IF
1553 274800 : END DO
1554 : END IF
1555 123660 : END DO
1556 : offset_inode = offset_inode + nnodes_proc
1557 : END DO
1558 : CALL timestop(handle)
1559 : END SUBROUTINE count_remote_dbcsr_elements
1560 12 :
1561 : ! **************************************************************************************************
1562 4 : !> \brief Construct list of non-zero matrix elements' indices in SIESTA format.
1563 4 : !> \param n_nonzero_cols number of non-zero matrix elements on each matrix row local to the current
1564 : !> MPI process
1565 : !> \param row_offset offset of the first non-zero matrix elements for each locally-stores row
1566 : !> \param col_index sorted list of column indices of non-zero matrix element
1567 : !> \param packed_index original order of non-sorted column indices
1568 : !> \param nl_repl replicated neighbour list
1569 : !> \param matrix_dbcsr_kp DBCSR matrix
1570 : !> \param symmetric whether the DBCSR matrix is a symmetric one
1571 : !> \param para_env parallel environment
1572 : !> \param gather_root if >=0, gather all non-zero matrix element on the MPI process with
1573 : !> gather_root rank (useful for bulk transport calculation).
1574 : !> If <0, distribute non-zero matrix element across all MPI processes
1575 : ! **************************************************************************************************
1576 : SUBROUTINE get_nonzero_element_indices(n_nonzero_cols, row_offset, col_index, packed_index, &
1577 : nl_repl, matrix_dbcsr_kp, symmetric, para_env, gather_root)
1578 : INTEGER, DIMENSION(:), INTENT(out) :: n_nonzero_cols, row_offset, col_index, &
1579 : packed_index
1580 8 : INTEGER, DIMENSION(:, :), INTENT(in) :: nl_repl
1581 4 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
1582 : LOGICAL, INTENT(in) :: symmetric
1583 : TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
1584 : INTEGER, INTENT(in) :: gather_root
1585 :
1586 : CHARACTER(len=*), PARAMETER :: routineN = 'get_nonzero_element_indices'
1587 :
1588 : INTEGER :: first_col_minus_one, first_row_minus_one, handle, icol_blk, icol_local, &
1589 : icol_offset, image, image_transp, inode, irow_blk, irow_local, irow_proc, mepos, &
1590 : ncols_blk, ncols_local, ncols_total, nnodes, nprocs, nrows_blk, nrows_local, nrows_total
1591 : INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
1592 : row_blk_offset, row_blk_size
1593 : LOGICAL :: do_distribute, is_root_rank
1594 :
1595 4 : CALL timeset(routineN, handle)
1596 4 : n_nonzero_cols(:) = 0
1597 : mepos = para_env%mepos
1598 : nprocs = para_env%num_pe
1599 4 : do_distribute = gather_root < 0
1600 618 : is_root_rank = gather_root == mepos
1601 4 :
1602 4 : CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
1603 4 : nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
1604 4 : nfullrows_total=nrows_total, nfullcols_total=ncols_total, &
1605 : row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
1606 : row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
1607 :
1608 : nnodes = SIZE(nl_repl, 2)
1609 : DO inode = 1, nnodes
1610 4 : irow_blk = nl_repl(neighbor_list_iatom_index, inode)
1611 : icol_blk = nl_repl(neighbor_list_jatom_index, inode)
1612 4 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
1613 29780 : image = nl_repl(neighbor_list_siesta_image_index, inode)
1614 29776 : image_transp = nl_repl(neighbor_list_siesta_transp_image_index, inode)
1615 29776 :
1616 29776 : IF (image > 0) THEN
1617 29776 : nrows_local = row_blk_size(irow_blk)
1618 29776 : first_row_minus_one = row_blk_offset(irow_blk) - 1
1619 : ncols_local = col_blk_size(icol_blk)
1620 29776 : DO irow_local = 1, nrows_local
1621 29776 : IF (do_distribute) THEN
1622 29776 : #if defined(__SMEAGOL)
1623 29776 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
1624 297760 : #else
1625 267984 : CALL cp_abort(__LOCATION__, &
1626 : "CP2K was compiled with no SMEAGOL support.")
1627 0 : #endif
1628 : ELSE
1629 : IF (is_root_rank) THEN
1630 : irow_proc = irow_local + first_row_minus_one
1631 : ELSE
1632 : irow_proc = 0
1633 267984 : END IF
1634 133992 : END IF
1635 : IF (irow_proc > 0) THEN
1636 133992 : n_nonzero_cols(irow_proc) = n_nonzero_cols(irow_proc) + ncols_local
1637 : END IF
1638 : END DO
1639 297760 : END IF
1640 133992 :
1641 : ! transposed block
1642 : IF (image_transp > 0) THEN
1643 : nrows_local = col_blk_size(icol_blk)
1644 : first_row_minus_one = col_blk_offset(icol_blk) - 1
1645 : ncols_local = row_blk_size(irow_blk)
1646 29780 : DO irow_local = 1, nrows_local
1647 27480 : IF (do_distribute) THEN
1648 27480 : #if defined(__SMEAGOL)
1649 27480 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
1650 274800 : #else
1651 247320 : CALL cp_abort(__LOCATION__, &
1652 : "CP2K was compiled with no SMEAGOL support.")
1653 0 : #endif
1654 : ELSE
1655 : IF (is_root_rank) THEN
1656 : irow_proc = irow_local + first_row_minus_one
1657 : ELSE
1658 : irow_proc = 0
1659 247320 : END IF
1660 123660 : END IF
1661 : IF (irow_proc > 0) THEN
1662 123660 : n_nonzero_cols(irow_proc) = n_nonzero_cols(irow_proc) + ncols_local
1663 : END IF
1664 : END DO
1665 274800 : END IF
1666 123660 : END DO
1667 :
1668 : row_offset(1) = 0
1669 : DO irow_local = 1, SIZE(n_nonzero_cols) - 1
1670 : row_offset(irow_local + 1) = row_offset(irow_local) + n_nonzero_cols(irow_local)
1671 : END DO
1672 4 :
1673 614 : n_nonzero_cols(:) = 0
1674 614 : col_index(:) = 0
1675 : DO inode = 1, nnodes
1676 : irow_blk = nl_repl(neighbor_list_iatom_index, inode)
1677 618 : icol_blk = nl_repl(neighbor_list_jatom_index, inode)
1678 2318874 : CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
1679 29780 : image = nl_repl(neighbor_list_siesta_image_index, inode)
1680 29776 : image_transp = nl_repl(neighbor_list_siesta_transp_image_index, inode)
1681 29776 :
1682 29776 : IF (image > 0) THEN
1683 29776 : nrows_local = row_blk_size(irow_blk)
1684 29776 : first_row_minus_one = row_blk_offset(irow_blk) - 1
1685 : ncols_local = col_blk_size(icol_blk)
1686 29776 : first_col_minus_one = col_blk_offset(icol_blk) + (image - 1)*ncols_total - 1
1687 29776 : DO irow_local = 1, nrows_local
1688 29776 : IF (do_distribute) THEN
1689 29776 : #if defined(__SMEAGOL)
1690 29776 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
1691 297760 : #else
1692 267984 : CALL cp_abort(__LOCATION__, &
1693 : "CP2K was compiled with no SMEAGOL support.")
1694 0 : #endif
1695 : ELSE
1696 : IF (is_root_rank) THEN
1697 : irow_proc = irow_local + first_row_minus_one
1698 : ELSE
1699 : irow_proc = 0
1700 267984 : END IF
1701 133992 : END IF
1702 : IF (irow_proc > 0) THEN
1703 133992 : icol_offset = row_offset(irow_proc) + n_nonzero_cols(irow_proc)
1704 : DO icol_local = 1, ncols_local
1705 : col_index(icol_offset + icol_local) = first_col_minus_one + icol_local
1706 297760 : END DO
1707 133992 : n_nonzero_cols(irow_proc) = n_nonzero_cols(irow_proc) + ncols_local
1708 1339920 : END IF
1709 1339920 : END DO
1710 : END IF
1711 133992 :
1712 : ! transposed block
1713 : IF (image_transp > 0) THEN
1714 : nrows_local = col_blk_size(icol_blk)
1715 : first_row_minus_one = col_blk_offset(icol_blk) - 1
1716 : ncols_local = row_blk_size(irow_blk)
1717 29780 : first_col_minus_one = row_blk_offset(irow_blk) + (image_transp - 1)*nrows_total - 1
1718 27480 : DO irow_local = 1, nrows_local
1719 27480 : IF (do_distribute) THEN
1720 27480 : #if defined(__SMEAGOL)
1721 27480 : CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
1722 274800 : #else
1723 247320 : CALL cp_abort(__LOCATION__, &
1724 : "CP2K was compiled with no SMEAGOL support.")
1725 0 : #endif
1726 : ELSE
1727 : IF (is_root_rank) THEN
1728 : irow_proc = irow_local + first_row_minus_one
1729 : ELSE
1730 : irow_proc = 0
1731 247320 : END IF
1732 123660 : END IF
1733 : IF (irow_proc > 0) THEN
1734 123660 : icol_offset = row_offset(irow_proc) + n_nonzero_cols(irow_proc)
1735 : DO icol_local = 1, ncols_local
1736 : col_index(icol_offset + icol_local) = first_col_minus_one + icol_local
1737 274800 : END DO
1738 123660 : n_nonzero_cols(irow_proc) = n_nonzero_cols(irow_proc) + ncols_local
1739 1236600 : END IF
1740 1236600 : END DO
1741 : END IF
1742 123660 : END DO
1743 :
1744 : IF (SIZE(n_nonzero_cols) > 0) THEN
1745 : DO irow_local = 1, SIZE(n_nonzero_cols)
1746 : CALL sort(col_index(row_offset(irow_local) + 1:row_offset(irow_local) + n_nonzero_cols(irow_local)), &
1747 : n_nonzero_cols(irow_local), &
1748 4 : packed_index(row_offset(irow_local) + 1:row_offset(irow_local) + n_nonzero_cols(irow_local)))
1749 618 : END DO
1750 : END IF
1751 :
1752 618 : CALL timestop(handle)
1753 : END SUBROUTINE get_nonzero_element_indices
1754 :
1755 : ! **************************************************************************************************
1756 4 : !> \brief Get absolute i, j, and k indices of cell image for the given DBCSR matrix block.
1757 4 : !> Effective cell image recorded in the neighbour-list (where the matrix block is actually
1758 : !> stored) depends on atomic coordinates can be significantly different.
1759 : !> \param cell_ijk array with 3 indices along the cell's vectors
1760 : !> \param r_ij actual interatomic distance (vector R_j - r_i), where R_j is the coordinates
1761 : !> of the j-th atom in the supercell
1762 : !> \param r_i coordinates of the i-th atom in the primary unit cell
1763 : !> \param r_j coordinates of the j-th atom in the primary unit cell
1764 : !> \param cell unit cell
1765 : ! **************************************************************************************************
1766 : SUBROUTINE get_negf_cell_ijk(cell_ijk, r_ij, r_i, r_j, cell)
1767 : INTEGER, DIMENSION(3), INTENT(out) :: cell_ijk
1768 : REAL(kind=dp), DIMENSION(3), INTENT(in) :: r_ij, r_i, r_j
1769 : TYPE(cell_type), INTENT(in), POINTER :: cell
1770 30320 :
1771 : REAL(kind=dp), DIMENSION(3) :: coords_scaled, r
1772 :
1773 : r(:) = r_ij(:) + r_i(:) - r_j(:)
1774 : CALL real_to_scaled(coords_scaled, r, cell)
1775 : cell_ijk(:) = NINT(coords_scaled(:))
1776 : END SUBROUTINE get_negf_cell_ijk
1777 121280 :
1778 30320 : ! **************************************************************************************************
1779 121280 : !> \brief Return the index of an integer number in the sequence 0, 1, -1, ..., n, -n, ...
1780 30320 : !> (canonical enumeration of integers, oeis.org/A001057).
1781 : !> \param inum integer number [in]
1782 : !> \return index of 'inum' in A001057
1783 : !> \note Cell images in SMEAGOL / SIESTA are ordered according to A001057. Therefore this
1784 : !> function converts the absolute index of a cell replica along some (x/y/z) dimension
1785 : !> into its corresponding SIESTA's index.
1786 : ! **************************************************************************************************
1787 : ELEMENTAL FUNCTION index_in_canonical_enumeration(inum) RESULT(ind)
1788 : INTEGER, INTENT(in) :: inum
1789 : INTEGER :: ind
1790 :
1791 85884 : INTEGER :: inum_abs, is_non_positive
1792 :
1793 : inum_abs = ABS(inum)
1794 : !IF (inum <= 0) THEN; is_non_positive = 1; ELSE; is_non_positive = 0; END IF
1795 : is_non_positive = MERGE(1, 0, inum <= 0)
1796 :
1797 85884 : ! inum = 0 -> inum_abs = 0, is_non_positive = 1 -> ind = 1
1798 : ! inum = 1 -> inum_abs = 1, is_non_positive = 0 -> ind = 2
1799 85884 : ! inum = -1 -> inum_abs = 1, is_non_positive = 1 -> ind = 3
1800 : ind = 2*inum_abs + is_non_positive
1801 : END FUNCTION index_in_canonical_enumeration
1802 :
1803 : ! **************************************************************************************************
1804 85884 : !> \brief Return an integer number according to its index in the sequence 0, 1, -1, ..., n, -n, ...
1805 85884 : !> (canonical enumeration of integers, oeis.org/A001057)
1806 : !> \param ind index in A001057 starting from 1
1807 : !> \return integer number according to its position 'ind' in A001057
1808 : !> \note Cell images in SMEAGOL / SIESTA are ordered according to A001057. Therefore this
1809 : !> function converts SIESTA's index of a cell replica along some (x/y/z) dimension
1810 : !> into the corresponding absolute index.
1811 : ! **************************************************************************************************
1812 : ELEMENTAL FUNCTION number_from_canonical_enumeration(ind) RESULT(inum)
1813 : INTEGER, INTENT(in) :: ind
1814 : INTEGER :: inum
1815 :
1816 372 : ! ind < 1 is invalid
1817 : ! ind = 1 -> SIGN(0, -1) = 0
1818 : ! ind = 2 -> SIGN(1, 0) = 1
1819 : ! ind = 3 -> SIGN(1, -1) = -1
1820 : inum = SIGN(ind/2, -MOD(ind, 2))
1821 : END FUNCTION number_from_canonical_enumeration
1822 :
1823 : ! **************************************************************************************************
1824 372 : !> \brief Apply periodic boundary conditions defined by a simulation cell to a position
1825 372 : !> vector r. Similar to pbc1 from cell_types.F but returns unscaled coordinates from
1826 : !> the scaled range [0, 1) instead of [-0.5, 0.5)
1827 : !> \param r_pbc position vector subject to the periodic boundary conditions [out]
1828 : !> \param r initial position vector [in]
1829 : !> \param cell simulation unit cell [in]
1830 : ! **************************************************************************************************
1831 : PURE SUBROUTINE pbc_0_1(r_pbc, r, cell)
1832 : REAL(KIND=dp), DIMENSION(3), INTENT(out) :: r_pbc
1833 : REAL(KIND=dp), DIMENSION(3), INTENT(in) :: r
1834 : TYPE(cell_type), INTENT(in), POINTER :: cell
1835 136 :
1836 : REAL(KIND=dp), DIMENSION(3) :: s
1837 :
1838 : IF (cell%orthorhombic) THEN
1839 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*REAL(FLOOR(cell%h_inv(1, 1)*r(1)), dp)
1840 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*REAL(FLOOR(cell%h_inv(2, 2)*r(2)), dp)
1841 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*REAL(FLOOR(cell%h_inv(3, 3)*r(3)), dp)
1842 136 : ELSE
1843 136 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
1844 136 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
1845 136 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
1846 : s(1) = s(1) - cell%perd(1)*REAL(FLOOR(s(1)), dp)
1847 0 : s(2) = s(2) - cell%perd(2)*REAL(FLOOR(s(2)), dp)
1848 0 : s(3) = s(3) - cell%perd(3)*REAL(FLOOR(s(3)), dp)
1849 0 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
1850 0 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
1851 0 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
1852 0 : END IF
1853 0 : END SUBROUTINE pbc_0_1
1854 0 :
1855 0 : ! **************************************************************************************************
1856 : !> \brief Computes the number of send requests from this MPI process to all the other processes.
1857 136 : !> Alternatively computes the number of recv requests per MPI process that the given process
1858 : !> expects.
1859 : !> \param mepos MPI rank of the given process
1860 : !> \param nelements_per_proc number of element to send / receive
1861 : !> \param max_nelements_per_packet maximum number of elements per single MPI request
1862 : !> \return number of MPI requests
1863 : ! **************************************************************************************************
1864 : PURE FUNCTION get_number_of_mpi_sendrecv_requests(mepos, nelements_per_proc, max_nelements_per_packet) RESULT(nrequests)
1865 : INTEGER, INTENT(in) :: mepos
1866 : INTEGER(kind=int_8), DIMENSION(0:), INTENT(in) :: nelements_per_proc
1867 : INTEGER(kind=int_8), INTENT(in) :: max_nelements_per_packet
1868 32 : INTEGER :: nrequests
1869 :
1870 : INTEGER :: iproc
1871 :
1872 : nrequests = 0
1873 : DO iproc = LBOUND(nelements_per_proc, 1), UBOUND(nelements_per_proc, 1)
1874 : ! there is no need to send data to the same MPI process
1875 : IF (iproc /= mepos) THEN
1876 32 : nrequests = nrequests + INT(nelements_per_proc(iproc)/max_nelements_per_packet)
1877 128 : IF (MOD(nelements_per_proc(iproc), max_nelements_per_packet) > 0) THEN
1878 : nrequests = nrequests + 1
1879 96 : END IF
1880 32 : END IF
1881 32 : END DO
1882 16 : END FUNCTION get_number_of_mpi_sendrecv_requests
1883 :
1884 : ! **************************************************************************************************
1885 : !> \brief Map non-zero matrix elements on to MPI requests.
1886 32 : !> \param element_offset offset (index-1) of the first element [out]
1887 : !> \param nelements_per_request number of element for each request [out]
1888 : !> \param peer_rank rank of a peering MPI process
1889 : !> \param tag MPI tag
1890 : !> \param mepos MPI rank of a given MPI process
1891 : !> \param nelements_per_proc number of element to send / receive by the current MPI process
1892 : !> \param max_nelements_per_packet maximum number of elements per single MPI request
1893 : ! **************************************************************************************************
1894 : SUBROUTINE assign_nonzero_elements_to_requests(element_offset, nelements_per_request, peer_rank, tag, &
1895 : mepos, nelements_per_proc, max_nelements_per_packet)
1896 : INTEGER(kind=int_8), DIMENSION(:), INTENT(out) :: element_offset, nelements_per_request
1897 : INTEGER, DIMENSION(:), INTENT(out) :: peer_rank, tag
1898 16 : INTEGER, INTENT(in) :: mepos
1899 16 : INTEGER(kind=int_8), DIMENSION(0:), INTENT(in) :: nelements_per_proc
1900 : INTEGER(kind=int_8), INTENT(in) :: max_nelements_per_packet
1901 :
1902 : INTEGER :: iproc, irequest, nrequests, &
1903 : request_offset
1904 : INTEGER(kind=int_8) :: element_offset_tmp, nelements
1905 :
1906 : request_offset = 0
1907 : element_offset_tmp = 0
1908 : DO iproc = LBOUND(nelements_per_proc, 1), UBOUND(nelements_per_proc, 1)
1909 : IF (iproc /= mepos) THEN
1910 16 : nrequests = INT(nelements_per_proc(iproc)/max_nelements_per_packet)
1911 16 : IF (MOD(nelements_per_proc(iproc), max_nelements_per_packet) > 0) nrequests = nrequests + 1
1912 64 : CPASSERT(nrequests <= max_mpi_rank + 1)
1913 32 : IF (nrequests > 0) THEN
1914 16 : nelements = nelements_per_proc(iproc)/nrequests
1915 16 : IF (nelements_per_proc(iproc) - nelements*nrequests > 0) nelements = nelements + 1
1916 16 : CPASSERT(nelements <= max_nelements_per_packet)
1917 16 :
1918 16 : DO irequest = 1, nrequests
1919 16 : element_offset(request_offset + irequest) = (irequest - 1)*nelements + element_offset_tmp
1920 16 : IF (irequest < nrequests) THEN
1921 : nelements_per_request(request_offset + irequest) = nelements
1922 32 : ELSE
1923 16 : nelements_per_request(request_offset + irequest) = nelements_per_proc(iproc) - nelements*(nrequests - 1)
1924 16 : END IF
1925 0 : peer_rank(request_offset + irequest) = iproc
1926 : tag(request_offset + irequest) = irequest - 1
1927 16 : END DO
1928 : END IF
1929 16 : request_offset = request_offset + nrequests
1930 32 : END IF
1931 : element_offset_tmp = element_offset_tmp + nelements_per_proc(iproc)
1932 : END DO
1933 16 :
1934 : IF (debug_this_module) THEN
1935 48 : CPASSERT(SIZE(element_offset) == request_offset)
1936 : CPASSERT(SIZE(nelements_per_request) == request_offset)
1937 : CPASSERT(SIZE(peer_rank) == request_offset)
1938 : CPASSERT(SIZE(tag) == request_offset)
1939 : END IF
1940 : END SUBROUTINE assign_nonzero_elements_to_requests
1941 :
1942 : END MODULE smeagol_matrix_utils
1943 : /* /opt/cp2k/src/smeagol_matrix_utils.F not long enough */
1944 16 : /* END: function "__smeagol_matrix_utils_MOD_assign_nonzero_elements_to_requests" */
1945 : /* ... */
1946 0 : /* END: function "__smeagol_matrix_utils_MOD___copy_26958B" */
|
