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