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