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