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