(git:ed6f26b)
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-2025 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! **************************************************************************************************
15
16 USE cp_dbcsr_api, ONLY: &
21 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
28 USE kinds, ONLY: dp
29 USE message_passing, ONLY: mp_comm_null,&
31#include "./base/base_uses.f90"
32
33 IMPLICIT NONE
34
35 PRIVATE
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'domain_submatrix_methods'
38
45
47 MODULE PROCEDURE init_submatrices_0d
48 MODULE PROCEDURE init_submatrices_1d
49 MODULE PROCEDURE init_submatrices_2d
50 END INTERFACE
51
53 MODULE PROCEDURE set_submatrix_array
54 MODULE PROCEDURE set_submatrix
55 END INTERFACE
56
58 MODULE PROCEDURE copy_submatrix_array
59 MODULE PROCEDURE copy_submatrix
60 END INTERFACE
61
63 MODULE PROCEDURE release_submatrix_array
64 MODULE PROCEDURE release_submatrix
65 END INTERFACE
66
68 MODULE PROCEDURE multiply_submatrices_once
69 MODULE PROCEDURE multiply_submatrices_array
70 END INTERFACE
71
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 .EQ. ndomainsb)
158 copy(:)%nnodes = original(:)%nnodes
159 copy(:)%group = original(:)%group
160 DO idomain = 1, ndomains
161 IF (original(idomain)%domain .GT. 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 .GT. 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) .NE. 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) .NE. 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) .NE. 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) .NE. 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 .GT. 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 .NE. ms1) .OR. (ds2 .NE. 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 .GT. 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 .GT. 0)
335 cpassert(submatrix%nrows .GT. 0)
336 cpassert(submatrix%ncols .GT. 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 .NE. ms1) .OR. (ds2 .NE. 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 .GT. 0)
452 cpassert(b%domain .GT. 0)
453 cpassert(c%domain .GT. 0)
454
455 lda = SIZE(a%mdata, 1)
456 ldb = SIZE(b%mdata, 1)
457
458 nota = (transa .EQ. 'N') .OR. (transa .EQ. 'n')
459 notb = (transb .EQ. 'N') .OR. (transb .EQ. '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 .EQ. 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 .EQ. 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 .NE. cs1) .OR. (c%ncols .NE. cs2)) THEN
534 !!! cannot deallocate data if beta is non-zero
535 cpassert(beta .EQ. 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 .EQ. ndomainsb)
584 cpassert(ndomainsb .EQ. ndomainsc)
585
586 DO idomain = 1, ndomains
587
588 idomaina = a(idomain)%domain
589 idomainb = b(idomain)%domain
590
591 cpassert(idomaina .EQ. idomainb)
592
593 c(idomain)%domain = idomaina
594
595 ! check if the submatrix exists
596 IF (idomaina .GT. 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 .GT. 0)
631 cpassert(b%domain .GT. 0)
632
633 r1 = a%nrows
634 c1 = a%ncols
635
636 notb = (transb .EQ. 'N') .OR. (transb .EQ. '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 .EQ. c2)
648 cpassert(r1 .EQ. 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 .EQ. ndomainsb)
693
694 DO idomain = 1, ndomains
695
696 idomaina = a(idomain)%domain
697 idomainb = b(idomain)%domain
698
699 cpassert(idomaina .EQ. idomainb)
700
701 ! check if the submatrix exists
702 IF (idomaina .GT. 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 .GT. 0) THEN
742 curr_norm = maxval(abs(submatrices(idomain)%mdata))
743 IF (curr_norm .GT. 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 .EQ. ndomainsb)
790
791 DO idomain = 1, ndomainsa
792
793 domaina = a(idomain)%domain
794 domainb = b(idomain)%domain
795
796 cpassert(domaina .EQ. domainb)
797
798 ! check if the submatrix is local
799 IF (domaina .GT. 0) THEN
800
801 cpassert(a(idomain)%nrows .EQ. b(idomain)%nrows)
802 cpassert(a(idomain)%ncols .EQ. 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 .GT. 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 .GT. 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 .EQ. 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 ! simple but quadratically scaling procedure
941 ! loop over local blocks
942 !CALL dbcsr_iterator_start(iter,matrix)
943 !DO WHILE (dbcsr_iterator_blocks_left(iter))
944 ! CALL dbcsr_iterator_next_block(iter,row,col,data_p,&
945 ! row_size=row_size,col_size=col_size)
946 ! DO idomain = 1, ndomains
947 ! IF (job_type==select_row_col) THEN
948 ! domain_needs_block=(qblk_exists(domain_map,col,idomain)&
949 ! .AND.qblk_exists(domain_map,row,idomain))
950 ! ELSE
951 ! domain_needs_block=(idomain==col&
952 ! .AND.qblk_exists(domain_map,row,idomain))
953 ! ENDIF
954 ! IF (domain_needs_block) THEN
955 ! transp=.FALSE.
956 ! dest_node=node_of_domain(idomain)
957 ! !CALL dbcsr_get_stored_coordinates(distr_pattern,&
958 ! ! idomain, idomain, transp, dest_node)
959 ! send_descriptor(1,dest_node+1)=send_descriptor(1,dest_node+1)+1
960 ! send_descriptor(2,dest_node+1)=send_descriptor(2,dest_node+1)+&
961 ! row_size*col_size
962 ! ENDIF
963 ! ENDDO
964 !ENDDO
965 !CALL dbcsr_iterator_stop(iter)
966
967 ! communicate number of blocks and their sizes to the other nodes
968 CALL group%alltoall(send_descriptor, recv_descriptor, ldesc)
969
970 ALLOCATE (send_size_cpu(nnodes), send_offset_cpu(nnodes))
971 send_offset_cpu(1) = 0
972 send_size_cpu(1) = send_descriptor(2, 1)
973 DO inode = 2, nnodes
974 send_size_cpu(inode) = send_descriptor(2, inode)
975 send_offset_cpu(inode) = send_offset_cpu(inode - 1) + &
976 send_size_cpu(inode - 1)
977 END DO
978 send_size_total = send_offset_cpu(nnodes) + send_size_cpu(nnodes)
979
980 ALLOCATE (recv_size_cpu(nnodes), recv_offset_cpu(nnodes))
981 recv_offset_cpu(1) = 0
982 recv_size_cpu(1) = recv_descriptor(2, 1)
983 DO inode = 2, nnodes
984 recv_size_cpu(inode) = recv_descriptor(2, inode)
985 recv_offset_cpu(inode) = recv_offset_cpu(inode - 1) + &
986 recv_size_cpu(inode - 1)
987 END DO
988 recv_size_total = recv_offset_cpu(nnodes) + recv_size_cpu(nnodes)
989
990 ALLOCATE (send_size2_cpu(nnodes), send_offset2_cpu(nnodes))
991 send_offset2_cpu(1) = 0
992 send_size2_cpu(1) = 2*send_descriptor(1, 1)
993 DO inode = 2, nnodes
994 send_size2_cpu(inode) = 2*send_descriptor(1, inode)
995 send_offset2_cpu(inode) = send_offset2_cpu(inode - 1) + &
996 send_size2_cpu(inode - 1)
997 END DO
998 send_size2_total = send_offset2_cpu(nnodes) + send_size2_cpu(nnodes)
999
1000 ALLOCATE (recv_size2_cpu(nnodes), recv_offset2_cpu(nnodes))
1001 recv_offset2_cpu(1) = 0
1002 recv_size2_cpu(1) = 2*recv_descriptor(1, 1)
1003 DO inode = 2, nnodes
1004 recv_size2_cpu(inode) = 2*recv_descriptor(1, inode)
1005 recv_offset2_cpu(inode) = recv_offset2_cpu(inode - 1) + &
1006 recv_size2_cpu(inode - 1)
1007 END DO
1008 recv_size2_total = recv_offset2_cpu(nnodes) + recv_size2_cpu(nnodes)
1009
1010 DEALLOCATE (send_descriptor)
1011 DEALLOCATE (recv_descriptor)
1012
1013 ! collect data from the matrix into send_data
1014 ALLOCATE (send_data(send_size_total))
1015 ALLOCATE (recv_data(recv_size_total))
1016 ALLOCATE (send_data2(send_size2_total))
1017 ALLOCATE (recv_data2(recv_size2_total))
1018 ALLOCATE (offset_block(nnodes))
1019 ALLOCATE (offset2_block(nnodes))
1020 offset_block(:) = 0
1021 offset2_block(:) = 0
1022 ! loop over all domains
1023 DO idomain = 1, ndomains
1024
1025 dest_node = node_of_domain(idomain)
1026
1027 ! loop over those rows that have non-zero quencher
1028 index_sr = 1 ! index start row
1029 IF (idomain .GT. 1) index_sr = domain_map%index1(idomain - 1)
1030 index_er = domain_map%index1(idomain) - 1 ! index end row
1031
1032 DO index_row = index_sr, index_er
1033
1034 row = domain_map%pairs(index_row, 1)
1035
1036 IF (job_type == select_row_col) THEN
1037 ! loop over those columns that have non-zero quencher
1038 index_sc = 1 ! index start col
1039 IF (idomain .GT. 1) index_sc = domain_map%index1(idomain - 1)
1040 index_ec = domain_map%index1(idomain) - 1 ! index end col
1041 ELSE
1042 ! fake loop
1043 index_sc = 1 ! index start col
1044 index_ec = 1 ! index end col
1045 END IF
1046
1047 DO index_col = index_sc, index_ec
1048
1049 IF (job_type == select_row_col) THEN
1050 col = domain_map%pairs(index_col, 1)
1051 ELSE
1052 col = idomain
1053 END IF
1054
1055 transp = .false.
1056 CALL dbcsr_get_stored_coordinates(matrix, &
1057 row, col, block_node)
1058 IF (block_node .EQ. mynode) THEN
1059 CALL dbcsr_get_block_p(matrix, row, col, block_p, found, row_size, col_size)
1060 IF (found) THEN
1061 !col_offset=0
1062 !DO icol=1,col_size
1063 ! start_data=send_offset_cpu(dest_node+1)+&
1064 ! offset_block(dest_node+1)+&
1065 ! col_offset
1066 ! send_data(start_data+1:start_data+row_size)=&
1067 ! data_p(1:row_size,icol)
1068 ! col_offset=col_offset+row_size
1069 !ENDDO
1070 col_offset = row_size*col_size
1071 start_data = send_offset_cpu(dest_node + 1) + &
1072 offset_block(dest_node + 1)
1073 send_data(start_data + 1:start_data + col_offset) = reshape(block_p, [col_offset])
1074 offset_block(dest_node + 1) = offset_block(dest_node + 1) + col_offset
1075 ! fill out row,col information
1076 send_data2(send_offset2_cpu(dest_node + 1) + &
1077 offset2_block(dest_node + 1) + 1) = row
1078 send_data2(send_offset2_cpu(dest_node + 1) + &
1079 offset2_block(dest_node + 1) + 2) = col
1080 offset2_block(dest_node + 1) = offset2_block(dest_node + 1) + 2
1081 END IF
1082 END IF
1083
1084 END DO ! loop over columns
1085
1086 END DO ! loop over rows
1087
1088 END DO
1089 ! more simple but quadratically scaling version
1090 !CALL dbcsr_iterator_start(iter,matrix)
1091 !DO WHILE (dbcsr_iterator_blocks_left(iter))
1092 ! CALL dbcsr_iterator_next_block(iter,row,col,data_p,&
1093 ! row_size=row_size,col_size=col_size)
1094 ! DO idomain = 1, ndomains
1095 ! IF (job_type==select_row_col) THEN
1096 ! domain_needs_block=(qblk_exists(domain_map,col,idomain)&
1097 ! .AND.qblk_exists(domain_map,row,idomain))
1098 ! ELSE
1099 ! domain_needs_block=(idomain==col&
1100 ! .AND.qblk_exists(domain_map,row,idomain))
1101 ! ENDIF
1102 ! IF (domain_needs_block) THEN
1103 ! transp=.FALSE.
1104 ! dest_node=node_of_domain(idomain)
1105 ! !CALL dbcsr_get_stored_coordinates(distr_pattern,&
1106 ! ! idomain, idomain, transp, dest_node)
1107 ! ! place the data appropriately
1108 ! col_offset=0
1109 ! DO icol=1,col_size
1110 ! start_data=send_offset_cpu(dest_node+1)+&
1111 ! offset_block(dest_node+1)+&
1112 ! col_offset
1113 ! send_data(start_data+1:start_data+row_size)=&
1114 ! data_p(1:row_size,icol)
1115 ! col_offset=col_offset+row_size
1116 ! ENDDO
1117 ! offset_block(dest_node+1)=offset_block(dest_node+1)+col_size*row_size
1118 ! ! fill out row,col information
1119 ! send_data2(send_offset2_cpu(dest_node+1)+&
1120 ! offset2_block(dest_node+1)+1)=row
1121 ! send_data2(send_offset2_cpu(dest_node+1)+&
1122 ! offset2_block(dest_node+1)+2)=col
1123 ! offset2_block(dest_node+1)=offset2_block(dest_node+1)+2
1124 ! ENDIF
1125 ! ENDDO
1126 !ENDDO
1127 !CALL dbcsr_iterator_stop(iter)
1128
1129 ! send-receive all blocks
1130 CALL group%alltoall(send_data, send_size_cpu, send_offset_cpu, &
1131 recv_data, recv_size_cpu, recv_offset_cpu)
1132 ! send-receive rows and cols of the blocks
1133 CALL group%alltoall(send_data2, send_size2_cpu, send_offset2_cpu, &
1134 recv_data2, recv_size2_cpu, recv_offset2_cpu)
1135
1136 DEALLOCATE (send_size_cpu, send_offset_cpu)
1137 DEALLOCATE (send_size2_cpu, send_offset2_cpu)
1138 DEALLOCATE (send_data)
1139 DEALLOCATE (send_data2)
1140 DEALLOCATE (offset_block)
1141 DEALLOCATE (offset2_block)
1142
1143 ! copy blocks into submatrices
1144 CALL dbcsr_get_info(matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
1145! ALLOCATE(subm_row_size(ndomains),subm_col_size(ndomains))
1146! subm_row_size(:)=0
1147! subm_col_size(:)=0
1148 ndomains2 = SIZE(submatrix)
1149 IF (ndomains2 .NE. ndomains) THEN
1150 cpabort("wrong submatrix size")
1151 END IF
1152 CALL release_submatrices(submatrix)
1153 submatrix(:)%nnodes = nnodes
1154 submatrix(:)%group = group
1155 submatrix(:)%nrows = 0
1156 submatrix(:)%ncols = 0
1157
1158 ALLOCATE (first_row(nblkrows_tot), first_col(nblkcols_tot))
1159 submatrix(:)%domain = -1
1160 DO idomain = 1, ndomains
1161 dest_node = node_of_domain(idomain)
1162 !transp=.FALSE.
1163 !CALL dbcsr_get_stored_coordinates(distr_pattern,&
1164 ! idomain, idomain, transp, dest_node)
1165 IF (dest_node .EQ. mynode) THEN
1166 submatrix(idomain)%domain = idomain
1167 submatrix(idomain)%nbrows = 0
1168 submatrix(idomain)%nbcols = 0
1169
1170 ! loop over those rows that have non-zero quencher
1171 first_row(:) = -1
1172 index_sr = 1 ! index start row
1173 IF (idomain .GT. 1) index_sr = domain_map%index1(idomain - 1)
1174 index_er = domain_map%index1(idomain) - 1 ! index end row
1175 DO index_row = index_sr, index_er
1176 row = domain_map%pairs(index_row, 1)
1177 !DO row = 1, nblkrows_tot
1178 ! IF (qblk_exists(domain_map,row,idomain)) THEN
1179 first_row(row) = submatrix(idomain)%nrows + 1
1180 submatrix(idomain)%nrows = submatrix(idomain)%nrows + row_blk_size(row)
1181 submatrix(idomain)%nbrows = submatrix(idomain)%nbrows + 1
1182 ! ENDIF
1183 END DO
1184 ALLOCATE (submatrix(idomain)%dbcsr_row(submatrix(idomain)%nbrows))
1185 ALLOCATE (submatrix(idomain)%size_brow(submatrix(idomain)%nbrows))
1186 smrow = 1
1187 ! again loop over those rows that have non-zero quencher
1188 index_sr = 1 ! index start row
1189 IF (idomain .GT. 1) index_sr = domain_map%index1(idomain - 1)
1190 index_er = domain_map%index1(idomain) - 1 ! index end row
1191 DO index_row = index_sr, index_er
1192 row = domain_map%pairs(index_row, 1)
1193 !DO row = 1, nblkrows_tot
1194 ! IF (first_row(row).ne.-1) THEN
1195 submatrix(idomain)%dbcsr_row(smrow) = row
1196 submatrix(idomain)%size_brow(smrow) = row_blk_size(row)
1197 smrow = smrow + 1
1198 ! ENDIF
1199 END DO
1200
1201 ! loop over the necessary columns
1202 first_col(:) = -1
1203 IF (job_type == select_row_col) THEN
1204 ! loop over those columns that have non-zero quencher
1205 index_sc = 1 ! index start col
1206 IF (idomain .GT. 1) index_sc = domain_map%index1(idomain - 1)
1207 index_ec = domain_map%index1(idomain) - 1 ! index end col
1208 ELSE
1209 ! fake loop
1210 index_sc = 1 ! index start col
1211 index_ec = 1 ! index end col
1212 END IF
1213 DO index_col = index_sc, index_ec
1214 IF (job_type == select_row_col) THEN
1215 col = domain_map%pairs(index_col, 1)
1216 ELSE
1217 col = idomain
1218 END IF
1219 !DO col = 1, nblkcols_tot
1220 ! IF (job_type==select_row_col) THEN
1221 ! domain_needs_block=(qblk_exists(domain_map,col,idomain))
1222 ! ELSE
1223 ! domain_needs_block=(col==idomain) ! RZK-warning col belongs to the domain
1224 ! ENDIF
1225 ! IF (domain_needs_block) THEN
1226 first_col(col) = submatrix(idomain)%ncols + 1
1227 submatrix(idomain)%ncols = submatrix(idomain)%ncols + col_blk_size(col)
1228 submatrix(idomain)%nbcols = submatrix(idomain)%nbcols + 1
1229 ! ENDIF
1230 END DO
1231
1232 ALLOCATE (submatrix(idomain)%dbcsr_col(submatrix(idomain)%nbcols))
1233 ALLOCATE (submatrix(idomain)%size_bcol(submatrix(idomain)%nbcols))
1234
1235 ! loop over the necessary columns again
1236 smcol = 1
1237 IF (job_type == select_row_col) THEN
1238 ! loop over those columns that have non-zero quencher
1239 index_sc = 1 ! index start col
1240 IF (idomain .GT. 1) index_sc = domain_map%index1(idomain - 1)
1241 index_ec = domain_map%index1(idomain) - 1 ! index end col
1242 ELSE
1243 ! fake loop
1244 index_sc = 1 ! index start col
1245 index_ec = 1 ! index end col
1246 END IF
1247 DO index_col = index_sc, index_ec
1248 IF (job_type == select_row_col) THEN
1249 col = domain_map%pairs(index_col, 1)
1250 ELSE
1251 col = idomain
1252 END IF
1253 !DO col = 1, nblkcols_tot
1254 ! IF (first_col(col).ne.-1) THEN
1255 submatrix(idomain)%dbcsr_col(smcol) = col
1256 submatrix(idomain)%size_bcol(smcol) = col_blk_size(col)
1257 smcol = smcol + 1
1258 ! ENDIF
1259 END DO
1260
1261 ALLOCATE (submatrix(idomain)%mdata( &
1262 submatrix(idomain)%nrows, &
1263 submatrix(idomain)%ncols))
1264 submatrix(idomain)%mdata(:, :) = 0.0_dp
1265 DO inode = 1, nnodes
1266 block_offset = 0
1267 DO iblock = 1, recv_size2_cpu(inode)/2
1268 ! read the (row,col) of the block
1269 row = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 1)
1270 col = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 2)
1271 ! check if this block should be in the submatrix of this domain
1272 IF ((first_col(col) .NE. -1) .AND. (first_row(row) .NE. -1)) THEN
1273 ! copy data from the received array into submatrix
1274 start_data = recv_offset_cpu(inode) + block_offset + 1
1275 DO icol = 0, col_blk_size(col) - 1
1276 submatrix(idomain)%mdata(first_row(row): &
1277 first_row(row) + row_blk_size(row) - 1, &
1278 first_col(col) + icol) = &
1279 recv_data(start_data:start_data + row_blk_size(row) - 1)
1280 start_data = start_data + row_blk_size(row)
1281 END DO
1282 IF (job_type == select_row_col) THEN
1283 IF (matrix_type == dbcsr_type_symmetric .OR. &
1284 matrix_type == dbcsr_type_antisymmetric) THEN
1285 ! copy data into the transposed block as well
1286 antifactor = 1.0_dp
1287 IF (matrix_type == dbcsr_type_antisymmetric) THEN
1288 antifactor = -1.0_dp
1289 END IF
1290 start_data = recv_offset_cpu(inode) + block_offset + 1
1291 DO icol = 0, col_blk_size(col) - 1
1292 submatrix(idomain)%mdata(first_row(col) + icol, &
1293 first_col(row): &
1294 first_col(row) + row_blk_size(row) - 1) = &
1295 antifactor*recv_data(start_data: &
1296 start_data + row_blk_size(row) - 1)
1297 start_data = start_data + row_blk_size(row)
1298 END DO
1299 ELSE IF (matrix_type == dbcsr_type_no_symmetry) THEN
1300 ELSE
1301 cpabort("matrix type is NYI")
1302 END IF
1303 END IF
1304 END IF
1305 block_offset = block_offset + col_blk_size(col)*row_blk_size(row)
1306 END DO
1307 END DO
1308 END IF ! mynode.eq.dest_node
1309 END DO ! loop over domains
1310
1311 DEALLOCATE (recv_size_cpu, recv_offset_cpu)
1312 DEALLOCATE (recv_size2_cpu, recv_offset2_cpu)
1313 DEALLOCATE (recv_data)
1314 DEALLOCATE (recv_data2)
1315! DEALLOCATE(subm_row_size,subm_col_size)
1316 DEALLOCATE (first_row, first_col)
1317
1318 CALL timestop(handle)
1319
1320 END SUBROUTINE construct_submatrices
1321
1322! **************************************************************************************************
1323!> \brief Constructs a DBCSR matrix from submatrices
1324!> \param matrix ...
1325!> \param submatrix ...
1326!> \param distr_pattern ...
1327!> \par History
1328!> 2013.01 created [Rustam Z. Khaliullin]
1329!> \author Rustam Z. Khaliullin
1330! **************************************************************************************************
1331 SUBROUTINE construct_dbcsr_from_submatrices(matrix, submatrix, distr_pattern)
1332
1333 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1334 TYPE(domain_submatrix_type), DIMENSION(:), &
1335 INTENT(IN) :: submatrix
1336 TYPE(dbcsr_type), INTENT(IN) :: distr_pattern
1337
1338 CHARACTER(len=*), PARAMETER :: routinen = 'construct_dbcsr_from_submatrices'
1339
1340 CHARACTER :: matrix_type
1341 INTEGER :: block_offset, col, col_offset, colsize, dest_node, groupid, handle, iblock, icol, &
1342 idomain, inode, irow_subm, ldesc, mynode, nblkcols_tot, nblkrows_tot, ndomains, &
1343 ndomains2, nnodes, recv_size2_total, recv_size_total, row, rowsize, send_size2_total, &
1344 send_size_total, smroff, start_data, unit_nr
1345 INTEGER, ALLOCATABLE, DIMENSION(:) :: offset2_block, offset_block, recv_data2, &
1346 recv_offset2_cpu, recv_offset_cpu, recv_size2_cpu, recv_size_cpu, send_data2, &
1347 send_offset2_cpu, send_offset_cpu, send_size2_cpu, send_size_cpu
1348 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: recv_descriptor, send_descriptor
1349 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
1350 LOGICAL :: transp
1351 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_data, send_data
1352 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
1353 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_p
1354 TYPE(cp_logger_type), POINTER :: logger
1355 TYPE(dbcsr_distribution_type) :: pattern_dist
1356 TYPE(dbcsr_iterator_type) :: iter
1357 TYPE(mp_comm_type) :: group
1358
1359 CALL timeset(routinen, handle)
1360
1361 ! get a useful output_unit
1362 logger => cp_get_default_logger()
1363 IF (logger%para_env%is_source()) THEN
1364 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1365 ELSE
1366 unit_nr = -1
1367 END IF
1368
1369 CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
1370 ndomains = nblkcols_tot ! RZK-warning not true for atomic distributions
1371 ndomains2 = SIZE(submatrix)
1372
1373 IF (ndomains .NE. ndomains2) THEN
1374 cpabort("domain mismatch")
1375 END IF
1376
1377 CALL dbcsr_get_info(distr_pattern, distribution=pattern_dist)
1378 CALL dbcsr_distribution_get(pattern_dist, numnodes=nnodes, group=groupid, mynode=mynode)
1379
1380 CALL group%set_handle(groupid)
1381
1382 matrix_type = dbcsr_get_matrix_type(matrix)
1383 IF (matrix_type .NE. dbcsr_type_no_symmetry) THEN
1384 cpabort("only non-symmetric matrices so far")
1385 END IF
1386
1387 ! remove all blocks from the dbcsr matrix
1388 CALL dbcsr_iterator_start(iter, matrix)
1389 DO WHILE (dbcsr_iterator_blocks_left(iter))
1390 CALL dbcsr_iterator_next_block(iter, row, col, block_p)
1391 block_p(:, :) = 0.0_dp
1392 END DO
1393 CALL dbcsr_iterator_stop(iter)
1394 CALL dbcsr_filter(matrix, 0.1_dp)
1395
1396 CALL dbcsr_work_create(matrix, work_mutable=.true.)
1397
1398 ldesc = 2
1399 ALLOCATE (send_descriptor(ldesc, nnodes))
1400 ALLOCATE (recv_descriptor(ldesc, nnodes))
1401 send_descriptor(:, :) = 0
1402
1403 ! loop over domains - find how much data to send
1404 DO idomain = 1, ndomains
1405
1406 IF (submatrix(idomain)%domain .GT. 0) THEN
1407
1408 DO irow_subm = 1, submatrix(idomain)%nbrows
1409
1410 IF (submatrix(idomain)%nbcols .NE. 1) THEN
1411 cpabort("corrupt submatrix structure")
1412 END IF
1413
1414 row = submatrix(idomain)%dbcsr_row(irow_subm)
1415 col = submatrix(idomain)%dbcsr_col(1)
1416
1417 IF (col .NE. idomain) THEN
1418 cpabort("corrupt submatrix structure")
1419 END IF
1420
1421 transp = .false.
1422 CALL dbcsr_get_stored_coordinates(distr_pattern, &
1423 row, idomain, dest_node)
1424
1425 send_descriptor(1, dest_node + 1) = send_descriptor(1, dest_node + 1) + 1
1426 send_descriptor(2, dest_node + 1) = send_descriptor(2, dest_node + 1) + &
1427 submatrix(idomain)%size_brow(irow_subm)* &
1428 submatrix(idomain)%size_bcol(1)
1429
1430 END DO ! loop over submatrix blocks
1431
1432 END IF
1433
1434 END DO ! loop over domains
1435
1436 ! communicate number of blocks and their sizes to the other nodes
1437 CALL group%alltoall(send_descriptor, recv_descriptor, ldesc)
1438
1439 ALLOCATE (send_size_cpu(nnodes), send_offset_cpu(nnodes))
1440 send_offset_cpu(1) = 0
1441 send_size_cpu(1) = send_descriptor(2, 1)
1442 DO inode = 2, nnodes
1443 send_size_cpu(inode) = send_descriptor(2, inode)
1444 send_offset_cpu(inode) = send_offset_cpu(inode - 1) + &
1445 send_size_cpu(inode - 1)
1446 END DO
1447 send_size_total = send_offset_cpu(nnodes) + send_size_cpu(nnodes)
1448
1449 ALLOCATE (recv_size_cpu(nnodes), recv_offset_cpu(nnodes))
1450 recv_offset_cpu(1) = 0
1451 recv_size_cpu(1) = recv_descriptor(2, 1)
1452 DO inode = 2, nnodes
1453 recv_size_cpu(inode) = recv_descriptor(2, inode)
1454 recv_offset_cpu(inode) = recv_offset_cpu(inode - 1) + &
1455 recv_size_cpu(inode - 1)
1456 END DO
1457 recv_size_total = recv_offset_cpu(nnodes) + recv_size_cpu(nnodes)
1458
1459 ALLOCATE (send_size2_cpu(nnodes), send_offset2_cpu(nnodes))
1460 send_offset2_cpu(1) = 0
1461 send_size2_cpu(1) = 2*send_descriptor(1, 1)
1462 DO inode = 2, nnodes
1463 send_size2_cpu(inode) = 2*send_descriptor(1, inode)
1464 send_offset2_cpu(inode) = send_offset2_cpu(inode - 1) + &
1465 send_size2_cpu(inode - 1)
1466 END DO
1467 send_size2_total = send_offset2_cpu(nnodes) + send_size2_cpu(nnodes)
1468
1469 ALLOCATE (recv_size2_cpu(nnodes), recv_offset2_cpu(nnodes))
1470 recv_offset2_cpu(1) = 0
1471 recv_size2_cpu(1) = 2*recv_descriptor(1, 1)
1472 DO inode = 2, nnodes
1473 recv_size2_cpu(inode) = 2*recv_descriptor(1, inode)
1474 recv_offset2_cpu(inode) = recv_offset2_cpu(inode - 1) + &
1475 recv_size2_cpu(inode - 1)
1476 END DO
1477 recv_size2_total = recv_offset2_cpu(nnodes) + recv_size2_cpu(nnodes)
1478
1479 DEALLOCATE (send_descriptor)
1480 DEALLOCATE (recv_descriptor)
1481
1482 ! collect data from the matrix into send_data
1483 ALLOCATE (send_data(send_size_total))
1484 ALLOCATE (recv_data(recv_size_total))
1485 ALLOCATE (send_data2(send_size2_total))
1486 ALLOCATE (recv_data2(recv_size2_total))
1487 ALLOCATE (offset_block(nnodes))
1488 ALLOCATE (offset2_block(nnodes))
1489 offset_block(:) = 0
1490 offset2_block(:) = 0
1491 ! loop over domains - collect data to send
1492 DO idomain = 1, ndomains
1493
1494 IF (submatrix(idomain)%domain .GT. 0) THEN
1495
1496 smroff = 0
1497 DO irow_subm = 1, submatrix(idomain)%nbrows
1498
1499 row = submatrix(idomain)%dbcsr_row(irow_subm)
1500 col = submatrix(idomain)%dbcsr_col(1)
1501
1502 rowsize = submatrix(idomain)%size_brow(irow_subm)
1503 colsize = submatrix(idomain)%size_bcol(1)
1504
1505 transp = .false.
1506 CALL dbcsr_get_stored_coordinates(distr_pattern, &
1507 row, idomain, dest_node)
1508
1509 ! place the data appropriately
1510 col_offset = 0
1511 DO icol = 1, colsize
1512 start_data = send_offset_cpu(dest_node + 1) + &
1513 offset_block(dest_node + 1) + &
1514 col_offset
1515 send_data(start_data + 1:start_data + rowsize) = &
1516 submatrix(idomain)%mdata(smroff + 1:smroff + rowsize, icol)
1517 col_offset = col_offset + rowsize
1518 END DO
1519 offset_block(dest_node + 1) = offset_block(dest_node + 1) + &
1520 colsize*rowsize
1521 ! fill out row,col information
1522 send_data2(send_offset2_cpu(dest_node + 1) + &
1523 offset2_block(dest_node + 1) + 1) = row
1524 send_data2(send_offset2_cpu(dest_node + 1) + &
1525 offset2_block(dest_node + 1) + 2) = col
1526 offset2_block(dest_node + 1) = offset2_block(dest_node + 1) + 2
1527
1528 smroff = smroff + rowsize
1529
1530 END DO
1531
1532 END IF
1533
1534 END DO ! loop over domains
1535
1536 ! send-receive all blocks
1537 CALL group%alltoall(send_data, send_size_cpu, send_offset_cpu, &
1538 recv_data, recv_size_cpu, recv_offset_cpu)
1539 ! send-receive rows and cols of the blocks
1540 CALL group%alltoall(send_data2, send_size2_cpu, send_offset2_cpu, &
1541 recv_data2, recv_size2_cpu, recv_offset2_cpu)
1542
1543 DEALLOCATE (send_size_cpu, send_offset_cpu)
1544 DEALLOCATE (send_size2_cpu, send_offset2_cpu)
1545 DEALLOCATE (send_data)
1546 DEALLOCATE (send_data2)
1547 DEALLOCATE (offset_block)
1548 DEALLOCATE (offset2_block)
1549
1550 ! copy received data into dbcsr matrix
1551 CALL dbcsr_get_info(matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
1552 DO inode = 1, nnodes
1553 block_offset = 0
1554 DO iblock = 1, recv_size2_cpu(inode)/2
1555 ! read the (row,col) of the block
1556 row = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 1)
1557 col = recv_data2(recv_offset2_cpu(inode) + (iblock - 1)*2 + 2)
1558 ! copy data from the received array into the matrix block
1559 start_data = recv_offset_cpu(inode) + block_offset + 1
1560 ALLOCATE (new_block(row_blk_size(row), col_blk_size(col)))
1561 DO icol = 1, col_blk_size(col)
1562 new_block(:, icol) = &
1563 recv_data(start_data:start_data + row_blk_size(row) - 1)
1564 start_data = start_data + row_blk_size(row)
1565 END DO
1566 CALL dbcsr_put_block(matrix, row, col, new_block)
1567 DEALLOCATE (new_block)
1568 block_offset = block_offset + col_blk_size(col)*row_blk_size(row)
1569 END DO
1570 END DO
1571
1572 DEALLOCATE (recv_size_cpu, recv_offset_cpu)
1573 DEALLOCATE (recv_size2_cpu, recv_offset2_cpu)
1574 DEALLOCATE (recv_data)
1575 DEALLOCATE (recv_data2)
1576
1577 CALL dbcsr_finalize(matrix)
1578
1579 CALL timestop(handle)
1580
1582
1583! **************************************************************************************************
1584!> \brief ...
1585!> \param submatrices ...
1586!> \param mpgroup ...
1587! **************************************************************************************************
1588 SUBROUTINE print_submatrices(submatrices, mpgroup)
1589
1590 TYPE(domain_submatrix_type), DIMENSION(:), &
1591 INTENT(IN) :: submatrices
1592 TYPE(mp_comm_type), INTENT(IN) :: mpgroup
1593
1594 CHARACTER(len=*), PARAMETER :: routinen = 'print_submatrices'
1595
1596 CHARACTER(len=30) :: colstr, formatstr
1597 INTEGER :: handle, i, irow, n, ncols, nrows
1598
1599 CALL timeset(routinen, handle)
1600
1601 n = SIZE(submatrices)
1602
1603 DO i = 1, n
1604 nrows = SIZE(submatrices(i)%mdata, 1)
1605 ncols = SIZE(submatrices(i)%mdata, 2)
1606 WRITE (colstr, *) ncols
1607 formatstr = '('//trim(adjustl(colstr))//'F16.9)'
1608 IF (submatrices(i)%domain .GT. 0) THEN
1609 WRITE (*, *) "SUBMATRIX: ", i, nrows, 'x', ncols
1610 nrows = SIZE(submatrices(i)%mdata, 1)
1611 DO irow = 1, nrows
1612 WRITE (*, formatstr) submatrices(i)%mdata(irow, :)
1613 END DO
1614 END IF
1615 CALL mpgroup%sync()
1616 END DO
1617
1618 CALL timestop(handle)
1619
1620 END SUBROUTINE print_submatrices
1621
1622! **************************************************************************************************
1623!> \brief Reports whether the DBCSR block (row,col) exists in the quencher
1624!> \param map ...
1625!> \param row ...
1626!> \param col ...
1627!> \return ...
1628!> \par History
1629!> 2013.01 created [Rustam Z. Khaliullin]
1630!> \author Rustam Z. Khaliullin
1631! **************************************************************************************************
1632 FUNCTION qblk_exists(map, row, col)
1633
1634 TYPE(domain_map_type), INTENT(IN) :: map
1635 INTEGER, INTENT(IN) :: row, col
1636 LOGICAL :: qblk_exists
1637
1638 INTEGER :: first, last, mid, ndomains
1639
1640!CALL timeset(routineN,handle)
1641
1642 ndomains = SIZE(map%index1)
1643
1644 qblk_exists = .false.
1645 IF (col .LT. 1 .OR. col .GT. ndomains) RETURN
1646 first = 1
1647 IF (col .GT. 1) first = map%index1(col - 1)
1648 last = map%index1(col) - 1
1649
1650 ! perform binary search within first-last
1651 DO WHILE (last .GE. first)
1652 mid = first + (last - first)/2
1653 IF (map%pairs(mid, 1) .GT. row) THEN
1654 last = mid - 1
1655 ELSE IF (map%pairs(mid, 1) .LT. row) THEN
1656 first = mid + 1
1657 ELSE
1658 qblk_exists = .true. ! SUCCESS!!
1659 EXIT
1660 END IF
1661 END DO
1662
1663 !CALL timestop(handle)
1664
1665 RETURN
1666
1667 END FUNCTION qblk_exists
1668
1669END MODULE domain_submatrix_methods
1670
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.
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
character function, public dbcsr_get_matrix_type(matrix)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
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)
...
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
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)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Subroutines to handle submatrices.
subroutine, public copy_submatrix_data(array, copy)
...
subroutine, public print_submatrices(submatrices, mpgroup)
...
subroutine, public construct_dbcsr_from_submatrices(matrix, submatrix, distr_pattern)
Constructs a DBCSR matrix from submatrices.
subroutine, public maxnorm_submatrices(submatrices, norm)
Computes the max norm of the collection of 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.
Types to handle submatrices.
integer, parameter, public select_row_col
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
type(mp_comm_type), parameter, public mp_comm_null
type of a logger, at the moment it contains just a print level starting at which level it should be l...