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 1867020 : 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 1867020 : CALL timeset(routineN, handle)
179 :
180 1867020 : IF (PRESENT(nrow) .OR. PRESENT(ncol)) THEN
181 1512 : CALL cp_fm_struct_get(matrix_struct, nrow_global=nrow_global, ncol_global=ncol_global)
182 1512 : IF (PRESENT(nrow)) nrow_global = nrow
183 1512 : IF (PRESENT(ncol)) ncol_global = ncol
184 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=matrix_struct, &
185 1512 : nrow_global=nrow_global, ncol_global=ncol_global)
186 :
187 1512 : context => fm_struct%context
188 1512 : matrix%matrix_struct => fm_struct
189 1512 : CALL cp_fm_struct_retain(matrix%matrix_struct)
190 :
191 1512 : nrow_local = fm_struct%local_leading_dimension
192 1512 : ncol_local = MAX(1, fm_struct%ncol_locals(context%mepos(2)))
193 :
194 1512 : CALL cp_fm_struct_release(fm_struct)
195 : ELSE
196 :
197 1865508 : context => matrix_struct%context
198 1865508 : matrix%matrix_struct => matrix_struct
199 1865508 : 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 1865508 : nrow_local = matrix_struct%local_leading_dimension
205 1865508 : ncol_local = MAX(1, matrix_struct%ncol_locals(context%mepos(2)))
206 : END IF
207 :
208 1867020 : NULLIFY (matrix%local_data)
209 :
210 7468080 : ALLOCATE (matrix%local_data(nrow_local, ncol_local))
211 :
212 1867020 : IF (PRESENT(set_zero)) THEN
213 229691 : IF (set_zero) THEN
214 69097314 : matrix%local_data(1:nrow_local, 1:ncol_local) = 0.0_dp
215 : END IF
216 : END IF
217 :
218 1867020 : IF (PRESENT(name)) THEN
219 862072 : matrix%name = name
220 : ELSE
221 1004948 : matrix%name = 'full matrix'
222 : END IF
223 :
224 1867020 : CALL timestop(handle)
225 :
226 1867020 : 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 1887321 : SUBROUTINE cp_fm_release_aa0(matrix)
236 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
237 :
238 1887321 : IF (ASSOCIATED(matrix%local_data)) THEN
239 1866296 : DEALLOCATE (matrix%local_data)
240 : NULLIFY (matrix%local_data)
241 : END IF
242 1887321 : matrix%name = ""
243 1887321 : CALL cp_fm_struct_release(matrix%matrix_struct)
244 :
245 1887321 : END SUBROUTINE cp_fm_release_aa0
246 :
247 : ! **************************************************************************************************
248 : !> \brief ...
249 : !> \param matrices ...
250 : ! **************************************************************************************************
251 81706 : SUBROUTINE cp_fm_release_aa1(matrices)
252 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: matrices
253 :
254 : INTEGER :: i
255 :
256 81706 : IF (ALLOCATED(matrices)) THEN
257 217601 : DO i = 1, SIZE(matrices)
258 217601 : CALL cp_fm_release(matrices(i))
259 : END DO
260 80260 : DEALLOCATE (matrices)
261 : END IF
262 81706 : END SUBROUTINE cp_fm_release_aa1
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param matrices ...
267 : ! **************************************************************************************************
268 9142 : SUBROUTINE cp_fm_release_aa2(matrices)
269 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: matrices
270 :
271 : INTEGER :: i, j
272 :
273 9142 : IF (ALLOCATED(matrices)) THEN
274 18012 : DO i = 1, SIZE(matrices, 1)
275 49968 : DO j = 1, SIZE(matrices, 2)
276 43820 : CALL cp_fm_release(matrices(i, j))
277 : END DO
278 : END DO
279 6148 : DEALLOCATE (matrices)
280 : END IF
281 9142 : END SUBROUTINE cp_fm_release_aa2
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param matrices ...
286 : ! **************************************************************************************************
287 50 : SUBROUTINE cp_fm_release_aa3(matrices)
288 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: matrices
289 :
290 : INTEGER :: i, j, k
291 :
292 50 : IF (ALLOCATED(matrices)) THEN
293 572 : DO i = 1, SIZE(matrices, 1)
294 1644 : DO j = 1, SIZE(matrices, 2)
295 2786 : DO k = 1, SIZE(matrices, 3)
296 2264 : CALL cp_fm_release(matrices(i, j, k))
297 : END DO
298 : END DO
299 : END DO
300 50 : DEALLOCATE (matrices)
301 : END IF
302 50 : END SUBROUTINE cp_fm_release_aa3
303 :
304 : ! **************************************************************************************************
305 : !> \brief ...
306 : !> \param matrices ...
307 : ! **************************************************************************************************
308 191898 : SUBROUTINE cp_fm_release_pa1(matrices)
309 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: matrices
310 :
311 : INTEGER :: i
312 :
313 191898 : IF (ASSOCIATED(matrices)) THEN
314 148737 : DO i = 1, SIZE(matrices)
315 148737 : CALL cp_fm_release(matrices(i))
316 : END DO
317 61518 : DEALLOCATE (matrices)
318 : NULLIFY (matrices)
319 : END IF
320 191898 : END SUBROUTINE cp_fm_release_pa1
321 :
322 : ! **************************************************************************************************
323 : !> \brief ...
324 : !> \param matrices ...
325 : ! **************************************************************************************************
326 72662 : SUBROUTINE cp_fm_release_pa2(matrices)
327 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: matrices
328 :
329 : INTEGER :: i, j
330 :
331 72662 : IF (ASSOCIATED(matrices)) THEN
332 82526 : DO i = 1, SIZE(matrices, 1)
333 153494 : DO j = 1, SIZE(matrices, 2)
334 130402 : CALL cp_fm_release(matrices(i, j))
335 : END DO
336 : END DO
337 23092 : DEALLOCATE (matrices)
338 : NULLIFY (matrices)
339 : END IF
340 72662 : 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 2926 : 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 5852 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
460 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buff
461 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
462 2926 : 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 2926 : CALL timeset(routineN, handle)
468 :
469 : ! guarantee same seed over all tasks
470 2926 : CALL matrix%matrix_struct%para_env%bcast(seed, 0)
471 :
472 : rng = rng_stream_type("cp_fm_init_random_stream", distribution_type=UNIFORM, &
473 2926 : 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 2926 : row_indices=row_indices, col_indices=col_indices)
479 :
480 2926 : my_start_col = 1
481 2926 : IF (PRESENT(start_col)) my_start_col = start_col
482 2926 : my_ncol = matrix%matrix_struct%ncol_global
483 2926 : IF (PRESENT(ncol)) my_ncol = ncol
484 :
485 2926 : 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 8778 : 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 2926 : icol_global = 0
495 22879 : DO icol_local = 1, ncol_local
496 19953 : CPASSERT(col_indices(icol_local) > icol_global)
497 : DO
498 19953 : CALL rng%reset_to_next_substream()
499 19953 : icol_global = icol_global + 1
500 19953 : IF (icol_global == col_indices(icol_local)) EXIT
501 : END DO
502 19953 : CALL rng%fill(buff)
503 722599 : DO irow_local = 1, nrow_local
504 719673 : local_data(irow_local, icol_local) = buff(row_indices(irow_local))
505 : END DO
506 : END DO
507 :
508 2926 : 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 2926 : CALL rng%get(ig=seed)
518 :
519 2926 : CALL timestop(handle)
520 :
521 81928 : 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 414079 : 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 414079 : CALL timeset(routineN, handle)
544 :
545 231898538 : matrix%local_data(:, :) = alpha
546 :
547 414079 : IF (PRESENT(beta)) THEN
548 66720 : n = MIN(matrix%matrix_struct%nrow_global, matrix%matrix_struct%ncol_global)
549 515216 : DO i = 1, n
550 515216 : CALL cp_fm_set_element(matrix, i, i, beta)
551 : END DO
552 : END IF
553 :
554 414079 : CALL timestop(handle)
555 :
556 414079 : END SUBROUTINE cp_fm_set_all
557 :
558 : ! **************************************************************************************************
559 : !> \brief returns the diagonal elements of a fm
560 : !> \param matrix ...
561 : !> \param diag ...
562 : ! **************************************************************************************************
563 16116 : SUBROUTINE cp_fm_get_diag(matrix, diag)
564 :
565 : ! arguments
566 : TYPE(cp_fm_type), INTENT(IN) :: matrix
567 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
568 :
569 : ! locals
570 : INTEGER :: i, nrow_global
571 :
572 : #if defined(__parallel)
573 : INTEGER, DIMENSION(9) :: desca
574 : TYPE(cp_blacs_env_type), POINTER :: context
575 : INTEGER :: icol_local, ipcol, iprow, irow_local, mypcol, myprow, npcol, &
576 : nprow
577 16116 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
578 : #endif
579 :
580 16116 : CALL cp_fm_get_info(matrix, nrow_global=nrow_global)
581 :
582 : #if defined(__parallel)
583 182016 : diag = 0.0_dp
584 16116 : context => matrix%matrix_struct%context
585 16116 : myprow = context%mepos(1)
586 16116 : mypcol = context%mepos(2)
587 16116 : nprow = context%num_pe(1)
588 16116 : npcol = context%num_pe(2)
589 :
590 16116 : a => matrix%local_data
591 161160 : desca(:) = matrix%matrix_struct%descriptor(:)
592 :
593 182016 : DO i = 1, nrow_global
594 : CALL infog2l(i, i, desca, nprow, npcol, myprow, mypcol, &
595 165900 : irow_local, icol_local, iprow, ipcol)
596 182016 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
597 83034 : diag(i) = a(irow_local, icol_local)
598 : END IF
599 : END DO
600 : #else
601 : DO i = 1, nrow_global
602 : diag(i) = matrix%local_data(i, i)
603 : END DO
604 : #endif
605 347916 : CALL matrix%matrix_struct%para_env%sum(diag)
606 :
607 16116 : END SUBROUTINE cp_fm_get_diag
608 :
609 : ! **************************************************************************************************
610 : !> \brief returns an element of a fm
611 : !> this value is valid on every cpu
612 : !> using this call is expensive
613 : !> \param matrix the matrix to read
614 : !> \param irow_global the row
615 : !> \param icol_global the col
616 : !> \param alpha the value of matrix(irow_global, icol_global)
617 : !> \param local true if the element is on this cpu, false otherwise
618 : !> \note
619 : !> - modified semantics. now this function always returns the value
620 : !> previously the value was zero on cpus that didn't own the relevant
621 : !> part of the matrix (Joost VandeVondele, May 2003)
622 : !> - usage of the function should be avoided, as it is likely to rather slow
623 : !> using row_indices/col_indices/local_data + some smart scheme normally
624 : !> yields a real parallel code
625 : ! **************************************************************************************************
626 1193594 : SUBROUTINE cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
627 :
628 : ! arguments
629 : TYPE(cp_fm_type), INTENT(IN) :: matrix
630 : REAL(KIND=dp), INTENT(OUT) :: alpha
631 : INTEGER, INTENT(IN) :: icol_global, &
632 : irow_global
633 : LOGICAL, INTENT(OUT), OPTIONAL :: local
634 :
635 : ! locals
636 : #if defined(__parallel)
637 : INTEGER, DIMENSION(9) :: desca
638 : TYPE(cp_blacs_env_type), POINTER :: context
639 : INTEGER :: icol_local, ipcol, iprow, irow_local, mypcol, myprow, npcol, &
640 : nprow
641 1193594 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
642 : #endif
643 :
644 : #if defined(__parallel)
645 1193594 : context => matrix%matrix_struct%context
646 1193594 : myprow = context%mepos(1)
647 1193594 : mypcol = context%mepos(2)
648 1193594 : nprow = context%num_pe(1)
649 1193594 : npcol = context%num_pe(2)
650 :
651 1193594 : a => matrix%local_data
652 11935940 : desca(:) = matrix%matrix_struct%descriptor(:)
653 :
654 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
655 1193594 : irow_local, icol_local, iprow, ipcol)
656 :
657 1193594 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
658 596847 : alpha = a(irow_local, icol_local)
659 596847 : CALL context%dgebs2d('All', ' ', 1, 1, alpha, 1)
660 596847 : IF (PRESENT(local)) local = .TRUE.
661 : ELSE
662 596747 : CALL context%dgebr2d('All', ' ', 1, 1, alpha, 1, iprow, ipcol)
663 596747 : IF (PRESENT(local)) local = .FALSE.
664 : END IF
665 :
666 : #else
667 : IF (PRESENT(local)) local = .TRUE.
668 : alpha = matrix%local_data(irow_global, icol_global)
669 : #endif
670 :
671 1193594 : END SUBROUTINE cp_fm_get_element
672 :
673 : ! **************************************************************************************************
674 : !> \brief sets an element of a matrix
675 : !> \param matrix ...
676 : !> \param irow_global ...
677 : !> \param icol_global ...
678 : !> \param alpha ...
679 : !> \note
680 : !> we expect all cpus to have the same arguments in the call to this function
681 : !> (otherwise one should use local_data tricks)
682 : ! **************************************************************************************************
683 707330 : SUBROUTINE cp_fm_set_element(matrix, irow_global, icol_global, alpha)
684 : TYPE(cp_fm_type), INTENT(IN) :: matrix
685 : INTEGER, INTENT(IN) :: irow_global, icol_global
686 : REAL(KIND=dp), INTENT(IN) :: alpha
687 :
688 : INTEGER :: mypcol, myprow, npcol, nprow
689 : TYPE(cp_blacs_env_type), POINTER :: context
690 : #if defined(__parallel)
691 : INTEGER :: icol_local, ipcol, iprow, &
692 : irow_local
693 : INTEGER, DIMENSION(9) :: desca
694 707330 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
695 : #endif
696 :
697 707330 : context => matrix%matrix_struct%context
698 707330 : myprow = context%mepos(1)
699 707330 : mypcol = context%mepos(2)
700 707330 : nprow = context%num_pe(1)
701 707330 : npcol = context%num_pe(2)
702 :
703 : #if defined(__parallel)
704 :
705 707330 : a => matrix%local_data
706 :
707 7073300 : desca(:) = matrix%matrix_struct%descriptor(:)
708 :
709 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
710 707330 : irow_local, icol_local, iprow, ipcol)
711 :
712 707330 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
713 355570 : a(irow_local, icol_local) = alpha
714 : END IF
715 :
716 : #else
717 :
718 : matrix%local_data(irow_global, icol_global) = alpha
719 :
720 : #endif
721 707330 : END SUBROUTINE cp_fm_set_element
722 :
723 : ! **************************************************************************************************
724 : !> \brief sets a submatrix of a full matrix
725 : !> fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
726 : !> = alpha*op(new_values)(1:n_rows,1:n_cols)+ beta
727 : !> * fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
728 : !> \param fm the full to change
729 : !> \param new_values a replicated full matrix with the new values
730 : !> \param start_row the starting row of b_matrix (defaults to 1)
731 : !> \param start_col the starting col of b_matrix (defaults to 1)
732 : !> \param n_rows the number of row to change in b (defaults to
733 : !> size(op(new_values),1))
734 : !> \param n_cols the number of columns to change in b (defaults to
735 : !> size(op(new_values),2))
736 : !> \param alpha rescaling factor for the new values (defaults to 1.0)
737 : !> \param beta rescaling factor for the old values (defaults to 0.0)
738 : !> \param transpose if new_values should be transposed: if true
739 : !> op(new_values)=new_values^T, else op(new_values)=new_values
740 : !> (defaults to false)
741 : !> \par History
742 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
743 : !> \author Fawzi Mohamed
744 : !> \note
745 : !> optimized for full column updates and alpha=1.0, beta=0.0
746 : !> the new_values need to be valid on all cpus
747 : ! **************************************************************************************************
748 70643 : SUBROUTINE cp_fm_set_submatrix(fm, new_values, start_row, &
749 : start_col, n_rows, n_cols, alpha, beta, transpose)
750 : TYPE(cp_fm_type), INTENT(IN) :: fm
751 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: new_values
752 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
753 : REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha, beta
754 : LOGICAL, INTENT(in), OPTIONAL :: transpose
755 :
756 : INTEGER :: i, i0, j, j0, ncol, ncol_block, &
757 : ncol_global, ncol_local, nrow, &
758 : nrow_block, nrow_global, nrow_local, &
759 : this_col, this_row
760 70643 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
761 : LOGICAL :: tr_a
762 : REAL(KIND=dp) :: al, be
763 70643 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
764 :
765 70643 : al = 1.0_dp; be = 0.0_dp; i0 = 1; j0 = 1; tr_a = .FALSE.
766 :
767 19356 : IF (PRESENT(alpha)) al = alpha
768 70643 : IF (PRESENT(beta)) be = beta
769 70643 : IF (PRESENT(start_row)) i0 = start_row
770 70643 : IF (PRESENT(start_col)) j0 = start_col
771 70643 : IF (PRESENT(transpose)) tr_a = transpose
772 16153 : IF (tr_a) THEN
773 16061 : nrow = SIZE(new_values, 2)
774 16061 : ncol = SIZE(new_values, 1)
775 : ELSE
776 54582 : nrow = SIZE(new_values, 1)
777 54582 : ncol = SIZE(new_values, 2)
778 : END IF
779 70643 : IF (PRESENT(n_rows)) nrow = n_rows
780 70643 : IF (PRESENT(n_cols)) ncol = n_cols
781 :
782 70643 : full_block => fm%local_data
783 :
784 : CALL cp_fm_get_info(matrix=fm, &
785 : nrow_global=nrow_global, ncol_global=ncol_global, &
786 : nrow_block=nrow_block, ncol_block=ncol_block, &
787 : nrow_local=nrow_local, ncol_local=ncol_local, &
788 70643 : row_indices=row_indices, col_indices=col_indices)
789 :
790 70643 : IF (al == 1.0 .AND. be == 0.0) THEN
791 1290994 : DO j = 1, ncol_local
792 1229937 : this_col = col_indices(j) - j0 + 1
793 1290994 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
794 296459 : IF (tr_a) THEN
795 6475 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
796 158370 : DO i = 1, nrow_local
797 158370 : full_block(i, j) = new_values(this_col, row_indices(i))
798 : END DO
799 : ELSE
800 594 : DO i = 1, nrow_local
801 510 : this_row = row_indices(i) - i0 + 1
802 594 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
803 255 : full_block(i, j) = new_values(this_col, this_row)
804 : END IF
805 : END DO
806 : END IF
807 : ELSE
808 289984 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
809 6955169 : DO i = 1, nrow_local
810 6955169 : full_block(i, j) = new_values(row_indices(i), this_col)
811 : END DO
812 : ELSE
813 538491 : DO i = 1, nrow_local
814 528241 : this_row = row_indices(i) - i0 + 1
815 538491 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
816 58493 : full_block(i, j) = new_values(this_row, this_col)
817 : END IF
818 : END DO
819 : END IF
820 : END IF
821 : END IF
822 : END DO
823 : ELSE
824 838496 : DO j = 1, ncol_local
825 828910 : this_col = col_indices(j) - j0 + 1
826 838496 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
827 828910 : IF (tr_a) THEN
828 88287891 : DO i = 1, nrow_local
829 87458981 : this_row = row_indices(i) - i0 + 1
830 88287891 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
831 : full_block(i, j) = al*new_values(this_col, this_row) + &
832 414455 : be*full_block(i, j)
833 : END IF
834 : END DO
835 : ELSE
836 0 : DO i = 1, nrow_local
837 0 : this_row = row_indices(i) - i0 + 1
838 0 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
839 : full_block(i, j) = al*new_values(this_row, this_col) + &
840 0 : be*full_block(i, j)
841 : END IF
842 : END DO
843 : END IF
844 : END IF
845 : END DO
846 : END IF
847 :
848 70643 : END SUBROUTINE cp_fm_set_submatrix
849 :
850 : ! **************************************************************************************************
851 : !> \brief sets a submatrix of a full matrix to a given value
852 : !> fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = value
853 : !> \param fm the full to change
854 : !> \param new_value ...
855 : !> \param start_row the starting row of matrix
856 : !> \param start_col the starting col of matrix
857 : !> \param n_rows the number of rows to change
858 : !> \param n_cols the number of columns to change
859 : !> \par History
860 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
861 : !> 12.2025 created from cp_fm_set_submatrix
862 : !> \author JGH
863 : ! **************************************************************************************************
864 1102292 : SUBROUTINE cp_fm_set_all_submatrix(fm, new_value, start_row, start_col, n_rows, n_cols)
865 : TYPE(cp_fm_type), INTENT(IN) :: fm
866 : REAL(KIND=dp), INTENT(in) :: new_value
867 : INTEGER, INTENT(in) :: start_row, start_col, n_rows, n_cols
868 :
869 : INTEGER :: i, i0, j, j0, ncol_global, ncol_local, &
870 : nrow_global, nrow_local, this_col, &
871 : this_row
872 1102292 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
873 1102292 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
874 :
875 1102292 : full_block => fm%local_data
876 :
877 : CALL cp_fm_get_info(matrix=fm, &
878 : nrow_global=nrow_global, ncol_global=ncol_global, &
879 : nrow_local=nrow_local, ncol_local=ncol_local, &
880 1102292 : row_indices=row_indices, col_indices=col_indices)
881 :
882 1102292 : i0 = start_row
883 1102292 : j0 = start_col
884 20761024 : DO j = 1, ncol_local
885 19658732 : this_col = col_indices(j) - j0 + 1
886 20761024 : IF (this_col >= 1 .AND. this_col <= n_cols) THEN
887 639770447 : DO i = 1, nrow_local
888 620171677 : this_row = row_indices(i) - i0 + 1
889 639770447 : IF (this_row >= 1 .AND. this_row <= n_rows) THEN
890 612983463 : full_block(i, j) = new_value
891 : END IF
892 : END DO
893 : END IF
894 : END DO
895 :
896 1102292 : END SUBROUTINE cp_fm_set_all_submatrix
897 :
898 : ! **************************************************************************************************
899 : !> \brief gets a submatrix of a full matrix
900 : !> op(target_m)(1:n_rows,1:n_cols)
901 : !> =fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
902 : !> target_m is replicated on all cpus
903 : !> using this call is expensive
904 : !> \param fm the full you want to get the info from
905 : !> \param target_m a replicated full matrix that will contain the result
906 : !> \param start_row the starting row of b_matrix (defaults to 1)
907 : !> \param start_col the starting col of b_matrix (defaults to 1)
908 : !> \param n_rows the number of row to change in b (defaults to
909 : !> size(op(new_values),1))
910 : !> \param n_cols the number of columns to change in b (defaults to
911 : !> size(op(new_values),2))
912 : !> \param transpose if target_m should be transposed: if true
913 : !> op(target_m)=target_m^T, else op(target_m)=target_m
914 : !> (defaults to false)
915 : !> \par History
916 : !> 07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
917 : !> \author Fawzi Mohamed
918 : !> \note
919 : !> optimized for full column updates. Zeros out a little too much
920 : !> of target_m
921 : !> the target_m is replicated and valid on all cpus
922 : ! **************************************************************************************************
923 100816 : SUBROUTINE cp_fm_get_submatrix(fm, target_m, start_row, &
924 : start_col, n_rows, n_cols, transpose)
925 : TYPE(cp_fm_type), INTENT(IN) :: fm
926 : REAL(KIND=dp), DIMENSION(:, :), INTENT(out) :: target_m
927 : INTEGER, INTENT(in), OPTIONAL :: start_row, start_col, n_rows, n_cols
928 : LOGICAL, INTENT(in), OPTIONAL :: transpose
929 :
930 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_get_submatrix'
931 :
932 : INTEGER :: handle, i, i0, j, j0, ncol, ncol_global, &
933 : ncol_local, nrow, nrow_global, &
934 : nrow_local, this_col, this_row
935 100816 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
936 : LOGICAL :: tr_a
937 100816 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: full_block
938 : TYPE(mp_para_env_type), POINTER :: para_env
939 :
940 100816 : CALL timeset(routineN, handle)
941 :
942 100816 : i0 = 1; j0 = 1; tr_a = .FALSE.
943 :
944 100816 : IF (PRESENT(start_row)) i0 = start_row
945 100816 : IF (PRESENT(start_col)) j0 = start_col
946 100816 : IF (PRESENT(transpose)) tr_a = transpose
947 6152 : IF (tr_a) THEN
948 2206 : nrow = SIZE(target_m, 2)
949 2206 : ncol = SIZE(target_m, 1)
950 : ELSE
951 98610 : nrow = SIZE(target_m, 1)
952 98610 : ncol = SIZE(target_m, 2)
953 : END IF
954 100816 : IF (PRESENT(n_rows)) nrow = n_rows
955 100816 : IF (PRESENT(n_cols)) ncol = n_cols
956 :
957 100816 : para_env => fm%matrix_struct%para_env
958 :
959 100816 : full_block => fm%local_data
960 : #if defined(__parallel)
961 : ! zero-out whole target_m
962 100816 : IF (SIZE(target_m, 1)*SIZE(target_m, 2) /= 0) THEN
963 119008 : CALL dcopy(SIZE(target_m, 1)*SIZE(target_m, 2), [0.0_dp], 0, target_m, 1)
964 : END IF
965 : #endif
966 :
967 : CALL cp_fm_get_info(matrix=fm, &
968 : nrow_global=nrow_global, ncol_global=ncol_global, &
969 : nrow_local=nrow_local, ncol_local=ncol_local, &
970 100816 : row_indices=row_indices, col_indices=col_indices)
971 :
972 616874 : DO j = 1, ncol_local
973 516058 : this_col = col_indices(j) - j0 + 1
974 616874 : IF (this_col >= 1 .AND. this_col <= ncol) THEN
975 377070 : IF (tr_a) THEN
976 2206 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
977 80516 : DO i = 1, nrow_local
978 80516 : target_m(this_col, row_indices(i)) = full_block(i, j)
979 : END DO
980 : ELSE
981 0 : DO i = 1, nrow_local
982 0 : this_row = row_indices(i) - i0 + 1
983 0 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
984 0 : target_m(this_col, this_row) = full_block(i, j)
985 : END IF
986 : END DO
987 : END IF
988 : ELSE
989 374864 : IF (i0 == 1 .AND. nrow_global == nrow) THEN
990 7989133 : DO i = 1, nrow_local
991 7989133 : target_m(row_indices(i), this_col) = full_block(i, j)
992 : END DO
993 : ELSE
994 1735300 : DO i = 1, nrow_local
995 1694848 : this_row = row_indices(i) - i0 + 1
996 1735300 : IF (this_row >= 1 .AND. this_row <= nrow) THEN
997 85570 : target_m(this_row, this_col) = full_block(i, j)
998 : END IF
999 : END DO
1000 : END IF
1001 : END IF
1002 : END IF
1003 : END DO
1004 :
1005 23827844 : CALL para_env%sum(target_m)
1006 :
1007 100816 : CALL timestop(handle)
1008 :
1009 100816 : END SUBROUTINE cp_fm_get_submatrix
1010 :
1011 : ! **************************************************************************************************
1012 : !> \brief returns all kind of information about the full matrix
1013 : !> \param matrix ...
1014 : !> \param name ...
1015 : !> \param nrow_global ...
1016 : !> \param ncol_global ...
1017 : !> \param nrow_block ...
1018 : !> \param ncol_block ...
1019 : !> \param nrow_local ...
1020 : !> \param ncol_local ...
1021 : !> \param row_indices ...
1022 : !> \param col_indices ...
1023 : !> \param local_data ...
1024 : !> \param context ...
1025 : !> \param nrow_locals ...
1026 : !> \param ncol_locals ...
1027 : !> \param matrix_struct ...
1028 : !> \param para_env ...
1029 : !> \note
1030 : !> see also cp_fm_struct for explanation
1031 : !> - nrow_local, ncol_local, row_indices, col_indices, local_data are hooks for efficient
1032 : !> access to the local blacs block
1033 : ! **************************************************************************************************
1034 6380303 : SUBROUTINE cp_fm_get_info(matrix, name, nrow_global, ncol_global, &
1035 : nrow_block, ncol_block, nrow_local, ncol_local, &
1036 : row_indices, col_indices, local_data, context, &
1037 : nrow_locals, ncol_locals, matrix_struct, para_env)
1038 :
1039 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1040 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: name
1041 : INTEGER, INTENT(OUT), OPTIONAL :: nrow_global, ncol_global, nrow_block, &
1042 : ncol_block, nrow_local, ncol_local
1043 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: row_indices, col_indices
1044 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1045 : OPTIONAL, POINTER :: local_data
1046 : TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: context
1047 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nrow_locals, ncol_locals
1048 : TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: matrix_struct
1049 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
1050 :
1051 8 : IF (PRESENT(name)) name = matrix%name
1052 6380303 : IF (PRESENT(matrix_struct)) matrix_struct => matrix%matrix_struct
1053 6380303 : IF (PRESENT(local_data)) local_data => matrix%local_data ! not hiding things anymore :-(
1054 :
1055 : CALL cp_fm_struct_get(matrix%matrix_struct, nrow_local=nrow_local, &
1056 : ncol_local=ncol_local, nrow_global=nrow_global, &
1057 : ncol_global=ncol_global, nrow_block=nrow_block, &
1058 : ncol_block=ncol_block, row_indices=row_indices, &
1059 : col_indices=col_indices, nrow_locals=nrow_locals, &
1060 6380303 : ncol_locals=ncol_locals, context=context, para_env=para_env)
1061 :
1062 6380303 : END SUBROUTINE cp_fm_get_info
1063 :
1064 : ! **************************************************************************************************
1065 : !> \brief Write nicely formatted info about the FM to the given I/O unit (including the underlying FM struct)
1066 : !> \param matrix a cp_fm_type instance
1067 : !> \param io_unit the I/O unit to use for writing
1068 : ! **************************************************************************************************
1069 3 : SUBROUTINE cp_fm_write_info(matrix, io_unit)
1070 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1071 : INTEGER, INTENT(IN) :: io_unit
1072 :
1073 3 : WRITE (io_unit, '(/,A,A12)') "CP_FM | Name: ", matrix%name
1074 3 : CALL cp_fm_struct_write_info(matrix%matrix_struct, io_unit)
1075 3 : END SUBROUTINE cp_fm_write_info
1076 :
1077 : ! **************************************************************************************************
1078 : !> \brief find the maximum absolute value of the matrix element
1079 : !> maxval(abs(matrix))
1080 : !> \param matrix ...
1081 : !> \param a_max ...
1082 : !> \param ir_max ...
1083 : !> \param ic_max ...
1084 : ! **************************************************************************************************
1085 127407 : SUBROUTINE cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
1086 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1087 : REAL(KIND=dp), INTENT(OUT) :: a_max
1088 : INTEGER, INTENT(OUT), OPTIONAL :: ir_max, ic_max
1089 :
1090 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_maxabsval'
1091 :
1092 : INTEGER :: handle, i, ic_max_local, ir_max_local, &
1093 : j, mepos, ncol_local, nrow_local, &
1094 : num_pe
1095 127407 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ic_max_vec, ir_max_vec
1096 127407 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1097 : REAL(dp) :: my_max
1098 127407 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_max_vec
1099 127407 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1100 :
1101 127407 : CALL timeset(routineN, handle)
1102 :
1103 127407 : my_block => matrix%local_data
1104 :
1105 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
1106 127407 : row_indices=row_indices, col_indices=col_indices)
1107 :
1108 77958727 : a_max = MAXVAL(ABS(my_block(1:nrow_local, 1:ncol_local)))
1109 :
1110 127407 : IF (PRESENT(ir_max)) THEN
1111 0 : num_pe = matrix%matrix_struct%para_env%num_pe
1112 0 : mepos = matrix%matrix_struct%para_env%mepos
1113 0 : ALLOCATE (ir_max_vec(0:num_pe - 1))
1114 0 : ir_max_vec(0:num_pe - 1) = 0
1115 0 : ALLOCATE (ic_max_vec(0:num_pe - 1))
1116 0 : ic_max_vec(0:num_pe - 1) = 0
1117 0 : ALLOCATE (a_max_vec(0:num_pe - 1))
1118 0 : a_max_vec(0:num_pe - 1) = 0.0_dp
1119 0 : my_max = 0.0_dp
1120 :
1121 0 : IF ((ncol_local > 0) .AND. (nrow_local > 0)) THEN
1122 0 : DO i = 1, ncol_local
1123 0 : DO j = 1, nrow_local
1124 0 : IF (ABS(my_block(j, i)) > my_max) THEN
1125 0 : my_max = my_block(j, i)
1126 0 : ir_max_local = j
1127 0 : ic_max_local = i
1128 : END IF
1129 : END DO
1130 : END DO
1131 :
1132 0 : a_max_vec(mepos) = my_max
1133 0 : ir_max_vec(mepos) = row_indices(ir_max_local)
1134 0 : ic_max_vec(mepos) = col_indices(ic_max_local)
1135 :
1136 : END IF
1137 :
1138 0 : CALL matrix%matrix_struct%para_env%sum(a_max_vec)
1139 0 : CALL matrix%matrix_struct%para_env%sum(ir_max_vec)
1140 0 : CALL matrix%matrix_struct%para_env%sum(ic_max_vec)
1141 :
1142 0 : my_max = 0.0_dp
1143 0 : DO i = 0, num_pe - 1
1144 0 : IF (a_max_vec(i) > my_max) THEN
1145 0 : ir_max = ir_max_vec(i)
1146 0 : ic_max = ic_max_vec(i)
1147 : END IF
1148 : END DO
1149 :
1150 0 : DEALLOCATE (ir_max_vec, ic_max_vec, a_max_vec)
1151 0 : CPASSERT(ic_max > 0)
1152 0 : CPASSERT(ir_max > 0)
1153 :
1154 : END IF
1155 :
1156 127407 : CALL matrix%matrix_struct%para_env%max(a_max)
1157 :
1158 127407 : CALL timestop(handle)
1159 :
1160 254814 : END SUBROUTINE cp_fm_maxabsval
1161 :
1162 : ! **************************************************************************************************
1163 : !> \brief find the maximum over the rows of the sum of the absolute values of the elements of a given row
1164 : !> = || A ||_infinity
1165 : !> \param matrix ...
1166 : !> \param a_max ...
1167 : !> \note
1168 : !> for a real symmetric matrix it holds that || A ||_2 = |lambda_max| < || A ||_infinity
1169 : !> Hence this can be used to estimate an upper bound for the eigenvalues of a matrix
1170 : !> http://mathworld.wolfram.com/MatrixNorm.html
1171 : !> (but the bound is not so tight in the general case)
1172 : ! **************************************************************************************************
1173 4986 : SUBROUTINE cp_fm_maxabsrownorm(matrix, a_max)
1174 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1175 : REAL(KIND=dp), INTENT(OUT) :: a_max
1176 :
1177 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_maxabsrownorm'
1178 :
1179 : INTEGER :: handle, i, j, ncol_local, nrow_global, &
1180 : nrow_local
1181 4986 : INTEGER, DIMENSION(:), POINTER :: row_indices
1182 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: values
1183 4986 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1184 :
1185 4986 : CALL timeset(routineN, handle)
1186 :
1187 4986 : my_block => matrix%local_data
1188 :
1189 : CALL cp_fm_get_info(matrix, row_indices=row_indices, nrow_global=nrow_global, &
1190 4986 : nrow_local=nrow_local, ncol_local=ncol_local)
1191 :
1192 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1193 14958 : ALLOCATE (values(nrow_global))
1194 4986 : values = 0.0_dp
1195 65834 : DO j = 1, ncol_local
1196 535250 : DO i = 1, nrow_local
1197 530264 : values(row_indices(i)) = values(row_indices(i)) + ABS(my_block(i, j))
1198 : END DO
1199 : END DO
1200 4986 : CALL matrix%matrix_struct%para_env%sum(values)
1201 65834 : a_max = MAXVAL(values)
1202 4986 : DEALLOCATE (values)
1203 :
1204 4986 : CALL timestop(handle)
1205 4986 : END SUBROUTINE cp_fm_maxabsrownorm
1206 :
1207 : ! **************************************************************************************************
1208 : !> \brief find the inorm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
1209 : !> \param matrix ...
1210 : !> \param norm_array ...
1211 : ! **************************************************************************************************
1212 1316 : SUBROUTINE cp_fm_vectorsnorm(matrix, norm_array)
1213 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1214 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: norm_array
1215 :
1216 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_vectorsnorm'
1217 :
1218 : INTEGER :: handle, i, j, ncol_global, ncol_local, &
1219 : nrow_local
1220 1316 : INTEGER, DIMENSION(:), POINTER :: col_indices
1221 1316 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1222 :
1223 1316 : CALL timeset(routineN, handle)
1224 :
1225 1316 : my_block => matrix%local_data
1226 :
1227 : CALL cp_fm_get_info(matrix, col_indices=col_indices, ncol_global=ncol_global, &
1228 1316 : nrow_local=nrow_local, ncol_local=ncol_local)
1229 :
1230 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1231 31004 : norm_array = 0.0_dp
1232 31004 : DO j = 1, ncol_local
1233 1325922 : DO i = 1, nrow_local
1234 1324606 : norm_array(col_indices(j)) = norm_array(col_indices(j)) + my_block(i, j)*my_block(i, j)
1235 : END DO
1236 : END DO
1237 60692 : CALL matrix%matrix_struct%para_env%sum(norm_array)
1238 31004 : norm_array = SQRT(norm_array)
1239 :
1240 1316 : CALL timestop(handle)
1241 1316 : END SUBROUTINE cp_fm_vectorsnorm
1242 :
1243 : ! **************************************************************************************************
1244 : !> \brief summing up all the elements along the matrix's i-th index
1245 : !> \f$ \mathrm{sum}_{j} = \sum_{i} A_{ij} \f$
1246 : !> or
1247 : !> \f$ \mathrm{sum}_{i} = \sum_{j} A_{ij} \f$
1248 : !> \param matrix an input matrix A
1249 : !> \param sum_array sums of elements in each column/row
1250 : !> \param dir ...
1251 : !> \note forked from cp_fm_vectorsnorm() to be used with
1252 : !> the maximum overlap method
1253 : !> added row variation
1254 : ! **************************************************************************************************
1255 11782 : SUBROUTINE cp_fm_vectorssum(matrix, sum_array, dir)
1256 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1257 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: sum_array
1258 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: dir
1259 :
1260 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_vectorssum'
1261 :
1262 : INTEGER :: handle, i, j, ncol_local, nrow_local
1263 11782 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1264 : LOGICAL :: docol
1265 11782 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_block
1266 :
1267 11782 : CALL timeset(routineN, handle)
1268 :
1269 11782 : IF (PRESENT(dir)) THEN
1270 11742 : IF (dir == 'c' .OR. dir == 'C') THEN
1271 : docol = .TRUE.
1272 : ELSE IF (dir == 'r' .OR. dir == 'R') THEN
1273 : docol = .FALSE.
1274 : ELSE
1275 0 : CPABORT('Wrong argument DIR')
1276 : END IF
1277 : ELSE
1278 : docol = .TRUE.
1279 : END IF
1280 :
1281 11782 : my_block => matrix%local_data
1282 :
1283 : CALL cp_fm_get_info(matrix, col_indices=col_indices, row_indices=row_indices, &
1284 11782 : nrow_local=nrow_local, ncol_local=ncol_local)
1285 :
1286 : ! the efficiency could be improved by making use of the row-col distribution of scalapack
1287 286428 : sum_array(:) = 0.0_dp
1288 11782 : IF (docol) THEN
1289 448 : DO j = 1, ncol_local
1290 3628 : DO i = 1, nrow_local
1291 3588 : sum_array(col_indices(j)) = sum_array(col_indices(j)) + my_block(i, j)
1292 : END DO
1293 : END DO
1294 : ELSE
1295 102950 : DO j = 1, ncol_local
1296 6528294 : DO i = 1, nrow_local
1297 6516552 : sum_array(row_indices(i)) = sum_array(row_indices(i)) + my_block(i, j)
1298 : END DO
1299 : END DO
1300 : END IF
1301 561074 : CALL matrix%matrix_struct%para_env%sum(sum_array)
1302 :
1303 11782 : CALL timestop(handle)
1304 11782 : END SUBROUTINE cp_fm_vectorssum
1305 :
1306 : ! **************************************************************************************************
1307 : !> \brief copy one identically sized matrix in the other
1308 : !> \param source ...
1309 : !> \param destination ...
1310 : !> \note
1311 : !> see also cp_fm_to_fm_columns
1312 : ! **************************************************************************************************
1313 584285 : SUBROUTINE cp_fm_to_fm_matrix(source, destination)
1314 :
1315 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1316 :
1317 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_matrix'
1318 :
1319 : INTEGER :: handle, npcol, nprow
1320 :
1321 584285 : CALL timeset(routineN, handle)
1322 :
1323 584285 : nprow = source%matrix_struct%context%num_pe(1)
1324 584285 : npcol = source%matrix_struct%context%num_pe(2)
1325 :
1326 584285 : IF ((.NOT. cp2k_is_parallel) .OR. &
1327 : cp_fm_struct_equivalent(source%matrix_struct, &
1328 : destination%matrix_struct)) THEN
1329 584285 : IF (SIZE(source%local_data, 1) /= SIZE(destination%local_data, 1) .OR. &
1330 : SIZE(source%local_data, 2) /= SIZE(destination%local_data, 2)) THEN
1331 : CALL cp_abort(__LOCATION__, &
1332 : "Cannot copy full matrix <"//TRIM(source%name)// &
1333 : "> to full matrix <"//TRIM(destination%name)// &
1334 0 : ">. The local_data blocks have different sizes.")
1335 : END IF
1336 : CALL dcopy(SIZE(source%local_data, 1)*SIZE(source%local_data, 2), &
1337 584285 : source%local_data, 1, destination%local_data, 1)
1338 : ELSE
1339 0 : CPABORT("Data structures of source and target full matrix are not equivalent")
1340 : END IF
1341 :
1342 584285 : CALL timestop(handle)
1343 :
1344 584285 : END SUBROUTINE cp_fm_to_fm_matrix
1345 :
1346 : ! **************************************************************************************************
1347 : !> \brief copy just a subset of columns of a fm to a fm
1348 : !> \param msource ...
1349 : !> \param mtarget ...
1350 : !> \param ncol ...
1351 : !> \param source_start ...
1352 : !> \param target_start ...
1353 : ! **************************************************************************************************
1354 173124 : SUBROUTINE cp_fm_to_fm_columns(msource, mtarget, ncol, source_start, &
1355 : target_start)
1356 :
1357 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1358 : INTEGER, INTENT(IN) :: ncol
1359 : INTEGER, INTENT(IN), OPTIONAL :: source_start, target_start
1360 :
1361 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_columns'
1362 :
1363 : INTEGER :: handle, n, ss, ts
1364 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1365 : #if defined(__parallel)
1366 : INTEGER :: i
1367 : INTEGER, DIMENSION(9) :: desca, descb
1368 : #endif
1369 :
1370 173124 : CALL timeset(routineN, handle)
1371 :
1372 173124 : ss = 1
1373 173124 : ts = 1
1374 :
1375 173124 : IF (PRESENT(source_start)) ss = source_start
1376 173124 : IF (PRESENT(target_start)) ts = target_start
1377 :
1378 173124 : n = msource%matrix_struct%nrow_global
1379 :
1380 173124 : a => msource%local_data
1381 173124 : b => mtarget%local_data
1382 :
1383 : #if defined(__parallel)
1384 1731240 : desca(:) = msource%matrix_struct%descriptor(:)
1385 1731240 : descb(:) = mtarget%matrix_struct%descriptor(:)
1386 744908 : DO i = 0, ncol - 1
1387 744908 : CALL pdcopy(n, a, 1, ss + i, desca, 1, b, 1, ts + i, descb, 1)
1388 : END DO
1389 : #else
1390 : IF (ss <= SIZE(a, 2) .AND. ts <= SIZE(b, 2)) THEN
1391 : CALL dcopy(ncol*n, a(:, ss), 1, b(:, ts), 1)
1392 : END IF
1393 : #endif
1394 :
1395 173124 : CALL timestop(handle)
1396 :
1397 173124 : END SUBROUTINE cp_fm_to_fm_columns
1398 :
1399 : ! **************************************************************************************************
1400 : !> \brief copy just a triangular matrix
1401 : !> \param msource ...
1402 : !> \param mtarget ...
1403 : !> \param uplo ...
1404 : ! **************************************************************************************************
1405 58 : SUBROUTINE cp_fm_to_fm_triangular(msource, mtarget, uplo)
1406 :
1407 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1408 : CHARACTER(LEN=1), OPTIONAL, INTENT(IN) :: uplo
1409 :
1410 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_triangular'
1411 :
1412 : CHARACTER(LEN=1) :: myuplo
1413 : INTEGER :: handle, ncol, nrow
1414 58 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1415 : #if defined(__parallel)
1416 : INTEGER, DIMENSION(9) :: desca, descb
1417 : #endif
1418 :
1419 58 : CALL timeset(routineN, handle)
1420 :
1421 58 : myuplo = 'U'
1422 58 : IF (PRESENT(uplo)) myuplo = uplo
1423 :
1424 58 : nrow = msource%matrix_struct%nrow_global
1425 58 : ncol = msource%matrix_struct%ncol_global
1426 :
1427 58 : a => msource%local_data
1428 58 : b => mtarget%local_data
1429 :
1430 : #if defined(__parallel)
1431 580 : desca(:) = msource%matrix_struct%descriptor(:)
1432 580 : descb(:) = mtarget%matrix_struct%descriptor(:)
1433 58 : CALL pdlacpy(myuplo, nrow, ncol, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb)
1434 : #else
1435 : CALL dlacpy(myuplo, nrow, ncol, a(1, 1), nrow, b(1, 1), nrow)
1436 : #endif
1437 :
1438 58 : CALL timestop(handle)
1439 :
1440 58 : END SUBROUTINE cp_fm_to_fm_triangular
1441 :
1442 : ! **************************************************************************************************
1443 : !> \brief copy just a part ot the matrix
1444 : !> \param msource ...
1445 : !> \param mtarget ...
1446 : !> \param nrow ...
1447 : !> \param ncol ...
1448 : !> \param s_firstrow ...
1449 : !> \param s_firstcol ...
1450 : !> \param t_firstrow ...
1451 : !> \param t_firstcol ...
1452 : ! **************************************************************************************************
1453 :
1454 16396 : SUBROUTINE cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
1455 :
1456 : TYPE(cp_fm_type), INTENT(IN) :: msource, mtarget
1457 : INTEGER, INTENT(IN) :: nrow, ncol, s_firstrow, &
1458 : s_firstcol, t_firstrow, &
1459 : t_firstcol
1460 :
1461 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat'
1462 :
1463 : INTEGER :: handle, i, na, nb, ss, ts
1464 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
1465 : #if defined(__parallel)
1466 : INTEGER, DIMENSION(9) :: desca, descb
1467 : #endif
1468 :
1469 16396 : CALL timeset(routineN, handle)
1470 :
1471 16396 : a => msource%local_data
1472 16396 : b => mtarget%local_data
1473 :
1474 16396 : na = msource%matrix_struct%nrow_global
1475 16396 : nb = mtarget%matrix_struct%nrow_global
1476 : ! nrow must be <= na and nb
1477 16396 : IF (nrow > na) THEN
1478 0 : CPABORT("cannot copy because nrow > number of rows of source matrix")
1479 : END IF
1480 16396 : IF (nrow > nb) THEN
1481 0 : CPABORT("cannot copy because nrow > number of rows of target matrix")
1482 : END IF
1483 16396 : na = msource%matrix_struct%ncol_global
1484 16396 : nb = mtarget%matrix_struct%ncol_global
1485 : ! ncol must be <= na_col and nb_col
1486 16396 : IF (ncol > na) THEN
1487 0 : CPABORT("cannot copy because nrow > number of rows of source matrix")
1488 : END IF
1489 16396 : IF (ncol > nb) THEN
1490 0 : CPABORT("cannot copy because nrow > number of rows of target matrix")
1491 : END IF
1492 :
1493 : #if defined(__parallel)
1494 163960 : desca(:) = msource%matrix_struct%descriptor(:)
1495 163960 : descb(:) = mtarget%matrix_struct%descriptor(:)
1496 184512 : DO i = 0, ncol - 1
1497 168116 : ss = s_firstcol + i
1498 168116 : ts = t_firstcol + i
1499 184512 : CALL pdcopy(nrow, a, s_firstrow, ss, desca, 1, b, t_firstrow, ts, descb, 1)
1500 : END DO
1501 : #else
1502 : DO i = 0, ncol - 1
1503 : ss = s_firstcol + i
1504 : ts = t_firstcol + i
1505 : CALL dcopy(nrow, a(s_firstrow:, ss), 1, b(t_firstrow:, ts), 1)
1506 : END DO
1507 : #endif
1508 :
1509 16396 : CALL timestop(handle)
1510 16396 : END SUBROUTINE cp_fm_to_fm_submat
1511 :
1512 : ! **************************************************************************************************
1513 : !> \brief General copy of a fm matrix to another fm matrix.
1514 : !> Uses non-blocking MPI rather than ScaLAPACK.
1515 : !>
1516 : !> \param source input fm matrix
1517 : !> \param destination output fm matrix
1518 : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
1519 : !> of the input and output matrices
1520 : !> \par History
1521 : !> 31-Jan-2017 : Re-implemented using non-blocking MPI [IainB, MarkT]
1522 : ! **************************************************************************************************
1523 19008 : SUBROUTINE cp_fm_copy_general(source, destination, para_env)
1524 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1525 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1526 :
1527 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_copy_general'
1528 :
1529 : INTEGER :: handle
1530 171072 : TYPE(copy_info_type) :: info
1531 :
1532 19008 : CALL timeset(routineN, handle)
1533 :
1534 19008 : CALL cp_fm_start_copy_general(source, destination, para_env, info)
1535 19008 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1536 18994 : CALL cp_fm_finish_copy_general(destination, info)
1537 : END IF
1538 19008 : IF (ASSOCIATED(source%matrix_struct)) THEN
1539 18803 : CALL cp_fm_cleanup_copy_general(info)
1540 : END IF
1541 :
1542 19008 : CALL timestop(handle)
1543 19008 : END SUBROUTINE cp_fm_copy_general
1544 :
1545 : ! **************************************************************************************************
1546 : !> \brief Initiates the copy operation: get distribution data, post MPI isend and irecvs
1547 : !> \param source input fm matrix
1548 : !> \param destination output fm matrix
1549 : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
1550 : !> of the input and output matrices
1551 : !> \param info all of the data that will be needed to complete the copy operation
1552 : ! **************************************************************************************************
1553 7476320 : SUBROUTINE cp_fm_start_copy_general(source, destination, para_env, info)
1554 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1555 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1556 : TYPE(copy_info_type), INTENT(OUT) :: info
1557 :
1558 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_start_copy_general'
1559 :
1560 : INTEGER :: dest_p_i, dest_q_j, global_rank, global_size, handle, i, j, k, mpi_rank, &
1561 : ncol_block_dest, ncol_block_src, ncol_local_recv, ncol_local_send, ncols, &
1562 : nrow_block_dest, nrow_block_src, nrow_local_recv, nrow_local_send, nrows, p, q, &
1563 : recv_rank, recv_size, send_rank, send_size
1564 747632 : INTEGER, ALLOCATABLE, DIMENSION(:) :: all_ranks, dest2global, dest_p, dest_q, &
1565 1495264 : recv_count, send_count, send_disp, &
1566 747632 : source2global, src_p, src_q
1567 747632 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: dest_blacs2mpi
1568 : INTEGER, DIMENSION(2) :: dest_block, dest_block_tmp, dest_num_pe, &
1569 : src_block, src_block_tmp, src_num_pe
1570 1495264 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices, &
1571 1495264 : send_col_indices, send_row_indices
1572 : TYPE(cp_fm_struct_type), POINTER :: recv_dist, send_dist
1573 10466848 : TYPE(mp_request_type), DIMENSION(6) :: recv_req, send_req
1574 :
1575 747632 : CALL timeset(routineN, handle)
1576 :
1577 : IF (.NOT. cp2k_is_parallel) THEN
1578 : ! Just copy all of the matrix data into a 'send buffer', to be unpacked later
1579 : nrow_local_send = SIZE(source%local_data, 1)
1580 : ncol_local_send = SIZE(source%local_data, 2)
1581 : ALLOCATE (info%send_buf(nrow_local_send*ncol_local_send))
1582 : k = 0
1583 : DO j = 1, ncol_local_send
1584 : DO i = 1, nrow_local_send
1585 : k = k + 1
1586 : info%send_buf(k) = source%local_data(i, j)
1587 : END DO
1588 : END DO
1589 : ELSE
1590 747632 : NULLIFY (recv_dist, send_dist)
1591 747632 : NULLIFY (recv_col_indices, recv_row_indices, send_col_indices, send_row_indices)
1592 :
1593 : ! The 'global' communicator contains both the source and destination decompositions
1594 747632 : global_size = para_env%num_pe
1595 747632 : global_rank = para_env%mepos
1596 :
1597 : ! The source/send decomposition and destination/recv decompositions may only exist on
1598 : ! on a subset of the processes involved in the communication
1599 : ! Check if the source and/or destination arguments are .not. ASSOCIATED():
1600 : ! if so, skip the send / recv parts (since these processes do not participate in the sending/receiving distribution)
1601 747632 : IF (ASSOCIATED(destination%matrix_struct)) THEN
1602 566770 : recv_dist => destination%matrix_struct
1603 566770 : recv_rank = recv_dist%para_env%mepos
1604 : ELSE
1605 180862 : recv_rank = mp_proc_null
1606 : END IF
1607 :
1608 747632 : IF (ASSOCIATED(source%matrix_struct)) THEN
1609 656671 : send_dist => source%matrix_struct
1610 656671 : send_rank = send_dist%para_env%mepos
1611 : ELSE
1612 90961 : send_rank = mp_proc_null
1613 : END IF
1614 :
1615 : ! Map the rank in the source/dest communicator to the global rank
1616 2242896 : ALLOCATE (all_ranks(0:global_size - 1))
1617 :
1618 747632 : CALL para_env%allgather(send_rank, all_ranks)
1619 747632 : IF (ASSOCIATED(recv_dist)) THEN
1620 2833850 : ALLOCATE (source2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
1621 1700310 : DO i = 0, global_size - 1
1622 1700310 : IF (all_ranks(i) /= mp_proc_null) THEN
1623 951618 : source2global(all_ranks(i)) = i
1624 : END IF
1625 : END DO
1626 : END IF
1627 :
1628 747632 : CALL para_env%allgather(recv_rank, all_ranks)
1629 747632 : IF (ASSOCIATED(send_dist)) THEN
1630 3283355 : ALLOCATE (dest2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
1631 1970013 : DO i = 0, global_size - 1
1632 1970013 : IF (all_ranks(i) /= mp_proc_null) THEN
1633 951618 : dest2global(all_ranks(i)) = i
1634 : END IF
1635 : END DO
1636 : END IF
1637 747632 : DEALLOCATE (all_ranks)
1638 :
1639 : ! Some data from the two decompositions will be needed by all processes in the global group :
1640 : ! process grid shape, block size, and the BLACS-to-MPI mapping
1641 :
1642 : ! The global root process will receive the data (from the root process in each decomposition)
1643 5233424 : send_req(:) = mp_request_null
1644 747632 : IF (global_rank == 0) THEN
1645 2616712 : recv_req(:) = mp_request_null
1646 373816 : CALL para_env%irecv(src_block, mp_any_source, recv_req(1), tag=src_tag)
1647 373816 : CALL para_env%irecv(dest_block, mp_any_source, recv_req(2), tag=dest_tag)
1648 373816 : CALL para_env%irecv(src_num_pe, mp_any_source, recv_req(3), tag=src_tag)
1649 373816 : CALL para_env%irecv(dest_num_pe, mp_any_source, recv_req(4), tag=dest_tag)
1650 : END IF
1651 :
1652 747632 : IF (ASSOCIATED(send_dist)) THEN
1653 656671 : IF ((send_rank == 0)) THEN
1654 : ! need to use separate buffers here in case this is actually global rank 0
1655 1121448 : src_block_tmp = [send_dist%nrow_block, send_dist%ncol_block]
1656 373816 : CALL para_env%isend(src_block_tmp, 0, send_req(1), tag=src_tag)
1657 373816 : CALL para_env%isend(send_dist%context%num_pe, 0, send_req(2), tag=src_tag)
1658 : END IF
1659 : END IF
1660 :
1661 747632 : IF (ASSOCIATED(recv_dist)) THEN
1662 566770 : IF ((recv_rank == 0)) THEN
1663 1121448 : dest_block_tmp = [recv_dist%nrow_block, recv_dist%ncol_block]
1664 373816 : CALL para_env%isend(dest_block_tmp, 0, send_req(3), tag=dest_tag)
1665 373816 : CALL para_env%isend(recv_dist%context%num_pe, 0, send_req(4), tag=dest_tag)
1666 : END IF
1667 : END IF
1668 :
1669 747632 : IF (global_rank == 0) THEN
1670 373816 : CALL mp_waitall(recv_req(1:4))
1671 : ! Now we know the process decomposition, we can allocate the arrays to hold the blacs2mpi mapping
1672 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1673 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
1674 2616712 : )
1675 373816 : CALL para_env%irecv(info%src_blacs2mpi, mp_any_source, recv_req(5), tag=src_tag)
1676 373816 : CALL para_env%irecv(dest_blacs2mpi, mp_any_source, recv_req(6), tag=dest_tag)
1677 : END IF
1678 :
1679 747632 : IF (ASSOCIATED(send_dist)) THEN
1680 656671 : IF ((send_rank == 0)) THEN
1681 373816 : CALL para_env%isend(send_dist%context%blacs2mpi(:, :), 0, send_req(5), tag=src_tag)
1682 : END IF
1683 : END IF
1684 :
1685 747632 : IF (ASSOCIATED(recv_dist)) THEN
1686 566770 : IF ((recv_rank == 0)) THEN
1687 373816 : CALL para_env%isend(recv_dist%context%blacs2mpi(:, :), 0, send_req(6), tag=dest_tag)
1688 : END IF
1689 : END IF
1690 :
1691 747632 : IF (global_rank == 0) THEN
1692 373816 : CALL mp_waitall(recv_req(5:6))
1693 : END IF
1694 :
1695 : ! Finally, broadcast the data to all processes in the global communicator
1696 747632 : CALL para_env%bcast(src_block, 0)
1697 747632 : CALL para_env%bcast(dest_block, 0)
1698 747632 : CALL para_env%bcast(src_num_pe, 0)
1699 747632 : CALL para_env%bcast(dest_num_pe, 0)
1700 2242896 : info%src_num_pe(1:2) = src_num_pe(1:2)
1701 2242896 : info%nblock_src(1:2) = src_block(1:2)
1702 747632 : IF (global_rank /= 0) THEN
1703 0 : ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
1704 0 : dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
1705 2616712 : )
1706 : END IF
1707 747632 : CALL para_env%bcast(info%src_blacs2mpi, 0)
1708 747632 : CALL para_env%bcast(dest_blacs2mpi, 0)
1709 :
1710 747632 : recv_size = dest_num_pe(1)*dest_num_pe(2)
1711 747632 : send_size = src_num_pe(1)*src_num_pe(2)
1712 747632 : info%send_size = send_size
1713 747632 : CALL mp_waitall(send_req(:))
1714 :
1715 : ! Setup is now complete, we can start the actual communication here.
1716 : ! The order implemented here is:
1717 : ! DEST_1
1718 : ! compute recv sizes
1719 : ! call irecv
1720 : ! SRC_1
1721 : ! compute send sizes
1722 : ! pack send buffers
1723 : ! call isend
1724 : ! DEST_2
1725 : ! wait for the recvs and unpack buffers (this part eventually will go into another
1726 : ! routine to allow comms to run concurrently)
1727 : ! SRC_2
1728 : ! wait for the sends
1729 :
1730 : ! DEST_1
1731 747632 : IF (ASSOCIATED(recv_dist)) THEN
1732 : CALL cp_fm_struct_get(recv_dist, row_indices=recv_row_indices, &
1733 : col_indices=recv_col_indices &
1734 566770 : )
1735 566770 : info%recv_col_indices => recv_col_indices
1736 566770 : info%recv_row_indices => recv_row_indices
1737 566770 : nrow_block_src = src_block(1)
1738 566770 : ncol_block_src = src_block(2)
1739 3785468 : ALLOCATE (recv_count(0:send_size - 1), info%recv_disp(0:send_size - 1), info%recv_request(0:send_size - 1))
1740 :
1741 : ! Determine the recv counts, allocate the receive buffers, call mpi_irecv for all the non-zero sized receives
1742 566770 : nrow_local_recv = recv_dist%nrow_locals(recv_dist%context%mepos(1))
1743 566770 : ncol_local_recv = recv_dist%ncol_locals(recv_dist%context%mepos(2))
1744 566770 : info%nlocal_recv(1) = nrow_local_recv
1745 566770 : info%nlocal_recv(2) = ncol_local_recv
1746 : ! Initialise src_p, src_q arrays (sized using number of rows/cols in the receiving distribution)
1747 2833850 : ALLOCATE (src_p(nrow_local_recv), src_q(ncol_local_recv))
1748 10243499 : DO i = 1, nrow_local_recv
1749 : ! For each local row we will receive, we look up its global row (in recv_row_indices),
1750 : ! then work out which row block it comes from, and which process row that row block comes from.
1751 10243499 : src_p(i) = MOD(((recv_row_indices(i) - 1)/nrow_block_src), src_num_pe(1))
1752 : END DO
1753 15081274 : DO j = 1, ncol_local_recv
1754 : ! Similarly for the columns
1755 15081274 : src_q(j) = MOD(((recv_col_indices(j) - 1)/ncol_block_src), src_num_pe(2))
1756 : END DO
1757 : ! src_p/q now contains the process row/column ID that will send data to that row/column
1758 :
1759 1133540 : DO q = 0, src_num_pe(2) - 1
1760 15081274 : ncols = COUNT(src_q == q)
1761 2085158 : DO p = 0, src_num_pe(1) - 1
1762 18055250 : nrows = COUNT(src_p == p)
1763 : ! Use the send_dist here as we are looking up the processes where the data comes from
1764 1518388 : recv_count(info%src_blacs2mpi(p, q)) = nrows*ncols
1765 : END DO
1766 : END DO
1767 566770 : DEALLOCATE (src_p, src_q)
1768 :
1769 : ! Use one long buffer (and displacements into that buffer)
1770 : ! this prevents the need for a rectangular array where not all elements will be populated
1771 2651928 : ALLOCATE (info%recv_buf(SUM(recv_count(:))))
1772 566770 : info%recv_disp(0) = 0
1773 951618 : DO i = 1, send_size - 1
1774 951618 : info%recv_disp(i) = info%recv_disp(i - 1) + recv_count(i - 1)
1775 : END DO
1776 :
1777 : ! Issue receive calls on ranks which expect data
1778 1518388 : DO k = 0, send_size - 1
1779 1518388 : IF (recv_count(k) > 0) THEN
1780 : CALL para_env%irecv(info%recv_buf(info%recv_disp(k) + 1:info%recv_disp(k) + recv_count(k)), &
1781 748198 : source2global(k), info%recv_request(k))
1782 : END IF
1783 : END DO
1784 566770 : DEALLOCATE (source2global)
1785 : END IF ! ASSOCIATED(recv_dist)
1786 :
1787 : ! SRC_1
1788 747632 : IF (ASSOCIATED(send_dist)) THEN
1789 : CALL cp_fm_struct_get(send_dist, row_indices=send_row_indices, &
1790 : col_indices=send_col_indices &
1791 656671 : )
1792 656671 : nrow_block_dest = dest_block(1)
1793 656671 : ncol_block_dest = dest_block(2)
1794 4234973 : ALLOCATE (send_count(0:recv_size - 1), send_disp(0:recv_size - 1), info%send_request(0:recv_size - 1))
1795 :
1796 : ! Determine the send counts, allocate the send buffers
1797 656671 : nrow_local_send = send_dist%nrow_locals(send_dist%context%mepos(1))
1798 656671 : ncol_local_send = send_dist%ncol_locals(send_dist%context%mepos(2))
1799 :
1800 : ! Initialise dest_p, dest_q arrays (sized nrow_local, ncol_local)
1801 : ! i.e. number of rows,cols in the sending distribution
1802 3283355 : ALLOCATE (dest_p(nrow_local_send), dest_q(ncol_local_send))
1803 :
1804 10333400 : DO i = 1, nrow_local_send
1805 : ! Use the send_dist%row_indices() here (we are looping over the local rows we will send)
1806 10333400 : dest_p(i) = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1807 : END DO
1808 17426607 : DO j = 1, ncol_local_send
1809 17426607 : dest_q(j) = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1810 : END DO
1811 : ! dest_p/q now contain the process row/column ID that will receive data from that row/column
1812 :
1813 1313342 : DO q = 0, dest_num_pe(2) - 1
1814 17426607 : ncols = COUNT(dest_q == q)
1815 2264960 : DO p = 0, dest_num_pe(1) - 1
1816 15800320 : nrows = COUNT(dest_p == p)
1817 1608289 : send_count(dest_blacs2mpi(p, q)) = nrows*ncols
1818 : END DO
1819 : END DO
1820 656671 : DEALLOCATE (dest_p, dest_q)
1821 :
1822 : ! Allocate the send buffer using send_count -- and calculate the offset into the buffer for each process
1823 2921631 : ALLOCATE (info%send_buf(SUM(send_count(:))))
1824 656671 : send_disp(0) = 0
1825 951618 : DO k = 1, recv_size - 1
1826 951618 : send_disp(k) = send_disp(k - 1) + send_count(k - 1)
1827 : END DO
1828 :
1829 : ! Loop over the smat, pack the send buffers
1830 656671 : send_count(:) = 0
1831 17426607 : DO j = 1, ncol_local_send
1832 : ! Use send_col_indices and row_indices here, as we are looking up the global row/column number of local rows.
1833 16769936 : dest_q_j = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
1834 529383379 : DO i = 1, nrow_local_send
1835 511956772 : dest_p_i = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
1836 511956772 : mpi_rank = dest_blacs2mpi(dest_p_i, dest_q_j)
1837 511956772 : send_count(mpi_rank) = send_count(mpi_rank) + 1
1838 528726708 : info%send_buf(send_disp(mpi_rank) + send_count(mpi_rank)) = source%local_data(i, j)
1839 : END DO
1840 : END DO
1841 :
1842 : ! For each non-zero send_count, call mpi_isend
1843 1608289 : DO k = 0, recv_size - 1
1844 1608289 : IF (send_count(k) > 0) THEN
1845 : CALL para_env%isend(info%send_buf(send_disp(k) + 1:send_disp(k) + send_count(k)), &
1846 748198 : dest2global(k), info%send_request(k))
1847 : END IF
1848 : END DO
1849 656671 : DEALLOCATE (send_count, send_disp, dest2global)
1850 : END IF ! ASSOCIATED(send_dist)
1851 747632 : DEALLOCATE (dest_blacs2mpi)
1852 :
1853 : END IF !IF (.NOT. cp2k_is_parallel)
1854 :
1855 747632 : CALL timestop(handle)
1856 :
1857 2990528 : END SUBROUTINE cp_fm_start_copy_general
1858 :
1859 : ! **************************************************************************************************
1860 : !> \brief Completes the copy operation: wait for comms, unpack, clean up MPI state
1861 : !> \param destination output fm matrix
1862 : !> \param info all of the data that will be needed to complete the copy operation
1863 : ! **************************************************************************************************
1864 566770 : SUBROUTINE cp_fm_finish_copy_general(destination, info)
1865 : TYPE(cp_fm_type), INTENT(IN) :: destination
1866 : TYPE(copy_info_type), INTENT(INOUT) :: info
1867 :
1868 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_finish_copy_general'
1869 :
1870 : INTEGER :: handle, i, j, k, mpi_rank, send_size, &
1871 : src_p_i, src_q_j
1872 566770 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_count
1873 : INTEGER, DIMENSION(2) :: nblock_src, nlocal_recv, src_num_pe
1874 566770 : INTEGER, DIMENSION(:), POINTER :: recv_col_indices, recv_row_indices
1875 :
1876 566770 : CALL timeset(routineN, handle)
1877 :
1878 : IF (.NOT. cp2k_is_parallel) THEN
1879 : ! Now unpack the data from the 'send buffer'
1880 : k = 0
1881 : DO j = 1, SIZE(destination%local_data, 2)
1882 : DO i = 1, SIZE(destination%local_data, 1)
1883 : k = k + 1
1884 : destination%local_data(i, j) = info%send_buf(k)
1885 : END DO
1886 : END DO
1887 : DEALLOCATE (info%send_buf)
1888 : ELSE
1889 : ! Set up local variables ...
1890 566770 : send_size = info%send_size
1891 1700310 : nlocal_recv(1:2) = info%nlocal_recv(:)
1892 1700310 : nblock_src(1:2) = info%nblock_src(:)
1893 1700310 : src_num_pe(1:2) = info%src_num_pe(:)
1894 566770 : recv_col_indices => info%recv_col_indices
1895 566770 : recv_row_indices => info%recv_row_indices
1896 :
1897 : ! ... use the local variables to do the work
1898 : ! DEST_2
1899 566770 : CALL mp_waitall(info%recv_request(:))
1900 1700310 : ALLOCATE (recv_count(0:send_size - 1))
1901 : ! Loop over the rmat, filling it in with data from the recv buffers
1902 : ! (here the block sizes, num_pes refer to the distribution of the source matrix)
1903 566770 : recv_count(:) = 0
1904 15081274 : DO j = 1, nlocal_recv(2)
1905 14514504 : src_q_j = MOD(((recv_col_indices(j) - 1)/nblock_src(2)), src_num_pe(2))
1906 527038046 : DO i = 1, nlocal_recv(1)
1907 511956772 : src_p_i = MOD(((recv_row_indices(i) - 1)/nblock_src(1)), src_num_pe(1))
1908 511956772 : mpi_rank = info%src_blacs2mpi(src_p_i, src_q_j)
1909 511956772 : recv_count(mpi_rank) = recv_count(mpi_rank) + 1
1910 526471276 : destination%local_data(i, j) = info%recv_buf(info%recv_disp(mpi_rank) + recv_count(mpi_rank))
1911 : END DO
1912 : END DO
1913 566770 : DEALLOCATE (recv_count, info%recv_disp, info%recv_request, info%recv_buf, info%src_blacs2mpi)
1914 : ! Invalidate the stored state
1915 : NULLIFY (info%recv_col_indices, &
1916 566770 : info%recv_row_indices)
1917 :
1918 : END IF
1919 :
1920 566770 : CALL timestop(handle)
1921 :
1922 566770 : END SUBROUTINE cp_fm_finish_copy_general
1923 :
1924 : ! **************************************************************************************************
1925 : !> \brief Completes the copy operation: wait for comms clean up MPI state
1926 : !> \param info all of the data that will be needed to complete the copy operation
1927 : ! **************************************************************************************************
1928 655975 : SUBROUTINE cp_fm_cleanup_copy_general(info)
1929 : TYPE(copy_info_type), INTENT(INOUT) :: info
1930 :
1931 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_cleanup_copy_general'
1932 :
1933 : INTEGER :: handle
1934 :
1935 655975 : CALL timeset(routineN, handle)
1936 :
1937 : IF (.NOT. cp2k_is_parallel) THEN
1938 : ! Don't do anything - no MPI state for the serial case
1939 : ELSE
1940 : ! SRC_2
1941 : ! If this process is also in the destination decomposition, this deallocate
1942 : ! Was already done in cp_fm_finish_copy_general
1943 655975 : IF (ALLOCATED(info%src_blacs2mpi)) THEN
1944 180166 : DEALLOCATE (info%src_blacs2mpi)
1945 : END IF
1946 655975 : CALL mp_waitall(info%send_request)
1947 655975 : DEALLOCATE (info%send_request, info%send_buf)
1948 :
1949 : END IF
1950 :
1951 655975 : CALL timestop(handle)
1952 :
1953 655975 : END SUBROUTINE cp_fm_cleanup_copy_general
1954 :
1955 : ! **************************************************************************************************
1956 : !> \brief General copy of a submatrix of fm matrix to a submatrix of another fm matrix.
1957 : !> The two matrices can have different contexts.
1958 : !>
1959 : !> Summary of distribution routines for dense matrices
1960 : !> The following will copy A(iA:iA+M-1,jA:jA+N-1) to B(iB:iB+M-1,jB:jB+N-1):
1961 : !>
1962 : !> call pdgemr2d(M,N,Aloc,iA,jA,descA,Bloc,iB,jB,descB,context)
1963 : !>
1964 : !> A process that is not a part of the context of A should set descA(2)
1965 : !> to -1, and similarly for B.
1966 : !>
1967 : !> \param source input fm matrix
1968 : !> \param destination output fm matrix
1969 : !> \param nrows number of rows of sub matrix to be copied
1970 : !> \param ncols number of cols of sub matrix to be copied
1971 : !> \param s_firstrow starting global row index of sub matrix in source
1972 : !> \param s_firstcol starting global col index of sub matrix in source
1973 : !> \param d_firstrow starting global row index of sub matrix in destination
1974 : !> \param d_firstcol starting global col index of sub matrix in destination
1975 : !> \param global_context process grid that covers all parts of either A or B.
1976 : ! **************************************************************************************************
1977 11579 : SUBROUTINE cp_fm_to_fm_submat_general(source, &
1978 : destination, &
1979 : nrows, &
1980 : ncols, &
1981 : s_firstrow, &
1982 : s_firstcol, &
1983 : d_firstrow, &
1984 : d_firstcol, &
1985 : global_context)
1986 :
1987 : TYPE(cp_fm_type), INTENT(IN) :: source, destination
1988 : INTEGER, INTENT(IN) :: nrows, ncols, s_firstrow, s_firstcol, &
1989 : d_firstrow, d_firstcol
1990 :
1991 : CLASS(cp_blacs_type), INTENT(IN) :: global_context
1992 :
1993 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat_general'
1994 :
1995 : LOGICAL :: debug
1996 : INTEGER :: handle
1997 : #if defined(__parallel)
1998 : INTEGER, DIMENSION(9) :: desca, descb
1999 : REAL(KIND=dp), DIMENSION(1, 1), TARGET :: dummy
2000 11579 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: smat, dmat
2001 : #endif
2002 :
2003 11579 : CALL timeset(routineN, handle)
2004 :
2005 11579 : debug = debug_this_module
2006 :
2007 : IF (.NOT. cp2k_is_parallel) THEN
2008 : CALL cp_fm_to_fm_submat(source, &
2009 : destination, &
2010 : nrows, &
2011 : ncols, &
2012 : s_firstrow, &
2013 : s_firstcol, &
2014 : d_firstrow, &
2015 : d_firstcol)
2016 : ELSE
2017 : #ifdef __parallel
2018 : NULLIFY (smat, dmat)
2019 : ! check whether source is available on this process
2020 11579 : IF (ASSOCIATED(source%matrix_struct)) THEN
2021 115790 : desca = source%matrix_struct%descriptor
2022 11579 : IF (nrows > source%matrix_struct%nrow_global) THEN
2023 0 : CPABORT("nrows is greater than nrow_global of source")
2024 : END IF
2025 11579 : IF (ncols > source%matrix_struct%ncol_global) THEN
2026 0 : CPABORT("ncols is greater than ncol_global of source")
2027 : END IF
2028 11579 : smat => source%local_data
2029 : ELSE
2030 0 : desca = -1
2031 0 : smat => dummy
2032 : END IF
2033 : ! check destination is available on this process
2034 11579 : IF (ASSOCIATED(destination%matrix_struct)) THEN
2035 115790 : descb = destination%matrix_struct%descriptor
2036 11579 : IF (nrows > destination%matrix_struct%nrow_global) THEN
2037 0 : CPABORT("nrows is greater than nrow_global of destination")
2038 : END IF
2039 11579 : IF (ncols > destination%matrix_struct%ncol_global) THEN
2040 0 : CPABORT("ncols is greater than ncol_global of destination")
2041 : END IF
2042 11579 : dmat => destination%local_data
2043 : ELSE
2044 0 : descb = -1
2045 0 : dmat => dummy
2046 : END IF
2047 : ! do copy
2048 :
2049 : CALL pdgemr2d(nrows, &
2050 : ncols, &
2051 : smat, &
2052 : s_firstrow, &
2053 : s_firstcol, &
2054 : desca, &
2055 : dmat, &
2056 : d_firstrow, &
2057 : d_firstcol, &
2058 : descb, &
2059 11579 : global_context%get_handle())
2060 : #else
2061 : MARK_USED(global_context)
2062 : CPABORT("this subroutine only supports SCALAPACK")
2063 : #endif
2064 : END IF
2065 :
2066 11579 : CALL timestop(handle)
2067 :
2068 11579 : END SUBROUTINE cp_fm_to_fm_submat_general
2069 :
2070 : ! **************************************************************************************************
2071 : !> \brief ...
2072 : !> \param matrix ...
2073 : !> \param irow_global ...
2074 : !> \param icol_global ...
2075 : !> \param alpha ...
2076 : ! **************************************************************************************************
2077 240 : SUBROUTINE cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
2078 :
2079 : ! Add alpha to the matrix element specified by the global indices
2080 : ! irow_global and icol_global
2081 :
2082 : ! - Creation (05.05.06,MK)
2083 :
2084 : TYPE(cp_fm_type), INTENT(IN) :: matrix
2085 : INTEGER, INTENT(IN) :: irow_global, icol_global
2086 : REAL(KIND=dp), INTENT(IN) :: alpha
2087 :
2088 : INTEGER :: mypcol, myprow, npcol, nprow
2089 240 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
2090 : TYPE(cp_blacs_env_type), POINTER :: context
2091 : #if defined(__parallel)
2092 : INTEGER :: icol_local, ipcol, iprow, &
2093 : irow_local
2094 : INTEGER, DIMENSION(9) :: desca
2095 : #endif
2096 :
2097 240 : context => matrix%matrix_struct%context
2098 :
2099 240 : myprow = context%mepos(1)
2100 240 : mypcol = context%mepos(2)
2101 :
2102 240 : nprow = context%num_pe(1)
2103 240 : npcol = context%num_pe(2)
2104 :
2105 240 : a => matrix%local_data
2106 :
2107 : #if defined(__parallel)
2108 :
2109 2400 : desca(:) = matrix%matrix_struct%descriptor(:)
2110 :
2111 : CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
2112 240 : irow_local, icol_local, iprow, ipcol)
2113 :
2114 240 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
2115 120 : a(irow_local, icol_local) = a(irow_local, icol_local) + alpha
2116 : END IF
2117 :
2118 : #else
2119 :
2120 : a(irow_global, icol_global) = a(irow_global, icol_global) + alpha
2121 :
2122 : #endif
2123 :
2124 240 : END SUBROUTINE cp_fm_add_to_element
2125 :
2126 : ! **************************************************************************************************
2127 : !> \brief ...
2128 : !> \param fm ...
2129 : !> \param unit ...
2130 : ! **************************************************************************************************
2131 82496 : SUBROUTINE cp_fm_write_unformatted(fm, unit)
2132 : TYPE(cp_fm_type), INTENT(IN) :: fm
2133 : INTEGER, INTENT(IN) :: unit
2134 :
2135 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_unformatted'
2136 :
2137 : INTEGER :: handle, j, max_block, &
2138 : ncol_global, nrow_global
2139 : TYPE(mp_para_env_type), POINTER :: para_env
2140 : #if defined(__parallel)
2141 : INTEGER :: i, i_block, icol_local, &
2142 : in, info, ipcol, &
2143 : iprow, irow_local, &
2144 : mepos, &
2145 : num_pe, rb, tag
2146 : INTEGER, DIMENSION(9) :: desc
2147 82496 : REAL(KIND=dp), DIMENSION(:), POINTER :: vecbuf
2148 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: newdat
2149 : TYPE(cp_blacs_type) :: ictxt_loc
2150 : INTEGER, EXTERNAL :: numroc
2151 : #endif
2152 :
2153 82496 : CALL timeset(routineN, handle)
2154 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2155 82496 : para_env=para_env)
2156 :
2157 : #if defined(__parallel)
2158 82496 : num_pe = para_env%num_pe
2159 82496 : mepos = para_env%mepos
2160 82496 : rb = nrow_global
2161 82496 : tag = 0
2162 : ! get a new context
2163 82496 : CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
2164 82496 : CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
2165 82496 : CPASSERT(info == 0)
2166 : ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
2167 : myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
2168 164992 : in = numroc(ncol_global, max_block, mypcol, 0, npcol)
2169 :
2170 329984 : ALLOCATE (newdat(nrow_global, MAX(1, in)))
2171 :
2172 : ! do the actual scalapack to cols reordering
2173 : CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
2174 : fm%matrix_struct%descriptor, &
2175 82496 : newdat, 1, 1, desc, ictxt_loc%get_handle())
2176 :
2177 247488 : ALLOCATE (vecbuf(nrow_global*max_block))
2178 30451326 : vecbuf = HUGE(1.0_dp) ! init for valgrind
2179 :
2180 359851 : DO i = 1, ncol_global, MAX(max_block, 1)
2181 194859 : i_block = MIN(max_block, ncol_global - i + 1)
2182 : CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
2183 194859 : irow_local, icol_local, iprow, ipcol)
2184 194859 : IF (ipcol == mypcol) THEN
2185 716301 : DO j = 1, i_block
2186 77838719 : vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
2187 : END DO
2188 : END IF
2189 :
2190 194859 : IF (ipcol == 0) THEN
2191 : ! do nothing
2192 : ELSE
2193 66134 : IF (ipcol == mypcol) THEN
2194 17111747 : CALL para_env%send(vecbuf(:), 0, tag)
2195 : END IF
2196 66134 : IF (mypcol == 0) THEN
2197 34190427 : CALL para_env%recv(vecbuf(:), ipcol, tag)
2198 : END IF
2199 : END IF
2200 :
2201 277355 : IF (unit > 0) THEN
2202 713589 : DO j = 1, i_block
2203 39247782 : WRITE (unit) vecbuf((j - 1)*nrow_global + 1:nrow_global*j)
2204 : END DO
2205 : END IF
2206 :
2207 : END DO
2208 : END ASSOCIATE
2209 82496 : DEALLOCATE (vecbuf)
2210 :
2211 82496 : CALL ictxt_loc%gridexit()
2212 :
2213 82496 : DEALLOCATE (newdat)
2214 :
2215 : #else
2216 :
2217 : IF (unit > 0) THEN
2218 : DO j = 1, ncol_global
2219 : WRITE (unit) fm%local_data(:, j)
2220 : END DO
2221 : END IF
2222 :
2223 : #endif
2224 82496 : CALL timestop(handle)
2225 :
2226 412480 : END SUBROUTINE cp_fm_write_unformatted
2227 :
2228 : ! **************************************************************************************************
2229 : !> \brief Write out a full matrix in plain text.
2230 : !> \param fm the matrix to be outputted
2231 : !> \param unit the unit number for I/O
2232 : !> \param header optional header
2233 : !> \param value_format ...
2234 : ! **************************************************************************************************
2235 1378 : SUBROUTINE cp_fm_write_formatted(fm, unit, header, value_format)
2236 : TYPE(cp_fm_type), INTENT(IN) :: fm
2237 : INTEGER, INTENT(IN) :: unit
2238 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: header, value_format
2239 :
2240 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_formatted'
2241 :
2242 : CHARACTER(LEN=21) :: my_value_format
2243 : INTEGER :: handle, i, j, max_block, &
2244 : ncol_global, nrow_global
2245 : TYPE(mp_para_env_type), POINTER :: para_env
2246 : #if defined(__parallel)
2247 : INTEGER :: i_block, icol_local, &
2248 : in, info, ipcol, &
2249 : iprow, irow_local, &
2250 : mepos, num_pe, rb, tag, k, &
2251 : icol, irow
2252 : INTEGER, DIMENSION(9) :: desc
2253 1378 : REAL(KIND=dp), DIMENSION(:), POINTER :: vecbuf
2254 1378 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: newdat
2255 : TYPE(cp_blacs_type) :: ictxt_loc
2256 : INTEGER, EXTERNAL :: numroc
2257 : #endif
2258 :
2259 1378 : CALL timeset(routineN, handle)
2260 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2261 1378 : para_env=para_env)
2262 :
2263 1378 : IF (PRESENT(value_format)) THEN
2264 0 : CPASSERT(LEN_TRIM(ADJUSTL(value_format)) < 11)
2265 0 : my_value_format = "(I10, I10, "//TRIM(ADJUSTL(value_format))//")"
2266 : ELSE
2267 1378 : my_value_format = "(I10, I10, ES24.12)"
2268 : END IF
2269 :
2270 1378 : IF (unit > 0) THEN
2271 11 : IF (PRESENT(header)) WRITE (unit, *) header
2272 11 : WRITE (unit, "(A2, A8, A10, A24)") "#", "Row", "Column", ADJUSTL("Value")
2273 : END IF
2274 :
2275 : #if defined(__parallel)
2276 1378 : num_pe = para_env%num_pe
2277 1378 : mepos = para_env%mepos
2278 1378 : rb = nrow_global
2279 1378 : tag = 0
2280 : ! get a new context
2281 1378 : CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
2282 1378 : CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
2283 1378 : CPASSERT(info == 0)
2284 : ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
2285 : myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
2286 2756 : in = numroc(ncol_global, max_block, mypcol, 0, npcol)
2287 :
2288 5512 : ALLOCATE (newdat(nrow_global, MAX(1, in)))
2289 :
2290 : ! do the actual scalapack to cols reordering
2291 : CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
2292 : fm%matrix_struct%descriptor, &
2293 1378 : newdat, 1, 1, desc, ictxt_loc%get_handle())
2294 :
2295 4134 : ALLOCATE (vecbuf(nrow_global*max_block))
2296 14136 : vecbuf = HUGE(1.0_dp) ! init for valgrind
2297 1378 : irow = 1
2298 1378 : icol = 1
2299 :
2300 5510 : DO i = 1, ncol_global, MAX(max_block, 1)
2301 2754 : i_block = MIN(max_block, ncol_global - i + 1)
2302 : CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
2303 2754 : irow_local, icol_local, iprow, ipcol)
2304 2754 : IF (ipcol == mypcol) THEN
2305 2922 : DO j = 1, i_block
2306 20802 : vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
2307 : END DO
2308 : END IF
2309 :
2310 2754 : IF (ipcol == 0) THEN
2311 : ! do nothing
2312 : ELSE
2313 1370 : IF (ipcol == mypcol) THEN
2314 5019 : CALL para_env%send(vecbuf(:), 0, tag)
2315 : END IF
2316 1370 : IF (mypcol == 0) THEN
2317 9353 : CALL para_env%recv(vecbuf(:), ipcol, tag)
2318 : END IF
2319 : END IF
2320 :
2321 4132 : IF (unit > 0) THEN
2322 210 : DO j = 1, i_block
2323 6438 : DO k = (j - 1)*nrow_global + 1, nrow_global*j
2324 6228 : WRITE (UNIT=unit, FMT=my_value_format) irow, icol, vecbuf(k)
2325 6228 : irow = irow + 1
2326 6417 : IF (irow > nrow_global) THEN
2327 189 : irow = 1
2328 189 : icol = icol + 1
2329 : END IF
2330 : END DO
2331 : END DO
2332 : END IF
2333 :
2334 : END DO
2335 : END ASSOCIATE
2336 1378 : DEALLOCATE (vecbuf)
2337 :
2338 1378 : CALL ictxt_loc%gridexit()
2339 :
2340 1378 : DEALLOCATE (newdat)
2341 :
2342 : #else
2343 :
2344 : IF (unit > 0) THEN
2345 : DO j = 1, ncol_global
2346 : DO i = 1, nrow_global
2347 : WRITE (UNIT=unit, FMT=my_value_format) i, j, fm%local_data(i, j)
2348 : END DO
2349 : END DO
2350 : END IF
2351 :
2352 : #endif
2353 1378 : CALL timestop(handle)
2354 :
2355 6890 : END SUBROUTINE cp_fm_write_formatted
2356 :
2357 : ! **************************************************************************************************
2358 : !> \brief ...
2359 : !> \param fm ...
2360 : !> \param unit ...
2361 : ! **************************************************************************************************
2362 1974 : SUBROUTINE cp_fm_read_unformatted(fm, unit)
2363 : TYPE(cp_fm_type), INTENT(INOUT) :: fm
2364 : INTEGER, INTENT(IN) :: unit
2365 :
2366 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_read_unformatted'
2367 :
2368 : INTEGER :: handle, j, max_block, &
2369 : ncol_global, nrow_global
2370 : TYPE(mp_para_env_type), POINTER :: para_env
2371 : #if defined(__parallel)
2372 : INTEGER :: k, n_cols
2373 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuf
2374 : #endif
2375 :
2376 1974 : CALL timeset(routineN, handle)
2377 :
2378 : CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
2379 1974 : para_env=para_env)
2380 :
2381 : #if defined(__parallel)
2382 :
2383 : ! the parallel case could be made more efficient (see cp_fm_write_unformatted)
2384 :
2385 7896 : ALLOCATE (vecbuf(nrow_global, max_block))
2386 :
2387 6152 : DO j = 1, ncol_global, max_block
2388 :
2389 4178 : n_cols = MIN(max_block, ncol_global - j + 1)
2390 4178 : IF (para_env%mepos == 0) THEN
2391 16918 : DO k = 1, n_cols
2392 368077 : READ (unit) vecbuf(:, k)
2393 : END DO
2394 : END IF
2395 1482498 : CALL para_env%bcast(vecbuf, 0)
2396 6152 : CALL cp_fm_set_submatrix(fm, vecbuf, start_row=1, start_col=j, n_cols=n_cols)
2397 :
2398 : END DO
2399 :
2400 1974 : DEALLOCATE (vecbuf)
2401 :
2402 : #else
2403 :
2404 : DO j = 1, ncol_global
2405 : READ (unit) fm%local_data(:, j)
2406 : END DO
2407 :
2408 : #endif
2409 :
2410 1974 : CALL timestop(handle)
2411 :
2412 1974 : END SUBROUTINE cp_fm_read_unformatted
2413 :
2414 : ! **************************************************************************************************
2415 : !> \brief ...
2416 : !> \param mm_type ...
2417 : ! **************************************************************************************************
2418 11087 : SUBROUTINE cp_fm_setup(mm_type)
2419 : INTEGER, INTENT(IN) :: mm_type
2420 :
2421 11087 : cp_fm_mm_type = mm_type
2422 11087 : END SUBROUTINE cp_fm_setup
2423 :
2424 : ! **************************************************************************************************
2425 : !> \brief ...
2426 : !> \return ...
2427 : ! **************************************************************************************************
2428 2284050 : FUNCTION cp_fm_get_mm_type() RESULT(res)
2429 : INTEGER :: res
2430 :
2431 2284050 : res = cp_fm_mm_type
2432 2284050 : END FUNCTION cp_fm_get_mm_type
2433 :
2434 : ! **************************************************************************************************
2435 : !> \brief ...
2436 : !> \param ictxt ...
2437 : !> \param prec ...
2438 : !> \return ...
2439 : ! **************************************************************************************************
2440 10 : FUNCTION cp_fm_pilaenv(ictxt, prec) RESULT(res)
2441 : INTEGER :: ictxt
2442 : CHARACTER(LEN=1) :: prec
2443 : INTEGER :: res
2444 : #if defined(__parallel)
2445 : INTEGER :: pilaenv
2446 10 : res = pilaenv(ictxt, prec)
2447 : #else
2448 : MARK_USED(ictxt)
2449 : MARK_USED(prec)
2450 : res = -1
2451 : #endif
2452 :
2453 10 : END FUNCTION cp_fm_pilaenv
2454 :
2455 0 : END MODULE cp_fm_types
|
