(git:98357aa)
Loading...
Searching...
No Matches
domain_submatrix_methods.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 Subroutines to handle submatrices
10!> \par History
11!> 2013.01 created [Rustam Z Khaliullin]
12!> \author Rustam Z Khaliullin
13! **************************************************************************************************
14MODULE domain_submatrix_methods
15
16 USE cp_dbcsr_api, ONLY: &
17 dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
18 dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_get_stored_coordinates, &
19 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
20 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_put_block, dbcsr_type, &
21 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
22 USE cp_log_handling, ONLY: cp_get_default_logger,&
23 cp_logger_get_default_unit_nr,&
24 cp_logger_type
25 USE domain_submatrix_types, ONLY: domain_map_type,&
26 domain_submatrix_type,&
27 select_row_col
28 USE kinds, ONLY: dp
29 USE message_passing, ONLY: mp_comm_null,&
30 mp_comm_type
31#include "./base/base_uses.f90"
32
33 IMPLICIT NONE
34
35 PRIVATE
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'domain_submatrix_methods'
38
39 PUBLIC :: copy_submatrices, copy_submatrix_data, &
40 release_submatrices, multiply_submatrices, add_submatrices, &
41 construct_submatrices, init_submatrices, &
42 construct_dbcsr_from_submatrices, &
43 set_submatrices, &
44 print_submatrices, maxnorm_submatrices
45
46 INTERFACE init_submatrices
47 MODULE PROCEDURE init_submatrices_0d
48 MODULE PROCEDURE init_submatrices_1d
49 MODULE PROCEDURE init_submatrices_2d
50 END INTERFACE
51
52 INTERFACE set_submatrices
53 MODULE PROCEDURE set_submatrix_array
54 MODULE PROCEDURE set_submatrix
55 END INTERFACE
56
57 INTERFACE copy_submatrices
58 MODULE PROCEDURE copy_submatrix_array
59 MODULE PROCEDURE copy_submatrix
60 END INTERFACE
61
62 INTERFACE release_submatrices
63 MODULE PROCEDURE release_submatrix_array
64 MODULE PROCEDURE release_submatrix
65 END INTERFACE
66
67 INTERFACE multiply_submatrices
68 MODULE PROCEDURE multiply_submatrices_once
69 MODULE PROCEDURE multiply_submatrices_array
70 END INTERFACE
71
72 INTERFACE add_submatrices
73 MODULE PROCEDURE add_submatrices_once
74 MODULE PROCEDURE add_submatrices_array
75 END INTERFACE
76
77CONTAINS
78
79! **************************************************************************************************
80!> \brief ...
81!> \param subm ...
82! **************************************************************************************************
83 SUBROUTINE init_submatrices_0d(subm)
84
85 TYPE(domain_submatrix_type), INTENT(INOUT) :: subm
86
87 subm%domain = -1
88 subm%nbrows = -1
89 subm%nbcols = -1
90 subm%nrows = -1
91 subm%ncols = -1
92 subm%nnodes = -1
93 subm%group = mp_comm_null
94
95 END SUBROUTINE init_submatrices_0d
96
97! **************************************************************************************************
98!> \brief ...
99!> \param subm ...
100! **************************************************************************************************
101 SUBROUTINE init_submatrices_1d(subm)
102
103 TYPE(domain_submatrix_type), DIMENSION(:), &
104 INTENT(INOUT) :: subm
105
106 subm(:)%domain = -1
107 subm(:)%nbrows = -1
108 subm(:)%nbcols = -1
109 subm(:)%nrows = -1
110 subm(:)%ncols = -1
111 subm(:)%nnodes = -1
112 subm(:)%group = mp_comm_null
113
114 END SUBROUTINE init_submatrices_1d
115
116! **************************************************************************************************
117!> \brief ...
118!> \param subm ...
119! **************************************************************************************************
120 SUBROUTINE init_submatrices_2d(subm)
121
122 TYPE(domain_submatrix_type), DIMENSION(:, :), &
123 INTENT(INOUT) :: subm
124
125 subm(:, :)%domain = -1
126 subm(:, :)%nbrows = -1
127 subm(:, :)%nbcols = -1
128 subm(:, :)%nrows = -1
129 subm(:, :)%ncols = -1
130 subm(:, :)%nnodes = -1
131 subm(:, :)%group = mp_comm_null
132
133 END SUBROUTINE init_submatrices_2d
134
135! **************************************************************************************************
136!> \brief ...
137!> \param original ...
138!> \param copy ...
139!> \param copy_data ...
140! **************************************************************************************************
141 SUBROUTINE copy_submatrix_array(original, copy, copy_data)
142
143 TYPE(domain_submatrix_type), DIMENSION(:), &
144 INTENT(IN) :: original
145 TYPE(domain_submatrix_type), DIMENSION(:), &
146 INTENT(INOUT) :: copy
147 LOGICAL, INTENT(IN) :: copy_data
148
149 CHARACTER(len=*), PARAMETER :: routineN = 'copy_submatrix_array'
150
151 INTEGER :: handle, idomain, ndomains, ndomainsB
152
153 CALL timeset(routinen, handle)
154
155 ndomains = SIZE(original)
156 ndomainsb = SIZE(copy)
157 cpassert(ndomains == ndomainsb)
158 copy(:)%nnodes = original(:)%nnodes
159 copy(:)%group = original(:)%group
160 DO idomain = 1, ndomains
161 IF (original(idomain)%domain > 0) THEN
162 CALL copy_submatrix(original(idomain), copy(idomain), copy_data)
163 END IF
164 END DO ! loop over domains
165
166 CALL timestop(handle)
167
168 END SUBROUTINE copy_submatrix_array
169
170! **************************************************************************************************
171!> \brief ...
172!> \param original ...
173!> \param copy ...
174!> \param copy_data ...
175! **************************************************************************************************
176 SUBROUTINE copy_submatrix(original, copy, copy_data)
177
178 TYPE(domain_submatrix_type), INTENT(IN) :: original
179 TYPE(domain_submatrix_type), INTENT(INOUT) :: copy
180 LOGICAL, INTENT(IN) :: copy_data
181
182 CHARACTER(len=*), PARAMETER :: routineN = 'copy_submatrix'
183
184 INTEGER :: handle, icol, irow
185
186 CALL timeset(routinen, handle)
187
188 copy%domain = original%domain
189 copy%nnodes = original%nnodes
190 copy%group = original%group
191
192 IF (original%domain > 0) THEN
193
194 copy%nbrows = original%nbrows
195 copy%nbcols = original%nbcols
196 copy%nrows = original%nrows
197 copy%ncols = original%ncols
198
199 IF (.NOT. ALLOCATED(copy%dbcsr_row)) THEN
200 ALLOCATE (copy%dbcsr_row(original%nbrows))
201 ELSE
202 IF (SIZE(copy%dbcsr_row) /= SIZE(original%dbcsr_row)) THEN
203 DEALLOCATE (copy%dbcsr_row)
204 ALLOCATE (copy%dbcsr_row(original%nbrows))
205 END IF
206 END IF
207 IF (.NOT. ALLOCATED(copy%dbcsr_col)) THEN
208 ALLOCATE (copy%dbcsr_col(original%nbcols))
209 ELSE
210 IF (SIZE(copy%dbcsr_col) /= SIZE(original%dbcsr_col)) THEN
211 DEALLOCATE (copy%dbcsr_col)
212 ALLOCATE (copy%dbcsr_col(original%nbcols))
213 END IF
214 END IF
215 IF (.NOT. ALLOCATED(copy%size_brow)) THEN
216 ALLOCATE (copy%size_brow(original%nbrows))
217 ELSE
218 IF (SIZE(copy%size_brow) /= SIZE(original%size_brow)) THEN
219 DEALLOCATE (copy%size_brow)
220 ALLOCATE (copy%size_brow(original%nbrows))
221 END IF
222 END IF
223 IF (.NOT. ALLOCATED(copy%size_bcol)) THEN
224 ALLOCATE (copy%size_bcol(original%nbcols))
225 ELSE
226 IF (SIZE(copy%size_bcol) /= SIZE(original%size_bcol)) THEN
227 DEALLOCATE (copy%size_bcol)
228 ALLOCATE (copy%size_bcol(original%nbcols))
229 END IF
230 END IF
231
232 DO irow = 1, original%nbrows
233 copy%dbcsr_row(irow) = original%dbcsr_row(irow)
234 copy%size_brow(irow) = original%size_brow(irow)
235 END DO
236
237 DO icol = 1, original%nbcols
238 copy%dbcsr_col(icol) = original%dbcsr_col(icol)
239 copy%size_bcol(icol) = original%size_bcol(icol)
240 END DO
241
242 IF (copy_data) THEN
243 CALL copy_submatrix_data(original%mdata, copy)
244 END IF
245
246 END IF ! do not copy empty submatrix
247
248 CALL timestop(handle)
249
250 END SUBROUTINE copy_submatrix
251
252! **************************************************************************************************
253!> \brief ...
254!> \param array ...
255!> \param copy ...
256! **************************************************************************************************
257 SUBROUTINE copy_submatrix_data(array, copy)
258
259 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: array
260 TYPE(domain_submatrix_type), INTENT(INOUT) :: copy
261
262 CHARACTER(len=*), PARAMETER :: routinen = 'copy_submatrix_data'
263
264 INTEGER :: ds1, ds2, handle, ms1, ms2
265
266 CALL timeset(routinen, handle)
267
268 cpassert(copy%domain > 0)
269
270 ds1 = SIZE(array, 1)
271 ds2 = SIZE(array, 2)
272
273 IF (.NOT. ALLOCATED(copy%mdata)) THEN
274 ALLOCATE (copy%mdata(ds1, ds2))
275 ELSE
276 ms1 = SIZE(copy%mdata, 1)
277 ms2 = SIZE(copy%mdata, 2)
278 IF ((ds1 /= ms1) .OR. (ds2 /= ms2)) THEN
279 DEALLOCATE (copy%mdata)
280 ALLOCATE (copy%mdata(ds1, ds2))
281 END IF
282 END IF
283
284 copy%mdata(:, :) = array(:, :)
285
286 CALL timestop(handle)
287
288 END SUBROUTINE copy_submatrix_data
289
290! **************************************************************************************************
291!> \brief ...
292!> \param submatrices ...
293!> \param scalar ...
294! **************************************************************************************************
295 SUBROUTINE set_submatrix_array(submatrices, scalar)
296
297 TYPE(domain_submatrix_type), DIMENSION(:), &
298 INTENT(INOUT) :: submatrices
299 REAL(KIND=dp), INTENT(IN) :: scalar
300
301 CHARACTER(len=*), PARAMETER :: routineN = 'set_submatrix_array'
302
303 INTEGER :: handle, idomain, ndomains
304
305 CALL timeset(routinen, handle)
306
307 ndomains = SIZE(submatrices)
308 DO idomain = 1, ndomains
309 IF (submatrices(idomain)%domain > 0) THEN
310 CALL set_submatrix(submatrices(idomain), scalar)
311 END IF
312 END DO ! loop over domains
313
314 CALL timestop(handle)
315
316 END SUBROUTINE set_submatrix_array
317
318! **************************************************************************************************
319!> \brief ...
320!> \param submatrix ...
321!> \param scalar ...
322! **************************************************************************************************
323 SUBROUTINE set_submatrix(submatrix, scalar)
324
325 TYPE(domain_submatrix_type), INTENT(INOUT) :: submatrix
326 REAL(KIND=dp), INTENT(IN) :: scalar
327
328 CHARACTER(len=*), PARAMETER :: routineN = 'set_submatrix'
329
330 INTEGER :: ds1, ds2, handle, ms1, ms2
331
332 CALL timeset(routinen, handle)
333
334 cpassert(submatrix%domain > 0)
335 cpassert(submatrix%nrows > 0)
336 cpassert(submatrix%ncols > 0)
337
338 ds1 = submatrix%nrows
339 ds2 = submatrix%ncols
340
341 IF (.NOT. ALLOCATED(submatrix%mdata)) THEN
342 ALLOCATE (submatrix%mdata(ds1, ds2))
343 ELSE
344 ms1 = SIZE(submatrix%mdata, 1)
345 ms2 = SIZE(submatrix%mdata, 2)
346 IF ((ds1 /= ms1) .OR. (ds2 /= ms2)) THEN
347 DEALLOCATE (submatrix%mdata)
348 ALLOCATE (submatrix%mdata(ds1, ds2))
349 END IF
350 END IF
351
352 submatrix%mdata(:, :) = scalar
353
354 CALL timestop(handle)
355
356 END SUBROUTINE set_submatrix
357
358! **************************************************************************************************
359!> \brief ...
360!> \param subm ...
361! **************************************************************************************************
362 SUBROUTINE release_submatrix_array(subm)
363
364 TYPE(domain_submatrix_type), DIMENSION(:), &
365 INTENT(INOUT) :: subm
366
367 CHARACTER(len=*), PARAMETER :: routineN = 'release_submatrix_array'
368
369 INTEGER :: handle, idomain, ndomains
370
371 CALL timeset(routinen, handle)
372
373 ndomains = SIZE(subm)
374 DO idomain = 1, ndomains
375 CALL release_submatrix(subm(idomain))
376 END DO ! loop over domains
377
378 CALL timestop(handle)
379
380 END SUBROUTINE release_submatrix_array
381
382! **************************************************************************************************
383!> \brief ...
384!> \param subm ...
385! **************************************************************************************************
386 SUBROUTINE release_submatrix(subm)
387
388 TYPE(domain_submatrix_type), INTENT(INOUT) :: subm
389
390 CHARACTER(len=*), PARAMETER :: routineN = 'release_submatrix'
391
392 INTEGER :: handle
393
394 CALL timeset(routinen, handle)
395
396 subm%domain = -1
397 subm%nbrows = -1
398 subm%nbcols = -1
399 subm%nrows = -1
400 subm%ncols = -1
401 subm%nnodes = -1
402 subm%group = mp_comm_null
403
404 IF (ALLOCATED(subm%dbcsr_row)) THEN
405 DEALLOCATE (subm%dbcsr_row)
406 END IF
407 IF (ALLOCATED(subm%dbcsr_col)) THEN
408 DEALLOCATE (subm%dbcsr_col)
409 END IF
410 IF (ALLOCATED(subm%size_brow)) THEN
411 DEALLOCATE (subm%size_brow)
412 END IF
413 IF (ALLOCATED(subm%size_bcol)) THEN
414 DEALLOCATE (subm%size_bcol)
415 END IF
416 IF (ALLOCATED(subm%mdata)) THEN
417 DEALLOCATE (subm%mdata)
418 END IF
419
420 CALL timestop(handle)
421
422 END SUBROUTINE release_submatrix
423
424 ! more complex routine might be necessary if submatrices are distributed
425! **************************************************************************************************
426!> \brief ...
427!> \param transA ...
428!> \param transB ...
429!> \param alpha ...
430!> \param A ...
431!> \param B ...
432!> \param beta ...
433!> \param C ...
434! **************************************************************************************************
435 SUBROUTINE multiply_submatrices_once(transA, transB, alpha, A, B, beta, C)
436
437 CHARACTER, INTENT(IN) :: transA, transB
438 REAL(KIND=dp), INTENT(IN) :: alpha
439 TYPE(domain_submatrix_type), INTENT(IN) :: a, b
440 REAL(kind=dp), INTENT(IN) :: beta
441 TYPE(domain_submatrix_type), INTENT(INOUT) :: c
442
443 CHARACTER(len=*), PARAMETER :: routinen = 'multiply_submatrices_once'
444
445 INTEGER :: cs1, cs2, handle, icol, irow, k, k1, &
446 lda, ldb, ldc, m, mblocks, n, nblocks
447 LOGICAL :: nota, notb
448
449 CALL timeset(routinen, handle)
450
451 cpassert(a%domain > 0)
452 cpassert(b%domain > 0)
453 cpassert(c%domain > 0)
454
455 lda = SIZE(a%mdata, 1)
456 ldb = SIZE(b%mdata, 1)
457
458 nota = (transa == 'N') .OR. (transa == 'n')
459 notb = (transb == 'N') .OR. (transb == 'n')
460
461 IF (nota) THEN
462 m = a%nrows
463 k = a%ncols
464 mblocks = a%nbrows
465 ELSE
466 m = a%ncols
467 k = a%nrows
468 mblocks = a%nbcols
469 END IF
470
471 IF (notb) THEN
472 k1 = b%nrows
473 n = b%ncols
474 nblocks = b%nbcols
475 ELSE
476 k1 = b%ncols
477 n = b%nrows
478 nblocks = b%nbrows
479 END IF
480
481 ! these checks are for debugging only
482 cpassert(k == k1)
483
484 ! conform C matrix
485 c%nrows = m
486 c%ncols = n
487 c%nbrows = mblocks
488 c%nbcols = nblocks
489 IF (ALLOCATED(c%dbcsr_row)) THEN
490 DEALLOCATE (c%dbcsr_row)
491 END IF
492 ALLOCATE (c%dbcsr_row(c%nbrows))
493 IF (ALLOCATED(c%dbcsr_col)) THEN
494 DEALLOCATE (c%dbcsr_col)
495 END IF
496 ALLOCATE (c%dbcsr_col(c%nbcols))
497 IF (ALLOCATED(c%size_brow)) THEN
498 DEALLOCATE (c%size_brow)
499 END IF
500 ALLOCATE (c%size_brow(c%nbrows))
501 IF (ALLOCATED(c%size_bcol)) THEN
502 DEALLOCATE (c%size_bcol)
503 END IF
504 ALLOCATE (c%size_bcol(c%nbcols))
505
506 DO irow = 1, c%nbrows
507 IF (nota) THEN
508 c%dbcsr_row(irow) = a%dbcsr_row(irow)
509 c%size_brow(irow) = a%size_brow(irow)
510 ELSE
511 c%dbcsr_row(irow) = a%dbcsr_col(irow)
512 c%size_brow(irow) = a%size_bcol(irow)
513 END IF
514 END DO
515
516 DO icol = 1, c%nbcols
517 IF (notb) THEN
518 c%dbcsr_col(icol) = b%dbcsr_col(icol)
519 c%size_bcol(icol) = b%size_bcol(icol)
520 ELSE
521 c%dbcsr_col(icol) = b%dbcsr_row(icol)
522 c%size_bcol(icol) = b%size_brow(icol)
523 END IF
524 END DO
525
526 IF (.NOT. ALLOCATED(c%mdata)) THEN
527 !!! cannot use non-zero beta if C is not allocated
528 cpassert(beta == 0.0_dp)
529 ALLOCATE (c%mdata(c%nrows, c%ncols))
530 ELSE
531 cs1 = SIZE(c%mdata, 1)
532 cs2 = SIZE(c%mdata, 2)
533 IF ((c%nrows /= cs1) .OR. (c%ncols /= cs2)) THEN
534 !!! cannot deallocate data if beta is non-zero
535 cpassert(beta == 0.0_dp)
536 DEALLOCATE (c%mdata)
537 ALLOCATE (c%mdata(c%nrows, c%ncols))
538 END IF
539 END IF
540
541 ldc = c%nrows
542
543 CALL dgemm(transa, transb, m, n, k, alpha, a%mdata, lda, b%mdata, ldb, beta, c%mdata, ldc)
544
545 c%nnodes = a%nnodes
546 c%group = a%group
547
548 CALL timestop(handle)
549
550 END SUBROUTINE multiply_submatrices_once
551
552! **************************************************************************************************
553!> \brief ...
554!> \param transA ...
555!> \param transB ...
556!> \param alpha ...
557!> \param A ...
558!> \param B ...
559!> \param beta ...
560!> \param C ...
561! **************************************************************************************************
562 SUBROUTINE multiply_submatrices_array(transA, transB, alpha, A, B, beta, C)
563
564 CHARACTER, INTENT(IN) :: transa, transb
565 REAL(kind=dp), INTENT(IN) :: alpha
566 TYPE(domain_submatrix_type), DIMENSION(:), &
567 INTENT(IN) :: a, b
568 REAL(kind=dp), INTENT(IN) :: beta
569 TYPE(domain_submatrix_type), DIMENSION(:), &
570 INTENT(INOUT) :: c
571
572 CHARACTER(len=*), PARAMETER :: routinen = 'multiply_submatrices_array'
573
574 INTEGER :: handle, idomain, idomaina, idomainb, &
575 ndomains, ndomainsb, ndomainsc
576
577 CALL timeset(routinen, handle)
578
579 ndomains = SIZE(a)
580 ndomainsb = SIZE(b)
581 ndomainsc = SIZE(c)
582
583 cpassert(ndomains == ndomainsb)
584 cpassert(ndomainsb == ndomainsc)
585
586 DO idomain = 1, ndomains
587
588 idomaina = a(idomain)%domain
589 idomainb = b(idomain)%domain
590
591 cpassert(idomaina == idomainb)
592
593 c(idomain)%domain = idomaina
594
595 ! check if the submatrix exists
596 IF (idomaina > 0) THEN
597 CALL multiply_submatrices_once(transa, transb, alpha, a(idomain), b(idomain), beta, c(idomain))
598 END IF ! submatrix for the domain exists
599
600 END DO ! loop over domains
601
602 CALL timestop(handle)
603
604 END SUBROUTINE multiply_submatrices_array
605
606 ! more complex routine might be necessary if submatrices are distributed
607! **************************************************************************************************
608!> \brief ...
609!> \param alpha ...
610!> \param A ...
611!> \param beta ...
612!> \param B ...
613!> \param transB ...
614! **************************************************************************************************
615 SUBROUTINE add_submatrices_once(alpha, A, beta, B, transB)
616
617 REAL(kind=dp), INTENT(IN) :: alpha
618 TYPE(domain_submatrix_type), INTENT(INOUT) :: a
619 REAL(kind=dp), INTENT(IN) :: beta
620 TYPE(domain_submatrix_type), INTENT(IN) :: b
621 CHARACTER, INTENT(IN) :: transb
622
623 CHARACTER(len=*), PARAMETER :: routinen = 'add_submatrices_once'
624
625 INTEGER :: c1, c2, handle, icol, r1, r2
626 LOGICAL :: notb
627
628 CALL timeset(routinen, handle)
629
630 cpassert(a%domain > 0)
631 cpassert(b%domain > 0)
632
633 r1 = a%nrows
634 c1 = a%ncols
635
636 notb = (transb == 'N') .OR. (transb == 'n')
637
638 IF (notb) THEN
639 r2 = b%nrows
640 c2 = b%ncols
641 ELSE
642 r2 = b%ncols
643 c2 = b%nrows
644 END IF
645
646 ! these checks are for debugging only
647 cpassert(c1 == c2)
648 cpassert(r1 == r2)
649
650 IF (notb) THEN
651 DO icol = 1, c1
652 a%mdata(:, icol) = alpha*a%mdata(:, icol) + beta*b%mdata(:, icol)
653 END DO
654 ELSE
655 DO icol = 1, c1
656 a%mdata(:, icol) = alpha*a%mdata(:, icol) + beta*b%mdata(icol, :)
657 END DO
658 END IF
659
660 CALL timestop(handle)
661
662 END SUBROUTINE add_submatrices_once
663
664! **************************************************************************************************
665!> \brief ...
666!> \param alpha ...
667!> \param A ...
668!> \param beta ...
669!> \param B ...
670!> \param transB ...
671! **************************************************************************************************
672 SUBROUTINE add_submatrices_array(alpha, A, beta, B, transB)
673
674 REAL(kind=dp), INTENT(IN) :: alpha
675 TYPE(domain_submatrix_type), DIMENSION(:), &
676 INTENT(INOUT) :: a
677 REAL(kind=dp), INTENT(IN) :: beta
678 TYPE(domain_submatrix_type), DIMENSION(:), &
679 INTENT(IN) :: b
680 CHARACTER, INTENT(IN) :: transb
681
682 CHARACTER(len=*), PARAMETER :: routinen = 'add_submatrices_array'
683
684 INTEGER :: handle, idomain, idomaina, idomainb, &
685 ndomains, ndomainsb
686
687 CALL timeset(routinen, handle)
688
689 ndomains = SIZE(a)
690 ndomainsb = SIZE(b)
691
692 cpassert(ndomains == ndomainsb)
693
694 DO idomain = 1, ndomains
695
696 idomaina = a(idomain)%domain
697 idomainb = b(idomain)%domain
698
699 cpassert(idomaina == idomainb)
700
701 ! check if the submatrix exists
702 IF (idomaina > 0) THEN
703 CALL add_submatrices_once(alpha, a(idomain), beta, b(idomain), transb)
704 END IF ! submatrix for the domain exists
705
706 END DO ! loop over domains
707
708 CALL timestop(handle)
709
710 END SUBROUTINE add_submatrices_array
711
712! **************************************************************************************************
713!> \brief Computes the max norm of the collection of submatrices
714!> \param submatrices ...
715!> \param norm ...
716!> \par History
717!> 2013.03 created [Rustam Z. Khaliullin]
718!> \author Rustam Z. Khaliullin
719! **************************************************************************************************
720 SUBROUTINE maxnorm_submatrices(submatrices, norm)
721
722 TYPE(domain_submatrix_type), DIMENSION(:), &
723 INTENT(IN) :: submatrices
724 REAL(kind=dp), INTENT(OUT) :: norm
725
726 CHARACTER(len=*), PARAMETER :: routinen = 'maxnorm_submatrices'
727
728 INTEGER :: handle, idomain, ndomains
729 REAL(kind=dp) :: curr_norm, send_norm
730 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_norm
731
732 CALL timeset(routinen, handle)
733
734 send_norm = 0.0_dp
735
736 ndomains = SIZE(submatrices)
737
738 DO idomain = 1, ndomains
739
740 ! check if the submatrix is local
741 IF (submatrices(idomain)%domain > 0) THEN
742 curr_norm = maxval(abs(submatrices(idomain)%mdata))
743 IF (curr_norm > send_norm) send_norm = curr_norm
744 END IF
745
746 END DO ! loop over domains
747
748 ! communicate local norm to the other nodes
749 ALLOCATE (recv_norm(submatrices(1)%nnodes))
750 CALL submatrices(1)%group%allgather(send_norm, recv_norm)
751
752 norm = maxval(recv_norm)
753
754 DEALLOCATE (recv_norm)
755
756 CALL timestop(handle)
757
758 END SUBROUTINE maxnorm_submatrices
759
760! **************************************************************************************************
761!> \brief Computes the sum of traces of the submatrix A.tr(B)
762!> \param A ...
763!> \param B ...
764!> \param trace ...
765!> \par History
766!> 2013.03 created [Rustam Z. Khaliullin]
767!> \author Rustam Z. Khaliullin
768! **************************************************************************************************
769 SUBROUTINE trace_submatrices(A, B, trace)
770
771 TYPE(domain_submatrix_type), DIMENSION(:), &
772 INTENT(IN) :: a, b
773 REAL(kind=dp), INTENT(OUT) :: trace
774
775 CHARACTER(len=*), PARAMETER :: routinen = 'trace_submatrices'
776
777 INTEGER :: domaina, domainb, handle, idomain, &
778 ndomainsa, ndomainsb
779 REAL(kind=dp) :: curr_trace, send_trace
780 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_trace
781
782 CALL timeset(routinen, handle)
783
784 send_trace = 0.0_dp
785
786 ndomainsa = SIZE(a)
787 ndomainsb = SIZE(b)
788
789 cpassert(ndomainsa == ndomainsb)
790
791 DO idomain = 1, ndomainsa
792
793 domaina = a(idomain)%domain
794 domainb = b(idomain)%domain
795
796 cpassert(domaina == domainb)
797
798 ! check if the submatrix is local
799 IF (domaina > 0) THEN
800
801 cpassert(a(idomain)%nrows == b(idomain)%nrows)
802 cpassert(a(idomain)%ncols == b(idomain)%ncols)
803
804 curr_trace = sum(a(idomain)%mdata(:, :)*b(idomain)%mdata(:, :))
805 send_trace = send_trace + curr_trace
806
807 END IF
808
809 END DO ! loop over domains
810
811 ! communicate local norm to the other nodes
812 ALLOCATE (recv_trace(a(1)%nnodes))
813 CALL a(1)%group%allgather(send_trace, recv_trace)
814
815 trace = sum(recv_trace)
816
817 DEALLOCATE (recv_trace)
818
819 CALL timestop(handle)
820
821 END SUBROUTINE trace_submatrices
822
823! **************************************************************************************************
824!> \brief Constructs submatrices for each ALMO domain by collecting distributed
825!> DBCSR blocks to local arrays
826!> \param matrix ...
827!> \param submatrix ...
828!> \param distr_pattern ...
829!> \param domain_map ...
830!> \param node_of_domain ...
831!> \param job_type ...
832!> \par History
833!> 2013.01 created [Rustam Z. Khaliullin]
834!> \author Rustam Z. Khaliullin
835! **************************************************************************************************
836 SUBROUTINE construct_submatrices(matrix, submatrix, distr_pattern, domain_map, &
837 node_of_domain, job_type)
838
839 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
840 TYPE(domain_submatrix_type), DIMENSION(:), &
841 INTENT(INOUT) :: submatrix
842 TYPE(dbcsr_type), INTENT(IN) :: distr_pattern
843 TYPE(domain_map_type), INTENT(IN) :: domain_map
844 INTEGER, DIMENSION(:), INTENT(IN) :: node_of_domain
845 INTEGER, INTENT(IN) :: job_type
846
847 CHARACTER(len=*), PARAMETER :: routinen = 'construct_submatrices'
848
849 CHARACTER :: matrix_type
850 INTEGER :: block_node, block_offset, col, col_offset, col_size, dest_node, groupid, handle, &
851 iblock, icol, idomain, index_col, index_ec, index_er, index_row, index_sc, index_sr, &
852 inode, ldesc, mynode, nblkcols_tot, nblkrows_tot, ndomains, ndomains2, nnodes, &
853 recv_size2_total, recv_size_total, row, row_size, send_size2_total, send_size_total, &
854 smcol, smrow, start_data
855 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, offset2_block, offset_block, &
856 recv_data2, recv_offset2_cpu, recv_offset_cpu, recv_size2_cpu, recv_size_cpu, send_data2, &
857 send_offset2_cpu, send_offset_cpu, send_size2_cpu, send_size_cpu
858 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: recv_descriptor, send_descriptor
859 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
860 LOGICAL :: found, transp
861 REAL(kind=dp) :: antifactor
862 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_data, send_data
863 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_p
864 TYPE(dbcsr_distribution_type) :: pattern_dist
865 TYPE(mp_comm_type) :: group
866
867!INTEGER, PARAMETER :: select_row_col = 1
868!INTEGER, PARAMETER :: select_row = 2
869! subm_row_size,&
870! subm_col_size,&
871
872 CALL timeset(routinen, handle)
873
874 CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
875 ndomains = nblkcols_tot ! RZK-warning not true for atomic distributions
876 CALL dbcsr_get_info(distr_pattern, distribution=pattern_dist)
877 CALL dbcsr_distribution_get(pattern_dist, numnodes=nnodes, group=groupid, mynode=mynode)
878
879 CALL group%set_handle(groupid)
880
881 matrix_type = dbcsr_get_matrix_type(matrix)
882
883 ldesc = 2
884 ALLOCATE (send_descriptor(ldesc, nnodes))
885 ALLOCATE (recv_descriptor(ldesc, nnodes))
886 send_descriptor(:, :) = 0
887
888 ! find: the number of blocks and their sizes that must be sent to each cpu
889 ! loop over all domains
890 DO idomain = 1, ndomains
891
892 dest_node = node_of_domain(idomain)
893
894 ! loop over those rows that have non-zero quencher
895 index_sr = 1 ! index start row
896 IF (idomain > 1) index_sr = domain_map%index1(idomain - 1)
897 index_er = domain_map%index1(idomain) - 1 ! index end row
898
899 DO index_row = index_sr, index_er
900
901 row = domain_map%pairs(index_row, 1)
902
903 IF (job_type == select_row_col) THEN
904 ! loop over those columns that have non-zero quencher
905 index_sc = 1 ! index start col
906 IF (idomain > 1) index_sc = domain_map%index1(idomain - 1)
907 index_ec = domain_map%index1(idomain) - 1 ! index end col
908 ELSE
909 ! fake loop
910 index_sc = 1 ! index start col
911 index_ec = 1 ! index end col
912 END IF
913
914 DO index_col = index_sc, index_ec
915
916 IF (job_type == select_row_col) THEN
917 col = domain_map%pairs(index_col, 1)
918 ELSE
919 col = idomain
920 END IF
921
922 transp = .false.
923 CALL dbcsr_get_stored_coordinates(matrix, &
924 row, col, block_node)
925 IF (block_node == mynode) THEN
926 CALL dbcsr_get_block_p(matrix, row, col, block_p, found, row_size, col_size)
927 IF (found) THEN
928 send_descriptor(1, dest_node + 1) = send_descriptor(1, dest_node + 1) + 1
929 send_descriptor(2, dest_node + 1) = send_descriptor(2, dest_node + 1) + &
930 row_size*col_size
931 END IF
932 END IF
933
934 END DO ! loop over columns
935
936 END DO ! loop over rows
937
938 END DO
939
940 ! communicate number of blocks and their sizes to the other nodes
941 CALL group%alltoall(send_descriptor, recv_descriptor, ldesc)
942
943 ALLOCATE (send_size_cpu(nnodes), send_offset_cpu(nnodes))
944 send_offset_cpu(1) = 0
945 send_size_cpu(1) = send_descriptor(2, 1)
946 DO inode = 2, nnodes
947 send_size_cpu(inode) = send_descriptor(2, inode)
948 send_offset_cpu(inode) = send_offset_cpu(inode - 1) + &
949 send_size_cpu(inode - 1)
950 END DO
951 send_size_total = send_offset_cpu(nnodes) + send_size_cpu(nnodes)
952
953 ALLOCATE (recv_size_cpu(nnodes), recv_offset_cpu(nnodes))
954 recv_offset_cpu(1) = 0
955 recv_size_cpu(1) = recv_descriptor(2, 1)
956 DO inode = 2, nnodes
957 recv_size_cpu(inode) = recv_descriptor(2, inode)
958 recv_offset_cpu(inode) = recv_offset_cpu(inode - 1) + &
959 recv_size_cpu(inode - 1)
960 END DO
961 recv_size_total = recv_offset_cpu(nnodes) + recv_size_cpu(nnodes)
962
963 ALLOCATE (send_size2_cpu(nnodes), send_offset2_cpu(nnodes))
964 send_offset2_cpu(1) = 0
965 send_size2_cpu(1) = 2*send_descriptor(1, 1)
966 DO inode = 2, nnodes
967 send_size2_cpu(inode) = 2*send_descriptor(1, inode)
968 send_offset2_cpu(inode) = send_offset2_cpu(inode - 1) + &
969 send_size2_cpu(inode - 1)
970 END DO
971 send_size2_total = send_offset2_cpu(nnodes) + send_size2_cpu(nnodes)
972
973 ALLOCATE (recv_size2_cpu(nnodes), recv_offset2_cpu(nnodes))
974 recv_offset2_cpu(1) = 0
975 recv_size2_cpu(1) = 2*recv_descriptor(1, 1)
976 DO inode = 2, nnodes
977 recv_size2_cpu(inode) = 2*recv_descriptor(1, inode)
978 recv_offset2_cpu(inode) = recv_offset2_cpu(inode - 1) + &
979 recv_size2_cpu(inode - 1)
980 END DO
981 recv_size2_total = recv_offset2_cpu(nnodes) + recv_size2_cpu(nnodes)
982
983 DEALLOCATE (send_descriptor)
984 DEALLOCATE (recv_descriptor)
985
986 ! collect data from the matrix into send_data
987 ALLOCATE (send_data(send_size_total))
988 ALLOCATE (recv_data(recv_size_total))
989 ALLOCATE (send_data2(send_size2_total))
990 ALLOCATE (recv_data2(recv_size2_total))
991 ALLOCATE (offset_block(nnodes))
992 ALLOCATE (offset2_block(nnodes))
993 offset_block(:) = 0
994 offset2_block(:) = 0
995 ! loop over all domains
996 DO idomain = 1, ndomains
997
998 dest_node = node_of_domain(idomain)
999
1000 ! loop over those rows that have non-zero quencher
1001 index_sr = 1 ! index start row
1002 IF (idomain > 1) index_sr = domain_map%index1(idomain - 1)
1003 index_er = domain_map%index1(idomain) - 1 ! index end row
1004
1005 DO index_row = index_sr, index_er
1006
1007 row = domain_map%pairs(index_row, 1)
1008
1009 IF (job_type == select_row_col) THEN
1010 ! loop over those columns that have non-zero quencher
1011 index_sc = 1 ! index start col
1012 IF (idomain > 1) index_sc = domain_map%index1(idomain - 1)
1013 index_ec = domain_map%index1(idomain) - 1 ! index end col
1014 ELSE
1015 ! fake loop
1016 index_sc = 1 ! index start col
1017 index_ec = 1 ! index end col
1018 END IF
1019
1020 DO index_col = index_sc, index_ec
1021
1022 IF (job_type == select_row_col) THEN
1023 col = domain_map%pairs(index_col, 1)
1024 ELSE
1025 col = idomain
1026 END IF
1027
1028 transp = .false.
1029 CALL dbcsr_get_stored_coordinates(matrix, &
1030 row, col, block_node)
1031 IF (block_node == mynode) THEN
1032 CALL dbcsr_get_block_p(matrix, row, col, block_p, found, row_size, col_size)
1033 IF (found) THEN
1034 col_offset = row_size*col_size
1035 start_data = send_offset_cpu(dest_node + 1) + &
1036 offset_block(dest_node + 1)
1037 send_data(start_data + 1:start_data + col_offset) = reshape(block_p, [col_offset])
1038 offset_block(dest_node + 1) = offset_block(dest_node + 1) + col_offset
1039 ! fill out row,col information
1040 send_data2(send_offset2_cpu(dest_node + 1) + &
1041 offset2_block(dest_node + 1) + 1) = row
1042 send_data2(send_offset2_cpu(dest_node + 1) + &
1043 offset2_block(dest_node + 1) + 2) = col
1044 offset2_block(dest_node + 1) = offset2_block(dest_node + 1) + 2
1045 END IF
1046 END IF
1047
1048 END DO ! loop over columns
1049
1050 END DO ! loop over rows
1051
1052 END DO
1053
1054 ! send-receive all blocks
1055 CALL group%alltoall(send_data, send_size_cpu, send_offset_cpu, &
1056 recv_data, recv_size_cpu, recv_offset_cpu)
1057 ! send-receive rows and cols of the blocks
1058 CALL group%alltoall(send_data2, send_size2_cpu, send_offset2_cpu, &
1059 recv_data2, recv_size2_cpu, recv_offset2_cpu)
1060
1061 DEALLOCATE (send_size_cpu, send_offset_cpu)
1062 DEALLOCATE (send_size2_cpu, send_offset2_cpu)
1063 DEALLOCATE (send_data)
1064 DEALLOCATE (send_data2)
1065 DEALLOCATE (offset_block)
1066 DEALLOCATE (offset2_block)
1067
1068 ! copy blocks into submatrices
1069 CALL dbcsr_get_info(matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
1070 ndomains2 = SIZE(submatrix)
1071 IF (ndomains2 /= ndomains) THEN
1072 cpabort("wrong submatrix size")
1073 END IF
1074 CALL release_submatrices(submatrix)
1075 submatrix(:)%nnodes = nnodes
1076 submatrix(:)%group = group
1077 submatrix(:)%nrows = 0
1078 submatrix(:)%ncols = 0
1079
1080 ALLOCATE (first_row(nblkrows_tot), first_col(nblkcols_tot))
1081 submatrix(:)%domain = -1
1082 DO idomain = 1, ndomains
1083 dest_node = node_of_domain(idomain)
1084 IF (dest_node == mynode) THEN
1085 submatrix(idomain)%domain = idomain
1086 submatrix(idomain)%nbrows = 0
1087 submatrix(idomain)%nbcols = 0
1088
1089 ! loop over those rows that have non-zero quencher
1090 first_row(:) = -1
1091 index_sr = 1 ! index start row
1092 IF (idomain > 1) index_sr = domain_map%index1(idomain - 1)
1093 index_er = domain_map%index1(idomain) - 1 ! index end row
1094 DO index_row = index_sr, index_er
1095 row = domain_map%pairs(index_row, 1)
1096 first_row(row) = submatrix(idomain)%nrows + 1
1097 submatrix(idomain)%nrows = submatrix(idomain)%nrows + row_blk_size(row)
1098 submatrix(idomain)%nbrows = submatrix(idomain)%nbrows + 1
1099 END DO
1100 ALLOCATE (submatrix(idomain)%dbcsr_row(submatrix(idomain)%nbrows))
1101 ALLOCATE (submatrix(idomain)%size_brow(submatrix(idomain)%nbrows))
1102 smrow = 1
1103 ! again loop over those rows that have non-zero quencher
1104 index_sr = 1 ! index start row
1105 IF (idomain > 1) index_sr = domain_map%index1(idomain - 1)
1106 index_er = domain_map%index1(idomain) - 1 ! index end row
1107 DO index_row = index_sr, index_er
1108 row = domain_map%pairs(index_row, 1)
1109 submatrix(idomain)%dbcsr_row(smrow) = row
1110 submatrix(idomain)%size_brow(smrow) = row_blk_size(row)
1111 smrow = smrow + 1
1112 END DO
1113
1114 ! loop over the necessary columns
1115 first_col(:) = -1
1116 IF (job_type == select_row_col) THEN
1117 ! loop over those columns that have non-zero quencher
1118 index_sc = 1 ! index start col
1119 IF (idomain > 1) index_sc = domain_map%index1(idomain - 1)
1120 index_ec = domain_map%index1(idomain) - 1 ! index end col
1121 ELSE
1122 ! fake loop
1123 index_sc = 1 ! index start col
1124 index_ec = 1 ! index end col
1125 END IF
1126 DO index_col = index_sc, index_ec
1127 IF (job_type == select_row_col) THEN
1128 col = domain_map%pairs(index_col, 1)
1129 ELSE
1130 col = idomain
1131 END IF
1132 first_col(col) = submatrix(idomain)%ncols + 1
1133 submatrix(idomain)%ncols = submatrix(idomain)%ncols + col_blk_size(col)
1134 submatrix(idomain)%nbcols = submatrix(idomain)%nbcols + 1
1135 END DO
1136
1137 ALLOCATE (submatrix(idomain)%dbcsr_col(submatrix(idomain)%nbcols))
1138 ALLOCATE (submatrix(idomain)%size_bcol(submatrix(idomain)%nbcols))
1139
1140 ! loop over the necessary columns again
1141 smcol = 1
1142 IF (job_type == select_row_col) THEN
1143 ! loop over those columns that have non-zero quencher
1144 index_sc = 1 ! index start col
1145 IF (idomain > 1) index_sc = domain_map%index1(idomain - 1)
1146 index_ec = domain_map%index1(idomain) - 1 ! index end col
1147 ELSE
1148 ! fake loop
1149 index_sc = 1 ! index start col
1150 index_ec = 1 ! index end col
1151 END IF
1152 DO index_col = index_sc, index_ec
1153 IF (job_type == select_row_col) THEN
1154 col = domain_map%pairs(index_col, 1)
1155 ELSE
1156 col = idomain
1157 END IF
1158 submatrix(idomain)%dbcsr_col(smcol) = col
1159 submatrix(idomain)%size_bcol(smcol) = col_blk_size(col)
1160 smcol = smcol + 1
1161 END DO
1162
1163 ALLOCATE (submatrix(idomain)%mdata( &
1164 submatrix(idomain)%nrows, &
1165 submatrix(idomain)%ncols))
1166 submatrix(idomain)%mdata(:, :) = 0.0_dp
1167 DO inode = 1, nnodes
1168 block_offset = 0
1169 DO iblock = 1, recv_size2_cpu(inode)/2
1170 ! read the (row,col) of the block
1171 row = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 1)
1172 col = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 2)
1173 ! check if this block should be in the submatrix of this domain
1174 IF ((first_col(col) /= -1) .AND. (first_row(row) /= -1)) THEN
1175 ! copy data from the received array into submatrix
1176 start_data = recv_offset_cpu(inode) + block_offset + 1
1177 DO icol = 0, col_blk_size(col) - 1
1178 submatrix(idomain)%mdata(first_row(row): &
1179 first_row(row) + row_blk_size(row) - 1, &
1180 first_col(col) + icol) = &
1181 recv_data(start_data:start_data + row_blk_size(row) - 1)
1182 start_data = start_data + row_blk_size(row)
1183 END DO
1184 IF (job_type == select_row_col) THEN
1185 IF (matrix_type == dbcsr_type_symmetric .OR. &
1186 matrix_type == dbcsr_type_antisymmetric) THEN
1187 ! copy data into the transposed block as well
1188 antifactor = 1.0_dp
1189 IF (matrix_type == dbcsr_type_antisymmetric) THEN
1190 antifactor = -1.0_dp
1191 END IF
1192 start_data = recv_offset_cpu(inode) + block_offset + 1
1193 DO icol = 0, col_blk_size(col) - 1
1194 submatrix(idomain)%mdata(first_row(col) + icol, &
1195 first_col(row): &
1196 first_col(row) + row_blk_size(row) - 1) = &
1197 antifactor*recv_data(start_data: &
1198 start_data + row_blk_size(row) - 1)
1199 start_data = start_data + row_blk_size(row)
1200 END DO
1201 ELSE IF (matrix_type == dbcsr_type_no_symmetry) THEN
1202 ELSE
1203 cpabort("matrix type is NYI")
1204 END IF
1205 END IF
1206 END IF
1207 block_offset = block_offset + col_blk_size(col)*row_blk_size(row)
1208 END DO
1209 END DO
1210 END IF ! mynode.eq.dest_node
1211 END DO ! loop over domains
1212
1213 DEALLOCATE (recv_size_cpu, recv_offset_cpu)
1214 DEALLOCATE (recv_size2_cpu, recv_offset2_cpu)
1215 DEALLOCATE (recv_data)
1216 DEALLOCATE (recv_data2)
1217! DEALLOCATE(subm_row_size,subm_col_size)
1218 DEALLOCATE (first_row, first_col)
1219
1220 CALL timestop(handle)
1221
1222 END SUBROUTINE construct_submatrices
1223
1224! **************************************************************************************************
1225!> \brief Constructs a DBCSR matrix from submatrices
1226!> \param matrix ...
1227!> \param submatrix ...
1228!> \param distr_pattern ...
1229!> \par History
1230!> 2013.01 created [Rustam Z. Khaliullin]
1231!> \author Rustam Z. Khaliullin
1232! **************************************************************************************************
1233 SUBROUTINE construct_dbcsr_from_submatrices(matrix, submatrix, distr_pattern)
1234
1235 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1236 TYPE(domain_submatrix_type), DIMENSION(:), &
1237 INTENT(IN) :: submatrix
1238 TYPE(dbcsr_type), INTENT(IN) :: distr_pattern
1239
1240 CHARACTER(len=*), PARAMETER :: routinen = 'construct_dbcsr_from_submatrices'
1241
1242 CHARACTER :: matrix_type
1243 INTEGER :: block_offset, col, col_offset, colsize, dest_node, groupid, handle, iblock, icol, &
1244 idomain, inode, irow_subm, ldesc, mynode, nblkcols_tot, nblkrows_tot, ndomains, &
1245 ndomains2, nnodes, recv_size2_total, recv_size_total, row, rowsize, send_size2_total, &
1246 send_size_total, smroff, start_data, unit_nr
1247 INTEGER, ALLOCATABLE, DIMENSION(:) :: offset2_block, offset_block, recv_data2, &
1248 recv_offset2_cpu, recv_offset_cpu, recv_size2_cpu, recv_size_cpu, send_data2, &
1249 send_offset2_cpu, send_offset_cpu, send_size2_cpu, send_size_cpu
1250 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: recv_descriptor, send_descriptor
1251 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
1252 LOGICAL :: transp
1253 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_data, send_data
1254 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
1255 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_p
1256 TYPE(cp_logger_type), POINTER :: logger
1257 TYPE(dbcsr_distribution_type) :: pattern_dist
1258 TYPE(dbcsr_iterator_type) :: iter
1259 TYPE(mp_comm_type) :: group
1260
1261 CALL timeset(routinen, handle)
1262
1263 ! get a useful output_unit
1264 logger => cp_get_default_logger()
1265 IF (logger%para_env%is_source()) THEN
1266 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1267 ELSE
1268 unit_nr = -1
1269 END IF
1270
1271 CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
1272 ndomains = nblkcols_tot ! RZK-warning not true for atomic distributions
1273 ndomains2 = SIZE(submatrix)
1274
1275 IF (ndomains /= ndomains2) THEN
1276 cpabort("domain mismatch")
1277 END IF
1278
1279 CALL dbcsr_get_info(distr_pattern, distribution=pattern_dist)
1280 CALL dbcsr_distribution_get(pattern_dist, numnodes=nnodes, group=groupid, mynode=mynode)
1281
1282 CALL group%set_handle(groupid)
1283
1284 matrix_type = dbcsr_get_matrix_type(matrix)
1285 IF (matrix_type /= dbcsr_type_no_symmetry) THEN
1286 cpabort("only non-symmetric matrices so far")
1287 END IF
1288
1289 ! remove all blocks from the dbcsr matrix
1290 CALL dbcsr_iterator_start(iter, matrix)
1291 DO WHILE (dbcsr_iterator_blocks_left(iter))
1292 CALL dbcsr_iterator_next_block(iter, row, col, block_p)
1293 block_p(:, :) = 0.0_dp
1294 END DO
1295 CALL dbcsr_iterator_stop(iter)
1296 CALL dbcsr_filter(matrix, 0.1_dp)
1297
1298 CALL dbcsr_work_create(matrix, work_mutable=.true.)
1299
1300 ldesc = 2
1301 ALLOCATE (send_descriptor(ldesc, nnodes))
1302 ALLOCATE (recv_descriptor(ldesc, nnodes))
1303 send_descriptor(:, :) = 0
1304
1305 ! loop over domains - find how much data to send
1306 DO idomain = 1, ndomains
1307
1308 IF (submatrix(idomain)%domain > 0) THEN
1309
1310 DO irow_subm = 1, submatrix(idomain)%nbrows
1311
1312 IF (submatrix(idomain)%nbcols /= 1) THEN
1313 cpabort("corrupt submatrix structure")
1314 END IF
1315
1316 row = submatrix(idomain)%dbcsr_row(irow_subm)
1317 col = submatrix(idomain)%dbcsr_col(1)
1318
1319 IF (col /= idomain) THEN
1320 cpabort("corrupt submatrix structure")
1321 END IF
1322
1323 transp = .false.
1324 CALL dbcsr_get_stored_coordinates(distr_pattern, &
1325 row, idomain, dest_node)
1326
1327 send_descriptor(1, dest_node + 1) = send_descriptor(1, dest_node + 1) + 1
1328 send_descriptor(2, dest_node + 1) = send_descriptor(2, dest_node + 1) + &
1329 submatrix(idomain)%size_brow(irow_subm)* &
1330 submatrix(idomain)%size_bcol(1)
1331
1332 END DO ! loop over submatrix blocks
1333
1334 END IF
1335
1336 END DO ! loop over domains
1337
1338 ! communicate number of blocks and their sizes to the other nodes
1339 CALL group%alltoall(send_descriptor, recv_descriptor, ldesc)
1340
1341 ALLOCATE (send_size_cpu(nnodes), send_offset_cpu(nnodes))
1342 send_offset_cpu(1) = 0
1343 send_size_cpu(1) = send_descriptor(2, 1)
1344 DO inode = 2, nnodes
1345 send_size_cpu(inode) = send_descriptor(2, inode)
1346 send_offset_cpu(inode) = send_offset_cpu(inode - 1) + &
1347 send_size_cpu(inode - 1)
1348 END DO
1349 send_size_total = send_offset_cpu(nnodes) + send_size_cpu(nnodes)
1350
1351 ALLOCATE (recv_size_cpu(nnodes), recv_offset_cpu(nnodes))
1352 recv_offset_cpu(1) = 0
1353 recv_size_cpu(1) = recv_descriptor(2, 1)
1354 DO inode = 2, nnodes
1355 recv_size_cpu(inode) = recv_descriptor(2, inode)
1356 recv_offset_cpu(inode) = recv_offset_cpu(inode - 1) + &
1357 recv_size_cpu(inode - 1)
1358 END DO
1359 recv_size_total = recv_offset_cpu(nnodes) + recv_size_cpu(nnodes)
1360
1361 ALLOCATE (send_size2_cpu(nnodes), send_offset2_cpu(nnodes))
1362 send_offset2_cpu(1) = 0
1363 send_size2_cpu(1) = 2*send_descriptor(1, 1)
1364 DO inode = 2, nnodes
1365 send_size2_cpu(inode) = 2*send_descriptor(1, inode)
1366 send_offset2_cpu(inode) = send_offset2_cpu(inode - 1) + &
1367 send_size2_cpu(inode - 1)
1368 END DO
1369 send_size2_total = send_offset2_cpu(nnodes) + send_size2_cpu(nnodes)
1370
1371 ALLOCATE (recv_size2_cpu(nnodes), recv_offset2_cpu(nnodes))
1372 recv_offset2_cpu(1) = 0
1373 recv_size2_cpu(1) = 2*recv_descriptor(1, 1)
1374 DO inode = 2, nnodes
1375 recv_size2_cpu(inode) = 2*recv_descriptor(1, inode)
1376 recv_offset2_cpu(inode) = recv_offset2_cpu(inode - 1) + &
1377 recv_size2_cpu(inode - 1)
1378 END DO
1379 recv_size2_total = recv_offset2_cpu(nnodes) + recv_size2_cpu(nnodes)
1380
1381 DEALLOCATE (send_descriptor)
1382 DEALLOCATE (recv_descriptor)
1383
1384 ! collect data from the matrix into send_data
1385 ALLOCATE (send_data(send_size_total))
1386 ALLOCATE (recv_data(recv_size_total))
1387 ALLOCATE (send_data2(send_size2_total))
1388 ALLOCATE (recv_data2(recv_size2_total))
1389 ALLOCATE (offset_block(nnodes))
1390 ALLOCATE (offset2_block(nnodes))
1391 offset_block(:) = 0
1392 offset2_block(:) = 0
1393 ! loop over domains - collect data to send
1394 DO idomain = 1, ndomains
1395
1396 IF (submatrix(idomain)%domain > 0) THEN
1397
1398 smroff = 0
1399 DO irow_subm = 1, submatrix(idomain)%nbrows
1400
1401 row = submatrix(idomain)%dbcsr_row(irow_subm)
1402 col = submatrix(idomain)%dbcsr_col(1)
1403
1404 rowsize = submatrix(idomain)%size_brow(irow_subm)
1405 colsize = submatrix(idomain)%size_bcol(1)
1406
1407 transp = .false.
1408 CALL dbcsr_get_stored_coordinates(distr_pattern, &
1409 row, idomain, dest_node)
1410
1411 ! place the data appropriately
1412 col_offset = 0
1413 DO icol = 1, colsize
1414 start_data = send_offset_cpu(dest_node + 1) + &
1415 offset_block(dest_node + 1) + &
1416 col_offset
1417 send_data(start_data + 1:start_data + rowsize) = &
1418 submatrix(idomain)%mdata(smroff + 1:smroff + rowsize, icol)
1419 col_offset = col_offset + rowsize
1420 END DO
1421 offset_block(dest_node + 1) = offset_block(dest_node + 1) + &
1422 colsize*rowsize
1423 ! fill out row,col information
1424 send_data2(send_offset2_cpu(dest_node + 1) + &
1425 offset2_block(dest_node + 1) + 1) = row
1426 send_data2(send_offset2_cpu(dest_node + 1) + &
1427 offset2_block(dest_node + 1) + 2) = col
1428 offset2_block(dest_node + 1) = offset2_block(dest_node + 1) + 2
1429
1430 smroff = smroff + rowsize
1431
1432 END DO
1433
1434 END IF
1435
1436 END DO ! loop over domains
1437
1438 ! send-receive all blocks
1439 CALL group%alltoall(send_data, send_size_cpu, send_offset_cpu, &
1440 recv_data, recv_size_cpu, recv_offset_cpu)
1441 ! send-receive rows and cols of the blocks
1442 CALL group%alltoall(send_data2, send_size2_cpu, send_offset2_cpu, &
1443 recv_data2, recv_size2_cpu, recv_offset2_cpu)
1444
1445 DEALLOCATE (send_size_cpu, send_offset_cpu)
1446 DEALLOCATE (send_size2_cpu, send_offset2_cpu)
1447 DEALLOCATE (send_data)
1448 DEALLOCATE (send_data2)
1449 DEALLOCATE (offset_block)
1450 DEALLOCATE (offset2_block)
1451
1452 ! copy received data into dbcsr matrix
1453 CALL dbcsr_get_info(matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
1454 DO inode = 1, nnodes
1455 block_offset = 0
1456 DO iblock = 1, recv_size2_cpu(inode)/2
1457 ! read the (row,col) of the block
1458 row = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 1)
1459 col = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 2)
1460 ! copy data from the received array into the matrix block
1461 start_data = recv_offset_cpu(inode) + block_offset + 1
1462 ALLOCATE (new_block(row_blk_size(row), col_blk_size(col)))
1463 DO icol = 1, col_blk_size(col)
1464 new_block(:, icol) = &
1465 recv_data(start_data:start_data + row_blk_size(row) - 1)
1466 start_data = start_data + row_blk_size(row)
1467 END DO
1468 CALL dbcsr_put_block(matrix, row, col, new_block)
1469 DEALLOCATE (new_block)
1470 block_offset = block_offset + col_blk_size(col)*row_blk_size(row)
1471 END DO
1472 END DO
1473
1474 DEALLOCATE (recv_size_cpu, recv_offset_cpu)
1475 DEALLOCATE (recv_size2_cpu, recv_offset2_cpu)
1476 DEALLOCATE (recv_data)
1477 DEALLOCATE (recv_data2)
1478
1479 CALL dbcsr_finalize(matrix)
1480
1481 CALL timestop(handle)
1482
1483 END SUBROUTINE construct_dbcsr_from_submatrices
1484
1485! **************************************************************************************************
1486!> \brief ...
1487!> \param submatrices ...
1488!> \param mpgroup ...
1489! **************************************************************************************************
1490 SUBROUTINE print_submatrices(submatrices, mpgroup)
1491
1492 TYPE(domain_submatrix_type), DIMENSION(:), &
1493 INTENT(IN) :: submatrices
1494 TYPE(mp_comm_type), INTENT(IN) :: mpgroup
1495
1496 CHARACTER(len=*), PARAMETER :: routinen = 'print_submatrices'
1497
1498 CHARACTER(len=30) :: colstr, formatstr
1499 INTEGER :: handle, i, irow, n, ncols, nrows
1500
1501 CALL timeset(routinen, handle)
1502
1503 n = SIZE(submatrices)
1504
1505 DO i = 1, n
1506 nrows = SIZE(submatrices(i)%mdata, 1)
1507 ncols = SIZE(submatrices(i)%mdata, 2)
1508 WRITE (colstr, *) ncols
1509 formatstr = '('//trim(adjustl(colstr))//'F16.9)'
1510 IF (submatrices(i)%domain > 0) THEN
1511 WRITE (*, *) "SUBMATRIX: ", i, nrows, 'x', ncols
1512 nrows = SIZE(submatrices(i)%mdata, 1)
1513 DO irow = 1, nrows
1514 WRITE (*, formatstr) submatrices(i)%mdata(irow, :)
1515 END DO
1516 END IF
1517 CALL mpgroup%sync()
1518 END DO
1519
1520 CALL timestop(handle)
1521
1522 END SUBROUTINE print_submatrices
1523
1524! **************************************************************************************************
1525!> \brief Reports whether the DBCSR block (row,col) exists in the quencher
1526!> \param map ...
1527!> \param row ...
1528!> \param col ...
1529!> \return ...
1530!> \par History
1531!> 2013.01 created [Rustam Z. Khaliullin]
1532!> \author Rustam Z. Khaliullin
1533! **************************************************************************************************
1534 FUNCTION qblk_exists(map, row, col)
1535
1536 TYPE(domain_map_type), INTENT(IN) :: map
1537 INTEGER, INTENT(IN) :: row, col
1538 LOGICAL :: qblk_exists
1539
1540 INTEGER :: first, last, mid, ndomains
1541
1542!CALL timeset(routineN,handle)
1543
1544 ndomains = SIZE(map%index1)
1545
1546 qblk_exists = .false.
1547 IF (col < 1 .OR. col > ndomains) RETURN
1548 first = 1
1549 IF (col > 1) first = map%index1(col - 1)
1550 last = map%index1(col) - 1
1551
1552 ! perform binary search within first-last
1553 DO WHILE (last >= first)
1554 mid = first + (last - first)/2
1555 IF (map%pairs(mid, 1) > row) THEN
1556 last = mid - 1
1557 ELSE IF (map%pairs(mid, 1) < row) THEN
1558 first = mid + 1
1559 ELSE
1560 qblk_exists = .true. ! SUCCESS!!
1561 EXIT
1562 END IF
1563 END DO
1564
1565 !CALL timestop(handle)
1566
1567 RETURN
1568
1569 END FUNCTION qblk_exists
1570
1571END MODULE domain_submatrix_methods
1572
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
domain_submatrix_methods::add_submatrices
Definition domain_submatrix_methods.F:72
domain_submatrix_methods::copy_submatrices
Definition domain_submatrix_methods.F:57
domain_submatrix_methods::init_submatrices
Definition domain_submatrix_methods.F:46
domain_submatrix_methods::multiply_submatrices
Definition domain_submatrix_methods.F:67
domain_submatrix_methods::release_submatrices
Definition domain_submatrix_methods.F:62
domain_submatrix_methods::set_submatrices
Definition domain_submatrix_methods.F:52
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_get_stored_coordinates
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
Definition cp_dbcsr_api.F:911
cp_dbcsr_api::dbcsr_work_create
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
Definition cp_dbcsr_api.F:1307
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
domain_submatrix_methods
Subroutines to handle submatrices.
Definition domain_submatrix_methods.F:14
domain_submatrix_methods::copy_submatrix_data
subroutine, public copy_submatrix_data(array, copy)
...
Definition domain_submatrix_methods.F:258
domain_submatrix_methods::print_submatrices
subroutine, public print_submatrices(submatrices, mpgroup)
...
Definition domain_submatrix_methods.F:1491
domain_submatrix_methods::construct_dbcsr_from_submatrices
subroutine, public construct_dbcsr_from_submatrices(matrix, submatrix, distr_pattern)
Constructs a DBCSR matrix from submatrices.
Definition domain_submatrix_methods.F:1234
domain_submatrix_methods::maxnorm_submatrices
subroutine, public maxnorm_submatrices(submatrices, norm)
Computes the max norm of the collection of submatrices.
Definition domain_submatrix_methods.F:721
domain_submatrix_methods::construct_submatrices
subroutine, public construct_submatrices(matrix, submatrix, distr_pattern, domain_map, node_of_domain, job_type)
Constructs submatrices for each ALMO domain by collecting distributed DBCSR blocks to local arrays.
Definition domain_submatrix_methods.F:838
domain_submatrix_types
Types to handle submatrices.
Definition domain_submatrix_types.F:14
domain_submatrix_types::select_row_col
integer, parameter, public select_row_col
Definition domain_submatrix_types.F:25
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:23
message_passing::mp_comm_null
type(mp_comm_type), parameter, public mp_comm_null
Definition message_passing.F:798
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
domain_submatrix_types::domain_map_type
Definition domain_submatrix_types.F:47
domain_submatrix_types::domain_submatrix_type
Definition domain_submatrix_types.F:32
message_passing::mp_comm_type
Definition message_passing.F:161