(git:0de0cc2)
cp_dbcsr_operations.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 DBCSR operations in CP2K
10 !> \author Urban Borstnik
11 !> \date 2009-05-12
12 !> \version 0.8
13 !>
14 !> <b>Modification history:</b>
15 !> - Created 2009-05-12
16 !> - Generalized sm_fm_mulitply for matrices w/ different row/col block size (A. Bussy, 11.2018)
17 ! **************************************************************************************************
19  USE mathlib, ONLY: lcm, gcd
20  USE cp_blacs_env, ONLY: cp_blacs_env_type
21  USE cp_cfm_types, ONLY: cp_cfm_type
22  USE dbcsr_api, ONLY: dbcsr_distribution_get, &
23  dbcsr_convert_sizes_to_offsets, dbcsr_add, &
24  dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, &
25  dbcsr_deallocate_matrix, &
26  dbcsr_desymmetrize, dbcsr_distribution_new, &
27  dbcsr_get_data_type, dbcsr_get_info, dbcsr_get_matrix_type, &
28  dbcsr_iterator_type, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
29  dbcsr_iterator_start, dbcsr_iterator_stop, &
30  dbcsr_multiply, dbcsr_norm, dbcsr_p_type, dbcsr_release, &
31  dbcsr_reserve_all_blocks, dbcsr_scale, dbcsr_type, &
32  dbcsr_valid_index, dbcsr_verify_matrix, &
33  dbcsr_distribution_type, dbcsr_distribution_release, &
34  dbcsr_norm_frobenius, &
35  dbcsr_type_antisymmetric, dbcsr_type_complex_8, dbcsr_type_no_symmetry, dbcsr_type_real_8, &
36  dbcsr_type_symmetric
38  USE cp_fm_struct, ONLY: cp_fm_struct_create, &
40  cp_fm_struct_type
41  USE cp_fm_types, ONLY: cp_fm_create, &
43  cp_fm_release, &
44  cp_fm_to_fm, &
45  cp_fm_type
46  USE message_passing, ONLY: mp_para_env_type
48  distribution_2d_type
49  USE kinds, ONLY: dp, default_string_length
50  USE message_passing, ONLY: mp_comm_type
51 
52 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
53 #include "base/base_uses.f90"
54 
55  IMPLICIT NONE
56 
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_operations'
58  LOGICAL, PARAMETER :: debug_mod = .false.
59 
60  INTEGER, SAVE, PUBLIC :: max_elements_per_block = 32
61 
62  PUBLIC :: dbcsr_multiply_local
63 
64  ! CP2K API emulation
71 
72  ! distribution_2d_type compatibility
73  PUBLIC :: cp_dbcsr_dist2d_to_dist
74 
75  PUBLIC :: dbcsr_copy_columns_hack
76 
77  ! matrix set
80 
82  MODULE PROCEDURE allocate_dbcsr_matrix_set_1d
83  MODULE PROCEDURE allocate_dbcsr_matrix_set_2d
84  MODULE PROCEDURE allocate_dbcsr_matrix_set_3d
85  MODULE PROCEDURE allocate_dbcsr_matrix_set_4d
86  MODULE PROCEDURE allocate_dbcsr_matrix_set_5d
87  END INTERFACE
88 
90  MODULE PROCEDURE deallocate_dbcsr_matrix_set_1d
91  MODULE PROCEDURE deallocate_dbcsr_matrix_set_2d
92  MODULE PROCEDURE deallocate_dbcsr_matrix_set_3d
93  MODULE PROCEDURE deallocate_dbcsr_matrix_set_4d
94  MODULE PROCEDURE deallocate_dbcsr_matrix_set_5d
95  END INTERFACE
96 
97  PRIVATE
98 
99 CONTAINS
100 
101 
102 ! **************************************************************************************************
103 !> \brief Copy a BLACS matrix to a dbcsr matrix.
104 !>
105 !> real_matrix=beta*real_matrix+alpha*fm
106 !> beta defaults to 0, alpha to 1
107 !> \param[in] fm full matrix
108 !> \param[out] matrix DBCSR matrix
109 !> \param[in] keep_sparsity (optional) retains the sparsity of the input
110 !> matrix
111 !> \date 2009-10-13
112 !> \par History
113 !> 2009-10-13 rewritten based on copy_dbcsr_to_fm
114 !> \author Urban Borstnik
115 !> \version 2.0
116 ! **************************************************************************************************
117  SUBROUTINE copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
118  TYPE(cp_fm_type), INTENT(IN) :: fm
119  TYPE(dbcsr_type), INTENT(INOUT) :: matrix
120  LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
121 
122  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_fm_to_dbcsr'
123 
124  TYPE(dbcsr_type) :: bc_mat
125  INTEGER :: handle
126 
127  CALL timeset(routinen, handle)
128 
129  CALL copy_fm_to_dbcsr_bc(fm, bc_mat)
130  CALL dbcsr_complete_redistribute(bc_mat, matrix, keep_sparsity=keep_sparsity)
131  CALL dbcsr_release(bc_mat)
132 
133  CALL timestop(handle)
134  END SUBROUTINE copy_fm_to_dbcsr
135 
136 ! **************************************************************************************************
137 !> \brief Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution,
138 !> which requires no complete redistribution.
139 !> \param fm ...
140 !> \param bc_mat ...
141 ! **************************************************************************************************
142  SUBROUTINE copy_fm_to_dbcsr_bc(fm, bc_mat)
143  TYPE(cp_fm_type), INTENT(IN) :: fm
144  TYPE(dbcsr_type), INTENT(INOUT) :: bc_mat
145 
146  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_fm_to_dbcsr_bc'
147 
148  INTEGER :: col, handle, ncol_block, ncol_global, nrow_block, nrow_global, row
149  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
150  INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
151  REAL (kind=dp), DIMENSION(:, :), POINTER :: fm_block, dbcsr_block
152  TYPE(dbcsr_distribution_type) :: bc_dist
153  TYPE(dbcsr_iterator_type) :: iter
154  INTEGER, DIMENSION(:, :), POINTER :: pgrid
155 
156  CALL timeset(routinen, handle)
157 
158  IF (fm%use_sp) cpabort("copy_fm_to_dbcsr_bc: single precision not supported")
159 
160  ! Create processor grid
161  pgrid => fm%matrix_struct%context%blacs2mpi
162 
163  ! Create a block-cyclic distribution compatible with the FM matrix.
164  nrow_block = fm%matrix_struct%nrow_block
165  ncol_block = fm%matrix_struct%ncol_block
166  nrow_global = fm%matrix_struct%nrow_global
167  ncol_global = fm%matrix_struct%ncol_global
168  NULLIFY (col_blk_size, row_blk_size)
169  CALL dbcsr_create_dist_block_cyclic(bc_dist, &
170  nrows=nrow_global, ncolumns=ncol_global, & ! Actual full matrix size
171  nrow_block=nrow_block, ncol_block=ncol_block, & ! BLACS parameters
172  group_handle=fm%matrix_struct%para_env%get_handle(), pgrid=pgrid, &
173  row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size) ! block-cyclic row/col sizes
174 
175  ! Create the block-cyclic DBCSR matrix
176  CALL dbcsr_create(bc_mat, "Block-cyclic ", bc_dist, &
177  dbcsr_type_no_symmetry, row_blk_size, col_blk_size, nze=0, &
178  reuse_arrays=.true., data_type=dbcsr_type_real_8)
179  CALL dbcsr_distribution_release(bc_dist)
180 
181  ! allocate all blocks
182  CALL dbcsr_reserve_all_blocks(bc_mat)
183 
184  CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
185 
186  ! Copy the FM data to the block-cyclic DBCSR matrix. This step
187  ! could be skipped with appropriate DBCSR index manipulation.
188  fm_block => fm%local_data
189 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
190 !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
191  CALL dbcsr_iterator_start(iter, bc_mat)
192  DO WHILE (dbcsr_iterator_blocks_left(iter))
193  CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
194  dbcsr_block(:, :) = fm_block(first_row(row):last_row(row), first_col(col):last_col(col))
195  END DO
196  CALL dbcsr_iterator_stop(iter)
197 !$OMP END PARALLEL
198 
199  CALL timestop(handle)
200  END SUBROUTINE copy_fm_to_dbcsr_bc
201 
202 ! **************************************************************************************************
203 !> \brief Copy a DBCSR matrix to a BLACS matrix
204 !> \param[in] matrix DBCSR matrix
205 !> \param[out] fm full matrix
206 ! **************************************************************************************************
207  SUBROUTINE copy_dbcsr_to_fm (matrix, fm)
208  TYPE(dbcsr_type), INTENT(IN) :: matrix
209  TYPE(cp_fm_type), INTENT(IN) :: fm
210 
211  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_fm'
212 
213  INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
214  INTEGER :: handle, ncol_block, nfullcols_total, group_handle, &
215  nfullrows_total, nrow_block
216  TYPE(dbcsr_type) :: bc_mat, matrix_nosym
217  TYPE(dbcsr_distribution_type) :: dist, bc_dist
218  CHARACTER(len=default_string_length) :: name
219  INTEGER, DIMENSION(:, :), POINTER :: pgrid
220 
221  CALL timeset(routinen, handle)
222 
223  ! check compatibility
224  CALL dbcsr_get_info(matrix, &
225  name=name, &
226  distribution=dist, &
227  nfullrows_total=nfullrows_total, &
228  nfullcols_total=nfullcols_total)
229 
230  cpassert(fm%matrix_struct%nrow_global == nfullrows_total)
231  cpassert(fm%matrix_struct%ncol_global == nfullcols_total)
232 
233  ! info about the full matrix
234  nrow_block = fm%matrix_struct%nrow_block
235  ncol_block = fm%matrix_struct%ncol_block
236 
237  ! Convert DBCSR to a block-cyclic
238  NULLIFY (col_blk_size, row_blk_size)
239  CALL dbcsr_distribution_get(dist, group=group_handle, pgrid=pgrid)
240  CALL dbcsr_create_dist_block_cyclic(bc_dist, &
241  nrows=nfullrows_total, ncolumns=nfullcols_total, &
242  nrow_block=nrow_block, ncol_block=ncol_block, &
243  group_handle=group_handle, pgrid=pgrid, &
244  row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size)
245 
246  CALL dbcsr_create(bc_mat, "Block-cyclic"//name, bc_dist, &
247  dbcsr_type_no_symmetry, row_blk_size, col_blk_size, &
248  nze=0, data_type=dbcsr_get_data_type(matrix), &
249  reuse_arrays=.true.)
250  CALL dbcsr_distribution_release(bc_dist)
251 
252  CALL dbcsr_create(matrix_nosym, template=matrix, matrix_type="N")
253  CALL dbcsr_desymmetrize(matrix, matrix_nosym)
254  CALL dbcsr_complete_redistribute(matrix_nosym, bc_mat)
255  CALL dbcsr_release(matrix_nosym)
256 
257  CALL copy_dbcsr_to_fm_bc(bc_mat, fm)
258 
259  CALL dbcsr_release(bc_mat)
260 
261  CALL timestop(handle)
262  END SUBROUTINE copy_dbcsr_to_fm
263 
264 ! **************************************************************************************************
265 !> \brief Copy a DBCSR_BLACS matrix to a BLACS matrix
266 !> \param bc_mat DBCSR matrix
267 !> \param[out] fm full matrix
268 ! **************************************************************************************************
269  SUBROUTINE copy_dbcsr_to_fm_bc(bc_mat, fm)
270  TYPE(dbcsr_type), INTENT(IN) :: bc_mat
271  TYPE(cp_fm_type), INTENT(IN) :: fm
272 
273  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_fm_bc'
274 
275  INTEGER :: col, handle, row
276  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
277  REAL (kind=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
278  TYPE(dbcsr_iterator_type) :: iter
279 
280  CALL timeset(routinen, handle)
281 
282  IF (fm%use_sp) cpabort("copy_dbcsr_to_fm_bc: single precision not supported")
283 
284  CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
285 
286  ! Now copy data to the FM matrix
287  fm_block => fm%local_data
288  fm_block = real(0.0, kind=dp)
289 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
290 !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
291  CALL dbcsr_iterator_start(iter, bc_mat)
292  DO WHILE (dbcsr_iterator_blocks_left(iter))
293  CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
294  fm_block(first_row(row):last_row(row), first_col(col):last_col(col)) = dbcsr_block(:, :)
295  END DO
296  CALL dbcsr_iterator_stop(iter)
297 !$OMP END PARALLEL
298 
299  CALL timestop(handle)
300  END SUBROUTINE copy_dbcsr_to_fm_bc
301 
302 
303 ! **************************************************************************************************
304 !> \brief Copy a BLACS matrix to a dbcsr matrix.
305 !>
306 !> real_matrix=beta*real_matrix+alpha*fm
307 !> beta defaults to 0, alpha to 1
308 !> \param[in] fm full matrix
309 !> \param[out] matrix DBCSR matrix
310 !> \param[in] keep_sparsity (optional) retains the sparsity of the input
311 !> matrix
312 !> \date 2009-10-13
313 !> \par History
314 !> 2009-10-13 rewritten based on copy_dbcsr_to_fm
315 !> \author Urban Borstnik
316 !> \version 2.0
317 ! **************************************************************************************************
318  SUBROUTINE copy_cfm_to_dbcsr(fm, matrix, keep_sparsity)
319  TYPE(cp_cfm_type), INTENT(IN) :: fm
320  TYPE(dbcsr_type), INTENT(INOUT) :: matrix
321  LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
322 
323  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_cfm_to_dbcsr'
324 
325  TYPE(dbcsr_type) :: bc_mat
326  INTEGER :: handle
327 
328  CALL timeset(routinen, handle)
329 
330  CALL copy_cfm_to_dbcsr_bc(fm, bc_mat)
331  CALL dbcsr_complete_redistribute(bc_mat, matrix, keep_sparsity=keep_sparsity)
332  CALL dbcsr_release(bc_mat)
333 
334  CALL timestop(handle)
335  END SUBROUTINE copy_cfm_to_dbcsr
336 
337 ! **************************************************************************************************
338 !> \brief Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution,
339 !> which requires no complete redistribution.
340 !> \param fm ...
341 !> \param bc_mat ...
342 ! **************************************************************************************************
343  SUBROUTINE copy_cfm_to_dbcsr_bc(fm, bc_mat)
344  TYPE(cp_cfm_type), INTENT(IN) :: fm
345  TYPE(dbcsr_type), INTENT(INOUT) :: bc_mat
346 
347  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_cfm_to_dbcsr_bc'
348 
349  INTEGER :: col, handle, ncol_block, ncol_global, nrow_block, nrow_global, row
350  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
351  INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
352  COMPLEX (KIND=dp), DIMENSION(:, :), POINTER :: fm_block, dbcsr_block
353  TYPE(dbcsr_distribution_type) :: bc_dist
354  TYPE(dbcsr_iterator_type) :: iter
355  INTEGER, DIMENSION(:, :), POINTER :: pgrid
356 
357  CALL timeset(routinen, handle)
358 
359 
360  ! Create processor grid
361  pgrid => fm%matrix_struct%context%blacs2mpi
362 
363  ! Create a block-cyclic distribution compatible with the FM matrix.
364  nrow_block = fm%matrix_struct%nrow_block
365  ncol_block = fm%matrix_struct%ncol_block
366  nrow_global = fm%matrix_struct%nrow_global
367  ncol_global = fm%matrix_struct%ncol_global
368  NULLIFY (col_blk_size, row_blk_size)
369  CALL dbcsr_create_dist_block_cyclic(bc_dist, &
370  nrows=nrow_global, ncolumns=ncol_global, & ! Actual full matrix size
371  nrow_block=nrow_block, ncol_block=ncol_block, & ! BLACS parameters
372  group_handle=fm%matrix_struct%para_env%get_handle(), pgrid=pgrid, &
373  row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size) ! block-cyclic row/col sizes
374 
375  ! Create the block-cyclic DBCSR matrix
376  CALL dbcsr_create(bc_mat, "Block-cyclic ", bc_dist, &
377  dbcsr_type_no_symmetry, row_blk_size, col_blk_size, nze=0, &
378  reuse_arrays=.true., data_type=dbcsr_type_complex_8)
379  CALL dbcsr_distribution_release(bc_dist)
380 
381  ! allocate all blocks
382  CALL dbcsr_reserve_all_blocks(bc_mat)
383 
384  CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
385 
386  ! Copy the FM data to the block-cyclic DBCSR matrix. This step
387  ! could be skipped with appropriate DBCSR index manipulation.
388  fm_block => fm%local_data
389 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
390 !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
391  CALL dbcsr_iterator_start(iter, bc_mat)
392  DO WHILE (dbcsr_iterator_blocks_left(iter))
393  CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
394  dbcsr_block(:, :) = fm_block(first_row(row):last_row(row), first_col(col):last_col(col))
395  END DO
396  CALL dbcsr_iterator_stop(iter)
397 !$OMP END PARALLEL
398 
399  CALL timestop(handle)
400  END SUBROUTINE copy_cfm_to_dbcsr_bc
401 
402 ! **************************************************************************************************
403 !> \brief Copy a DBCSR matrix to a BLACS matrix
404 !> \param[in] matrix DBCSR matrix
405 !> \param[out] fm full matrix
406 ! **************************************************************************************************
407  SUBROUTINE copy_dbcsr_to_cfm (matrix, fm)
408  TYPE(dbcsr_type), INTENT(IN) :: matrix
409  TYPE(cp_cfm_type), INTENT(IN) :: fm
410 
411  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_cfm'
412 
413  INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
414  INTEGER :: handle, ncol_block, nfullcols_total, group_handle, &
415  nfullrows_total, nrow_block
416  TYPE(dbcsr_type) :: bc_mat, matrix_nosym
417  TYPE(dbcsr_distribution_type) :: dist, bc_dist
418  CHARACTER(len=default_string_length) :: name
419  INTEGER, DIMENSION(:, :), POINTER :: pgrid
420 
421  CALL timeset(routinen, handle)
422 
423  ! check compatibility
424  CALL dbcsr_get_info(matrix, &
425  name=name, &
426  distribution=dist, &
427  nfullrows_total=nfullrows_total, &
428  nfullcols_total=nfullcols_total)
429 
430  cpassert(fm%matrix_struct%nrow_global == nfullrows_total)
431  cpassert(fm%matrix_struct%ncol_global == nfullcols_total)
432 
433  ! info about the full matrix
434  nrow_block = fm%matrix_struct%nrow_block
435  ncol_block = fm%matrix_struct%ncol_block
436 
437  ! Convert DBCSR to a block-cyclic
438  NULLIFY (col_blk_size, row_blk_size)
439  CALL dbcsr_distribution_get(dist, group=group_handle, pgrid=pgrid)
440  CALL dbcsr_create_dist_block_cyclic(bc_dist, &
441  nrows=nfullrows_total, ncolumns=nfullcols_total, &
442  nrow_block=nrow_block, ncol_block=ncol_block, &
443  group_handle=group_handle, pgrid=pgrid, &
444  row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size)
445 
446  CALL dbcsr_create(bc_mat, "Block-cyclic"//name, bc_dist, &
447  dbcsr_type_no_symmetry, row_blk_size, col_blk_size, &
448  nze=0, data_type=dbcsr_get_data_type(matrix), &
449  reuse_arrays=.true.)
450  CALL dbcsr_distribution_release(bc_dist)
451 
452  CALL dbcsr_create(matrix_nosym, template=matrix, matrix_type="N")
453  CALL dbcsr_desymmetrize(matrix, matrix_nosym)
454  CALL dbcsr_complete_redistribute(matrix_nosym, bc_mat)
455  CALL dbcsr_release(matrix_nosym)
456 
457  CALL copy_dbcsr_to_cfm_bc(bc_mat, fm)
458 
459  CALL dbcsr_release(bc_mat)
460 
461  CALL timestop(handle)
462  END SUBROUTINE copy_dbcsr_to_cfm
463 
464 ! **************************************************************************************************
465 !> \brief Copy a DBCSR_BLACS matrix to a BLACS matrix
466 !> \param bc_mat DBCSR matrix
467 !> \param[out] fm full matrix
468 ! **************************************************************************************************
469  SUBROUTINE copy_dbcsr_to_cfm_bc(bc_mat, fm)
470  TYPE(dbcsr_type), INTENT(IN) :: bc_mat
471  TYPE(cp_cfm_type), INTENT(IN) :: fm
472 
473  CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_cfm_bc'
474 
475  INTEGER :: col, handle, row
476  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
477  COMPLEX (KIND=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
478  TYPE(dbcsr_iterator_type) :: iter
479 
480  CALL timeset(routinen, handle)
481 
482 
483  CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
484 
485  ! Now copy data to the FM matrix
486  fm_block => fm%local_data
487  fm_block = cmplx(0.0, kind=dp)
488 !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
489 !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
490  CALL dbcsr_iterator_start(iter, bc_mat)
491  DO WHILE (dbcsr_iterator_blocks_left(iter))
492  CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
493  fm_block(first_row(row):last_row(row), first_col(col):last_col(col)) = dbcsr_block(:, :)
494  END DO
495  CALL dbcsr_iterator_stop(iter)
496 !$OMP END PARALLEL
497 
498  CALL timestop(handle)
499  END SUBROUTINE copy_dbcsr_to_cfm_bc
500 
501 
502 ! **************************************************************************************************
503 !> \brief Helper routine used to copy blocks from DBCSR into FM matrices and vice versa
504 !> \param bc_mat ...
505 !> \param first_row ...
506 !> \param last_row ...
507 !> \param first_col ...
508 !> \param last_col ...
509 !> \author Ole Schuett
510 ! **************************************************************************************************
511  SUBROUTINE calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
512  TYPE(dbcsr_type), INTENT(IN) :: bc_mat
513  INTEGER :: col, nblkcols_local, nblkcols_total, nblkrows_local, nblkrows_total, row
514  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, &
515  last_row, local_col_sizes, &
516  local_row_sizes
517  INTEGER, DIMENSION(:), POINTER :: col_blk_size, local_cols, local_rows, &
518  row_blk_size
519 
520  CALL dbcsr_get_info(bc_mat, &
521  nblkrows_total=nblkrows_total, &
522  nblkcols_total=nblkcols_total, &
523  nblkrows_local=nblkrows_local, &
524  nblkcols_local=nblkcols_local, &
525  local_rows=local_rows, &
526  local_cols=local_cols, &
527  row_blk_size=row_blk_size, &
528  col_blk_size=col_blk_size)
529 
530  ! calculate first_row and last_row
531  ALLOCATE (local_row_sizes(nblkrows_total))
532  local_row_sizes(:) = 0
533  IF (nblkrows_local .GE. 1) THEN
534  DO row = 1, nblkrows_local
535  local_row_sizes(local_rows(row)) = row_blk_size(local_rows(row))
536  END DO
537  END IF
538  ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total))
539  CALL dbcsr_convert_sizes_to_offsets(local_row_sizes, first_row, last_row)
540 
541  ! calculate first_col and last_col
542  ALLOCATE (local_col_sizes(nblkcols_total))
543  local_col_sizes(:) = 0
544  IF (nblkcols_local .GE. 1) THEN
545  DO col = 1, nblkcols_local
546  local_col_sizes(local_cols(col)) = col_blk_size(local_cols(col))
547  END DO
548  END IF
549  ALLOCATE (first_col(nblkcols_total), last_col(nblkcols_total))
550  CALL dbcsr_convert_sizes_to_offsets(local_col_sizes, first_col, last_col)
551 
552  END SUBROUTINE calculate_fm_block_ranges
553 
554 ! **************************************************************************************************
555 !> \brief hack for dbcsr_copy_columns
556 !> \param matrix_b ...
557 !> \param matrix_a ...
558 !> \param ncol ...
559 !> \param source_start ...
560 !> \param target_start ...
561 !> \param para_env ...
562 !> \param blacs_env ...
563 !> \author vw
564 ! **************************************************************************************************
565  SUBROUTINE dbcsr_copy_columns_hack(matrix_b, matrix_a, &
566  ncol, source_start, target_start, para_env, blacs_env)
567 
568  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b
569  TYPE(dbcsr_type), INTENT(IN) :: matrix_a
570  INTEGER, INTENT(IN) :: ncol, source_start, target_start
571  TYPE(mp_para_env_type), POINTER :: para_env
572  TYPE(cp_blacs_env_type), POINTER :: blacs_env
573 
574  INTEGER :: nfullcols_total, nfullrows_total
575  TYPE(cp_fm_struct_type), POINTER :: fm_struct
576  TYPE(cp_fm_type) :: fm_matrix_a, fm_matrix_b
577 
578  NULLIFY (fm_struct)
579  CALL dbcsr_get_info(matrix_a, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
580  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
581  ncol_global=nfullcols_total, para_env=para_env)
582  CALL cp_fm_create(fm_matrix_a, fm_struct, name="fm_matrix_a")
583  CALL cp_fm_struct_release(fm_struct)
584 
585  CALL dbcsr_get_info(matrix_b, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
586  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
587  ncol_global=nfullcols_total, para_env=para_env)
588  CALL cp_fm_create(fm_matrix_b, fm_struct, name="fm_matrix_b")
589  CALL cp_fm_struct_release(fm_struct)
590 
591  CALL copy_dbcsr_to_fm(matrix_a, fm_matrix_a)
592  CALL copy_dbcsr_to_fm(matrix_b, fm_matrix_b)
593 
594  CALL cp_fm_to_fm(fm_matrix_a, fm_matrix_b, ncol, source_start, target_start)
595 
596  CALL copy_fm_to_dbcsr(fm_matrix_b, matrix_b)
597 
598  CALL cp_fm_release(fm_matrix_a)
599  CALL cp_fm_release(fm_matrix_b)
600 
601  END SUBROUTINE dbcsr_copy_columns_hack
602 
603 ! **************************************************************************************************
604 !> \brief Creates a DBCSR distribution from a distribution_2d
605 !> \param[in] dist2d distribution_2d
606 !> \param[out] dist DBCSR distribution
607 !> \par History
608 !> move form dbcsr_operation 01.2010
609 ! **************************************************************************************************
610  SUBROUTINE cp_dbcsr_dist2d_to_dist(dist2d, dist)
611  TYPE(distribution_2d_type), INTENT(IN), TARGET :: dist2d
612  TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist
613 
614  INTEGER, DIMENSION(:, :), POINTER :: pgrid, col_dist_data, row_dist_data
615  TYPE(cp_blacs_env_type), POINTER :: blacs_env
616  TYPE(mp_para_env_type), POINTER :: para_env
617  TYPE(distribution_2d_type), POINTER :: dist2d_p
618  INTEGER, DIMENSION(:), POINTER :: row_dist, col_dist
619 
620  dist2d_p => dist2d
621  CALL distribution_2d_get(dist2d_p, &
622  row_distribution=row_dist_data, &
623  col_distribution=col_dist_data, &
624  blacs_env=blacs_env)
625  CALL blacs_env%get(para_env=para_env, blacs2mpi=pgrid)
626 
627  ! map to 1D arrays
628  row_dist => row_dist_data(:, 1)
629  col_dist => col_dist_data(:, 1)
630  !row_cluster => row_dist_data(:, 2)
631  !col_cluster => col_dist_data(:, 2)
632 
633  CALL dbcsr_distribution_new(dist, &
634  group=para_env%get_handle(), pgrid=pgrid, &
635  row_dist=row_dist, &
636  col_dist=col_dist)
637 
638  END SUBROUTINE cp_dbcsr_dist2d_to_dist
639 
640 ! **************************************************************************************************
641 !> \brief multiply a dbcsr with a replicated array
642 !> c = alpha_scalar * A (dbscr) * b + c
643 !> \param[in] matrix_a DBSCR matrxx
644 !> \param[in] vec_b vectors b
645 !> \param[inout] vec_c vectors c
646 !> \param[in] ncol nbr of columns
647 !> \param[in] alpha alpha
648 !>
649 ! **************************************************************************************************
650  SUBROUTINE dbcsr_multiply_local(matrix_a, vec_b, vec_c, ncol, alpha)
651  TYPE(dbcsr_type), INTENT(IN) :: matrix_a
652  REAL(dp), DIMENSION(:, :), INTENT(IN) :: vec_b
653  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vec_c
654  INTEGER, INTENT(in), OPTIONAL :: ncol
655  REAL(dp), INTENT(IN), OPTIONAL :: alpha
656 
657  CHARACTER(LEN=*), PARAMETER :: routinen = 'dbcsr_multiply_local'
658 
659  INTEGER :: blk, col, coloff, my_ncol, row, rowoff, &
660  timing_handle
661  LOGICAL :: has_symm
662  REAL(dp) :: my_alpha, my_alpha2
663  REAL(dp), DIMENSION(:, :), POINTER :: data_d
664  TYPE(dbcsr_iterator_type) :: iter
665 
666  CALL timeset(routinen, timing_handle)
667 
668  my_alpha = 1.0_dp
669  IF (PRESENT(alpha)) my_alpha = alpha
670 
671  my_ncol = SIZE(vec_b, 2)
672  IF (PRESENT(ncol)) my_ncol = ncol
673 
674  my_alpha2 = 0.0_dp
675  IF (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_symmetric) my_alpha2 = my_alpha
676  IF (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_antisymmetric) my_alpha2 = -my_alpha
677 
678  has_symm = (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_symmetric .OR. &
679  dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_antisymmetric)
680 
681 !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix_a,vec_b,vec_c,ncol,my_alpha2,my_alpha,my_ncol,has_symm) &
682 !$OMP PRIVATE(iter,row,col,data_d,blk,rowoff,coloff)
683  CALL dbcsr_iterator_start(iter, matrix_a, read_only=.true., dynamic=.true., dynamic_byrows=.true.)
684  DO WHILE (dbcsr_iterator_blocks_left(iter))
685  CALL dbcsr_iterator_next_block(iter, row, col, data_d, blk, row_offset=rowoff, col_offset=coloff)
686  IF (my_ncol .NE. 1) THEN
687  CALL dgemm('N', 'N', &
688  SIZE(data_d, 1), my_ncol, SIZE(data_d, 2), &
689  my_alpha, data_d(1, 1), SIZE(data_d, 1), &
690  vec_b(coloff, 1), SIZE(vec_b, 1), &
691  1.0_dp, vec_c(rowoff, 1), SIZE(vec_c, 1))
692  ELSE
693  CALL dgemv('N', SIZE(data_d, 1), SIZE(data_d, 2), &
694  my_alpha, data_d(1, 1), SIZE(data_d, 1), &
695  vec_b(coloff, 1), 1, &
696  1.0_dp, vec_c(rowoff, 1), 1)
697  END IF
698  END DO
699  CALL dbcsr_iterator_stop(iter)
700 !$OMP END PARALLEL
701 
702  ! FIXME ... in the symmetric case, the writes to vec_c depend on the column, not the row. This makes OMP-ing more difficult
703  ! needs e.g. a buffer for vec_c and a reduction of that buffer.
704  IF (has_symm) THEN
705  CALL dbcsr_iterator_start(iter, matrix_a)
706  DO WHILE (dbcsr_iterator_blocks_left(iter))
707  CALL dbcsr_iterator_next_block(iter, row, col, data_d, blk, row_offset=rowoff, col_offset=coloff)
708  IF (row .NE. col) THEN
709  IF (my_ncol .NE. 1) THEN
710  CALL dgemm('T', 'N', &
711  SIZE(data_d, 2), my_ncol, SIZE(data_d, 1), &
712  my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
713  vec_b(rowoff, 1), SIZE(vec_b, 1), &
714  1.0_dp, vec_c(coloff, 1), SIZE(vec_c, 1))
715  ELSE
716  CALL dgemv('T', SIZE(data_d, 1), SIZE(data_d, 2), &
717  my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
718  vec_b(rowoff, 1), 1, &
719  1.0_dp, vec_c(coloff, 1), 1)
720  END IF
721  END IF
722  END DO
723  CALL dbcsr_iterator_stop(iter)
724  END IF
725 
726  CALL timestop(timing_handle)
727  END SUBROUTINE dbcsr_multiply_local
728 
729 ! **************************************************************************************************
730 !> \brief multiply a dbcsr with a fm matrix
731 !>
732 !> For backwards compatibility with BLAS XGEMM, this routine supports
733 !> the multiplication of matrices with incompatible dimensions.
734 !>
735 !> \param[in] matrix DBCSR matrix
736 !> \param fm_in full matrix
737 !> \param fm_out full matrix
738 !> \param[in] ncol nbr of columns
739 !> \param[in] alpha alpha
740 !> \param[in] beta beta
741 !>
742 ! **************************************************************************************************
743  SUBROUTINE cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
744  TYPE(dbcsr_type), INTENT(IN) :: matrix
745  TYPE(cp_fm_type), INTENT(IN) :: fm_in, fm_out
746  INTEGER, INTENT(IN) :: ncol
747  REAL(dp), INTENT(IN), OPTIONAL :: alpha, beta
748 
749  CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_dbcsr_sm_fm_multiply'
750 
751  INTEGER :: k_in, k_out, timing_handle, &
752  timing_handle_mult, &
753  a_ncol, a_nrow, b_ncol, b_nrow, c_ncol, c_nrow
754  INTEGER, DIMENSION(:), POINTER :: col_blk_size_right_in, &
755  col_blk_size_right_out, row_blk_size, &
756  !row_cluster, col_cluster,&
757  row_dist, col_dist, col_blk_size
758  TYPE(dbcsr_type) :: in, out
759  REAL(dp) :: my_alpha, my_beta
760  TYPE(dbcsr_distribution_type) :: dist, dist_right_in, product_dist
761 
762  CALL timeset(routinen, timing_handle)
763 
764  my_alpha = 1.0_dp
765  my_beta = 0.0_dp
766  IF (PRESENT(alpha)) my_alpha = alpha
767  IF (PRESENT(beta)) my_beta = beta
768 
769  ! TODO
770  CALL cp_fm_get_info(fm_in, ncol_global=b_ncol, nrow_global=b_nrow)
771  CALL cp_fm_get_info(fm_out, ncol_global=c_ncol, nrow_global=c_nrow)
772  CALL dbcsr_get_info(matrix, nfullrows_total=a_nrow, nfullcols_total=a_ncol)
773  !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: A ", a_nrow, "x", a_ncol
774  !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: B ", b_nrow, "x", b_ncol
775  !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: C ", c_nrow, "x", c_ncol
776 
777  CALL cp_fm_get_info(fm_out, ncol_global=k_out)
778 
779  CALL cp_fm_get_info(fm_in, ncol_global=k_in)
780  !write(*,*)routineN//" -----------------------------------"
781  !IF (k_in .NE. k_out) &
782  ! WRITE(*,'(3(A,I5,1X),2(A,F5.2,1X))')&
783  ! routineN//" ncol", ncol,'k_in',k_in,'k_out',k_out,&
784  ! 'alpha',my_alpha,'beta',my_beta
785 
786  IF (ncol .GT. 0 .AND. k_out .GT. 0 .AND. k_in .GT. 0) THEN
787  CALL dbcsr_get_info(matrix, row_blk_size=row_blk_size, col_blk_size=col_blk_size, distribution=dist)
788  CALL dbcsr_create_dist_r_unrot(dist_right_in, dist, k_in, col_blk_size_right_in)
789 
790  CALL dbcsr_create(in, "D", dist_right_in, dbcsr_type_no_symmetry, &
791  col_blk_size, col_blk_size_right_in, nze=0)
792 
793  CALL dbcsr_distribution_get(dist, row_dist=row_dist)
794  CALL dbcsr_distribution_get(dist_right_in, col_dist=col_dist)
795  CALL dbcsr_distribution_new(product_dist, template=dist, &
796  row_dist=row_dist, col_dist=col_dist)
797  ALLOCATE (col_blk_size_right_out(SIZE(col_blk_size_right_in)))
798  col_blk_size_right_out = col_blk_size_right_in
799  CALL match_col_sizes(col_blk_size_right_out, col_blk_size_right_in, k_out)
800 
801  !if (k_in .ne. k_out) then
802  ! write(*,*)routineN//" in cs", col_blk_size_right_in
803  ! write(*,*)routineN//" out cs", col_blk_size_right_out
804  !endif
805 
806  CALL dbcsr_create(out, "D", product_dist, dbcsr_type_no_symmetry, &
807  row_blk_size, col_blk_size_right_out, nze=0)
808 
809  CALL copy_fm_to_dbcsr(fm_in, in)
810  IF (ncol .NE. k_out .OR. my_beta .NE. 0.0_dp) &
811  CALL copy_fm_to_dbcsr(fm_out, out)
812 
813  CALL timeset(routinen//'_core', timing_handle_mult)
814  CALL dbcsr_multiply("N", "N", my_alpha, matrix, in, my_beta, out, &
815  last_column=ncol)
816  CALL timestop(timing_handle_mult)
817 
818  CALL copy_dbcsr_to_fm(out, fm_out)
819 
820  CALL dbcsr_release(in)
821  CALL dbcsr_release(out)
822  DEALLOCATE (col_blk_size_right_in, col_blk_size_right_out)
823  CALL dbcsr_distribution_release(dist_right_in)
824  CALL dbcsr_distribution_release(product_dist)
825 
826  END IF
827 
828  CALL timestop(timing_handle)
829 
830  END SUBROUTINE cp_dbcsr_sm_fm_multiply
831 
832 ! **************************************************************************************************
833 !> \brief ...
834 !> \param sizes1 ...
835 !> \param sizes2 ...
836 !> \param full_num ...
837 ! **************************************************************************************************
838  SUBROUTINE match_col_sizes(sizes1, sizes2, full_num)
839  INTEGER, DIMENSION(:), INTENT(INOUT) :: sizes1
840  INTEGER, DIMENSION(:), INTENT(IN) :: sizes2
841  INTEGER, INTENT(IN) :: full_num
842 
843  INTEGER :: left, n1, n2, p, rm, used
844 
845  n1 = SIZE(sizes1)
846  n2 = SIZE(sizes2)
847  IF (n1 .NE. n2) &
848  cpabort("distributions must be equal!")
849  sizes1(1:n1) = sizes2(1:n1)
850  used = sum(sizes1(1:n1))
851  ! If sizes1 does not cover everything, then we increase the
852  ! size of the last block; otherwise we reduce the blocks
853  ! (from the end) until it is small enough.
854  IF (used .LT. full_num) THEN
855  sizes1(n1) = sizes1(n1) + full_num - used
856  ELSE
857  left = used - full_num
858  p = n1
859  DO WHILE (left .GT. 0 .AND. p .GT. 0)
860  rm = min(left, sizes1(p))
861  sizes1(p) = sizes1(p) - rm
862  left = left - rm
863  p = p - 1
864  END DO
865  END IF
866  END SUBROUTINE match_col_sizes
867 
868 ! **************************************************************************************************
869 !> \brief performs the multiplication sparse_matrix+dense_mat*dens_mat^T
870 !> if matrix_g is not explicitly given, matrix_v^T will be used
871 !> this can be important to save the necessary redistribute for a
872 !> different matrix_g and increase performance.
873 !> \param sparse_matrix ...
874 !> \param matrix_v ...
875 !> \param matrix_g ...
876 !> \param ncol ...
877 !> \param alpha ...
878 !> \param keep_sparsity Determines if the sparsity of sparse_matrix is retained
879 !> by default it is TRUE
880 !> \param symmetry_mode There are the following modes
881 !> 1: sparse_matrix += 0.5*alpha*(v*g^T+g^T*v) (symmetric update)
882 !> -1: sparse_matrix += 0.5*alpha*(v*g^T-g^T*v) (skewsymmetric update)
883 !> else: sparse_matrix += alpha*v*g^T (no symmetry, default)
884 !> saves some redistribution steps
885 ! **************************************************************************************************
886  SUBROUTINE cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
887  TYPE(dbcsr_type), INTENT(INOUT) :: sparse_matrix
888  TYPE(cp_fm_type), INTENT(IN) :: matrix_v
889  TYPE(cp_fm_type), OPTIONAL, INTENT(IN) :: matrix_g
890  INTEGER, INTENT(IN) :: ncol
891  REAL(kind=dp), INTENT(IN), OPTIONAL :: alpha
892  LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
893  INTEGER, INTENT(IN), OPTIONAL :: symmetry_mode
894 
895  CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_dbcsr_plus_fm_fm_t_native'
896 
897  INTEGER :: npcols, k, nao, timing_handle, data_type, my_symmetry_mode
898  INTEGER, DIMENSION(:), POINTER :: col_blk_size_left, &
899  col_dist_left, row_blk_size, row_dist
900  LOGICAL :: check_product, my_keep_sparsity
901  REAL(kind=dp) :: my_alpha, norm
902  TYPE(dbcsr_type) :: mat_g, mat_v, sparse_matrix2, &
903  sparse_matrix3
904  TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
905  TYPE(cp_fm_type) :: fm_matrix
906  TYPE(dbcsr_distribution_type) :: dist_left, sparse_dist
907 
908  check_product = .false.
909 
910  CALL timeset(routinen, timing_handle)
911 
912  my_keep_sparsity = .true.
913  IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
914 
915  my_symmetry_mode = 0
916  IF (PRESENT(symmetry_mode)) my_symmetry_mode = symmetry_mode
917 
918  NULLIFY (col_dist_left)
919 
920  IF (ncol .GT. 0) THEN
921  IF (.NOT. dbcsr_valid_index(sparse_matrix)) &
922  cpabort("sparse_matrix must pre-exist")
923  !
924  ! Setup matrix_v
925  CALL cp_fm_get_info(matrix_v, ncol_global=k)
926  !WRITE(*,*)routineN//'truncated mult k, ncol',k,ncol,' PRESENT (matrix_g)',PRESENT (matrix_g)
927  CALL dbcsr_get_info(sparse_matrix, distribution=sparse_dist)
928  CALL dbcsr_distribution_get(sparse_dist, npcols=npcols, row_dist=row_dist)
929  CALL create_bl_distribution(col_dist_left, col_blk_size_left, k, npcols)
930  CALL dbcsr_distribution_new(dist_left, template=sparse_dist, &
931  row_dist=row_dist, col_dist=col_dist_left)
932  DEALLOCATE (col_dist_left)
933  CALL dbcsr_get_info(sparse_matrix, row_blk_size=row_blk_size, data_type=data_type)
934  CALL dbcsr_create(mat_v, "DBCSR matrix_v", dist_left, dbcsr_type_no_symmetry, &
935  row_blk_size, col_blk_size_left, nze=0, data_type=data_type)
936  CALL copy_fm_to_dbcsr(matrix_v, mat_v)
937  CALL dbcsr_verify_matrix(mat_v)
938  !
939  ! Setup matrix_g
940  IF (PRESENT(matrix_g)) THEN
941  CALL dbcsr_create(mat_g, "DBCSR matrix_g", dist_left, dbcsr_type_no_symmetry, &
942  row_blk_size, col_blk_size_left, data_type=data_type)
943  CALL copy_fm_to_dbcsr(matrix_g, mat_g)
944  END IF
945  !
946  DEALLOCATE (col_blk_size_left)
947  CALL dbcsr_distribution_release(dist_left)
948  !
949  !
950  IF (check_product) THEN
951  CALL cp_fm_get_info(matrix_v, nrow_global=nao)
952  CALL cp_fm_struct_create(fm_struct_tmp, context=matrix_v%matrix_struct%context, nrow_global=nao, &
953  ncol_global=nao, para_env=matrix_v%matrix_struct%para_env)
954  CALL cp_fm_create(fm_matrix, fm_struct_tmp, name="fm matrix")
955  CALL cp_fm_struct_release(fm_struct_tmp)
956  CALL copy_dbcsr_to_fm(sparse_matrix, fm_matrix)
957  CALL dbcsr_copy(sparse_matrix3, sparse_matrix)
958  END IF
959  !
960  my_alpha = 1.0_dp
961  IF (PRESENT(alpha)) my_alpha = alpha
962  IF (PRESENT(matrix_g)) THEN
963  IF (my_symmetry_mode == 1) THEN
964  ! Symmetric mode
965  CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
966  1.0_dp, sparse_matrix, &
967  retain_sparsity=my_keep_sparsity, &
968  last_k=ncol)
969  CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_g, mat_v, &
970  1.0_dp, sparse_matrix, &
971  retain_sparsity=my_keep_sparsity, &
972  last_k=ncol)
973  ELSE IF (my_symmetry_mode == -1) THEN
974  ! Skewsymmetric mode
975  CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
976  1.0_dp, sparse_matrix, &
977  retain_sparsity=my_keep_sparsity, &
978  last_k=ncol)
979  CALL dbcsr_multiply("N", "T", -0.5_dp*my_alpha, mat_g, mat_v, &
980  1.0_dp, sparse_matrix, &
981  retain_sparsity=my_keep_sparsity, &
982  last_k=ncol)
983  ELSE
984  ! Normal mode
985  CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_g, &
986  1.0_dp, sparse_matrix, &
987  retain_sparsity=my_keep_sparsity, &
988  last_k=ncol)
989  END IF
990  ELSE
991  CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_v, &
992  1.0_dp, sparse_matrix, &
993  retain_sparsity=my_keep_sparsity, &
994  last_k=ncol)
995  END IF
996 
997  IF (check_product) THEN
998  IF (PRESENT(matrix_g)) THEN
999  IF (my_symmetry_mode == 1) THEN
1000  CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
1001  1.0_dp, fm_matrix)
1002  CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_g, matrix_v, &
1003  1.0_dp, fm_matrix)
1004  ELSE IF (my_symmetry_mode == -1) THEN
1005  CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
1006  1.0_dp, fm_matrix)
1007  CALL cp_fm_gemm("N", "T", nao, nao, ncol, -0.5_dp*my_alpha, matrix_g, matrix_v, &
1008  1.0_dp, fm_matrix)
1009  ELSE
1010  CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_g, &
1011  1.0_dp, fm_matrix)
1012  END IF
1013  ELSE
1014  CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_v, &
1015  1.0_dp, fm_matrix)
1016  END IF
1017 
1018  CALL dbcsr_copy(sparse_matrix2, sparse_matrix)
1019  CALL dbcsr_scale(sparse_matrix2, alpha_scalar=0.0_dp)
1020  CALL copy_fm_to_dbcsr(fm_matrix, sparse_matrix2, keep_sparsity=my_keep_sparsity)
1021  CALL dbcsr_add(sparse_matrix2, sparse_matrix, alpha_scalar=1.0_dp, &
1022  beta_scalar=-1.0_dp)
1023  CALL dbcsr_norm(sparse_matrix2, which_norm=dbcsr_norm_frobenius, &
1024  norm_scalar=norm)
1025  WRITE (*, *) 'nao=', nao, ' k=', k, ' ncol=', ncol, ' my_alpha=', my_alpha
1026  WRITE (*, *) 'PRESENT (matrix_g)', PRESENT(matrix_g)
1027  WRITE (*, *) 'matrix_type=', dbcsr_get_matrix_type(sparse_matrix)
1028  WRITE (*, *) 'norm(sm+alpha*v*g^t - fm+alpha*v*g^t)/n=', norm/real(nao, dp)
1029  IF (norm/real(nao, dp) .GT. 1e-12_dp) THEN
1030  !WRITE(*,*) 'fm_matrix'
1031  !DO j=1,SIZE(fm_matrix%local_data,2)
1032  ! DO i=1,SIZE(fm_matrix%local_data,1)
1033  ! WRITE(*,'(A,I3,A,I3,A,E26.16,A)') 'a(',i,',',j,')=',fm_matrix%local_data(i,j),';'
1034  ! ENDDO
1035  !ENDDO
1036  !WRITE(*,*) 'mat_v'
1037  !CALL dbcsr_print(mat_v,matlab_format=.TRUE.)
1038  !WRITE(*,*) 'mat_g'
1039  !CALL dbcsr_print(mat_g,matlab_format=.TRUE.)
1040  !WRITE(*,*) 'sparse_matrix'
1041  !CALL dbcsr_print(sparse_matrix,matlab_format=.TRUE.)
1042  !WRITE(*,*) 'sparse_matrix2 (-sm + sparse(fm))'
1043  !CALL dbcsr_print(sparse_matrix2,matlab_format=.TRUE.)
1044  !WRITE(*,*) 'sparse_matrix3 (copy of sm input)'
1045  !CALL dbcsr_print(sparse_matrix3,matlab_format=.TRUE.)
1046  !stop
1047  END IF
1048  CALL dbcsr_release(sparse_matrix2)
1049  CALL dbcsr_release(sparse_matrix3)
1050  CALL cp_fm_release(fm_matrix)
1051  END IF
1052  CALL dbcsr_release(mat_v)
1053  IF (PRESENT(matrix_g)) CALL dbcsr_release(mat_g)
1054  END IF
1055  CALL timestop(timing_handle)
1056 
1057  END SUBROUTINE cp_dbcsr_plus_fm_fm_t
1058 
1059 ! **************************************************************************************************
1060 !> \brief Utility function to copy a specially shaped fm to dbcsr_matrix
1061 !> The result matrix will be the matrix in dbcsr format
1062 !> with the row blocks sizes according to the block_sizes of the template
1063 !> and the col blocks sizes evenly blocked with the internal dbcsr conversion
1064 !> size (32 is the current default)
1065 !> \param matrix ...
1066 !> \param fm_in ...
1067 !> \param template ...
1068 ! **************************************************************************************************
1069  SUBROUTINE cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
1070  TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1071  TYPE(cp_fm_type), INTENT(IN) :: fm_in
1072  TYPE(dbcsr_type), INTENT(IN) :: template
1073 
1074  INTEGER :: k_in, data_type
1075  INTEGER, DIMENSION(:), POINTER :: col_blk_size_right_in, row_blk_size
1076  TYPE(dbcsr_distribution_type) :: tmpl_dist, dist_right_in
1077 
1078  CALL cp_fm_get_info(fm_in, ncol_global=k_in)
1079 
1080  CALL dbcsr_get_info(template, distribution=tmpl_dist)
1081  CALL dbcsr_create_dist_r_unrot(dist_right_in, tmpl_dist, k_in, col_blk_size_right_in)
1082  CALL dbcsr_get_info(template, row_blk_size=row_blk_size, data_type=data_type)
1083  CALL dbcsr_create(matrix, "D", dist_right_in, dbcsr_type_no_symmetry, &
1084  row_blk_size, col_blk_size_right_in, nze=0, data_type=data_type)
1085 
1086  CALL copy_fm_to_dbcsr(fm_in, matrix)
1087  DEALLOCATE (col_blk_size_right_in)
1088  CALL dbcsr_distribution_release(dist_right_in)
1089 
1090  END SUBROUTINE cp_fm_to_dbcsr_row_template
1091 
1092 ! **************************************************************************************************
1093 !> \brief Utility function to create an arbitrary shaped dbcsr matrix
1094 !> with the same processor grid as the template matrix
1095 !> both row sizes and col sizes are evenly blocked with the internal
1096 !> dbcsr_conversion size (32 is the current default)
1097 !> \param matrix dbcsr matrix to be created
1098 !> \param template template dbcsr matrix giving its mp_env
1099 !> \param m global row size of output matrix
1100 !> \param n global col size of output matrix
1101 !> \param sym ...
1102 !> \param data_type ...
1103 ! **************************************************************************************************
1104  SUBROUTINE cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
1105  TYPE(dbcsr_type), INTENT(INOUT) :: matrix, template
1106  INTEGER, INTENT(IN) :: m, n
1107  CHARACTER, OPTIONAL, INTENT(IN) :: sym
1108  INTEGER, OPTIONAL, INTENT(IN) :: data_type
1109 
1110  CHARACTER :: mysym
1111  INTEGER :: my_data_type, nprows, npcols
1112  INTEGER, DIMENSION(:), POINTER :: col_blk_size, &
1113  col_dist, row_blk_size, &
1114  row_dist
1115  TYPE(dbcsr_distribution_type) :: tmpl_dist, dist_m_n
1116 
1117  CALL dbcsr_get_info(template, &
1118  matrix_type=mysym, &
1119  data_type=my_data_type, &
1120  distribution=tmpl_dist)
1121 
1122  IF (PRESENT(sym)) mysym = sym
1123  IF (PRESENT(data_type)) my_data_type = data_type
1124 
1125  NULLIFY (row_dist, col_dist)
1126  NULLIFY (row_blk_size, col_blk_size)
1127  !NULLIFY (row_cluster, col_cluster)
1128 
1129  CALL dbcsr_distribution_get(tmpl_dist, nprows=nprows, npcols=npcols)
1130  CALL create_bl_distribution(row_dist, row_blk_size, m, nprows)
1131  CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
1132  CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
1133  row_dist=row_dist, col_dist=col_dist, &
1134  !row_cluster=row_cluster, col_cluster=col_cluster, &
1135  reuse_arrays=.true.)
1136 
1137  CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, &
1138  row_blk_size, col_blk_size, nze=0, data_type=my_data_type, &
1139  reuse_arrays=.true.)
1140  CALL dbcsr_distribution_release(dist_m_n)
1141 
1142  END SUBROUTINE cp_dbcsr_m_by_n_from_template
1143 
1144 ! **************************************************************************************************
1145 !> \brief Utility function to create dbcsr matrix, m x n matrix (n arbitrary)
1146 !> with the same processor grid and row distribution as the template matrix
1147 !> col sizes are evenly blocked with the internal
1148 !> dbcsr_conversion size (32 is the current default)
1149 !> \param matrix dbcsr matrix to be created
1150 !> \param template template dbcsr matrix giving its mp_env
1151 !> \param n global col size of output matrix
1152 !> \param sym ...
1153 !> \param data_type ...
1154 ! **************************************************************************************************
1155  SUBROUTINE cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym, data_type)
1156  TYPE(dbcsr_type), INTENT(INOUT) :: matrix, template
1157  INTEGER :: n
1158  CHARACTER, OPTIONAL :: sym
1159  INTEGER, OPTIONAL :: data_type
1160 
1161  CHARACTER :: mysym
1162  INTEGER :: my_data_type, npcols
1163  INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1164  row_dist
1165  TYPE(dbcsr_distribution_type) :: dist_m_n, tmpl_dist
1166 
1167  mysym = dbcsr_get_matrix_type(template)
1168  IF (PRESENT(sym)) mysym = sym
1169  my_data_type = dbcsr_get_data_type(template)
1170  IF (PRESENT(data_type)) my_data_type = data_type
1171 
1172  CALL dbcsr_get_info(template, distribution=tmpl_dist)
1173  CALL dbcsr_distribution_get(tmpl_dist, &
1174  npcols=npcols, &
1175  row_dist=row_dist)
1176 
1177  NULLIFY (col_dist, col_blk_size)
1178  CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
1179  CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
1180  row_dist=row_dist, col_dist=col_dist)
1181 
1182  CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
1183  CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, &
1184  row_blk_size, col_blk_size, nze=0, data_type=my_data_type)
1185 
1186  DEALLOCATE (col_dist, col_blk_size)
1187  CALL dbcsr_distribution_release(dist_m_n)
1188 
1189  END SUBROUTINE cp_dbcsr_m_by_n_from_row_template
1190 
1191 ! **************************************************************************************************
1192 !> \brief Distributes elements into blocks and into bins
1193 !>
1194 !> \param[out] block_distribution block distribution to bins
1195 !> \param[out] block_size sizes of blocks
1196 !> \param[in] nelements number of elements to bin
1197 !> \param[in] nbins number of bins
1198 !> \par Term clarification
1199 !> An example: blocks are atom blocks and bins are process rows/columns.
1200 ! **************************************************************************************************
1201  SUBROUTINE create_bl_distribution(block_distribution, &
1202  block_size, nelements, nbins)
1203  INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: block_distribution, block_size
1204  INTEGER, INTENT(IN) :: nelements, nbins
1205 
1206  CHARACTER(len=*), PARAMETER :: routinen = 'create_bl_distribution', &
1207  routinep = modulen//':'//routinen
1208 
1209  INTEGER :: bin, blk_layer, element_stack, els, &
1210  estimated_blocks, max_blocks_per_bin, &
1211  nblks, nblocks, stat
1212  INTEGER, DIMENSION(:), POINTER :: blk_dist, blk_sizes
1213 
1214 ! ---------------------------------------------------------------------------
1215 
1216  NULLIFY (block_distribution)
1217  NULLIFY (block_size)
1218  ! Define the sizes on which we build the distribution.
1219  IF (nelements .GT. 0) THEN
1220 
1221  nblocks = ceiling(real(nelements, kind=dp)/real(max_elements_per_block, kind=dp))
1222  max_blocks_per_bin = ceiling(real(nblocks, kind=dp)/real(nbins, kind=dp))
1223 
1224  IF (debug_mod) THEN
1225  WRITE (*, '(1X,A,1X,A,I7,A,I7,A)') routinep, "For", nelements, &
1226  " elements and", nbins, " bins"
1227  WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1228  max_elements_per_block, " max elements per block"
1229  WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1230  nblocks, " blocks"
1231  WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1232  max_blocks_per_bin, " max blocks/bin"
1233  END IF
1234 
1235  estimated_blocks = max_blocks_per_bin*nbins
1236  ALLOCATE (blk_dist(estimated_blocks), stat=stat)
1237  IF (stat /= 0) &
1238  cpabort("blk_dist")
1239  ALLOCATE (blk_sizes(estimated_blocks), stat=stat)
1240  IF (stat /= 0) &
1241  cpabort("blk_sizes")
1242  element_stack = 0
1243  nblks = 0
1244  DO blk_layer = 1, max_blocks_per_bin
1245  DO bin = 0, nbins - 1
1246  els = min(max_elements_per_block, nelements - element_stack)
1247  IF (els .GT. 0) THEN
1248  element_stack = element_stack + els
1249  nblks = nblks + 1
1250  blk_dist(nblks) = bin
1251  blk_sizes(nblks) = els
1252  IF (debug_mod) WRITE (*, '(1X,A,I5,A,I5,A,I5)') routinep//" Assigning", &
1253  els, " elements as block", nblks, " to bin", bin
1254  END IF
1255  END DO
1256  END DO
1257  ! Create the output arrays.
1258  IF (nblks .EQ. estimated_blocks) THEN
1259  block_distribution => blk_dist
1260  block_size => blk_sizes
1261  ELSE
1262  ALLOCATE (block_distribution(nblks), stat=stat)
1263  IF (stat /= 0) &
1264  cpabort("blk_dist")
1265  block_distribution(:) = blk_dist(1:nblks)
1266  DEALLOCATE (blk_dist)
1267  ALLOCATE (block_size(nblks), stat=stat)
1268  IF (stat /= 0) &
1269  cpabort("blk_sizes")
1270  block_size(:) = blk_sizes(1:nblks)
1271  DEALLOCATE (blk_sizes)
1272  END IF
1273  ELSE
1274  ALLOCATE (block_distribution(0), stat=stat)
1275  IF (stat /= 0) &
1276  cpabort("blk_dist")
1277  ALLOCATE (block_size(0), stat=stat)
1278  IF (stat /= 0) &
1279  cpabort("blk_sizes")
1280  END IF
1281 1579 FORMAT(i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5)
1282  IF (debug_mod) THEN
1283  WRITE (*, '(1X,A,A)') routinep//" Distribution"
1284  WRITE (*, 1579) block_distribution(:)
1285  WRITE (*, '(1X,A,A)') routinep//" Sizes"
1286  WRITE (*, 1579) block_size(:)
1287  END IF
1288  END SUBROUTINE create_bl_distribution
1289 
1290 ! **************************************************************************************************
1291 !> \brief Creates a new distribution for the right matrix in a matrix
1292 !> multiplication with unrotated grid.
1293 !> \param[out] dist_right new distribution for the right matrix
1294 !> \param[in] dist_left the distribution of the left matrix
1295 !> \param[in] ncolumns number of columns in right matrix
1296 !> \param[out] right_col_blk_sizes sizes of blocks in the created column
1297 !> \par The new row distribution for the right matrix is the same as the row
1298 !> distribution of the left matrix, while the column distribution is
1299 !> created so that it is appropriate to the parallel environment.
1300 ! **************************************************************************************************
1301  SUBROUTINE dbcsr_create_dist_r_unrot(dist_right, dist_left, ncolumns, &
1302  right_col_blk_sizes)
1303  TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist_right
1304  TYPE(dbcsr_distribution_type), INTENT(IN) :: dist_left
1305  INTEGER, INTENT(IN) :: ncolumns
1306  INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: right_col_blk_sizes
1307 
1308  INTEGER :: multiplicity, nimages, ncols, nprows, npcols
1309  INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_images
1310  INTEGER, DIMENSION(:), POINTER :: right_col_dist, right_row_dist, old_col_dist, dummy
1311 
1312  CALL dbcsr_distribution_get(dist_left, &
1313  ncols=ncols, &
1314  col_dist=old_col_dist, &
1315  nprows=nprows, &
1316  npcols=npcols)
1317 
1318  ! Create the column distribution
1319  CALL create_bl_distribution(right_col_dist, right_col_blk_sizes, ncolumns, npcols)
1320  ! Create an even row distribution.
1321  ALLOCATE (right_row_dist(ncols), tmp_images(ncols))
1322  nimages = lcm(nprows, npcols)/nprows
1323  multiplicity = nprows/gcd(nprows, npcols)
1324  CALL rebin_distribution(right_row_dist, tmp_images, old_col_dist, nprows, multiplicity, nimages)
1325 
1326  NULLIFY (dummy)
1327  CALL dbcsr_distribution_new(dist_right, &
1328  template=dist_left, &
1329  row_dist=right_row_dist, &
1330  col_dist=right_col_dist, &
1331  !row_cluster=dummy,&
1332  !col_cluster=dummy,&
1333  reuse_arrays=.true.)
1334  DEALLOCATE (tmp_images)
1335  END SUBROUTINE dbcsr_create_dist_r_unrot
1336 
1337 ! **************************************************************************************************
1338 !> \brief Makes new distribution with decimation and multiplicity
1339 !> \param[out] new_bins new real distribution
1340 !> \param[out] images new image distribution
1341 !> \param[in] source_bins Basis for the new distribution and images
1342 !> \param[in] nbins number of bins in the new real distribution
1343 !> \param[in] multiplicity multiplicity
1344 !> \param[in] nimages number of images in the new distribution
1345 !> \par Definition of multiplicity and nimages
1346 !> Multiplicity and decimation (number of images) are used to
1347 !> match process grid coordinates on non-square process
1348 !> grids. Given source_nbins and target_nbins, their relation is
1349 !> source_nbins * target_multiplicity
1350 !> = target_nbins * target_nimages.
1351 !> It is best when both multiplicity and nimages are small. To
1352 !> get these two factors, then, one can use the following formulas:
1353 !> nimages = lcm(source_nbins, target_nbins) / target_nbins
1354 !> multiplicity = target_nbins / gcd(source_nbins, target_nbins)
1355 !> from the target's point of view (nimages = target_nimages).
1356 !> \par Mapping
1357 !> The new distribution comprises of real bins and images within
1358 !> bins. These can be view as target_nbins*nimages virtual
1359 !> columns. These same virtual columns are also
1360 !> source_nbins*multiplicity in number. Therefore these virtual
1361 !> columns are mapped from source_nbins*multiplicity onto
1362 !> target_bins*nimages (each target bin has nimages images):
1363 !> Source 4: |1 2 3|4 5 6|7 8 9|A B C| (4*3)
1364 !> Target 6: |1 2|3 4|5 6|7 8|9 A|B C| (6*2)
1365 !> multiplicity=3, nimages=2, 12 virtual columns (1-C).
1366 !> Source bin elements are evenly mapped into one of multiplicity
1367 !> virtual columns. Other (non-even, block-size aware) mappings
1368 !> could be better.
1369 ! **************************************************************************************************
1370  SUBROUTINE rebin_distribution(new_bins, images, source_bins, &
1371  nbins, multiplicity, nimages)
1372  INTEGER, DIMENSION(:), INTENT(OUT) :: new_bins, images
1373  INTEGER, DIMENSION(:), INTENT(IN) :: source_bins
1374  INTEGER, INTENT(IN) :: nbins, multiplicity, nimages
1375 
1376  INTEGER :: bin, i, old_nbins, virtual_bin
1377  INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_multiplier
1378 
1379 ! ---------------------------------------------------------------------------
1380 
1381  IF (mod(nbins*nimages, multiplicity) .NE. 0) &
1382  cpwarn("mulitplicity is not divisor of new process grid coordinate")
1383  old_nbins = (nbins*nimages)/multiplicity
1384  ALLOCATE (bin_multiplier(0:old_nbins - 1))
1385  bin_multiplier(:) = 0
1386  DO i = 1, SIZE(new_bins)
1387  IF (i .LE. SIZE(source_bins)) THEN
1388  bin = source_bins(i)
1389  ELSE
1390  ! Fill remainder with a cyclic distribution
1391  bin = mod(i, old_nbins)
1392  END IF
1393  virtual_bin = bin*multiplicity + bin_multiplier(bin)
1394  new_bins(i) = virtual_bin/nimages
1395  images(i) = 1 + mod(virtual_bin, nimages)
1396  bin_multiplier(bin) = bin_multiplier(bin) + 1
1397  IF (bin_multiplier(bin) .GE. multiplicity) THEN
1398  bin_multiplier(bin) = 0
1399  END IF
1400  END DO
1401  END SUBROUTINE rebin_distribution
1402 
1403 ! **************************************************************************************************
1404 !> \brief Creates a block-cyclic compatible distribution
1405 !>
1406 !> All blocks in a dimension, except for possibly the last
1407 !> block, have the same size.
1408 !> \param[out] dist the elemental distribution
1409 !> \param[in] nrows number of full rows
1410 !> \param[in] ncolumns number of full columns
1411 !> \param[in] nrow_block size of row blocks
1412 !> \param[in] ncol_block size of column blocks
1413 !> \param[in] mp_env multiprocess environment
1414 !> \param[out] row_blk_sizes row block sizes
1415 !> \param[out] col_blk_sizes column block sizes
1416 ! **************************************************************************************************
1417  SUBROUTINE dbcsr_create_dist_block_cyclic(dist, nrows, ncolumns, &
1418  nrow_block, ncol_block, group_handle, pgrid, row_blk_sizes, col_blk_sizes)
1419  TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist
1420  INTEGER, INTENT(IN) :: nrows, ncolumns, nrow_block, ncol_block
1421  INTEGER, INTENT(IN) :: group_handle
1422  INTEGER, DIMENSION(:, :), POINTER :: pgrid
1423  INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: row_blk_sizes, col_blk_sizes
1424 
1425  CHARACTER(len=*), PARAMETER :: routinen = 'dbcsr_create_dist_block_cyclic'
1426 
1427  INTEGER :: nblkcols, nblkrows, npcols, nprows, &
1428  pdim, sz
1429  INTEGER, DIMENSION(:), POINTER :: cd_data, rd_data
1430 
1431  ! Row sizes
1432  IF (nrow_block .EQ. 0) THEN
1433  nblkrows = 0
1434  sz = 0
1435  ELSE
1436  nblkrows = nrows/nrow_block
1437  sz = mod(nrows, nrow_block)
1438  END IF
1439  IF (sz .GT. 0) nblkrows = nblkrows + 1
1440  ALLOCATE (row_blk_sizes(nblkrows), rd_data(nblkrows))
1441  row_blk_sizes = nrow_block
1442  IF (sz .NE. 0) row_blk_sizes(nblkrows) = sz
1443 
1444  ! Column sizes
1445  IF (ncol_block .EQ. 0) THEN
1446  nblkcols = 0
1447  sz = 0
1448  ELSE
1449  nblkcols = ncolumns/ncol_block
1450  sz = mod(ncolumns, ncol_block)
1451  END IF
1452  IF (sz .GT. 0) nblkcols = nblkcols + 1
1453  ALLOCATE (col_blk_sizes(nblkcols), cd_data(nblkcols))
1454  col_blk_sizes = ncol_block
1455  IF (sz .NE. 0) col_blk_sizes(nblkcols) = sz
1456  !
1457  IF (debug_mod) THEN
1458  WRITE (*, *) routinen//" nrows,nrow_block,nblkrows=", &
1459  nrows, nrow_block, nblkrows
1460  WRITE (*, *) routinen//" ncols,ncol_block,nblkcols=", &
1461  ncolumns, ncol_block, nblkcols
1462  END IF
1463  ! Calculate process row distribution
1464  nprows = SIZE(pgrid, 1)
1465  DO pdim = 0, min(nprows - 1, nblkrows - 1)
1466  rd_data(1 + pdim:nblkrows:nprows) = pdim
1467  END DO
1468  ! Calculate process column distribution
1469  npcols = SIZE(pgrid, 2)
1470  DO pdim = 0, min(npcols - 1, nblkcols - 1)
1471  cd_data(1 + pdim:nblkcols:npcols) = pdim
1472  END DO
1473  !
1474  IF (debug_mod) THEN
1475  WRITE (*, *) routinen//" row_dist", &
1476  rd_data
1477  WRITE (*, *) routinen//" col_dist", &
1478  cd_data
1479  END IF
1480  !
1481  CALL dbcsr_distribution_new(dist, &
1482  group=group_handle, pgrid=pgrid, &
1483  row_dist=rd_data, &
1484  col_dist=cd_data, &
1485  reuse_arrays=.true.)
1486 
1487  END SUBROUTINE dbcsr_create_dist_block_cyclic
1488 
1489 ! **************************************************************************************************
1490 !> \brief Allocate and initialize a real matrix 1-dimensional set.
1491 !> \param[in,out] matrix_set Set containing the DBCSR matrices
1492 !> \param[in] nmatrix Size of set
1493 !> \par History
1494 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1495 ! **************************************************************************************************
1496  SUBROUTINE allocate_dbcsr_matrix_set_1d(matrix_set, nmatrix)
1497  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_set
1498  INTEGER, INTENT(IN) :: nmatrix
1499 
1500  INTEGER :: imatrix
1501 
1502  IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1503  ALLOCATE (matrix_set(nmatrix))
1504  DO imatrix = 1, nmatrix
1505  NULLIFY (matrix_set(imatrix)%matrix)
1506  END DO
1507  END SUBROUTINE allocate_dbcsr_matrix_set_1d
1508 
1509 ! **************************************************************************************************
1510 !> \brief Allocate and initialize a real matrix 2-dimensional set.
1511 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1512 !> \param[in] nmatrix Size of set
1513 !> \param mmatrix ...
1514 !> \par History
1515 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1516 ! **************************************************************************************************
1517  SUBROUTINE allocate_dbcsr_matrix_set_2d(matrix_set, nmatrix, mmatrix)
1518  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_set
1519  INTEGER, INTENT(IN) :: nmatrix, mmatrix
1520 
1521  INTEGER :: imatrix, jmatrix
1522 
1523  IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1524  ALLOCATE (matrix_set(nmatrix, mmatrix))
1525  DO jmatrix = 1, mmatrix
1526  DO imatrix = 1, nmatrix
1527  NULLIFY (matrix_set(imatrix, jmatrix)%matrix)
1528  END DO
1529  END DO
1530  END SUBROUTINE allocate_dbcsr_matrix_set_2d
1531 
1532 ! **************************************************************************************************
1533 !> \brief Allocate and initialize a real matrix 3-dimensional set.
1534 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1535 !> \param[in] nmatrix Size of set
1536 !> \param mmatrix ...
1537 !> \param pmatrix ...
1538 !> \par History
1539 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1540 ! **************************************************************************************************
1541  SUBROUTINE allocate_dbcsr_matrix_set_3d(matrix_set, nmatrix, mmatrix, pmatrix)
1542  TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: matrix_set
1543  INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix
1544 
1545  INTEGER :: imatrix, jmatrix, kmatrix
1546 
1547  IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1548  ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix))
1549  DO kmatrix = 1, pmatrix
1550  DO jmatrix = 1, mmatrix
1551  DO imatrix = 1, nmatrix
1552  NULLIFY (matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1553  END DO
1554  END DO
1555  END DO
1556  END SUBROUTINE allocate_dbcsr_matrix_set_3d
1557 
1558 ! **************************************************************************************************
1559 !> \brief Allocate and initialize a real matrix 4-dimensional set.
1560 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1561 !> \param[in] nmatrix Size of set
1562 !> \param mmatrix ...
1563 !> \param pmatrix ...
1564 !> \par History
1565 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1566 ! **************************************************************************************************
1567  SUBROUTINE allocate_dbcsr_matrix_set_4d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix)
1568  TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
1569  INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix, qmatrix
1570 
1571  INTEGER :: imatrix, jmatrix, kmatrix, lmatrix
1572 
1573  IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1574  ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix))
1575  DO lmatrix = 1, qmatrix
1576  DO kmatrix = 1, pmatrix
1577  DO jmatrix = 1, mmatrix
1578  DO imatrix = 1, nmatrix
1579  NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1580  END DO
1581  END DO
1582  END DO
1583  END DO
1584  END SUBROUTINE allocate_dbcsr_matrix_set_4d
1585 
1586 ! **************************************************************************************************
1587 !> \brief Allocate and initialize a real matrix 5-dimensional set.
1588 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1589 !> \param[in] nmatrix Size of set
1590 !> \param mmatrix ...
1591 !> \param pmatrix ...
1592 !> \par History
1593 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1594 ! **************************************************************************************************
1595  SUBROUTINE allocate_dbcsr_matrix_set_5d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix, smatrix)
1596  TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), POINTER :: matrix_set
1597  INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix, qmatrix, smatrix
1598 
1599  INTEGER :: imatrix, jmatrix, kmatrix, lmatrix, hmatrix
1600 
1601  IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1602  ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix, smatrix))
1603  DO hmatrix = 1, smatrix
1604  DO lmatrix = 1, qmatrix
1605  DO kmatrix = 1, pmatrix
1606  DO jmatrix = 1, mmatrix
1607  DO imatrix = 1, nmatrix
1608  NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
1609  END DO
1610  END DO
1611  END DO
1612  END DO
1613  END DO
1614  END SUBROUTINE allocate_dbcsr_matrix_set_5d
1615 
1616  ! **************************************************************************************************
1617 !> \brief Deallocate a real matrix set and release all of the member matrices.
1618 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1619 !> \par History
1620 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1621 ! **************************************************************************************************
1622  SUBROUTINE deallocate_dbcsr_matrix_set_1d(matrix_set)
1623 
1624  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_set
1625 
1626  INTEGER :: imatrix
1627 
1628  IF (ASSOCIATED(matrix_set)) THEN
1629  DO imatrix = 1, SIZE(matrix_set)
1630  CALL dbcsr_deallocate_matrix(matrix_set(imatrix)%matrix)
1631  END DO
1632  DEALLOCATE (matrix_set)
1633  END IF
1634 
1635  END SUBROUTINE deallocate_dbcsr_matrix_set_1d
1636 
1637 ! **************************************************************************************************
1638 !> \brief Deallocate a real matrix set and release all of the member matrices.
1639 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1640 !> \par History
1641 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1642 ! **************************************************************************************************
1643  SUBROUTINE deallocate_dbcsr_matrix_set_2d(matrix_set)
1644 
1645  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_set
1646 
1647  INTEGER :: imatrix, jmatrix
1648 
1649  IF (ASSOCIATED(matrix_set)) THEN
1650  DO jmatrix = 1, SIZE(matrix_set, 2)
1651  DO imatrix = 1, SIZE(matrix_set, 1)
1652  CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix)%matrix)
1653  END DO
1654  END DO
1655  DEALLOCATE (matrix_set)
1656  END IF
1657  END SUBROUTINE deallocate_dbcsr_matrix_set_2d
1658 
1659 ! **************************************************************************************************
1660 !> \brief Deallocate a real matrix set and release all of the member matrices.
1661 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1662 !> \par History
1663 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1664 ! **************************************************************************************************
1665  SUBROUTINE deallocate_dbcsr_matrix_set_3d(matrix_set)
1666 
1667  TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: matrix_set
1668 
1669  INTEGER :: imatrix, jmatrix, kmatrix
1670 
1671  IF (ASSOCIATED(matrix_set)) THEN
1672  DO kmatrix = 1, SIZE(matrix_set, 3)
1673  DO jmatrix = 1, SIZE(matrix_set, 2)
1674  DO imatrix = 1, SIZE(matrix_set, 1)
1675  CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1676  END DO
1677  END DO
1678  END DO
1679  DEALLOCATE (matrix_set)
1680  END IF
1681  END SUBROUTINE deallocate_dbcsr_matrix_set_3d
1682 
1683 ! **************************************************************************************************
1684 !> \brief Deallocate a real matrix set and release all of the member matrices.
1685 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1686 !> \par History
1687 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1688 ! **************************************************************************************************
1689  SUBROUTINE deallocate_dbcsr_matrix_set_4d(matrix_set)
1690 
1691  TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
1692 
1693  INTEGER :: imatrix, jmatrix, kmatrix, lmatrix
1694 
1695  IF (ASSOCIATED(matrix_set)) THEN
1696  DO lmatrix = 1, SIZE(matrix_set, 4)
1697  DO kmatrix = 1, SIZE(matrix_set, 3)
1698  DO jmatrix = 1, SIZE(matrix_set, 2)
1699  DO imatrix = 1, SIZE(matrix_set, 1)
1700  CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1701  END DO
1702  END DO
1703  END DO
1704  END DO
1705  DEALLOCATE (matrix_set)
1706  END IF
1707  END SUBROUTINE deallocate_dbcsr_matrix_set_4d
1708 
1709 ! **************************************************************************************************
1710 !> \brief Deallocate a real matrix set and release all of the member matrices.
1711 !> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1712 !> \par History
1713 !> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1714 ! **************************************************************************************************
1715  SUBROUTINE deallocate_dbcsr_matrix_set_5d(matrix_set)
1716 
1717  TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), POINTER :: matrix_set
1718 
1719  INTEGER :: imatrix, jmatrix, kmatrix, hmatrix, lmatrix
1720 
1721  IF (ASSOCIATED(matrix_set)) THEN
1722  DO hmatrix = 1, SIZE(matrix_set, 5)
1723  DO lmatrix = 1, SIZE(matrix_set, 4)
1724  DO kmatrix = 1, SIZE(matrix_set, 3)
1725  DO jmatrix = 1, SIZE(matrix_set, 2)
1726  DO imatrix = 1, SIZE(matrix_set, 1)
1727  CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
1728  END DO
1729  END DO
1730  END DO
1731  END DO
1732  END DO
1733  DEALLOCATE (matrix_set)
1734  END IF
1735  END SUBROUTINE deallocate_dbcsr_matrix_set_5d
1736 
1737 END MODULE cp_dbcsr_operations
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.
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
integer, save, public max_elements_per_block
subroutine, public copy_cfm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
subroutine, public dbcsr_create_dist_r_unrot(dist_right, dist_left, ncolumns, right_col_blk_sizes)
Creates a new distribution for the right matrix in a matrix multiplication with unrotated grid.
subroutine, public dbcsr_multiply_local(matrix_a, vec_b, vec_c, ncol, alpha)
multiply a dbcsr with a replicated array c = alpha_scalar * A (dbscr) * b + c
logical, parameter debug_mod
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr_bc(fm, bc_mat)
Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution, which requires no com...
subroutine, public copy_dbcsr_to_fm_bc(bc_mat, fm)
Copy a DBCSR_BLACS matrix to a BLACS matrix.
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
subroutine, public copy_dbcsr_to_cfm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym, data_type)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
subroutine, public dbcsr_copy_columns_hack(matrix_b, matrix_a, ncol, source_start, target_start, para_env, blacs_env)
hack for dbcsr_copy_columns
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
computes matrix_c = beta * matrix_c + alpha * ( matrix_a ** transa ) * ( matrix_b ** transb )
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
subroutine, public distribution_2d_get(distribution_2d, row_distribution, col_distribution, n_row_distribution, n_col_distribution, n_local_rows, n_local_cols, local_rows, local_cols, flat_local_rows, flat_local_cols, n_flat_local_rows, n_flat_local_cols, blacs_env)
returns various attributes about the distribution_2d
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
Definition: mathlib.F:1334
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition: mathlib.F:1299
Interface to the message passing library MPI.