(git:374b731)
Loading...
Searching...
No Matches
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
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
41 USE cp_fm_types, ONLY: cp_fm_create, &
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
74
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
99CONTAINS
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
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
12811579 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
1737END 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
Represents a complex full matrix distributed on many processors.
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
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
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
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
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:1326
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1291
Interface to the message passing library MPI.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Represent a complex full matrix.
keeps the information about the structure of a full matrix
represent a full matrix
distributes pairs on a 2d grid of processors
stores all the informations relevant to an mpi environment