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