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 represent a full matrix distributed on many processors
10 : !> \par History
11 : !> 3) separated structure object, removed globenv, renamed to full matrix
12 : !> many changes (fawzi 08.2002)
13 : !> \author Matthias Krack (22.05.2001)
14 : ! **************************************************************************************************
15 : MODULE cp_fm_types
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_blacs_types, ONLY: cp_blacs_type
18 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
19 : cp_fm_struct_equivalent,&
20 : cp_fm_struct_get,&
21 : cp_fm_struct_release,&
22 : cp_fm_struct_retain,&
23 : cp_fm_struct_type,&
24 : cp_fm_struct_write_info
25 : USE kinds, ONLY: dp
26 : USE message_passing, ONLY: cp2k_is_parallel,&
27 : mp_any_source,&
28 : mp_para_env_type,&
29 : mp_proc_null,&
30 : mp_request_null,&
31 : mp_request_type,&
32 : mp_waitall
33 : USE parallel_rng_types, ONLY: UNIFORM,&
34 : rng_stream_type
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_types'
42 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
43 : INTEGER, PARAMETER :: src_tag = 3, dest_tag = 5, send_tag = 7, recv_tag = 11
44 :
45 : INTEGER, PRIVATE :: cp_fm_mm_type = 1
46 :
47 : PUBLIC :: cp_fm_type, &
48 : cp_fm_p_type, copy_info_type
49 :
50 : PUBLIC :: cp_fm_add_to_element, &
51 : cp_fm_create, &
52 : cp_fm_release, &
53 : cp_fm_get_info, &
54 : cp_fm_set_element, &
55 : cp_fm_get_element, &
56 : cp_fm_get_diag, & ! get diagonal
57 : cp_fm_set_all, & ! set all elements and diagonal
58 : cp_fm_set_all_submatrix, & ! set a submatrix to a given value
59 : cp_fm_set_submatrix, & ! set a submatrix to given values
60 : cp_fm_get_submatrix, & ! get a submatrix of given values
61 : cp_fm_init_random, &
62 : cp_fm_maxabsval, & ! find the maximum absolute value
63 : cp_fm_maxabsrownorm, & ! find the maximum of the sum of the abs of the elements of a row
64 : cp_fm_to_fm, & ! copy (parts of) a fm to a fm
65 : cp_fm_vectorsnorm, & ! compute the norm of the column-vectors
66 : cp_fm_vectorssum, & ! compute the sum of all elements of the column-vectors
67 : cp_fm_to_fm_submat, & ! copy (parts of) a fm to a fm
68 : cp_fm_to_fm_triangular, &
69 : cp_fm_copy_general, &
70 : cp_fm_start_copy_general, &
71 : cp_fm_finish_copy_general, &
72 : cp_fm_cleanup_copy_general, &
73 : cp_fm_write_unformatted, & ! writes a full matrix to an open unit
74 : cp_fm_write_formatted, & ! writes a full matrix to an open unit
75 : cp_fm_read_unformatted, & ! reads a full matrix from an open unit
76 : cp_fm_setup, & ! allows to set flags for fms
77 : cp_fm_get_mm_type, &
78 : cp_fm_write_info, &
79 : cp_fm_to_fm_submat_general ! copy matrix across different contexts
80 :
81 : PUBLIC :: cp_fm_pilaenv
82 :
83 : INTERFACE cp_fm_to_fm
84 : MODULE PROCEDURE cp_fm_to_fm_matrix, & ! a full matrix
85 : cp_fm_to_fm_columns ! just a number of columns
86 : END INTERFACE
87 :
88 : INTERFACE cp_fm_release
89 : MODULE PROCEDURE cp_fm_release_aa0, &
90 : cp_fm_release_aa1, &
91 : cp_fm_release_aa2, &
92 : cp_fm_release_aa3, &
93 : cp_fm_release_ap1, &
94 : cp_fm_release_ap2, &
95 : cp_fm_release_pa1, &
96 : cp_fm_release_pa2, &
97 : cp_fm_release_pa3, &
98 : cp_fm_release_pp1, &
99 : cp_fm_release_pp2
100 : END INTERFACE
101 :
102 : ! **************************************************************************************************
103 : !> \brief represent a full matrix
104 : !> \param name the name of the matrix, used for printing
105 : !> \param matrix_struct structure of this matrix
106 : !> \param local_data array with the data of the matrix (its contents
107 : !> depend on the matrix type used: in parallel runs it will be
108 : !> in scalapack format, in sequential, it will simply contain
109 : !> the matrix)
110 : !> \par History
111 : !> 08.2002 created [fawzi]
112 : !> \author fawzi
113 : ! **************************************************************************************************
114 : TYPE cp_fm_type
115 : ! PRIVATE
116 : CHARACTER(LEN=60) :: name = ""
117 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct => NULL()
118 : REAL(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: local_data => NULL()
119 : END TYPE cp_fm_type
120 :
121 : ! **************************************************************************************************
122 : !> \brief just to build arrays of pointers to matrices
123 : !> \param matrix the pointer to the matrix
124 : !> \par History
125 : !> 08.2002 created [fawzi]
126 : !> \author fawzi
127 : ! **************************************************************************************************
128 : TYPE cp_fm_p_type
129 : TYPE(cp_fm_type), POINTER :: matrix => NULL()
130 : END TYPE cp_fm_p_type
131 :
132 : ! **************************************************************************************************
133 : !> \brief Stores the state of a copy between cp_fm_start_copy_general
134 : !> and cp_fm_finish_copy_general
135 : !> \par History
136 : !> Jan 2017 [Mark T]
137 : ! **************************************************************************************************
138 : TYPE copy_info_type
139 : INTEGER :: send_size = -1
140 : INTEGER, DIMENSION(2) :: nlocal_recv = -1, nblock_src = -1, src_num_pe = -1 ! 1->row 2->col
141 : TYPE(mp_request_type), DIMENSION(:), ALLOCATABLE :: send_request, recv_request
142 : INTEGER, DIMENSION(:), ALLOCATABLE :: recv_disp
143 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices => NULL(), recv_row_indices => NULL()
144 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: src_blacs2mpi
145 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: recv_buf, send_buf
146 : END TYPE copy_info_type
147 :
148 : CONTAINS
149 :
150 : ! **************************************************************************************************
151 : !> \brief creates a new full matrix with the given structure
152 : !> \param matrix the matrix to be created
153 : !> \param matrix_struct the structure of matrix
154 : !> \param name ...
155 : !> \param nrow ...
156 : !> \param ncol ...
157 : !> \param set_zero ...
158 : !> \par History
159 : !> 08.2002 created [fawzi]
160 : !> \author Fawzi Mohamed
161 : !> \note
162 : !> preferred allocation routine
163 : ! **************************************************************************************************
164 2020374 : SUBROUTINE cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
165 : TYPE(cp_fm_type), INTENT(OUT) :: matrix
166 : TYPE(cp_fm_struct_type), INTENT(IN), TARGET :: matrix_struct
167 : CHARACTER(LEN=*), INTENT(in), OPTIONAL :: name
168 : INTEGER, INTENT(IN), OPTIONAL :: nrow, ncol
169 : LOGICAL, INTENT(in), OPTIONAL :: set_zero
170 :
171 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_create'
172 :
173 : INTEGER :: handle, ncol_global, ncol_local, &
174 : nrow_global, nrow_local
175 : TYPE(cp_blacs_env_type), POINTER :: context
176 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
177 :
178 2020374 : CALL timeset(routineN, handle)
179 :
180 2020374 : IF (PRESENT(nrow) .OR. PRESENT(ncol)) THEN
181 30196 : CALL cp_fm_struct_get(matrix_struct, nrow_global=nrow_global, ncol_global=ncol_global)
182 30196 : IF (PRESENT(nrow)) nrow_global = nrow
183 30196 : IF (PRESENT(ncol)) ncol_global = ncol
184 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=matrix_struct, &
185 30196 : nrow_global=nrow_global, ncol_global=ncol_global)
186 :
187 30196 : context => fm_struct%context
188 30196 : matrix%matrix_struct => fm_struct
189 30196 : CALL cp_fm_struct_retain(matrix%matrix_struct)
190 :
191 30196 : nrow_local = fm_struct%local_leading_dimension
192 30196 : ncol_local = MAX(1, fm_struct%ncol_locals(context%mepos(2)))
193 :
194 30196 : CALL cp_fm_struct_release(fm_struct)
195 : ELSE
196 :
197 1990178 : context => matrix_struct%context
198 1990178 : matrix%matrix_struct => matrix_struct
199 1990178 : CALL cp_fm_struct_retain(matrix%matrix_struct)
200 :
201 : ! OK, we allocate here at least a 1 x 1 matrix
202 : ! this must (and is) compatible with the descinit call
203 : ! in cp_fm_struct
204 1990178 : nrow_local = matrix_struct%local_leading_dimension
205 1990178 : ncol_local = MAX(1, matrix_struct%ncol_locals(context%mepos(2)))
206 : END IF
207 :
208 2020374 : NULLIFY (matrix%local_data)
209 :
210 8081496 : ALLOCATE (matrix%local_data(nrow_local, ncol_local))
211 :
212 2020374 : IF (PRESENT(set_zero)) THEN
213 253299 : IF (set_zero) THEN
214 84025874 : matrix%local_data(1:nrow_local, 1:ncol_local) = 0.0_dp
215 : END IF
216 : END IF
217 :
218 2020374 : IF (PRESENT(name)) THEN
219 980572 : matrix%name = name
220 : ELSE
221 1039802 : matrix%name = 'full matrix'
222 : END IF
223 :
224 2020374 : CALL timestop(handle)
225 :
226 2020374 : END SUBROUTINE cp_fm_create
227 :
228 : ! **************************************************************************************************
229 : !> \brief releases a full matrix
230 : !> \param matrix the matrix to release
231 : !> \par History
232 : !> 08.2002 created [fawzi]
233 : !> \author Fawzi Mohamed
234 : ! **************************************************************************************************
235 2076008 : SUBROUTINE cp_fm_release_aa0(matrix)
236 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
237 :
238 2076008 : IF (ASSOCIATED(matrix%local_data)) THEN
239 2019654 : DEALLOCATE (matrix%local_data)
240 : NULLIFY (matrix%local_data)
241 : END IF
242 2076008 : matrix%name = ""
243 2076008 : CALL cp_fm_struct_release(matrix%matrix_struct)
244 :
245 2076008 : END SUBROUTINE cp_fm_release_aa0
246 :
247 : ! **************************************************************************************************
248 : !> \brief ...
249 : !> \param matrices ...
250 : ! **************************************************************************************************
251 85190 : SUBROUTINE cp_fm_release_aa1(matrices)
252 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: matrices
253 :
254 : INTEGER :: i
255 :
256 85190 : IF (ALLOCATED(matrices)) THEN
257 224871 : DO i = 1, SIZE(matrices)
258 224871 : CALL cp_fm_release(matrices(i))
259 : END DO
260 81160 : DEALLOCATE (matrices)
261 : END IF
262 85190 : END SUBROUTINE cp_fm_release_aa1
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param matrices ...
267 : ! **************************************************************************************************
268 8884 : SUBROUTINE cp_fm_release_aa2(matrices)
269 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: matrices
270 :
271 : INTEGER :: i, j
272 :
273 8884 : IF (ALLOCATED(matrices)) THEN
274 16160 : DO i = 1, SIZE(matrices, 1)
275 44570 : DO j = 1, SIZE(matrices, 2)
276 38574 : CALL cp_fm_release(matrices(i, j))
277 : END DO
278 : END DO
279 5996 : DEALLOCATE (matrices)
280 : END IF
281 8884 : END SUBROUTINE cp_fm_release_aa2
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param matrices ...
286 : ! **************************************************************************************************
287 142 : SUBROUTINE cp_fm_release_aa3(matrices)
288 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: matrices
289 :
290 : INTEGER :: i, j, k
291 :
292 142 : IF (ALLOCATED(matrices)) THEN
293 2224 : DO i = 1, SIZE(matrices, 1)
294 6416 : DO j = 1, SIZE(matrices, 2)
295 11146 : DO k = 1, SIZE(matrices, 3)
296 9064 : CALL cp_fm_release(matrices(i, j, k))
297 : END DO
298 : END DO
299 : END DO
300 142 : DEALLOCATE (matrices)
301 : END IF
302 142 : END SUBROUTINE cp_fm_release_aa3
303 :
304 : ! **************************************************************************************************
305 : !> \brief ...
306 : !> \param matrices ...
307 : ! **************************************************************************************************
308 216125 : SUBROUTINE cp_fm_release_pa1(matrices)
309 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: matrices
310 :
311 : INTEGER :: i
312 :
313 216125 : IF (ASSOCIATED(matrices)) THEN
314 151898 : DO i = 1, SIZE(matrices)
315 151898 : CALL cp_fm_release(matrices(i))
316 : END DO
317 62111 : DEALLOCATE (matrices)
318 : NULLIFY (matrices)
319 : END IF
320 216125 : END SUBROUTINE cp_fm_release_pa1
321 :
322 : ! **************************************************************************************************
323 : !> \brief ...
324 : !> \param matrices ...
325 : ! **************************************************************************************************
326 69606 : SUBROUTINE cp_fm_release_pa2(matrices)
327 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: matrices
328 :
329 : INTEGER :: i, j
330 :
331 69606 : IF (ASSOCIATED(matrices)) THEN
332 86671 : DO i = 1, SIZE(matrices, 1)
333 161491 : DO j = 1, SIZE(matrices, 2)
334 137284 : CALL cp_fm_release(matrices(i, j))
335 : END DO
336 : END DO
337 24207 : DEALLOCATE (matrices)
338 : NULLIFY (matrices)
339 : END IF
340 69606 : END SUBROUTINE cp_fm_release_pa2
341 :
342 : ! **************************************************************************************************
343 : !> \brief ...
344 : !> \param matrices ...
345 : ! **************************************************************************************************
346 0 : SUBROUTINE cp_fm_release_pa3(matrices)
347 : TYPE(cp_fm_type), DIMENSION(:, :, :), POINTER :: matrices
348 :
349 : INTEGER :: i, j, k
350 :
351 0 : IF (ASSOCIATED(matrices)) THEN
352 0 : DO i = 1, SIZE(matrices, 1)
353 0 : DO j = 1, SIZE(matrices, 2)
354 0 : DO k = 1, SIZE(matrices, 3)
355 0 : CALL cp_fm_release(matrices(i, j, k))
356 : END DO
357 : END DO
358 : END DO
359 0 : DEALLOCATE (matrices)
360 : NULLIFY (matrices)
361 : END IF
362 0 : END SUBROUTINE cp_fm_release_pa3
363 :
364 : ! **************************************************************************************************
365 : !> \brief ...
366 : !> \param matrices ...
367 : ! **************************************************************************************************
368 0 : SUBROUTINE cp_fm_release_ap1(matrices)
369 : TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: matrices
370 :
371 : INTEGER :: i
372 :
373 0 : IF (ALLOCATED(matrices)) THEN
374 0 : DO i = 1, SIZE(matrices)
375 0 : CALL cp_fm_release(matrices(i)%matrix)
376 0 : DEALLOCATE (matrices(i)%matrix)
377 : END DO
378 0 : DEALLOCATE (matrices)
379 : END IF
380 0 : END SUBROUTINE cp_fm_release_ap1
381 :
382 : ! **************************************************************************************************
383 : !> \brief ...
384 : !> \param matrices ...
385 : ! **************************************************************************************************
386 0 : SUBROUTINE cp_fm_release_ap2(matrices)
387 : TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: matrices
388 :
389 : INTEGER :: i, j
390 :
391 0 : IF (ALLOCATED(matrices)) THEN
392 0 : DO i = 1, SIZE(matrices, 1)
393 0 : DO j = 1, SIZE(matrices, 2)
394 0 : CALL cp_fm_release(matrices(i, j)%matrix)
395 0 : DEALLOCATE (matrices(i, j)%matrix)
396 : END DO
397 : END DO
398 0 : DEALLOCATE (matrices)
399 : END IF
400 0 : END SUBROUTINE cp_fm_release_ap2
401 :
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param matrices ...
405 : ! **************************************************************************************************
406 0 : SUBROUTINE cp_fm_release_pp1(matrices)
407 : TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: matrices
408 :
409 : INTEGER :: i
410 :
411 0 : IF (ASSOCIATED(matrices)) THEN
412 0 : DO i = 1, SIZE(matrices)
413 0 : CALL cp_fm_release(matrices(i)%matrix)
414 0 : DEALLOCATE (matrices(i)%matrix)
415 : END DO
416 0 : DEALLOCATE (matrices)
417 : NULLIFY (matrices)
418 : END IF
419 0 : END SUBROUTINE cp_fm_release_pp1
420 :
421 : ! **************************************************************************************************
422 : !> \brief ...
423 : !> \param matrices ...
424 : ! **************************************************************************************************
425 0 : SUBROUTINE cp_fm_release_pp2(matrices)
426 : TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: matrices
427 :
428 : INTEGER :: i, j
429 :
430 0 : IF (ASSOCIATED(matrices)) THEN
431 0 : DO i = 1, SIZE(matrices, 1)
432 0 : DO j = 1, SIZE(matrices, 2)
433 0 : CALL cp_fm_release(matrices(i, j)%matrix)
434 0 : DEALLOCATE (matrices(i, j)%matrix)
435 : END DO
436 : END DO
437 0 : DEALLOCATE (matrices)
438 : NULLIFY (matrices)
439 : END IF
440 0 : END SUBROUTINE cp_fm_release_pp2
441 :
442 : ! **************************************************************************************************
443 : !> \brief fills a matrix with random numbers
444 : !> \param matrix : to be initialized
445 : !> \param ncol : numbers of cols to fill
446 : !> \param start_col : starting at coll number
447 : !> \author Joost VandeVondele
448 : !> \note
449 : !> the value of a_ij is independent of the number of cpus
450 : ! **************************************************************************************************
451 3240 : SUBROUTINE cp_fm_init_random(matrix, ncol, start_col)
452 : TYPE(cp_fm_type), INTENT(IN) :: matrix
453 : INTEGER, INTENT(IN), OPTIONAL :: ncol, start_col
454 :
455 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_init_random'
456 :
457 : INTEGER :: handle, icol_global, icol_local, irow_local, my_ncol, my_start_col, ncol_global, &
458 : ncol_local, nrow_global, nrow_local
459 6480 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
460 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buff
461 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
462 3240 : POINTER :: local_data
463 : REAL(KIND=dp), DIMENSION(3, 2), SAVE :: &
464 : seed = RESHAPE([1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp, 6.0_dp], [3, 2])
465 : TYPE(rng_stream_type) :: rng
466 :
467 3240 : CALL timeset(routineN, handle)
468 :
469 : ! guarantee same seed over all tasks
470 3240 : CALL matrix%matrix_struct%para_env%bcast(seed, 0)
471 :
472 : rng = rng_stream_type("cp_fm_init_random_stream", distribution_type=UNIFORM, &
473 3240 : extended_precision=.TRUE., seed=seed)
474 :
475 : CALL cp_fm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global, &
476 : nrow_local=nrow_local, ncol_local=ncol_local, &
477 : local_data=local_data, &
478 3240 : row_indices=row_indices, col_indices=col_indices)
479 :
480 3240 : my_start_col = 1
481 3240 : IF (PRESENT(start_col)) my_start_col = start_col
482 3240 : my_ncol = matrix%matrix_struct%ncol_global
483 3240 : IF (PRESENT(ncol)) my_ncol = ncol
484 :
485 3240 : IF (ncol_global < (my_start_col + my_ncol - 1)) THEN
486 0 : CPABORT("ncol_global>=(my_start_col+my_ncol-1)")
487 : END IF
488 :
489 9720 : ALLOCATE (buff(nrow_global))
490 :
491 : ! each global row has its own substream, in order to reach the stream for the local col,
492 : ! we just reset to the next substream
493 : ! following this, we fill the full buff with random numbers, and pick those we need
494 3240 : icol_global = 0
495 24571 : DO icol_local = 1, ncol_local
496 21331 : CPASSERT(col_indices(icol_local) > icol_global)
497 : DO
498 21331 : CALL rng%reset_to_next_substream()
499 21331 : icol_global = icol_global + 1
500 21331 : IF (icol_global == col_indices(icol_local)) EXIT
501 : END DO
502 21331 : CALL rng%fill(buff)
503 741834 : DO irow_local = 1, nrow_local
504 738594 : local_data(irow_local, icol_local) = buff(row_indices(irow_local))
505 : END DO
506 : END DO
507 :
508 3240 : DEALLOCATE (buff)
509 :
510 : ! store seed before deletion (unclear if this is the proper seed)
511 :
512 : ! Note that, the initial state (rng%ig) instead of the current state (rng%cg) is stored in the
513 : ! seed variable. As a consequence, each invocation of cp_fm_init_random uses exactly the same
514 : ! stream of random numbers. While this seems odd and contrary to the original design,
515 : ! it was probably introduced to improve reproducibility.
516 : ! See also https://github.com/cp2k/cp2k/pull/506
517 3240 : CALL rng%get(ig=seed)
518 :
519 3240 : CALL timestop(handle)
520 :
521 90720 : END SUBROUTINE cp_fm_init_random
522 :
523 : ! **************************************************************************************************
524 : !> \brief set all elements of a matrix to the same value,
525 : !> and optionally the diagonal to a different one
526 : !> \param matrix input matrix
527 : !> \param alpha scalar used to set all elements of the matrix
528 : !> \param beta scalar used to set diagonal of the matrix
529 : !> \note
530 : !> can be used to zero a matrix
531 : !> can be used to create a unit matrix (I-matrix) alpha=0.0_dp beta=1.0_dp
532 : ! **************************************************************************************************
533 404188 : SUBROUTINE cp_fm_set_all(matrix, alpha, beta)
534 :
535 : TYPE(cp_fm_type), INTENT(IN) :: matrix
536 : REAL(KIND=dp), INTENT(IN) :: alpha
537 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: beta
538 :
539 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_set_all'
540 :
541 : INTEGER :: handle, i, n
542 :
543 404188 : CALL timeset(routineN, handle)
544 :
545 226647393 : matrix%local_data(:, :) = alpha
546 :
547 404188 : IF (PRESENT(beta)) THEN
548 67228 : n = MIN(matrix%matrix_struct%nrow_global, matrix%matrix_struct%ncol_global)
549 520086 : DO i = 1, n
550 520086 : CALL cp_fm_set_element(matrix, i, i, beta)
551 : END DO
552 : END IF
553 :
554 404188 : CALL timestop(handle)
555 :
556 404188 : END SUBROUTINE cp_fm_set_all
557 :
558 : ! **************************************************************************************************
559 : !> \brief returns the diagonal elements of a fm
560 : !> \param matrix ...
561 : !> \param diag ...
562 : !> \param n_active optional number of leading diagonal elements to return
563 : ! **************************************************************************************************
564 35828 : SUBROUTINE cp_fm_get_diag(matrix, diag, n_active)
565 :
566 : TYPE(cp_fm_type), INTENT(IN) :: matrix
567 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
568 : INTEGER, INTENT(IN), OPTIONAL :: n_active
569 :
570 : INTEGER :: i_global, n_diag, ncol_global, nrow_global
571 : #if defined(__parallel)
572 : INTEGER :: icol_local, ipcol, iprow, irow, irow_local, &
573 : mypcol, myprow, npcol, nprow, nrow_local
574 : INTEGER, DIMENSION(9) :: desca
575 35828 : INTEGER, DIMENSION(:), POINTER :: row_indices
576 : TYPE(cp_blacs_env_type), POINTER :: context
577 : #endif
578 :
579 35828 : CALL cp_fm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
580 :
581 35828 : n_diag = MIN(nrow_global, ncol_global)
582 35828 : IF (PRESENT(n_active)) n_diag = n_active
583 35828 : CPASSERT(n_diag >= 0 .AND. n_diag <= MIN(nrow_global, ncol_global))
584 35828 : CPASSERT(SIZE(diag) >= n_diag)
585 :
586 2249574 : diag = 0.0_dp
587 : #if defined(__parallel)
588 35828 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, row_indices=row_indices)
589 35828 : context => matrix%matrix_struct%context
590 35828 : myprow = context%mepos(1)
591 35828 : mypcol = context%mepos(2)
592 35828 : nprow = context%num_pe(1)
593 35828 : npcol = context%num_pe(2)
594 :
595 358280 : desca(:) = matrix%matrix_struct%descriptor(:)
596 :
597 1196200 : DO irow = 1, nrow_local
598 1160372 : i_global = row_indices(irow)
599 1160372 : IF (i_global > n_diag) EXIT
600 : CALL infog2l(i_global, i_global, desca, nprow, npcol, myprow, mypcol, &
601 1160372 : irow_local, icol_local, iprow, ipcol)
602 1196200 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
603 1160372 : diag(i_global) = matrix%local_data(irow_local, icol_local)
604 : END IF
605 : END DO
606 : #else
607 : DO i_global = 1, n_diag
608 : diag(i_global) = matrix%local_data(i_global, i_global)
609 : END DO
610 : #endif
611 4463320 : CALL matrix%matrix_struct%para_env%sum(diag)
612 :
613 35828 : END SUBROUTINE cp_fm_get_diag
614 :
615 : ! **************************************************************************************************
616 : !> \brief returns an element of a fm
617 : !> this value is valid on every cpu
618 : !> using this call is expensive
619 : !> \param matrix the matrix to read
620 : !> \param irow_global the row
621 : !> \param icol_global the col
622 : !> \param alpha the value of matrix(irow_global, icol_global)
623 : !> \param local true if the element is on this cpu, false otherwise
624 : !> \note
625 : !> - modified semantics. now this function always returns the value
626 : !> previously the value was zero on cpus that didn't own the relevant
627 : !> part of the matrix (Joost VandeVondele, May 2003)
628 : !> - usage of the function should be avoided, as it is likely to rather slow
629 : !> using row_indices/col_indices/local_data + some smart scheme normally
630 : !> yields a real parallel code
631 : ! **************************************************************************************************
632 1095074 : SUBROUTINE cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
633 :
634 : ! arguments
635 : TYPE(cp_fm_type), INTENT(IN) :: matrix
636 : REAL(KIND=dp), INTENT(OUT) :: alpha
637 : INTEGER, INTENT(IN) :: icol_global, &
638 : irow_global
639 : LOGICAL, INTENT(OUT), OPTIONAL :: local
640 :
641 : ! locals
642 : #if defined(__parallel)
643 : INTEGER, DIMENSION(9) :: desca
644 : TYPE(cp_blacs_env_type), POINTER :: context
645 : INTEGER :: icol_local, ipcol, iprow, irow_local, mypcol, myprow, npcol, &
646 : nprow
647 1095074 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
648 : #endif
649 :
650 : #if defined(__parallel)
651 1095074 : context => matrix%matrix_struct%context
652 1095074 : myprow = context%mepos(1)
653 1095074 : mypcol = context%mepos(2)
654 1095074 : nprow = context%num_pe(1)
655 1095074 : npcol = context%num_pe(2)
656 :
657 1095074 : a => matrix%local_data
658 10950740 : desca(:) = matrix%matrix_struct%descriptor(:)
659 :
660 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
661 1095074 : irow_local, icol_local, iprow, ipcol)
662 :
663 1095074 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
664 547587 : alpha = a(irow_local, icol_local)
665 547587 : CALL context%dgebs2d('A', ' ', 1, 1, alpha, 1)
666 547587 : IF (PRESENT(local)) local = .TRUE.
667 : ELSE
668 547487 : CALL context%dgebr2d('A', ' ', 1, 1, alpha, 1, iprow, ipcol)
669 547487 : IF (PRESENT(local)) local = .FALSE.
670 : END IF
671 :
672 : #else
673 : IF (PRESENT(local)) local = .TRUE.
674 : alpha = matrix%local_data(irow_global, icol_global)
675 : #endif
676 :
677 1095074 : END SUBROUTINE cp_fm_get_element
678 :
679 : ! **************************************************************************************************
680 : !> \brief sets an element of a matrix
681 : !> \param matrix ...
682 : !> \param irow_global ...
683 : !> \param icol_global ...
684 : !> \param alpha ...
685 : !> \note
686 : !> we expect all cpus to have the same arguments in the call to this function
687 : !> (otherwise one should use local_data tricks)
688 : ! **************************************************************************************************
689 621884 : SUBROUTINE cp_fm_set_element(matrix, irow_global, icol_global, alpha)
690 : TYPE(cp_fm_type), INTENT(IN) :: matrix
691 : INTEGER, INTENT(IN) :: irow_global, icol_global
692 : REAL(KIND=dp), INTENT(IN) :: alpha
693 :
694 : INTEGER :: mypcol, myprow, npcol, nprow
695 : TYPE(cp_blacs_env_type), POINTER :: context
696 : #if defined(__parallel)
697 : INTEGER :: icol_local, ipcol, iprow, &
698 : irow_local
699 : INTEGER, DIMENSION(9) :: desca
700 621884 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
701 : #endif
702 :
703 621884 : context => matrix%matrix_struct%context
704 621884 : myprow = context%mepos(1)
705 621884 : mypcol = context%mepos(2)
706 621884 : nprow = context%num_pe(1)
707 621884 : npcol = context%num_pe(2)
708 :
709 : #if defined(__parallel)
710 :
711 621884 : a => matrix%local_data
712 :
713 6218840 : desca(:) = matrix%matrix_struct%descriptor(:)
714 :
715 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
716 621884 : irow_local, icol_local, iprow, ipcol)
717 :
718 621884 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
719 312847 : a(irow_local, icol_local) = alpha
720 : END IF
721 :
722 : #else
723 :
724 : matrix%local_data(irow_global, icol_global) = alpha
725 :
726 : #endif
727 621884 : END SUBROUTINE cp_fm_set_element
728 :
729 : ! **************************************************************************************************
730 : !> \brief sets a submatrix of a full matrix
731 : !> fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
732 : !> = alpha*op(new_values)(1:n_rows,1:n_cols)+ beta
733 : !> * fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
734 : !> \param fm the full to change
735 : !> \param new_values a replicated full matrix with the new values
736 : !> \param start_row the starting row of b_matrix (defaults to 1)
737 : !> \param start_col the starting col of b_matrix (defaults to 1)
738 : !> \param n_rows the number of row to change in b (defaults to
739 : !> size(op(new_values),1))
740 : !> \param n_cols the number of columns to change in b (defaults to
741 : !> size(op(new_values),2))
742 : !> \param alpha rescaling factor for the new values (defaults to 1.0)
743 : !> \param beta rescaling factor for the old values (defaults to 0.0)
744 : !> \param transpose if new_values should be transposed: if true
745 : !> op(new_values)=new_values^T, else op(new_values)=new_values
746 : !> (defaults to false)
747 : !> \par History
748 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
749 : !> \author Fawzi Mohamed
750 : !> \note
751 : !> optimized for full column updates and alpha=1.0, beta=0.0
752 : !> the new_values need to be valid on all cpus
753 : ! **************************************************************************************************
754 90087 : SUBROUTINE cp_fm_set_submatrix(fm, new_values, start_row, &
755 : start_col, n_rows, n_cols, alpha, beta, transpose)
756 : TYPE(cp_fm_type), INTENT(IN) :: fm
757 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: new_values
758 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
759 : REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha, beta
760 : LOGICAL, INTENT(in), OPTIONAL :: transpose
761 :
762 : INTEGER :: i, i0, j, j0, ncol, ncol_block, &
763 : ncol_global, ncol_local, nrow, &
764 : nrow_block, nrow_global, nrow_local, &
765 : this_col, this_row
766 90087 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
767 : LOGICAL :: tr_a
768 : REAL(KIND=dp) :: al, be
769 90087 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
770 :
771 90087 : al = 1.0_dp; be = 0.0_dp; i0 = 1; j0 = 1; tr_a = .FALSE.
772 :
773 25048 : IF (PRESENT(alpha)) al = alpha
774 90087 : IF (PRESENT(beta)) be = beta
775 90087 : IF (PRESENT(start_row)) i0 = start_row
776 90087 : IF (PRESENT(start_col)) j0 = start_col
777 90087 : IF (PRESENT(transpose)) tr_a = transpose
778 16223 : IF (tr_a) THEN
779 16131 : nrow = SIZE(new_values, 2)
780 16131 : ncol = SIZE(new_values, 1)
781 : ELSE
782 73956 : nrow = SIZE(new_values, 1)
783 73956 : ncol = SIZE(new_values, 2)
784 : END IF
785 90087 : IF (PRESENT(n_rows)) nrow = n_rows
786 90087 : IF (PRESENT(n_cols)) ncol = n_cols
787 :
788 90087 : full_block => fm%local_data
789 :
790 : CALL cp_fm_get_info(matrix=fm, &
791 : nrow_global=nrow_global, ncol_global=ncol_global, &
792 : nrow_block=nrow_block, ncol_block=ncol_block, &
793 : nrow_local=nrow_local, ncol_local=ncol_local, &
794 90087 : row_indices=row_indices, col_indices=col_indices)
795 :
796 90087 : IF (al == 1.0_dp .AND. be == 0.0_dp) THEN
797 1434130 : DO j = 1, ncol_local
798 1359643 : this_col = col_indices(j) - j0 + 1
799 1434130 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
800 325491 : IF (tr_a) THEN
801 7595 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
802 162320 : DO i = 1, nrow_local
803 162320 : full_block(i, j) = new_values(this_col, row_indices(i))
804 : END DO
805 : ELSE
806 11840 : DO i = 1, nrow_local
807 10968 : this_row = row_indices(i) - i0 + 1
808 11840 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
809 1750 : full_block(i, j) = new_values(this_col, this_row)
810 : END IF
811 : END DO
812 : END IF
813 : ELSE
814 317896 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
815 6875427 : DO i = 1, nrow_local
816 6875427 : full_block(i, j) = new_values(row_indices(i), this_col)
817 : END DO
818 : ELSE
819 9070 : DO i = 1, nrow_local
820 8154 : this_row = row_indices(i) - i0 + 1
821 9070 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
822 2703 : full_block(i, j) = new_values(this_row, this_col)
823 : END IF
824 : END DO
825 : END IF
826 : END IF
827 : END IF
828 : END DO
829 : ELSE
830 873320 : DO j = 1, ncol_local
831 857720 : this_col = col_indices(j) - j0 + 1
832 873320 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
833 857720 : IF (tr_a) THEN
834 88195316 : DO i = 1, nrow_local
835 87373812 : this_row = row_indices(i) - i0 + 1
836 88195316 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
837 : full_block(i, j) = al*new_values(this_col, this_row) + &
838 410752 : be*full_block(i, j)
839 : END IF
840 : END DO
841 : ELSE
842 140364 : DO i = 1, nrow_local
843 104148 : this_row = row_indices(i) - i0 + 1
844 140364 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
845 : full_block(i, j) = al*new_values(this_row, this_col) + &
846 104148 : be*full_block(i, j)
847 : END IF
848 : END DO
849 : END IF
850 : END IF
851 : END DO
852 : END IF
853 :
854 90087 : END SUBROUTINE cp_fm_set_submatrix
855 :
856 : ! **************************************************************************************************
857 : !> \brief sets a submatrix of a full matrix to a given value
858 : !> fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = value
859 : !> \param fm the full to change
860 : !> \param new_value ...
861 : !> \param start_row the starting row of matrix
862 : !> \param start_col the starting col of matrix
863 : !> \param n_rows the number of rows to change
864 : !> \param n_cols the number of columns to change
865 : !> \par History
866 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
867 : !> 12.2025 created from cp_fm_set_submatrix
868 : !> \author JGH
869 : ! **************************************************************************************************
870 1053455 : SUBROUTINE cp_fm_set_all_submatrix(fm, new_value, start_row, start_col, n_rows, n_cols)
871 : TYPE(cp_fm_type), INTENT(IN) :: fm
872 : REAL(KIND=dp), INTENT(in) :: new_value
873 : INTEGER, INTENT(in) :: start_row, start_col, n_rows, n_cols
874 :
875 : INTEGER :: i, i0, j, j0, ncol_global, ncol_local, &
876 : nrow_global, nrow_local, this_col, &
877 : this_row
878 1053455 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
879 1053455 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
880 :
881 1053455 : full_block => fm%local_data
882 :
883 : CALL cp_fm_get_info(matrix=fm, &
884 : nrow_global=nrow_global, ncol_global=ncol_global, &
885 : nrow_local=nrow_local, ncol_local=ncol_local, &
886 1053455 : row_indices=row_indices, col_indices=col_indices)
887 :
888 1053455 : i0 = start_row
889 1053455 : j0 = start_col
890 20913910 : DO j = 1, ncol_local
891 19860455 : this_col = col_indices(j) - j0 + 1
892 20913910 : IF (this_col >= 1 .AND. this_col <= n_cols) THEN
893 699343653 : DO i = 1, nrow_local
894 679543884 : this_row = row_indices(i) - i0 + 1
895 699343653 : IF (this_row >= 1 .AND. this_row <= n_rows) THEN
896 672474778 : full_block(i, j) = new_value
897 : END IF
898 : END DO
899 : END IF
900 : END DO
901 :
902 1053455 : END SUBROUTINE cp_fm_set_all_submatrix
903 :
904 : ! **************************************************************************************************
905 : !> \brief gets a submatrix of a full matrix
906 : !> op(target_m)(1:n_rows,1:n_cols)
907 : !> =fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
908 : !> target_m is replicated on all cpus
909 : !> using this call is expensive
910 : !> \param fm the full you want to get the info from
911 : !> \param target_m a replicated full matrix that will contain the result
912 : !> \param start_row the starting row of b_matrix (defaults to 1)
913 : !> \param start_col the starting col of b_matrix (defaults to 1)
914 : !> \param n_rows the number of row to change in b (defaults to
915 : !> size(op(new_values),1))
916 : !> \param n_cols the number of columns to change in b (defaults to
917 : !> size(op(new_values),2))
918 : !> \param transpose if target_m should be transposed: if true
919 : !> op(target_m)=target_m^T, else op(target_m)=target_m
920 : !> (defaults to false)
921 : !> \par History
922 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
923 : !> \author Fawzi Mohamed
924 : !> \note
925 : !> optimized for full column updates. Zeros out a little too much
926 : !> of target_m
927 : !> the target_m is replicated and valid on all cpus
928 : ! **************************************************************************************************
929 88184 : SUBROUTINE cp_fm_get_submatrix(fm, target_m, start_row, &
930 : start_col, n_rows, n_cols, transpose)
931 : TYPE(cp_fm_type), INTENT(IN) :: fm
932 : REAL(KIND=dp), DIMENSION(:, :), INTENT(out) :: target_m
933 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
934 : LOGICAL, INTENT(in), OPTIONAL :: transpose
935 :
936 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_get_submatrix'
937 :
938 : INTEGER :: handle, i, i0, j, j0, ncol, ncol_global, &
939 : ncol_local, nrow, nrow_global, &
940 : nrow_local, this_col, this_row
941 88184 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
942 : LOGICAL :: tr_a
943 88184 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
944 : TYPE(mp_para_env_type), POINTER :: para_env
945 :
946 88184 : CALL timeset(routineN, handle)
947 :
948 88184 : i0 = 1; j0 = 1; tr_a = .FALSE.
949 :
950 88184 : IF (PRESENT(start_row)) i0 = start_row
951 88184 : IF (PRESENT(start_col)) j0 = start_col
952 88184 : IF (PRESENT(transpose)) tr_a = transpose
953 6152 : IF (tr_a) THEN
954 2206 : nrow = SIZE(target_m, 2)
955 2206 : ncol = SIZE(target_m, 1)
956 : ELSE
957 85978 : nrow = SIZE(target_m, 1)
958 85978 : ncol = SIZE(target_m, 2)
959 : END IF
960 88184 : IF (PRESENT(n_rows)) nrow = n_rows
961 88184 : IF (PRESENT(n_cols)) ncol = n_cols
962 :
963 88184 : para_env => fm%matrix_struct%para_env
964 :
965 88184 : full_block => fm%local_data
966 : #if defined(__parallel)
967 : ! zero-out whole target_m
968 88184 : IF (SIZE(target_m, 1)*SIZE(target_m, 2) /= 0) THEN
969 207166 : CALL dcopy(SIZE(target_m, 1)*SIZE(target_m, 2), [0.0_dp], 0, target_m, 1)
970 : END IF
971 : #endif
972 :
973 : CALL cp_fm_get_info(matrix=fm, &
974 : nrow_global=nrow_global, ncol_global=ncol_global, &
975 : nrow_local=nrow_local, ncol_local=ncol_local, &
976 88184 : row_indices=row_indices, col_indices=col_indices)
977 :
978 573620 : DO j = 1, ncol_local
979 485436 : this_col = col_indices(j) - j0 + 1
980 573620 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
981 347266 : IF (tr_a) THEN
982 2206 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
983 80516 : DO i = 1, nrow_local
984 80516 : target_m(this_col, row_indices(i)) = full_block(i, j)
985 : END DO
986 : ELSE
987 0 : DO i = 1, nrow_local
988 0 : this_row = row_indices(i) - i0 + 1
989 0 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
990 0 : target_m(this_col, this_row) = full_block(i, j)
991 : END IF
992 : END DO
993 : END IF
994 : ELSE
995 345060 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
996 6449551 : DO i = 1, nrow_local
997 6449551 : target_m(row_indices(i), this_col) = full_block(i, j)
998 : END DO
999 : ELSE
1000 1591252 : DO i = 1, nrow_local
1001 1545930 : this_row = row_indices(i) - i0 + 1
1002 1591252 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
1003 54993 : target_m(this_row, this_col) = full_block(i, j)
1004 : END IF
1005 : END DO
1006 : END IF
1007 : END IF
1008 : END IF
1009 : END DO
1010 :
1011 17726744 : CALL para_env%sum(target_m)
1012 :
1013 88184 : CALL timestop(handle)
1014 :
1015 88184 : END SUBROUTINE cp_fm_get_submatrix
1016 :
1017 : ! **************************************************************************************************
1018 : !> \brief returns all kind of information about the full matrix
1019 : !> \param matrix ...
1020 : !> \param name ...
1021 : !> \param nrow_global ...
1022 : !> \param ncol_global ...
1023 : !> \param nrow_block ...
1024 : !> \param ncol_block ...
1025 : !> \param nrow_local ...
1026 : !> \param ncol_local ...
1027 : !> \param row_indices ...
1028 : !> \param col_indices ...
1029 : !> \param local_data ...
1030 : !> \param context ...
1031 : !> \param nrow_locals ...
1032 : !> \param ncol_locals ...
1033 : !> \param matrix_struct ...
1034 : !> \param para_env ...
1035 : !> \note
1036 : !> see also cp_fm_struct for explanation
1037 : !> - nrow_local, ncol_local, row_indices, col_indices, local_data are hooks for efficient
1038 : !> access to the local blacs block
1039 : ! **************************************************************************************************
1040 6700849 : SUBROUTINE cp_fm_get_info(matrix, name, nrow_global, ncol_global, &
1041 : nrow_block, ncol_block, nrow_local, ncol_local, &
1042 : row_indices, col_indices, local_data, context, &
1043 : nrow_locals, ncol_locals, matrix_struct, para_env)
1044 :
1045 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1046 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: name
1047 : INTEGER, INTENT(OUT), OPTIONAL :: nrow_global, ncol_global, nrow_block, &
1048 : ncol_block, nrow_local, ncol_local
1049 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: row_indices, col_indices
1050 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1051 : OPTIONAL, POINTER :: local_data
1052 : TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: context
1053 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nrow_locals, ncol_locals
1054 : TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: matrix_struct
1055 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
1056 :
1057 8 : IF (PRESENT(name)) name = matrix%name
1058 6700849 : IF (PRESENT(matrix_struct)) matrix_struct => matrix%matrix_struct
1059 6700849 : IF (PRESENT(local_data)) local_data => matrix%local_data ! not hiding things anymore :-(
1060 :
1061 : CALL cp_fm_struct_get(matrix%matrix_struct, nrow_local=nrow_local, &
1062 : ncol_local=ncol_local, nrow_global=nrow_global, &
1063 : ncol_global=ncol_global, nrow_block=nrow_block, &
1064 : ncol_block=ncol_block, row_indices=row_indices, &
1065 : col_indices=col_indices, nrow_locals=nrow_locals, &
1066 6700849 : ncol_locals=ncol_locals, context=context, para_env=para_env)
1067 :
1068 6700849 : END SUBROUTINE cp_fm_get_info
1069 :
1070 : ! **************************************************************************************************
1071 : !> \brief Write nicely formatted info about the FM to the given I/O unit (including the underlying FM struct)
1072 : !> \param matrix a cp_fm_type instance
1073 : !> \param io_unit the I/O unit to use for writing
1074 : ! **************************************************************************************************
1075 3 : SUBROUTINE cp_fm_write_info(matrix, io_unit)
1076 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1077 : INTEGER, INTENT(IN) :: io_unit
1078 :
1079 3 : WRITE (io_unit, '(/,A,A12)') "CP_FM | Name: ", matrix%name
1080 3 : CALL cp_fm_struct_write_info(matrix%matrix_struct, io_unit)
1081 3 : END SUBROUTINE cp_fm_write_info
1082 :
1083 : ! **************************************************************************************************
1084 : !> \brief find the maximum absolute value of the matrix element
1085 : !> maxval(abs(matrix))
1086 : !> \param matrix ...
1087 : !> \param a_max ...
1088 : !> \param ir_max ...
1089 : !> \param ic_max ...
1090 : ! **************************************************************************************************
1091 131657 : SUBROUTINE cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
1092 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1093 : REAL(KIND=dp), INTENT(OUT) :: a_max
1094 : INTEGER, INTENT(OUT), OPTIONAL :: ir_max, ic_max
1095 :
1096 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_maxabsval'
1097 :
1098 : INTEGER :: handle, i, ic_max_local, ir_max_local, &
1099 : j, mepos, ncol_local, nrow_local, &
1100 : num_pe
1101 131657 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ic_max_vec, ir_max_vec
1102 131657 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1103 : REAL(dp) :: my_max
1104 131657 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_max_vec
1105 131657 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1106 :
1107 131657 : CALL timeset(routineN, handle)
1108 :
1109 131657 : my_block => matrix%local_data
1110 :
1111 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
1112 131657 : row_indices=row_indices, col_indices=col_indices)
1113 :
1114 74491237 : a_max = MAXVAL(ABS(my_block(1:nrow_local, 1:ncol_local)))
1115 :
1116 131657 : IF (PRESENT(ir_max)) THEN
1117 0 : num_pe = matrix%matrix_struct%para_env%num_pe
1118 0 : mepos = matrix%matrix_struct%para_env%mepos
1119 0 : ALLOCATE (ir_max_vec(0:num_pe - 1))
1120 0 : ir_max_vec(0:num_pe - 1) = 0
1121 0 : ALLOCATE (ic_max_vec(0:num_pe - 1))
1122 0 : ic_max_vec(0:num_pe - 1) = 0
1123 0 : ALLOCATE (a_max_vec(0:num_pe - 1))
1124 0 : a_max_vec(0:num_pe - 1) = 0.0_dp
1125 0 : my_max = 0.0_dp
1126 :
1127 0 : IF ((ncol_local > 0) .AND. (nrow_local > 0)) THEN
1128 0 : DO i = 1, ncol_local
1129 0 : DO j = 1, nrow_local
1130 0 : IF (ABS(my_block(j, i)) > my_max) THEN
1131 0 : my_max = my_block(j, i)
1132 0 : ir_max_local = j
1133 0 : ic_max_local = i
1134 : END IF
1135 : END DO
1136 : END DO
1137 :
1138 0 : a_max_vec(mepos) = my_max
1139 0 : ir_max_vec(mepos) = row_indices(ir_max_local)
1140 0 : ic_max_vec(mepos) = col_indices(ic_max_local)
1141 :
1142 : END IF
1143 :
1144 0 : CALL matrix%matrix_struct%para_env%sum(a_max_vec)
1145 0 : CALL matrix%matrix_struct%para_env%sum(ir_max_vec)
1146 0 : CALL matrix%matrix_struct%para_env%sum(ic_max_vec)
1147 :
1148 0 : my_max = 0.0_dp
1149 0 : DO i = 0, num_pe - 1
1150 0 : IF (a_max_vec(i) > my_max) THEN
1151 0 : ir_max = ir_max_vec(i)
1152 0 : ic_max = ic_max_vec(i)
1153 : END IF
1154 : END DO
1155 :
1156 0 : DEALLOCATE (ir_max_vec, ic_max_vec, a_max_vec)
1157 0 : CPASSERT(ic_max > 0)
1158 0 : CPASSERT(ir_max > 0)
1159 :
1160 : END IF
1161 :
1162 131657 : CALL matrix%matrix_struct%para_env%max(a_max)
1163 :
1164 131657 : CALL timestop(handle)
1165 :
1166 263314 : END SUBROUTINE cp_fm_maxabsval
1167 :
1168 : ! **************************************************************************************************
1169 : !> \brief find the maximum over the rows of the sum of the absolute values of the elements of a given row
1170 : !> = || A ||_infinity
1171 : !> \param matrix ...
1172 : !> \param a_max ...
1173 : !> \note
1174 : !> for a real symmetric matrix it holds that || A ||_2 = |lambda_max| < || A ||_infinity
1175 : !> Hence this can be used to estimate an upper bound for the eigenvalues of a matrix
1176 : !> http://mathworld.wolfram.com/MatrixNorm.html
1177 : !> (but the bound is not so tight in the general case)
1178 : ! **************************************************************************************************
1179 4984 : SUBROUTINE cp_fm_maxabsrownorm(matrix, a_max)
1180 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1181 : REAL(KIND=dp), INTENT(OUT) :: a_max
1182 :
1183 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_maxabsrownorm'
1184 :
1185 : INTEGER :: handle, i, j, ncol_local, nrow_global, &
1186 : nrow_local
1187 4984 : INTEGER, DIMENSION(:), POINTER :: row_indices
1188 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: values
1189 4984 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1190 :
1191 4984 : CALL timeset(routineN, handle)
1192 :
1193 4984 : my_block => matrix%local_data
1194 :
1195 : CALL cp_fm_get_info(matrix, row_indices=row_indices, nrow_global=nrow_global, &
1196 4984 : nrow_local=nrow_local, ncol_local=ncol_local)
1197 :
1198 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1199 14952 : ALLOCATE (values(nrow_global))
1200 4984 : values = 0.0_dp
1201 65804 : DO j = 1, ncol_local
1202 535024 : DO i = 1, nrow_local
1203 530040 : values(row_indices(i)) = values(row_indices(i)) + ABS(my_block(i, j))
1204 : END DO
1205 : END DO
1206 4984 : CALL matrix%matrix_struct%para_env%sum(values)
1207 65804 : a_max = MAXVAL(values)
1208 4984 : DEALLOCATE (values)
1209 :
1210 4984 : CALL timestop(handle)
1211 4984 : END SUBROUTINE cp_fm_maxabsrownorm
1212 :
1213 : ! **************************************************************************************************
1214 : !> \brief find the inorm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
1215 : !> \param matrix ...
1216 : !> \param norm_array ...
1217 : ! **************************************************************************************************
1218 1504 : SUBROUTINE cp_fm_vectorsnorm(matrix, norm_array)
1219 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1220 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: norm_array
1221 :
1222 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_vectorsnorm'
1223 :
1224 : INTEGER :: handle, i, j, ncol_global, ncol_local, &
1225 : nrow_local
1226 1504 : INTEGER, DIMENSION(:), POINTER :: col_indices
1227 1504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1228 :
1229 1504 : CALL timeset(routineN, handle)
1230 :
1231 1504 : my_block => matrix%local_data
1232 :
1233 : CALL cp_fm_get_info(matrix, col_indices=col_indices, ncol_global=ncol_global, &
1234 1504 : nrow_local=nrow_local, ncol_local=ncol_local)
1235 :
1236 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1237 34252 : norm_array = 0.0_dp
1238 34252 : DO j = 1, ncol_local
1239 1401764 : DO i = 1, nrow_local
1240 1400260 : norm_array(col_indices(j)) = norm_array(col_indices(j)) + my_block(i, j)*my_block(i, j)
1241 : END DO
1242 : END DO
1243 67000 : CALL matrix%matrix_struct%para_env%sum(norm_array)
1244 34252 : norm_array = SQRT(norm_array)
1245 :
1246 1504 : CALL timestop(handle)
1247 1504 : END SUBROUTINE cp_fm_vectorsnorm
1248 :
1249 : ! **************************************************************************************************
1250 : !> \brief summing up all the elements along the matrix's i-th index
1251 : !> \f$ \mathrm{sum}_{j} = \sum_{i} A_{ij} \f$
1252 : !> or
1253 : !> \f$ \mathrm{sum}_{i} = \sum_{j} A_{ij} \f$
1254 : !> \param matrix an input matrix A
1255 : !> \param sum_array sums of elements in each column/row
1256 : !> \param dir ...
1257 : !> \note forked from cp_fm_vectorsnorm() to be used with
1258 : !> the maximum overlap method
1259 : !> added row variation
1260 : ! **************************************************************************************************
1261 11172 : SUBROUTINE cp_fm_vectorssum(matrix, sum_array, dir)
1262 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1263 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: sum_array
1264 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: dir
1265 :
1266 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_vectorssum'
1267 :
1268 : INTEGER :: handle, i, j, ncol_local, nrow_local
1269 11172 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1270 : LOGICAL :: docol
1271 11172 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1272 :
1273 11172 : CALL timeset(routineN, handle)
1274 :
1275 11172 : IF (PRESENT(dir)) THEN
1276 11132 : IF (dir == 'c' .OR. dir == 'C') THEN
1277 : docol = .TRUE.
1278 : ELSE IF (dir == 'r' .OR. dir == 'R') THEN
1279 : docol = .FALSE.
1280 : ELSE
1281 0 : CPABORT('Wrong argument DIR')
1282 : END IF
1283 : ELSE
1284 : docol = .TRUE.
1285 : END IF
1286 :
1287 11172 : my_block => matrix%local_data
1288 :
1289 : CALL cp_fm_get_info(matrix, col_indices=col_indices, row_indices=row_indices, &
1290 11172 : nrow_local=nrow_local, ncol_local=ncol_local)
1291 :
1292 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1293 279588 : sum_array(:) = 0.0_dp
1294 11172 : IF (docol) THEN
1295 448 : DO j = 1, ncol_local
1296 3628 : DO i = 1, nrow_local
1297 3588 : sum_array(col_indices(j)) = sum_array(col_indices(j)) + my_block(i, j)
1298 : END DO
1299 : END DO
1300 : ELSE
1301 100132 : DO j = 1, ncol_local
1302 6513944 : DO i = 1, nrow_local
1303 6502812 : sum_array(row_indices(i)) = sum_array(row_indices(i)) + my_block(i, j)
1304 : END DO
1305 : END DO
1306 : END IF
1307 548004 : CALL matrix%matrix_struct%para_env%sum(sum_array)
1308 :
1309 11172 : CALL timestop(handle)
1310 11172 : END SUBROUTINE cp_fm_vectorssum
1311 :
1312 : ! **************************************************************************************************
1313 : !> \brief copy one identically sized matrix in the other
1314 : !> \param source ...
1315 : !> \param destination ...
1316 : !> \note
1317 : !> see also cp_fm_to_fm_columns
1318 : ! **************************************************************************************************
1319 435071 : SUBROUTINE cp_fm_to_fm_matrix(source, destination)
1320 :
1321 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1322 :
1323 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_matrix'
1324 :
1325 : INTEGER :: handle, npcol, nprow
1326 :
1327 435071 : CALL timeset(routineN, handle)
1328 :
1329 435071 : nprow = source%matrix_struct%context%num_pe(1)
1330 435071 : npcol = source%matrix_struct%context%num_pe(2)
1331 :
1332 435071 : IF ((.NOT. cp2k_is_parallel) .OR. &
1333 : cp_fm_struct_equivalent(source%matrix_struct, &
1334 : destination%matrix_struct)) THEN
1335 435071 : IF (SIZE(source%local_data, 1) /= SIZE(destination%local_data, 1) .OR. &
1336 : SIZE(source%local_data, 2) /= SIZE(destination%local_data, 2)) THEN
1337 : CALL cp_abort(__LOCATION__, &
1338 : "Cannot copy full matrix <"//TRIM(source%name)// &
1339 : "> to full matrix <"//TRIM(destination%name)// &
1340 0 : ">. The local_data blocks have different sizes.")
1341 : END IF
1342 : CALL dcopy(SIZE(source%local_data, 1)*SIZE(source%local_data, 2), &
1343 435071 : source%local_data, 1, destination%local_data, 1)
1344 : ELSE
1345 0 : CPABORT("Data structures of source and target full matrix are not equivalent")
1346 : END IF
1347 :
1348 435071 : CALL timestop(handle)
1349 :
1350 435071 : END SUBROUTINE cp_fm_to_fm_matrix
1351 :
1352 : ! **************************************************************************************************
1353 : !> \brief copy just a subset of columns of a fm to a fm
1354 : !> \param msource ...
1355 : !> \param mtarget ...
1356 : !> \param ncol ...
1357 : !> \param source_start ...
1358 : !> \param target_start ...
1359 : ! **************************************************************************************************
1360 209488 : SUBROUTINE cp_fm_to_fm_columns(msource, mtarget, ncol, source_start, &
1361 : target_start)
1362 :
1363 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1364 : INTEGER, INTENT(IN) :: ncol
1365 : INTEGER, INTENT(IN), OPTIONAL :: source_start, target_start
1366 :
1367 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_columns'
1368 :
1369 : INTEGER :: handle, n, ss, ts
1370 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1371 : #if defined(__parallel)
1372 : INTEGER :: i
1373 : INTEGER, DIMENSION(9) :: desca, descb
1374 : #endif
1375 :
1376 209488 : CALL timeset(routineN, handle)
1377 :
1378 209488 : ss = 1
1379 209488 : ts = 1
1380 :
1381 209488 : IF (PRESENT(source_start)) ss = source_start
1382 209488 : IF (PRESENT(target_start)) ts = target_start
1383 :
1384 209488 : n = msource%matrix_struct%nrow_global
1385 :
1386 209488 : a => msource%local_data
1387 209488 : b => mtarget%local_data
1388 :
1389 : #if defined(__parallel)
1390 2094880 : desca(:) = msource%matrix_struct%descriptor(:)
1391 2094880 : descb(:) = mtarget%matrix_struct%descriptor(:)
1392 999112 : DO i = 0, ncol - 1
1393 999112 : CALL pdcopy(n, a, 1, ss + i, desca, 1, b, 1, ts + i, descb, 1)
1394 : END DO
1395 : #else
1396 : IF (ss <= SIZE(a, 2) .AND. ts <= SIZE(b, 2)) THEN
1397 : CALL dcopy(ncol*n, a(:, ss), 1, b(:, ts), 1)
1398 : END IF
1399 : #endif
1400 :
1401 209488 : CALL timestop(handle)
1402 :
1403 209488 : END SUBROUTINE cp_fm_to_fm_columns
1404 :
1405 : ! **************************************************************************************************
1406 : !> \brief copy just a triangular matrix
1407 : !> \param msource ...
1408 : !> \param mtarget ...
1409 : !> \param uplo ...
1410 : ! **************************************************************************************************
1411 58 : SUBROUTINE cp_fm_to_fm_triangular(msource, mtarget, uplo)
1412 :
1413 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1414 : CHARACTER(LEN=1), OPTIONAL, INTENT(IN) :: uplo
1415 :
1416 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_triangular'
1417 :
1418 : CHARACTER(LEN=1) :: myuplo
1419 : INTEGER :: handle, ncol, nrow
1420 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1421 : #if defined(__parallel)
1422 : INTEGER, DIMENSION(9) :: desca, descb
1423 : #endif
1424 :
1425 58 : CALL timeset(routineN, handle)
1426 :
1427 58 : myuplo = 'U'
1428 58 : IF (PRESENT(uplo)) myuplo = uplo
1429 :
1430 58 : nrow = msource%matrix_struct%nrow_global
1431 58 : ncol = msource%matrix_struct%ncol_global
1432 :
1433 58 : a => msource%local_data
1434 58 : b => mtarget%local_data
1435 :
1436 : #if defined(__parallel)
1437 580 : desca(:) = msource%matrix_struct%descriptor(:)
1438 580 : descb(:) = mtarget%matrix_struct%descriptor(:)
1439 58 : CALL pdlacpy(myuplo, nrow, ncol, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb)
1440 : #else
1441 : CALL dlacpy(myuplo, nrow, ncol, a(1, 1), nrow, b(1, 1), nrow)
1442 : #endif
1443 :
1444 58 : CALL timestop(handle)
1445 :
1446 58 : END SUBROUTINE cp_fm_to_fm_triangular
1447 :
1448 : ! **************************************************************************************************
1449 : !> \brief copy just a part ot the matrix
1450 : !> \param msource ...
1451 : !> \param mtarget ...
1452 : !> \param nrow ...
1453 : !> \param ncol ...
1454 : !> \param s_firstrow ...
1455 : !> \param s_firstcol ...
1456 : !> \param t_firstrow ...
1457 : !> \param t_firstcol ...
1458 : ! **************************************************************************************************
1459 :
1460 16846 : SUBROUTINE cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
1461 :
1462 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1463 : INTEGER, INTENT(IN) :: nrow, ncol, s_firstrow, &
1464 : s_firstcol, t_firstrow, &
1465 : t_firstcol
1466 :
1467 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat'
1468 :
1469 : INTEGER :: handle, i, na, nb, ss, ts
1470 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1471 : #if defined(__parallel)
1472 : INTEGER, DIMENSION(9) :: desca, descb
1473 : #endif
1474 :
1475 16846 : CALL timeset(routineN, handle)
1476 :
1477 16846 : a => msource%local_data
1478 16846 : b => mtarget%local_data
1479 :
1480 16846 : na = msource%matrix_struct%nrow_global
1481 16846 : nb = mtarget%matrix_struct%nrow_global
1482 : ! nrow must be <= na and nb
1483 16846 : IF (nrow > na) THEN
1484 0 : CPABORT("cannot copy because nrow > number of rows of source matrix")
1485 : END IF
1486 16846 : IF (nrow > nb) THEN
1487 0 : CPABORT("cannot copy because nrow > number of rows of target matrix")
1488 : END IF
1489 16846 : na = msource%matrix_struct%ncol_global
1490 16846 : nb = mtarget%matrix_struct%ncol_global
1491 : ! ncol must be <= na_col and nb_col
1492 16846 : IF (ncol > na) THEN
1493 0 : CPABORT("cannot copy because nrow > number of rows of source matrix")
1494 : END IF
1495 16846 : IF (ncol > nb) THEN
1496 0 : CPABORT("cannot copy because nrow > number of rows of target matrix")
1497 : END IF
1498 :
1499 : #if defined(__parallel)
1500 168460 : desca(:) = msource%matrix_struct%descriptor(:)
1501 168460 : descb(:) = mtarget%matrix_struct%descriptor(:)
1502 211064 : DO i = 0, ncol - 1
1503 194218 : ss = s_firstcol + i
1504 194218 : ts = t_firstcol + i
1505 211064 : CALL pdcopy(nrow, a, s_firstrow, ss, desca, 1, b, t_firstrow, ts, descb, 1)
1506 : END DO
1507 : #else
1508 : DO i = 0, ncol - 1
1509 : ss = s_firstcol + i
1510 : ts = t_firstcol + i
1511 : CALL dcopy(nrow, a(s_firstrow:, ss), 1, b(t_firstrow:, ts), 1)
1512 : END DO
1513 : #endif
1514 :
1515 16846 : CALL timestop(handle)
1516 16846 : END SUBROUTINE cp_fm_to_fm_submat
1517 :
1518 : ! **************************************************************************************************
1519 : !> \brief General copy of a fm matrix to another fm matrix.
1520 : !> Uses non-blocking MPI rather than ScaLAPACK.
1521 : !>
1522 : !> \param source input fm matrix
1523 : !> \param destination output fm matrix
1524 : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
1525 : !> of the input and output matrices
1526 : !> \par History
1527 : !> 31-Jan-2017 : Re-implemented using non-blocking MPI [IainB, MarkT]
1528 : ! **************************************************************************************************
1529 5184 : SUBROUTINE cp_fm_copy_general(source, destination, para_env)
1530 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1531 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1532 :
1533 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_copy_general'
1534 :
1535 : INTEGER :: handle
1536 46656 : TYPE(copy_info_type) :: info
1537 :
1538 5184 : CALL timeset(routineN, handle)
1539 :
1540 5184 : CALL cp_fm_start_copy_general(source, destination, para_env, info)
1541 5184 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1542 5170 : CALL cp_fm_finish_copy_general(destination, info)
1543 : END IF
1544 5184 : IF (ASSOCIATED(source%matrix_struct)) THEN
1545 5023 : CALL cp_fm_cleanup_copy_general(info)
1546 : END IF
1547 :
1548 5184 : CALL timestop(handle)
1549 5184 : END SUBROUTINE cp_fm_copy_general
1550 :
1551 : ! **************************************************************************************************
1552 : !> \brief Initiates the copy operation: get distribution data, post MPI isend and irecvs
1553 : !> \param source input fm matrix
1554 : !> \param destination output fm matrix
1555 : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
1556 : !> of the input and output matrices
1557 : !> \param info all of the data that will be needed to complete the copy operation
1558 : ! **************************************************************************************************
1559 198320 : SUBROUTINE cp_fm_start_copy_general(source, destination, para_env, info)
1560 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1561 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1562 : TYPE(copy_info_type), INTENT(OUT) :: info
1563 :
1564 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_start_copy_general'
1565 :
1566 : INTEGER :: dest_p_i, dest_q_j, global_rank, global_size, handle, i, j, k, mpi_rank, &
1567 : ncol_block_dest, ncol_block_src, ncol_local_recv, ncol_local_send, ncols, &
1568 : nrow_block_dest, nrow_block_src, nrow_local_recv, nrow_local_send, nrows, p, q, &
1569 : recv_rank, recv_size, send_rank, send_size
1570 19832 : INTEGER, ALLOCATABLE, DIMENSION(:) :: all_ranks, dest2global, dest_p, dest_q, &
1571 39664 : recv_count, send_count, send_disp, &
1572 19832 : source2global, src_p, src_q
1573 19832 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dest_blacs2mpi
1574 : INTEGER, DIMENSION(2) :: dest_block, dest_block_tmp, dest_num_pe, &
1575 : src_block, src_block_tmp, src_num_pe
1576 39664 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices, &
1577 39664 : send_col_indices, send_row_indices
1578 : TYPE(cp_fm_struct_type), POINTER :: recv_dist, send_dist
1579 277648 : TYPE(mp_request_type), DIMENSION(6) :: recv_req, send_req
1580 :
1581 19832 : CALL timeset(routineN, handle)
1582 :
1583 : IF (.NOT. cp2k_is_parallel) THEN
1584 : ! Just copy all of the matrix data into a 'send buffer', to be unpacked later
1585 : nrow_local_send = SIZE(source%local_data, 1)
1586 : ncol_local_send = SIZE(source%local_data, 2)
1587 : ALLOCATE (info%send_buf(nrow_local_send*ncol_local_send))
1588 : k = 0
1589 : DO j = 1, ncol_local_send
1590 : DO i = 1, nrow_local_send
1591 : k = k + 1
1592 : info%send_buf(k) = source%local_data(i, j)
1593 : END DO
1594 : END DO
1595 : ELSE
1596 19832 : NULLIFY (recv_dist, send_dist)
1597 19832 : NULLIFY (recv_col_indices, recv_row_indices, send_col_indices, send_row_indices)
1598 :
1599 : ! The 'global' communicator contains both the source and destination decompositions
1600 19832 : global_size = para_env%num_pe
1601 19832 : global_rank = para_env%mepos
1602 :
1603 : ! The source/send decomposition and destination/recv decompositions may only exist on
1604 : ! on a subset of the processes involved in the communication
1605 : ! Check if the source and/or destination arguments are .not. ASSOCIATED():
1606 : ! if so, skip the send / recv parts (since these processes do not participate in the sending/receiving distribution)
1607 19832 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1608 15160 : recv_dist => destination%matrix_struct
1609 15160 : recv_rank = recv_dist%para_env%mepos
1610 : ELSE
1611 4672 : recv_rank = mp_proc_null
1612 : END IF
1613 :
1614 19832 : IF (ASSOCIATED(source%matrix_struct)) THEN
1615 17079 : send_dist => source%matrix_struct
1616 17079 : send_rank = send_dist%para_env%mepos
1617 : ELSE
1618 2753 : send_rank = mp_proc_null
1619 : END IF
1620 :
1621 : ! Map the rank in the source/dest communicator to the global rank
1622 59496 : ALLOCATE (all_ranks(0:global_size - 1))
1623 :
1624 19832 : CALL para_env%allgather(send_rank, all_ranks)
1625 19832 : IF (ASSOCIATED(recv_dist)) THEN
1626 75800 : ALLOCATE (source2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
1627 45480 : DO i = 0, global_size - 1
1628 45480 : IF (all_ranks(i) /= mp_proc_null) THEN
1629 24814 : source2global(all_ranks(i)) = i
1630 : END IF
1631 : END DO
1632 : END IF
1633 :
1634 19832 : CALL para_env%allgather(recv_rank, all_ranks)
1635 19832 : IF (ASSOCIATED(send_dist)) THEN
1636 85395 : ALLOCATE (dest2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
1637 51237 : DO i = 0, global_size - 1
1638 51237 : IF (all_ranks(i) /= mp_proc_null) THEN
1639 24814 : dest2global(all_ranks(i)) = i
1640 : END IF
1641 : END DO
1642 : END IF
1643 19832 : DEALLOCATE (all_ranks)
1644 :
1645 : ! Some data from the two decompositions will be needed by all processes in the global group :
1646 : ! process grid shape, block size, and the BLACS-to-MPI mapping
1647 :
1648 : ! The global root process will receive the data (from the root process in each decomposition)
1649 138824 : send_req(:) = mp_request_null
1650 19832 : IF (global_rank == 0) THEN
1651 69412 : recv_req(:) = mp_request_null
1652 9916 : CALL para_env%irecv(src_block, mp_any_source, recv_req(1), tag=src_tag)
1653 9916 : CALL para_env%irecv(dest_block, mp_any_source, recv_req(2), tag=dest_tag)
1654 9916 : CALL para_env%irecv(src_num_pe, mp_any_source, recv_req(3), tag=src_tag)
1655 9916 : CALL para_env%irecv(dest_num_pe, mp_any_source, recv_req(4), tag=dest_tag)
1656 : END IF
1657 :
1658 19832 : IF (ASSOCIATED(send_dist)) THEN
1659 17079 : IF ((send_rank == 0)) THEN
1660 : ! need to use separate buffers here in case this is actually global rank 0
1661 29748 : src_block_tmp = [send_dist%nrow_block, send_dist%ncol_block]
1662 9916 : CALL para_env%isend(src_block_tmp, 0, send_req(1), tag=src_tag)
1663 9916 : CALL para_env%isend(send_dist%context%num_pe, 0, send_req(2), tag=src_tag)
1664 : END IF
1665 : END IF
1666 :
1667 19832 : IF (ASSOCIATED(recv_dist)) THEN
1668 15160 : IF ((recv_rank == 0)) THEN
1669 29748 : dest_block_tmp = [recv_dist%nrow_block, recv_dist%ncol_block]
1670 9916 : CALL para_env%isend(dest_block_tmp, 0, send_req(3), tag=dest_tag)
1671 9916 : CALL para_env%isend(recv_dist%context%num_pe, 0, send_req(4), tag=dest_tag)
1672 : END IF
1673 : END IF
1674 :
1675 19832 : IF (global_rank == 0) THEN
1676 9916 : CALL mp_waitall(recv_req(1:4))
1677 : ! Now we know the process decomposition, we can allocate the arrays to hold the blacs2mpi mapping
1678 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1679 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
1680 69412 : )
1681 9916 : CALL para_env%irecv(info%src_blacs2mpi, mp_any_source, recv_req(5), tag=src_tag)
1682 9916 : CALL para_env%irecv(dest_blacs2mpi, mp_any_source, recv_req(6), tag=dest_tag)
1683 : END IF
1684 :
1685 19832 : IF (ASSOCIATED(send_dist)) THEN
1686 17079 : IF ((send_rank == 0)) THEN
1687 9916 : CALL para_env%isend(send_dist%context%blacs2mpi(:, :), 0, send_req(5), tag=src_tag)
1688 : END IF
1689 : END IF
1690 :
1691 19832 : IF (ASSOCIATED(recv_dist)) THEN
1692 15160 : IF ((recv_rank == 0)) THEN
1693 9916 : CALL para_env%isend(recv_dist%context%blacs2mpi(:, :), 0, send_req(6), tag=dest_tag)
1694 : END IF
1695 : END IF
1696 :
1697 19832 : IF (global_rank == 0) THEN
1698 9916 : CALL mp_waitall(recv_req(5:6))
1699 : END IF
1700 :
1701 : ! Finally, broadcast the data to all processes in the global communicator
1702 19832 : CALL para_env%bcast(src_block, 0)
1703 19832 : CALL para_env%bcast(dest_block, 0)
1704 19832 : CALL para_env%bcast(src_num_pe, 0)
1705 19832 : CALL para_env%bcast(dest_num_pe, 0)
1706 59496 : info%src_num_pe(1:2) = src_num_pe(1:2)
1707 59496 : info%nblock_src(1:2) = src_block(1:2)
1708 19832 : IF (global_rank /= 0) THEN
1709 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1710 0 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
1711 69412 : )
1712 : END IF
1713 19832 : CALL para_env%bcast(info%src_blacs2mpi, 0)
1714 19832 : CALL para_env%bcast(dest_blacs2mpi, 0)
1715 :
1716 19832 : recv_size = dest_num_pe(1)*dest_num_pe(2)
1717 19832 : send_size = src_num_pe(1)*src_num_pe(2)
1718 19832 : info%send_size = send_size
1719 19832 : CALL mp_waitall(send_req(:))
1720 :
1721 : ! Setup is now complete, we can start the actual communication here.
1722 : ! The order implemented here is:
1723 : ! DEST_1
1724 : ! compute recv sizes
1725 : ! call irecv
1726 : ! SRC_1
1727 : ! compute send sizes
1728 : ! pack send buffers
1729 : ! call isend
1730 : ! DEST_2
1731 : ! wait for the recvs and unpack buffers (this part eventually will go into another
1732 : ! routine to allow comms to run concurrently)
1733 : ! SRC_2
1734 : ! wait for the sends
1735 :
1736 : ! DEST_1
1737 19832 : IF (ASSOCIATED(recv_dist)) THEN
1738 : CALL cp_fm_struct_get(recv_dist, row_indices=recv_row_indices, &
1739 : col_indices=recv_col_indices &
1740 15160 : )
1741 15160 : info%recv_col_indices => recv_col_indices
1742 15160 : info%recv_row_indices => recv_row_indices
1743 15160 : nrow_block_src = src_block(1)
1744 15160 : ncol_block_src = src_block(2)
1745 100614 : ALLOCATE (recv_count(0:send_size - 1), info%recv_disp(0:send_size - 1), info%recv_request(0:send_size - 1))
1746 :
1747 : ! Determine the recv counts, allocate the receive buffers, call mpi_irecv for all the non-zero sized receives
1748 15160 : nrow_local_recv = recv_dist%nrow_locals(recv_dist%context%mepos(1))
1749 15160 : ncol_local_recv = recv_dist%ncol_locals(recv_dist%context%mepos(2))
1750 15160 : info%nlocal_recv(1) = nrow_local_recv
1751 15160 : info%nlocal_recv(2) = ncol_local_recv
1752 : ! Initialise src_p, src_q arrays (sized using number of rows/cols in the receiving distribution)
1753 75800 : ALLOCATE (src_p(nrow_local_recv), src_q(ncol_local_recv))
1754 217833 : DO i = 1, nrow_local_recv
1755 : ! For each local row we will receive, we look up its global row (in recv_row_indices),
1756 : ! then work out which row block it comes from, and which process row that row block comes from.
1757 217833 : src_p(i) = MOD(((recv_row_indices(i) - 1)/nrow_block_src), src_num_pe(1))
1758 : END DO
1759 342584 : DO j = 1, ncol_local_recv
1760 : ! Similarly for the columns
1761 342584 : src_q(j) = MOD(((recv_col_indices(j) - 1)/ncol_block_src), src_num_pe(2))
1762 : END DO
1763 : ! src_p/q now contains the process row/column ID that will send data to that row/column
1764 :
1765 30320 : DO q = 0, src_num_pe(2) - 1
1766 342584 : ncols = COUNT(src_q == q)
1767 55134 : DO p = 0, src_num_pe(1) - 1
1768 418770 : nrows = COUNT(src_p == p)
1769 : ! Use the send_dist here as we are looking up the processes where the data comes from
1770 39974 : recv_count(info%src_blacs2mpi(p, q)) = nrows*ncols
1771 : END DO
1772 : END DO
1773 15160 : DEALLOCATE (src_p, src_q)
1774 :
1775 : ! Use one long buffer (and displacements into that buffer)
1776 : ! this prevents the need for a rectangular array where not all elements will be populated
1777 70294 : ALLOCATE (info%recv_buf(SUM(recv_count(:))))
1778 15160 : info%recv_disp(0) = 0
1779 24814 : DO i = 1, send_size - 1
1780 24814 : info%recv_disp(i) = info%recv_disp(i - 1) + recv_count(i - 1)
1781 : END DO
1782 :
1783 : ! Issue receive calls on ranks which expect data
1784 39974 : DO k = 0, send_size - 1
1785 39974 : IF (recv_count(k) > 0) THEN
1786 : CALL para_env%irecv(info%recv_buf(info%recv_disp(k) + 1:info%recv_disp(k) + recv_count(k)), &
1787 19834 : source2global(k), info%recv_request(k))
1788 : END IF
1789 : END DO
1790 15160 : DEALLOCATE (source2global)
1791 : END IF ! ASSOCIATED(recv_dist)
1792 :
1793 : ! SRC_1
1794 19832 : IF (ASSOCIATED(send_dist)) THEN
1795 : CALL cp_fm_struct_get(send_dist, row_indices=send_row_indices, &
1796 : col_indices=send_col_indices &
1797 17079 : )
1798 17079 : nrow_block_dest = dest_block(1)
1799 17079 : ncol_block_dest = dest_block(2)
1800 110209 : ALLOCATE (send_count(0:recv_size - 1), send_disp(0:recv_size - 1), info%send_request(0:recv_size - 1))
1801 :
1802 : ! Determine the send counts, allocate the send buffers
1803 17079 : nrow_local_send = send_dist%nrow_locals(send_dist%context%mepos(1))
1804 17079 : ncol_local_send = send_dist%ncol_locals(send_dist%context%mepos(2))
1805 :
1806 : ! Initialise dest_p, dest_q arrays (sized nrow_local, ncol_local)
1807 : ! i.e. number of rows,cols in the sending distribution
1808 85395 : ALLOCATE (dest_p(nrow_local_send), dest_q(ncol_local_send))
1809 :
1810 219752 : DO i = 1, nrow_local_send
1811 : ! Use the send_dist%row_indices() here (we are looping over the local rows we will send)
1812 219752 : dest_p(i) = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1813 : END DO
1814 399231 : DO j = 1, ncol_local_send
1815 399231 : dest_q(j) = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1816 : END DO
1817 : ! dest_p/q now contain the process row/column ID that will receive data from that row/column
1818 :
1819 34158 : DO q = 0, dest_num_pe(2) - 1
1820 399231 : ncols = COUNT(dest_q == q)
1821 58972 : DO p = 0, dest_num_pe(1) - 1
1822 367080 : nrows = COUNT(dest_p == p)
1823 41893 : send_count(dest_blacs2mpi(p, q)) = nrows*ncols
1824 : END DO
1825 : END DO
1826 17079 : DEALLOCATE (dest_p, dest_q)
1827 :
1828 : ! Allocate the send buffer using send_count -- and calculate the offset into the buffer for each process
1829 76051 : ALLOCATE (info%send_buf(SUM(send_count(:))))
1830 17079 : send_disp(0) = 0
1831 24814 : DO k = 1, recv_size - 1
1832 24814 : send_disp(k) = send_disp(k - 1) + send_count(k - 1)
1833 : END DO
1834 :
1835 : ! Loop over the smat, pack the send buffers
1836 17079 : send_count(:) = 0
1837 399231 : DO j = 1, ncol_local_send
1838 : ! Use send_col_indices and row_indices here, as we are looking up the global row/column number of local rows.
1839 382152 : dest_q_j = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1840 11553271 : DO i = 1, nrow_local_send
1841 11154040 : dest_p_i = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1842 11154040 : mpi_rank = dest_blacs2mpi(dest_p_i, dest_q_j)
1843 11154040 : send_count(mpi_rank) = send_count(mpi_rank) + 1
1844 11536192 : info%send_buf(send_disp(mpi_rank) + send_count(mpi_rank)) = source%local_data(i, j)
1845 : END DO
1846 : END DO
1847 :
1848 : ! For each non-zero send_count, call mpi_isend
1849 41893 : DO k = 0, recv_size - 1
1850 41893 : IF (send_count(k) > 0) THEN
1851 : CALL para_env%isend(info%send_buf(send_disp(k) + 1:send_disp(k) + send_count(k)), &
1852 19834 : dest2global(k), info%send_request(k))
1853 : END IF
1854 : END DO
1855 17079 : DEALLOCATE (send_count, send_disp, dest2global)
1856 : END IF ! ASSOCIATED(send_dist)
1857 19832 : DEALLOCATE (dest_blacs2mpi)
1858 :
1859 : END IF !IF (.NOT. cp2k_is_parallel)
1860 :
1861 19832 : CALL timestop(handle)
1862 :
1863 79328 : END SUBROUTINE cp_fm_start_copy_general
1864 :
1865 : ! **************************************************************************************************
1866 : !> \brief Completes the copy operation: wait for comms, unpack, clean up MPI state
1867 : !> \param destination output fm matrix
1868 : !> \param info all of the data that will be needed to complete the copy operation
1869 : ! **************************************************************************************************
1870 15160 : SUBROUTINE cp_fm_finish_copy_general(destination, info)
1871 : TYPE(cp_fm_type), INTENT(IN) :: destination
1872 : TYPE(copy_info_type), INTENT(INOUT) :: info
1873 :
1874 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_finish_copy_general'
1875 :
1876 : INTEGER :: handle, i, j, k, mpi_rank, send_size, &
1877 : src_p_i, src_q_j
1878 15160 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_count
1879 : INTEGER, DIMENSION(2) :: nblock_src, nlocal_recv, src_num_pe
1880 15160 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices
1881 :
1882 15160 : CALL timeset(routineN, handle)
1883 :
1884 : IF (.NOT. cp2k_is_parallel) THEN
1885 : ! Now unpack the data from the 'send buffer'
1886 : k = 0
1887 : DO j = 1, SIZE(destination%local_data, 2)
1888 : DO i = 1, SIZE(destination%local_data, 1)
1889 : k = k + 1
1890 : destination%local_data(i, j) = info%send_buf(k)
1891 : END DO
1892 : END DO
1893 : DEALLOCATE (info%send_buf)
1894 : ELSE
1895 : ! Set up local variables ...
1896 15160 : send_size = info%send_size
1897 45480 : nlocal_recv(1:2) = info%nlocal_recv(:)
1898 45480 : nblock_src(1:2) = info%nblock_src(:)
1899 45480 : src_num_pe(1:2) = info%src_num_pe(:)
1900 15160 : recv_col_indices => info%recv_col_indices
1901 15160 : recv_row_indices => info%recv_row_indices
1902 :
1903 : ! ... use the local variables to do the work
1904 : ! DEST_2
1905 15160 : CALL mp_waitall(info%recv_request(:))
1906 45480 : ALLOCATE (recv_count(0:send_size - 1))
1907 : ! Loop over the rmat, filling it in with data from the recv buffers
1908 : ! (here the block sizes, num_pes refer to the distribution of the source matrix)
1909 15160 : recv_count(:) = 0
1910 342584 : DO j = 1, nlocal_recv(2)
1911 327424 : src_q_j = MOD(((recv_col_indices(j) - 1)/nblock_src(2)), src_num_pe(2))
1912 11496624 : DO i = 1, nlocal_recv(1)
1913 11154040 : src_p_i = MOD(((recv_row_indices(i) - 1)/nblock_src(1)), src_num_pe(1))
1914 11154040 : mpi_rank = info%src_blacs2mpi(src_p_i, src_q_j)
1915 11154040 : recv_count(mpi_rank) = recv_count(mpi_rank) + 1
1916 11481464 : destination%local_data(i, j) = info%recv_buf(info%recv_disp(mpi_rank) + recv_count(mpi_rank))
1917 : END DO
1918 : END DO
1919 15160 : DEALLOCATE (recv_count, info%recv_disp, info%recv_request, info%recv_buf, info%src_blacs2mpi)
1920 : ! Invalidate the stored state
1921 : NULLIFY (info%recv_col_indices, &
1922 15160 : info%recv_row_indices)
1923 :
1924 : END IF
1925 :
1926 15160 : CALL timestop(handle)
1927 :
1928 15160 : END SUBROUTINE cp_fm_finish_copy_general
1929 :
1930 : ! **************************************************************************************************
1931 : !> \brief Completes the copy operation: wait for comms clean up MPI state
1932 : !> \param info all of the data that will be needed to complete the copy operation
1933 : ! **************************************************************************************************
1934 17079 : SUBROUTINE cp_fm_cleanup_copy_general(info)
1935 : TYPE(copy_info_type), INTENT(INOUT) :: info
1936 :
1937 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_cleanup_copy_general'
1938 :
1939 : INTEGER :: handle
1940 :
1941 17079 : CALL timeset(routineN, handle)
1942 :
1943 : IF (.NOT. cp2k_is_parallel) THEN
1944 : ! Don't do anything - no MPI state for the serial case
1945 : ELSE
1946 : ! SRC_2
1947 : ! If this process is also in the destination decomposition, this deallocate
1948 : ! Was already done in cp_fm_finish_copy_general
1949 17079 : IF (ALLOCATED(info%src_blacs2mpi)) THEN
1950 4672 : DEALLOCATE (info%src_blacs2mpi)
1951 : END IF
1952 17079 : CALL mp_waitall(info%send_request)
1953 17079 : DEALLOCATE (info%send_request, info%send_buf)
1954 :
1955 : END IF
1956 :
1957 17079 : CALL timestop(handle)
1958 :
1959 17079 : END SUBROUTINE cp_fm_cleanup_copy_general
1960 :
1961 : ! **************************************************************************************************
1962 : !> \brief General copy of a submatrix of fm matrix to a submatrix of another fm matrix.
1963 : !> The two matrices can have different contexts.
1964 : !>
1965 : !> Summary of distribution routines for dense matrices
1966 : !> The following will copy A(iA:iA+M-1,jA:jA+N-1) to B(iB:iB+M-1,jB:jB+N-1):
1967 : !>
1968 : !> call pdgemr2d(M,N,Aloc,iA,jA,descA,Bloc,iB,jB,descB,context)
1969 : !>
1970 : !> A process that is not a part of the context of A should set descA(2)
1971 : !> to -1, and similarly for B.
1972 : !>
1973 : !> \param source input fm matrix
1974 : !> \param destination output fm matrix
1975 : !> \param nrows number of rows of sub matrix to be copied
1976 : !> \param ncols number of cols of sub matrix to be copied
1977 : !> \param s_firstrow starting global row index of sub matrix in source
1978 : !> \param s_firstcol starting global col index of sub matrix in source
1979 : !> \param d_firstrow starting global row index of sub matrix in destination
1980 : !> \param d_firstcol starting global col index of sub matrix in destination
1981 : !> \param global_context process grid that covers all parts of either A or B.
1982 : ! **************************************************************************************************
1983 13707 : SUBROUTINE cp_fm_to_fm_submat_general(source, &
1984 : destination, &
1985 : nrows, &
1986 : ncols, &
1987 : s_firstrow, &
1988 : s_firstcol, &
1989 : d_firstrow, &
1990 : d_firstcol, &
1991 : global_context)
1992 :
1993 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1994 : INTEGER, INTENT(IN) :: nrows, ncols, s_firstrow, s_firstcol, &
1995 : d_firstrow, d_firstcol
1996 :
1997 : CLASS(cp_blacs_type), INTENT(IN) :: global_context
1998 :
1999 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat_general'
2000 :
2001 : LOGICAL :: debug
2002 : INTEGER :: handle
2003 : #if defined(__parallel)
2004 : INTEGER, DIMENSION(9) :: desca, descb
2005 : REAL(KIND=dp), DIMENSION(1, 1), TARGET :: dummy
2006 13707 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: smat, dmat
2007 : #endif
2008 :
2009 13707 : CALL timeset(routineN, handle)
2010 :
2011 13707 : debug = debug_this_module
2012 :
2013 : IF (.NOT. cp2k_is_parallel) THEN
2014 : CALL cp_fm_to_fm_submat(source, &
2015 : destination, &
2016 : nrows, &
2017 : ncols, &
2018 : s_firstrow, &
2019 : s_firstcol, &
2020 : d_firstrow, &
2021 : d_firstcol)
2022 : ELSE
2023 : #ifdef __parallel
2024 : NULLIFY (smat, dmat)
2025 : ! check whether source is available on this process
2026 13707 : IF (ASSOCIATED(source%matrix_struct)) THEN
2027 137070 : desca = source%matrix_struct%descriptor
2028 13707 : IF (nrows > source%matrix_struct%nrow_global) THEN
2029 0 : CPABORT("nrows is greater than nrow_global of source")
2030 : END IF
2031 13707 : IF (ncols > source%matrix_struct%ncol_global) THEN
2032 0 : CPABORT("ncols is greater than ncol_global of source")
2033 : END IF
2034 13707 : smat => source%local_data
2035 : ELSE
2036 0 : desca = -1
2037 0 : smat => dummy
2038 : END IF
2039 : ! check destination is available on this process
2040 13707 : IF (ASSOCIATED(destination%matrix_struct)) THEN
2041 137070 : descb = destination%matrix_struct%descriptor
2042 13707 : IF (nrows > destination%matrix_struct%nrow_global) THEN
2043 0 : CPABORT("nrows is greater than nrow_global of destination")
2044 : END IF
2045 13707 : IF (ncols > destination%matrix_struct%ncol_global) THEN
2046 0 : CPABORT("ncols is greater than ncol_global of destination")
2047 : END IF
2048 13707 : dmat => destination%local_data
2049 : ELSE
2050 0 : descb = -1
2051 0 : dmat => dummy
2052 : END IF
2053 : ! do copy
2054 :
2055 : CALL pdgemr2d(nrows, &
2056 : ncols, &
2057 : smat, &
2058 : s_firstrow, &
2059 : s_firstcol, &
2060 : desca, &
2061 : dmat, &
2062 : d_firstrow, &
2063 : d_firstcol, &
2064 : descb, &
2065 13707 : global_context%get_handle())
2066 : #else
2067 : MARK_USED(global_context)
2068 : CPABORT("this subroutine only supports SCALAPACK")
2069 : #endif
2070 : END IF
2071 :
2072 13707 : CALL timestop(handle)
2073 :
2074 13707 : END SUBROUTINE cp_fm_to_fm_submat_general
2075 :
2076 : ! **************************************************************************************************
2077 : !> \brief ...
2078 : !> \param matrix ...
2079 : !> \param irow_global ...
2080 : !> \param icol_global ...
2081 : !> \param alpha ...
2082 : ! **************************************************************************************************
2083 240 : SUBROUTINE cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
2084 :
2085 : ! Add alpha to the matrix element specified by the global indices
2086 : ! irow_global and icol_global
2087 :
2088 : ! - Creation (05.05.06,MK)
2089 :
2090 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2091 : INTEGER, INTENT(IN) :: irow_global, icol_global
2092 : REAL(KIND=dp), INTENT(IN) :: alpha
2093 :
2094 : INTEGER :: mypcol, myprow, npcol, nprow
2095 240 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2096 : TYPE(cp_blacs_env_type), POINTER :: context
2097 : #if defined(__parallel)
2098 : INTEGER :: icol_local, ipcol, iprow, &
2099 : irow_local
2100 : INTEGER, DIMENSION(9) :: desca
2101 : #endif
2102 :
2103 240 : context => matrix%matrix_struct%context
2104 :
2105 240 : myprow = context%mepos(1)
2106 240 : mypcol = context%mepos(2)
2107 :
2108 240 : nprow = context%num_pe(1)
2109 240 : npcol = context%num_pe(2)
2110 :
2111 240 : a => matrix%local_data
2112 :
2113 : #if defined(__parallel)
2114 :
2115 2400 : desca(:) = matrix%matrix_struct%descriptor(:)
2116 :
2117 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
2118 240 : irow_local, icol_local, iprow, ipcol)
2119 :
2120 240 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
2121 120 : a(irow_local, icol_local) = a(irow_local, icol_local) + alpha
2122 : END IF
2123 :
2124 : #else
2125 :
2126 : a(irow_global, icol_global) = a(irow_global, icol_global) + alpha
2127 :
2128 : #endif
2129 :
2130 240 : END SUBROUTINE cp_fm_add_to_element
2131 :
2132 : ! **************************************************************************************************
2133 : !> \brief ...
2134 : !> \param fm ...
2135 : !> \param unit ...
2136 : ! **************************************************************************************************
2137 141026 : SUBROUTINE cp_fm_write_unformatted(fm, unit)
2138 : TYPE(cp_fm_type), INTENT(IN) :: fm
2139 : INTEGER, INTENT(IN) :: unit
2140 :
2141 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_unformatted'
2142 :
2143 : INTEGER :: handle, j, max_block, &
2144 : ncol_global, nrow_global
2145 : TYPE(mp_para_env_type), POINTER :: para_env
2146 : #if defined(__parallel)
2147 : INTEGER :: i, i_block, icol_local, &
2148 : in, info, ipcol, &
2149 : iprow, irow_local, &
2150 : mepos, &
2151 : num_pe, rb, tag
2152 : INTEGER, DIMENSION(9) :: desc
2153 141026 : REAL(KIND=dp), DIMENSION(:), POINTER :: vecbuf
2154 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: newdat
2155 : TYPE(cp_blacs_type) :: ictxt_loc
2156 : INTEGER, EXTERNAL :: numroc
2157 : #endif
2158 :
2159 141026 : CALL timeset(routineN, handle)
2160 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2161 141026 : para_env=para_env)
2162 :
2163 : #if defined(__parallel)
2164 141026 : num_pe = para_env%num_pe
2165 141026 : mepos = para_env%mepos
2166 141026 : rb = nrow_global
2167 141026 : tag = 0
2168 : ! get a new context
2169 141026 : CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
2170 141026 : CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
2171 141026 : CPASSERT(info == 0)
2172 : ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
2173 : myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
2174 282052 : in = numroc(ncol_global, max_block, mypcol, 0, npcol)
2175 :
2176 564104 : ALLOCATE (newdat(nrow_global, MAX(1, in)))
2177 :
2178 : ! do the actual scalapack to cols reordering
2179 : CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
2180 : fm%matrix_struct%descriptor, &
2181 141026 : newdat, 1, 1, desc, ictxt_loc%get_handle())
2182 :
2183 423078 : ALLOCATE (vecbuf(nrow_global*max_block))
2184 50098786 : vecbuf = HUGE(1.0_dp) ! init for valgrind
2185 :
2186 622501 : DO i = 1, ncol_global, MAX(max_block, 1)
2187 340449 : i_block = MIN(max_block, ncol_global - i + 1)
2188 : CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
2189 340449 : irow_local, icol_local, iprow, ipcol)
2190 340449 : IF (ipcol == mypcol) THEN
2191 1413018 : DO j = 1, i_block
2192 120781760 : vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
2193 : END DO
2194 : END IF
2195 :
2196 340449 : IF (ipcol == 0) THEN
2197 : ! do nothing
2198 : ELSE
2199 135386 : IF (ipcol == mypcol) THEN
2200 27067761 : CALL para_env%send(vecbuf(:), 0, tag)
2201 : END IF
2202 135386 : IF (mypcol == 0) THEN
2203 54067829 : CALL para_env%recv(vecbuf(:), ipcol, tag)
2204 : END IF
2205 : END IF
2206 :
2207 481475 : IF (unit > 0) THEN
2208 1400923 : DO j = 1, i_block
2209 60980266 : WRITE (unit) vecbuf((j - 1)*nrow_global + 1:nrow_global*j)
2210 : END DO
2211 : END IF
2212 :
2213 : END DO
2214 : END ASSOCIATE
2215 141026 : DEALLOCATE (vecbuf)
2216 :
2217 141026 : CALL ictxt_loc%gridexit()
2218 :
2219 141026 : DEALLOCATE (newdat)
2220 :
2221 : #else
2222 :
2223 : IF (unit > 0) THEN
2224 : DO j = 1, ncol_global
2225 : WRITE (unit) fm%local_data(:, j)
2226 : END DO
2227 : END IF
2228 :
2229 : #endif
2230 141026 : CALL timestop(handle)
2231 :
2232 705130 : END SUBROUTINE cp_fm_write_unformatted
2233 :
2234 : ! **************************************************************************************************
2235 : !> \brief Write out a full matrix in plain text.
2236 : !> \param fm the matrix to be outputted
2237 : !> \param unit the unit number for I/O
2238 : !> \param header optional header
2239 : !> \param value_format ...
2240 : ! **************************************************************************************************
2241 3278 : SUBROUTINE cp_fm_write_formatted(fm, unit, header, value_format)
2242 : TYPE(cp_fm_type), INTENT(IN) :: fm
2243 : INTEGER, INTENT(IN) :: unit
2244 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: header, value_format
2245 :
2246 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_formatted'
2247 :
2248 : CHARACTER(LEN=21) :: my_value_format
2249 : INTEGER :: handle, i, j, max_block, &
2250 : ncol_global, nrow_global
2251 : TYPE(mp_para_env_type), POINTER :: para_env
2252 : #if defined(__parallel)
2253 : INTEGER :: i_block, icol_local, &
2254 : in, info, ipcol, &
2255 : iprow, irow_local, &
2256 : mepos, num_pe, rb, tag, k, &
2257 : icol, irow
2258 : INTEGER, DIMENSION(9) :: desc
2259 3278 : REAL(KIND=dp), DIMENSION(:), POINTER :: vecbuf
2260 3278 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: newdat
2261 : TYPE(cp_blacs_type) :: ictxt_loc
2262 : INTEGER, EXTERNAL :: numroc
2263 : #endif
2264 :
2265 3278 : CALL timeset(routineN, handle)
2266 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2267 3278 : para_env=para_env)
2268 :
2269 3278 : IF (PRESENT(value_format)) THEN
2270 0 : CPASSERT(LEN_TRIM(ADJUSTL(value_format)) < 11)
2271 0 : my_value_format = "(I10, I10, "//TRIM(ADJUSTL(value_format))//")"
2272 : ELSE
2273 3278 : my_value_format = "(I10, I10, ES24.12)"
2274 : END IF
2275 :
2276 3278 : IF (unit > 0) THEN
2277 11 : IF (PRESENT(header)) WRITE (unit, *) header
2278 11 : WRITE (unit, "(A2, A8, A10, A24)") "#", "Row", "Column", ADJUSTL("Value")
2279 : END IF
2280 :
2281 : #if defined(__parallel)
2282 3278 : num_pe = para_env%num_pe
2283 3278 : mepos = para_env%mepos
2284 3278 : rb = nrow_global
2285 3278 : tag = 0
2286 : ! get a new context
2287 3278 : CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
2288 3278 : CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
2289 3278 : CPASSERT(info == 0)
2290 : ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
2291 : myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
2292 6556 : in = numroc(ncol_global, max_block, mypcol, 0, npcol)
2293 :
2294 13112 : ALLOCATE (newdat(nrow_global, MAX(1, in)))
2295 :
2296 : ! do the actual scalapack to cols reordering
2297 : CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
2298 : fm%matrix_struct%descriptor, &
2299 3278 : newdat, 1, 1, desc, ictxt_loc%get_handle())
2300 :
2301 9834 : ALLOCATE (vecbuf(nrow_global*max_block))
2302 103876 : vecbuf = HUGE(1.0_dp) ! init for valgrind
2303 3278 : irow = 1
2304 3278 : icol = 1
2305 :
2306 15250 : DO i = 1, ncol_global, MAX(max_block, 1)
2307 8694 : i_block = MIN(max_block, ncol_global - i + 1)
2308 : CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
2309 8694 : irow_local, icol_local, iprow, ipcol)
2310 8694 : IF (ipcol == mypcol) THEN
2311 17052 : DO j = 1, i_block
2312 299892 : vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
2313 : END DO
2314 : END IF
2315 :
2316 8694 : IF (ipcol == 0) THEN
2317 : ! do nothing
2318 : ELSE
2319 3430 : IF (ipcol == mypcol) THEN
2320 50449 : CALL para_env%send(vecbuf(:), 0, tag)
2321 : END IF
2322 3430 : IF (mypcol == 0) THEN
2323 99183 : CALL para_env%recv(vecbuf(:), ipcol, tag)
2324 : END IF
2325 : END IF
2326 :
2327 11972 : IF (unit > 0) THEN
2328 210 : DO j = 1, i_block
2329 6438 : DO k = (j - 1)*nrow_global + 1, nrow_global*j
2330 6228 : WRITE (UNIT=unit, FMT=my_value_format) irow, icol, vecbuf(k)
2331 6228 : irow = irow + 1
2332 6417 : IF (irow > nrow_global) THEN
2333 189 : irow = 1
2334 189 : icol = icol + 1
2335 : END IF
2336 : END DO
2337 : END DO
2338 : END IF
2339 :
2340 : END DO
2341 : END ASSOCIATE
2342 3278 : DEALLOCATE (vecbuf)
2343 :
2344 3278 : CALL ictxt_loc%gridexit()
2345 :
2346 3278 : DEALLOCATE (newdat)
2347 :
2348 : #else
2349 :
2350 : IF (unit > 0) THEN
2351 : DO j = 1, ncol_global
2352 : DO i = 1, nrow_global
2353 : WRITE (UNIT=unit, FMT=my_value_format) i, j, fm%local_data(i, j)
2354 : END DO
2355 : END DO
2356 : END IF
2357 :
2358 : #endif
2359 3278 : CALL timestop(handle)
2360 :
2361 16390 : END SUBROUTINE cp_fm_write_formatted
2362 :
2363 : ! **************************************************************************************************
2364 : !> \brief ...
2365 : !> \param fm ...
2366 : !> \param unit ...
2367 : ! **************************************************************************************************
2368 5814 : SUBROUTINE cp_fm_read_unformatted(fm, unit)
2369 : TYPE(cp_fm_type), INTENT(INOUT) :: fm
2370 : INTEGER, INTENT(IN) :: unit
2371 :
2372 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_read_unformatted'
2373 :
2374 : INTEGER :: handle, j, max_block, &
2375 : ncol_global, nrow_global
2376 : TYPE(mp_para_env_type), POINTER :: para_env
2377 : #if defined(__parallel)
2378 : INTEGER :: k, n_cols
2379 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuf
2380 : #endif
2381 :
2382 5814 : CALL timeset(routineN, handle)
2383 :
2384 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2385 5814 : para_env=para_env)
2386 :
2387 : #if defined(__parallel)
2388 :
2389 : ! the parallel case could be made more efficient (see cp_fm_write_unformatted)
2390 :
2391 23256 : ALLOCATE (vecbuf(nrow_global, max_block))
2392 :
2393 23474 : DO j = 1, ncol_global, max_block
2394 :
2395 17660 : n_cols = MIN(max_block, ncol_global - j + 1)
2396 17660 : IF (para_env%mepos == 0) THEN
2397 40708 : DO k = 1, n_cols
2398 626116 : READ (unit) vecbuf(:, k)
2399 : END DO
2400 : END IF
2401 2775924 : CALL para_env%bcast(vecbuf, 0)
2402 23474 : CALL cp_fm_set_submatrix(fm, vecbuf, start_row=1, start_col=j, n_cols=n_cols)
2403 :
2404 : END DO
2405 :
2406 5814 : DEALLOCATE (vecbuf)
2407 :
2408 : #else
2409 :
2410 : DO j = 1, ncol_global
2411 : READ (unit) fm%local_data(:, j)
2412 : END DO
2413 :
2414 : #endif
2415 :
2416 5814 : CALL timestop(handle)
2417 :
2418 5814 : END SUBROUTINE cp_fm_read_unformatted
2419 :
2420 : ! **************************************************************************************************
2421 : !> \brief ...
2422 : !> \param mm_type ...
2423 : ! **************************************************************************************************
2424 11712 : SUBROUTINE cp_fm_setup(mm_type)
2425 : INTEGER, INTENT(IN) :: mm_type
2426 :
2427 11712 : cp_fm_mm_type = mm_type
2428 11712 : END SUBROUTINE cp_fm_setup
2429 :
2430 : ! **************************************************************************************************
2431 : !> \brief ...
2432 : !> \return ...
2433 : ! **************************************************************************************************
2434 2018518 : FUNCTION cp_fm_get_mm_type() RESULT(res)
2435 : INTEGER :: res
2436 :
2437 2018518 : res = cp_fm_mm_type
2438 2018518 : END FUNCTION cp_fm_get_mm_type
2439 :
2440 : ! **************************************************************************************************
2441 : !> \brief ...
2442 : !> \param ictxt ...
2443 : !> \param prec ...
2444 : !> \return ...
2445 : ! **************************************************************************************************
2446 10 : FUNCTION cp_fm_pilaenv(ictxt, prec) RESULT(res)
2447 : INTEGER :: ictxt
2448 : CHARACTER(LEN=1) :: prec
2449 : INTEGER :: res
2450 : #if defined(__parallel)
2451 : INTEGER :: pilaenv
2452 10 : res = pilaenv(ictxt, prec)
2453 : #else
2454 : MARK_USED(ictxt)
2455 : MARK_USED(prec)
2456 : res = -1
2457 : #endif
2458 :
2459 10 : END FUNCTION cp_fm_pilaenv
2460 :
2461 0 : END MODULE cp_fm_types
|
