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 :: 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 :
112 : CONTAINS
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 772 : 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 772 : TYPE(local_sum_type), ALLOCATABLE :: work(:)
131 :
132 772 : CPASSERT(SIZE(matrix_c) == SIZE(matrix_a, 1))
133 772 : CPASSERT(SIZE(alpha) == SIZE(matrix_a, 2))
134 772 : IF (SIZE(matrix_c) == 0) RETURN
135 772 : nworkers = 1
136 772 : !$ nworkers = MIN(SIZE(matrix_c), omp_get_max_threads())
137 3088 : ALLOCATE (work(nworkers))
138 4568 : DO first = 1, SIZE(matrix_c), nworkers
139 3796 : nbatch = MIN(nworkers, SIZE(matrix_c) - first + 1)
140 7592 : DO j = 1, nbatch
141 3796 : CALL dbcsr_get_info(matrix_c(first + j - 1)%matrix, nfullrows_total=n, nfullcols_total=k)
142 3796 : IF (ALLOCATED(work(j)%value)) THEN
143 3024 : IF (SIZE(work(j)%value, 1) /= n .OR. SIZE(work(j)%value, 2) /= k) THEN
144 0 : DEALLOCATE (work(j)%value, work(j)%buffer)
145 : END IF
146 : END IF
147 7656 : IF (.NOT. ALLOCATED(work(j)%value)) ALLOCATE (work(j)%value(n, k), work(j)%buffer(n, k))
148 699244 : work(j)%value(:, :) = 0.0_dp
149 : END DO
150 37124 : DO i = 1, SIZE(alpha)
151 66656 : DO j = 1, nbatch
152 66656 : CALL copy_dbcsr_to_local(matrix_a(first + j - 1, i)%matrix, work(j)%buffer)
153 : END DO
154 37124 : !$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 8364 : DO j = 1, nbatch
161 7592 : CALL copy_local_to_dbcsr(matrix_c(first + j - 1)%matrix, work(j)%value)
162 : END DO
163 : END DO
164 1544 : 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 99280 : 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 99280 : REAL(dp), POINTER :: block(:, :)
178 : TYPE(dbcsr_iterator_type) :: iter
179 :
180 99280 : CPASSERT(dbcsr_get_matrix_type(matrix) == dbcsr_type_no_symmetry)
181 20658348 : dense(:, :) = 0.0_dp
182 99280 : CALL dbcsr_iterator_readonly_start(iter, matrix)
183 274700 : DO WHILE (dbcsr_iterator_blocks_left(iter))
184 : CALL dbcsr_iterator_next_block(iter, block=block, transposed=transposed, &
185 175420 : row_size=nr, col_size=nc, row_offset=row, col_offset=col)
186 274700 : IF (transposed) THEN
187 0 : dense(row:row + nr - 1, col:col + nc - 1) = TRANSPOSE(block)
188 : ELSE
189 20765976 : dense(row:row + nr - 1, col:col + nc - 1) = block
190 : END IF
191 : END DO
192 99280 : CALL dbcsr_iterator_stop(iter)
193 99280 : 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 52206 : 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 52206 : REAL(dp), POINTER :: block(:, :)
207 : TYPE(dbcsr_iterator_type) :: iter
208 :
209 52206 : CPASSERT(dbcsr_get_matrix_type(matrix) == dbcsr_type_no_symmetry)
210 52206 : CALL dbcsr_reserve_all_blocks(matrix)
211 52206 : CALL dbcsr_iterator_start(iter, matrix)
212 131636 : DO WHILE (dbcsr_iterator_blocks_left(iter))
213 : CALL dbcsr_iterator_next_block(iter, block=block, transposed=transposed, &
214 79430 : row_size=nr, col_size=nc, row_offset=row, col_offset=col)
215 131636 : IF (transposed) THEN
216 0 : block(:, :) = TRANSPOSE(dense(row:row + nr - 1, col:col + nc - 1))
217 : ELSE
218 10245886 : block(:, :) = dense(row:row + nr - 1, col:col + nc - 1)
219 : END IF
220 : END DO
221 52206 : CALL dbcsr_iterator_stop(iter)
222 52206 : 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 1319433 : 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 1319433 : CALL timeset(routineN, handle)
251 :
252 1319433 : my_keep_sparsity = .FALSE.
253 1319433 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
254 :
255 1319433 : CALL copy_fm_to_dbcsr_bc(fm, bc_mat)
256 :
257 1319433 : IF (my_keep_sparsity) THEN
258 46742 : CALL dbcsr_create(redist_mat, template=matrix)
259 46742 : CALL dbcsr_complete_redistribute(bc_mat, redist_mat)
260 46742 : CALL dbcsr_copy(matrix, redist_mat, keep_sparsity=.TRUE.)
261 46742 : CALL dbcsr_release(redist_mat)
262 : ELSE
263 1272691 : CALL dbcsr_complete_redistribute(bc_mat, matrix)
264 : END IF
265 :
266 1319433 : CALL dbcsr_release(bc_mat)
267 :
268 1319433 : CALL timestop(handle)
269 1319433 : 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 1322409 : 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 1322409 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
286 1322409 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
287 1322409 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
288 1322409 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
289 : TYPE(dbcsr_distribution_type) :: bc_dist
290 : TYPE(dbcsr_iterator_type) :: iter
291 :
292 1322409 : CALL timeset(routineN, handle)
293 :
294 : ! Create processor grid
295 1322409 : pgrid => fm%matrix_struct%context%blacs2mpi
296 :
297 : ! Create a block-cyclic distribution compatible with the FM matrix.
298 1322409 : nrow_block = fm%matrix_struct%nrow_block
299 1322409 : ncol_block = fm%matrix_struct%ncol_block
300 1322409 : nrow_global = fm%matrix_struct%nrow_global
301 1322409 : ncol_global = fm%matrix_struct%ncol_global
302 1322409 : 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 1322409 : 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 1322409 : dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
312 1322409 : CALL dbcsr_distribution_release(bc_dist)
313 :
314 : ! allocate all blocks
315 1322409 : CALL dbcsr_reserve_all_blocks(bc_mat)
316 :
317 1322409 : 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 1322409 : fm_block => fm%local_data
322 : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
323 1322409 : !$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 1322409 : CALL timestop(handle)
333 2644818 : 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 1274793 : 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 1274793 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
357 1274793 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
358 : LOGICAL :: validate_plan
359 : TYPE(dbcsr_distribution_type) :: bc_dist, dist
360 1274793 : TYPE(dbcsr_fm_copy_plan_type) :: local_plan
361 : TYPE(dbcsr_type) :: bc_mat, matrix_nosym
362 :
363 1274793 : CALL timeset(routineN, handle)
364 1274793 : validate_plan = .FALSE.
365 1274793 : 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 1274793 : nfullcols_total=nfullcols_total)
373 :
374 1274793 : CPASSERT(fm%matrix_struct%nrow_global == nfullrows_total)
375 1274793 : CPASSERT(fm%matrix_struct%ncol_global == nfullcols_total)
376 :
377 419887 : SELECT CASE (dbcsr_get_matrix_type(matrix))
378 : CASE (dbcsr_type_symmetric, dbcsr_type_antisymmetric)
379 419887 : IF (PRESENT(plan)) THEN
380 243046 : CALL copy_dbcsr_to_fm_sym(matrix, fm, plan, validate_plan)
381 : ELSE
382 176841 : CALL copy_dbcsr_to_fm_sym(matrix, fm, local_plan, .FALSE.)
383 : END IF
384 419887 : CALL timestop(handle)
385 1694680 : RETURN
386 : END SELECT
387 :
388 : ! Other matrix types use the uncached redistribution path.
389 854906 : IF (PRESENT(plan)) plan = dbcsr_fm_copy_plan_type()
390 :
391 : ! info about the full matrix
392 854906 : nrow_block = fm%matrix_struct%nrow_block
393 854906 : ncol_block = fm%matrix_struct%ncol_block
394 :
395 : ! Convert DBCSR to a block-cyclic
396 854906 : NULLIFY (col_blk_size, row_blk_size)
397 854906 : 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 854906 : 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 854906 : dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
406 854906 : CALL dbcsr_distribution_release(bc_dist)
407 :
408 854906 : IF (dbcsr_get_matrix_type(matrix) == dbcsr_type_no_symmetry) THEN
409 854906 : CALL dbcsr_complete_redistribute(matrix, bc_mat)
410 : ELSE
411 0 : CALL dbcsr_create(matrix_nosym, template=matrix, matrix_type=dbcsr_type_no_symmetry)
412 0 : CALL dbcsr_desymmetrize(matrix, matrix_nosym)
413 0 : CALL dbcsr_complete_redistribute(matrix_nosym, bc_mat)
414 854906 : CALL dbcsr_release(matrix_nosym)
415 : END IF
416 :
417 854906 : CALL copy_dbcsr_to_fm_bc(bc_mat, fm)
418 :
419 854906 : CALL dbcsr_release(bc_mat)
420 :
421 854906 : CALL timestop(handle)
422 2726427 : 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 419887 : 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 419887 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
445 : TYPE(dbcsr_iterator_type) :: iter
446 : TYPE(mp_para_env_type), POINTER :: para_env
447 :
448 199835 : IF (.NOT. ALLOCATED(plan%thread)) CALL dbcsr_to_fm_plan_create(matrix, fm, plan)
449 419887 : para_env => fm%matrix_struct%context%para_env
450 419887 : symmetry_sign = 1.0_dp
451 419887 : IF (dbcsr_get_matrix_type(matrix) == dbcsr_type_antisymmetric) symmetry_sign = -1.0_dp
452 :
453 420071 : DO attempt = 1, 2
454 : layout_matches = ASSOCIATED(plan%matrix_struct, fm%matrix_struct) .AND. &
455 420071 : 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 420071 : !$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 420071 : IF (check_plan) CALL para_env%max(rebuild)
504 420071 : IF (rebuild == 0) EXIT
505 184 : CPASSERT(check_plan .AND. attempt == 1)
506 420071 : 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 419887 : plan%recv_values, plan%recv_count, plan%recv_disp)
512 1257252 : nrecv = SUM(plan%recv_count)
513 206789356 : fm%local_data = 0.0_dp
514 419887 : !$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 419887 : 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 200203 : 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 200019 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: local_pos
543 200019 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:, :) :: thread_pos
544 200019 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nblocks
545 200019 : 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 200019 : CALL timeset(routineN, handle)
552 200019 : fms => fm%matrix_struct
553 200019 : plan%matrix_struct => fms
554 200019 : plan%matrix_type = dbcsr_get_matrix_type(matrix)
555 200019 : para_env => fms%context%para_env
556 200019 : nproc = para_env%num_pe
557 200019 : nthreads = 1
558 200019 : !$ nthreads = omp_get_max_threads()
559 1200114 : ALLOCATE (thread_pos(0:nproc - 1, 0:nthreads - 1), nblocks(0:nthreads - 1))
560 0 : ALLOCATE (plan%send_count(0:nproc - 1), plan%send_disp(0:nproc - 1), &
561 1200114 : plan%recv_count(0:nproc - 1), plan%recv_disp(0:nproc - 1))
562 200019 : thread_pos(:, :) = 0
563 200019 : 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 200019 : !$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 200019 : DEALLOCATE (thread_pos, nblocks)
642 :
643 200019 : CALL para_env%alltoall(plan%send_count, plan%recv_count, 1)
644 599092 : CPASSERT(SUM(INT(plan%recv_count, int_8)) <= INT(HUGE(0), int_8)/2)
645 200019 : plan%recv_disp(0) = 0
646 399073 : DO dest = 1, nproc - 1
647 399073 : plan%recv_disp(dest) = plan%recv_disp(dest - 1) + plan%recv_count(dest - 1)
648 : END DO
649 0 : ALLOCATE (plan%recv_index(2, MAX(1, SUM(plan%recv_count))), &
650 1798241 : plan%recv_values(MAX(1, SUM(plan%recv_count))))
651 209375643 : plan%recv_index(:, :) = 0
652 69925227 : plan%recv_values(:) = 0.0_dp
653 : CALL para_env%alltoall(send_index, 2*plan%send_count, 2*plan%send_disp, &
654 1796311 : plan%recv_index, 2*plan%recv_count, 2*plan%recv_disp)
655 200019 : CALL timestop(handle)
656 :
657 600057 : 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 854906 : 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 854906 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
672 854906 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
673 : TYPE(dbcsr_iterator_type) :: iter
674 :
675 854906 : CALL timeset(routineN, handle)
676 :
677 854906 : 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 854906 : fm_block => fm%local_data
681 198092118 : fm_block = REAL(0.0, KIND=dp)
682 : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
683 854906 : !$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 854906 : CALL timestop(handle)
693 1709812 : 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 2177315 : 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 2177315 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, local_cols, local_rows, &
712 2177315 : 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 2177315 : col_blk_size=col_blk_size)
723 :
724 : ! calculate first_row and last_row
725 6530697 : ALLOCATE (local_row_sizes(nblkrows_total))
726 2177315 : local_row_sizes(:) = 0
727 2177315 : IF (nblkrows_local >= 1) THEN
728 5401182 : DO row = 1, nblkrows_local
729 5401182 : local_row_sizes(local_rows(row)) = row_blk_size(local_rows(row))
730 : END DO
731 : END IF
732 6529449 : ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total))
733 2177315 : CALL dbcsr_convert_sizes_to_offsets(local_row_sizes, first_row, last_row)
734 2177315 : DEALLOCATE (local_row_sizes)
735 :
736 : ! calculate first_col and last_col
737 6528329 : ALLOCATE (local_col_sizes(nblkcols_total))
738 2177315 : local_col_sizes(:) = 0
739 2177315 : IF (nblkcols_local >= 1) THEN
740 6209480 : DO col = 1, nblkcols_local
741 6209480 : local_col_sizes(local_cols(col)) = col_blk_size(local_cols(col))
742 : END DO
743 : END IF
744 6524713 : ALLOCATE (first_col(nblkcols_total), last_col(nblkcols_total))
745 2177315 : CALL dbcsr_convert_sizes_to_offsets(local_col_sizes, first_col, last_col)
746 2177315 : DEALLOCATE (local_col_sizes)
747 :
748 2177315 : 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 9232 : 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 2308 : NULLIFY (fm_struct)
775 2308 : 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 2308 : ncol_global=nfullcols_total, para_env=para_env)
778 2308 : CALL cp_fm_create(fm_matrix_a, fm_struct, name="fm_matrix_a")
779 2308 : CALL cp_fm_struct_release(fm_struct)
780 :
781 2308 : 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 2308 : ncol_global=nfullcols_total, para_env=para_env)
784 2308 : CALL cp_fm_create(fm_matrix_b, fm_struct, name="fm_matrix_b")
785 2308 : CALL cp_fm_struct_release(fm_struct)
786 :
787 2308 : CALL copy_dbcsr_to_fm(matrix_a, fm_matrix_a)
788 2308 : CALL copy_dbcsr_to_fm(matrix_b, fm_matrix_b)
789 :
790 2308 : CALL cp_fm_to_fm(fm_matrix_a, fm_matrix_b, ncol, source_start, target_start)
791 :
792 2308 : CALL copy_fm_to_dbcsr(fm_matrix_b, matrix_b)
793 :
794 2308 : CALL cp_fm_release(fm_matrix_a)
795 2308 : CALL cp_fm_release(fm_matrix_b)
796 :
797 2308 : 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 12596 : 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 12596 : INTEGER, DIMENSION(:), POINTER :: col_dist, row_dist
811 12596 : 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 12596 : dist2d_p => dist2d
817 : CALL distribution_2d_get(dist2d_p, &
818 : row_distribution=row_dist_data, &
819 : col_distribution=col_dist_data, &
820 12596 : blacs_env=blacs_env)
821 12596 : CALL blacs_env%get(para_env=para_env, blacs2mpi=pgrid)
822 :
823 : ! map to 1D arrays
824 12596 : row_dist => row_dist_data(:, 1)
825 12596 : 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 12596 : col_dist=col_dist)
833 :
834 12596 : 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 0 : 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 0 : REAL(dp), DIMENSION(:, :), POINTER :: data_d
860 : TYPE(dbcsr_iterator_type) :: iter
861 :
862 0 : CALL timeset(routineN, timing_handle)
863 :
864 0 : my_alpha = 1.0_dp
865 0 : IF (PRESENT(alpha)) my_alpha = alpha
866 :
867 0 : my_ncol = SIZE(vec_b, 2)
868 0 : IF (PRESENT(ncol)) my_ncol = ncol
869 :
870 0 : my_alpha2 = 0.0_dp
871 0 : IF (dbcsr_get_matrix_type(matrix_a) == dbcsr_type_symmetric) my_alpha2 = my_alpha
872 0 : 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 0 : 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 0 : !$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 0 : IF (has_symm) THEN
901 0 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
902 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
903 0 : CALL dbcsr_iterator_next_block(iter, row, col, data_d, row_offset=rowoff, col_offset=coloff)
904 0 : IF (row /= col) THEN
905 0 : 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 0 : 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 0 : 1.0_dp, vec_c(coloff, 1), 1)
916 : END IF
917 : END IF
918 : END DO
919 0 : CALL dbcsr_iterator_stop(iter)
920 : END IF
921 :
922 0 : CALL timestop(timing_handle)
923 0 : 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 3191070 : 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 531845 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_blk_size_right_in, &
952 531845 : col_blk_size_right_out, col_dist, &
953 531845 : 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 531845 : CALL timeset(routineN, timing_handle)
959 :
960 531845 : my_alpha = 1.0_dp
961 531845 : my_beta = 0.0_dp
962 531845 : IF (PRESENT(alpha)) my_alpha = alpha
963 531845 : IF (PRESENT(beta)) my_beta = beta
964 :
965 : ! TODO
966 531845 : CALL cp_fm_get_info(fm_in, ncol_global=b_ncol, nrow_global=b_nrow)
967 531845 : CALL cp_fm_get_info(fm_out, ncol_global=c_ncol, nrow_global=c_nrow)
968 531845 : 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 531845 : CALL cp_fm_get_info(fm_out, ncol_global=k_out)
974 :
975 531845 : 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 531845 : IF (ncol > 0 .AND. k_out > 0 .AND. k_in > 0) THEN
983 530565 : CALL dbcsr_get_info(matrix, row_blk_size=row_blk_size, col_blk_size=col_blk_size, distribution=dist)
984 530565 : 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 530565 : col_blk_size, col_blk_size_right_in)
988 :
989 530565 : CALL dbcsr_distribution_get(dist, row_dist=row_dist)
990 530565 : CALL dbcsr_distribution_get(dist_right_in, col_dist=col_dist)
991 : CALL dbcsr_distribution_new(product_dist, template=dist, &
992 530565 : row_dist=row_dist, col_dist=col_dist)
993 1591695 : ALLOCATE (col_blk_size_right_out(SIZE(col_blk_size_right_in)))
994 2159384 : col_blk_size_right_out = col_blk_size_right_in
995 530565 : 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 530565 : row_blk_size, col_blk_size_right_out)
1004 :
1005 530565 : CALL copy_fm_to_dbcsr(fm_in, in)
1006 530565 : IF (ncol /= k_out .OR. my_beta /= 0.0_dp) THEN
1007 119154 : CALL copy_fm_to_dbcsr(fm_out, out)
1008 : END IF
1009 :
1010 530565 : CALL timeset(routineN//'_core', timing_handle_mult)
1011 : CALL dbcsr_multiply("N", "N", my_alpha, matrix, in, my_beta, out, &
1012 530565 : last_column=ncol)
1013 530565 : CALL timestop(timing_handle_mult)
1014 :
1015 530565 : CALL copy_dbcsr_to_fm(out, fm_out)
1016 :
1017 530565 : CALL dbcsr_release(in)
1018 530565 : CALL dbcsr_release(out)
1019 530565 : DEALLOCATE (col_blk_size_right_in, col_blk_size_right_out)
1020 530565 : CALL dbcsr_distribution_release(dist_right_in)
1021 2122260 : CALL dbcsr_distribution_release(product_dist)
1022 :
1023 : END IF
1024 :
1025 531845 : CALL timestop(timing_handle)
1026 :
1027 531845 : END SUBROUTINE cp_dbcsr_sm_fm_multiply
1028 :
1029 : ! **************************************************************************************************
1030 : !> \brief ...
1031 : !> \param sizes1 ...
1032 : !> \param sizes2 ...
1033 : !> \param full_num ...
1034 : ! **************************************************************************************************
1035 530565 : 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 530565 : n1 = SIZE(sizes1)
1043 530565 : n2 = SIZE(sizes2)
1044 530565 : IF (n1 /= n2) THEN
1045 0 : CPABORT("distributions must be equal!")
1046 : END IF
1047 1079692 : sizes1(1:n1) = sizes2(1:n1)
1048 1079692 : 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 530565 : IF (used < full_num) THEN
1053 0 : sizes1(n1) = sizes1(n1) + full_num - used
1054 : ELSE
1055 530565 : left = used - full_num
1056 530565 : p = n1
1057 530565 : DO WHILE (left > 0 .AND. p > 0)
1058 0 : rm = MIN(left, sizes1(p))
1059 0 : sizes1(p) = sizes1(p) - rm
1060 0 : left = left - rm
1061 0 : p = p - 1
1062 : END DO
1063 : END IF
1064 530565 : 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 265034 : 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 265034 : INTEGER, DIMENSION(:), POINTER :: col_blk_size_left, col_dist_left, &
1098 265034 : 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 265034 : check_product = .FALSE.
1108 :
1109 265034 : CALL timeset(routineN, timing_handle)
1110 :
1111 265034 : my_keep_sparsity = .TRUE.
1112 265034 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
1113 :
1114 265034 : my_symmetry_mode = 0
1115 265034 : IF (PRESENT(symmetry_mode)) my_symmetry_mode = symmetry_mode
1116 :
1117 265034 : NULLIFY (col_dist_left)
1118 :
1119 265034 : IF (ncol > 0) THEN
1120 263044 : IF (.NOT. dbcsr_valid_index(sparse_matrix)) THEN
1121 0 : CPABORT("sparse_matrix must pre-exist")
1122 : END IF
1123 : !
1124 : ! Setup matrix_v
1125 263044 : 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 263044 : CALL dbcsr_get_info(sparse_matrix, distribution=sparse_dist)
1128 263044 : CALL dbcsr_distribution_get(sparse_dist, npcols=npcols, row_dist=row_dist)
1129 263044 : CALL create_bl_distribution(col_dist_left, col_blk_size_left, k, npcols)
1130 : CALL dbcsr_distribution_new(dist_left, template=sparse_dist, &
1131 263044 : row_dist=row_dist, col_dist=col_dist_left)
1132 263044 : DEALLOCATE (col_dist_left)
1133 263044 : 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 263044 : row_blk_size, col_blk_size_left)
1136 263044 : CALL copy_fm_to_dbcsr(matrix_v, mat_v)
1137 263044 : CALL dbcsr_verify_matrix(mat_v)
1138 : !
1139 : ! Setup matrix_g
1140 263044 : IF (PRESENT(matrix_g)) THEN
1141 : CALL dbcsr_create(mat_g, "DBCSR matrix_g", dist_left, dbcsr_type_no_symmetry, &
1142 128331 : row_blk_size, col_blk_size_left)
1143 128331 : CALL copy_fm_to_dbcsr(matrix_g, mat_g)
1144 : END IF
1145 : !
1146 263044 : DEALLOCATE (col_blk_size_left)
1147 263044 : 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 263044 : my_alpha = 1.0_dp
1161 263044 : IF (PRESENT(alpha)) my_alpha = alpha
1162 263044 : IF (PRESENT(matrix_g)) THEN
1163 128331 : 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 43940 : 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 43940 : last_k=ncol)
1173 84391 : 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 2588 : 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 2588 : 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 81803 : 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 134713 : 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 263044 : CALL dbcsr_release(mat_v)
1233 263044 : IF (PRESENT(matrix_g)) CALL dbcsr_release(mat_g)
1234 : END IF
1235 265034 : CALL timestop(timing_handle)
1236 :
1237 265034 : 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 17169 : 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 5723 : INTEGER, DIMENSION(:), POINTER :: col_blk_size_right_in, row_blk_size
1256 : TYPE(dbcsr_distribution_type) :: dist_right_in, tmpl_dist
1257 :
1258 5723 : CALL cp_fm_get_info(fm_in, ncol_global=k_in)
1259 :
1260 5723 : CALL dbcsr_get_info(template, distribution=tmpl_dist)
1261 5723 : CALL dbcsr_create_dist_r_unrot(dist_right_in, tmpl_dist, k_in, col_blk_size_right_in)
1262 5723 : 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 5723 : row_blk_size, col_blk_size_right_in)
1265 :
1266 5723 : CALL copy_fm_to_dbcsr(fm_in, matrix)
1267 5723 : DEALLOCATE (col_blk_size_right_in)
1268 5723 : CALL dbcsr_distribution_release(dist_right_in)
1269 :
1270 5723 : 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 354656 : 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 177328 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1291 177328 : row_dist
1292 : TYPE(dbcsr_distribution_type) :: dist_m_n, tmpl_dist
1293 :
1294 177328 : CALL dbcsr_get_info(template, matrix_type=mysym, distribution=tmpl_dist)
1295 :
1296 177328 : IF (PRESENT(sym)) mysym = sym
1297 :
1298 177328 : NULLIFY (row_dist, col_dist)
1299 177328 : NULLIFY (row_blk_size, col_blk_size)
1300 : !NULLIFY (row_cluster, col_cluster)
1301 :
1302 177328 : CALL dbcsr_distribution_get(tmpl_dist, nprows=nprows, npcols=npcols)
1303 177328 : CALL create_bl_distribution(row_dist, row_blk_size, m, nprows)
1304 177328 : 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 177328 : reuse_arrays=.TRUE.)
1309 :
1310 : CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, &
1311 177328 : row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
1312 177328 : CALL dbcsr_distribution_release(dist_m_n)
1313 :
1314 177328 : 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 703648 : 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 175912 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1334 175912 : row_dist
1335 : TYPE(dbcsr_distribution_type) :: dist_m_n, tmpl_dist
1336 :
1337 351824 : mysym = dbcsr_get_matrix_type(template)
1338 175912 : IF (PRESENT(sym)) mysym = sym
1339 :
1340 175912 : CALL dbcsr_get_info(template, distribution=tmpl_dist)
1341 : CALL dbcsr_distribution_get(tmpl_dist, &
1342 : npcols=npcols, &
1343 175912 : row_dist=row_dist)
1344 :
1345 175912 : NULLIFY (col_dist, col_blk_size)
1346 175912 : CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
1347 : CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
1348 175912 : row_dist=row_dist, col_dist=col_dist)
1349 :
1350 175912 : CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
1351 175912 : CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, row_blk_size, col_blk_size)
1352 :
1353 175912 : DEALLOCATE (col_dist, col_blk_size)
1354 175912 : CALL dbcsr_distribution_release(dist_m_n)
1355 :
1356 175912 : 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 1329900 : 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 1329900 : INTEGER, DIMENSION(:), POINTER :: blk_dist, blk_sizes
1380 :
1381 : ! ---------------------------------------------------------------------------
1382 :
1383 1329900 : NULLIFY (block_distribution)
1384 1329900 : NULLIFY (block_size)
1385 : ! Define the sizes on which we build the distribution.
1386 1329900 : IF (nelements > 0) THEN
1387 :
1388 1314502 : nblocks = CEILING(REAL(nelements, KIND=dp)/REAL(max_elements_per_block, KIND=dp))
1389 1314502 : 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 1314502 : estimated_blocks = max_blocks_per_bin*nbins
1403 3943506 : ALLOCATE (blk_dist(estimated_blocks), stat=stat)
1404 1314502 : IF (stat /= 0) THEN
1405 0 : CPABORT("blk_dist")
1406 : END IF
1407 2629004 : ALLOCATE (blk_sizes(estimated_blocks), stat=stat)
1408 : IF (stat /= 0) THEN
1409 0 : CPABORT("blk_sizes")
1410 : END IF
1411 1314502 : element_stack = 0
1412 1314502 : nblks = 0
1413 2706620 : DO blk_layer = 1, max_blocks_per_bin
1414 4246490 : DO bin = 0, nbins - 1
1415 1539870 : els = MIN(max_elements_per_block, nelements - element_stack)
1416 2931988 : IF (els > 0) THEN
1417 1403042 : element_stack = element_stack + els
1418 1403042 : nblks = nblks + 1
1419 1403042 : blk_dist(nblks) = bin
1420 1403042 : 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 1314502 : IF (nblks == estimated_blocks) THEN
1428 1177674 : block_distribution => blk_dist
1429 1177674 : block_size => blk_sizes
1430 : ELSE
1431 410484 : ALLOCATE (block_distribution(nblks), stat=stat)
1432 : IF (stat /= 0) THEN
1433 0 : CPABORT("blk_dist")
1434 : END IF
1435 550864 : block_distribution(:) = blk_dist(1:nblks)
1436 136828 : DEALLOCATE (blk_dist)
1437 273656 : ALLOCATE (block_size(nblks), stat=stat)
1438 : IF (stat /= 0) THEN
1439 0 : CPABORT("blk_sizes")
1440 : END IF
1441 550864 : block_size(:) = blk_sizes(1:nblks)
1442 136828 : DEALLOCATE (blk_sizes)
1443 : END IF
1444 : ELSE
1445 15398 : ALLOCATE (block_distribution(0), stat=stat)
1446 : IF (stat /= 0) THEN
1447 0 : CPABORT("blk_dist")
1448 : END IF
1449 15398 : ALLOCATE (block_size(0), stat=stat)
1450 : IF (stat /= 0) THEN
1451 0 : CPABORT("blk_sizes")
1452 : END IF
1453 : END IF
1454 : 1579 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 1329900 : 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 536288 : 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 536288 : INTEGER, DIMENSION(:), POINTER :: old_col_dist, right_col_dist, &
1485 536288 : right_row_dist
1486 :
1487 : CALL dbcsr_distribution_get(dist_left, &
1488 : ncols=ncols, &
1489 : col_dist=old_col_dist, &
1490 : nprows=nprows, &
1491 536288 : npcols=npcols)
1492 :
1493 : ! Create the column distribution
1494 536288 : CALL create_bl_distribution(right_col_dist, right_col_blk_sizes, ncolumns, npcols)
1495 : ! Create an even row distribution.
1496 2145152 : ALLOCATE (right_row_dist(ncols), tmp_images(ncols))
1497 536288 : nimages = lcm(nprows, npcols)/nprows
1498 536288 : multiplicity = nprows/gcd(nprows, npcols)
1499 536288 : 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 536288 : reuse_arrays=.TRUE.)
1508 536288 : DEALLOCATE (tmp_images)
1509 536288 : 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 536288 : 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 536288 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_multiplier
1552 :
1553 : ! ---------------------------------------------------------------------------
1554 :
1555 536288 : IF (MOD(nbins*nimages, multiplicity) /= 0) THEN
1556 0 : CPWARN("mulitplicity is not divisor of new process grid coordinate")
1557 : END IF
1558 536288 : old_nbins = (nbins*nimages)/multiplicity
1559 1608864 : ALLOCATE (bin_multiplier(0:old_nbins - 1))
1560 536288 : bin_multiplier(:) = 0
1561 2514177 : DO i = 1, SIZE(new_bins)
1562 1977889 : IF (i <= SIZE(source_bins)) THEN
1563 1977889 : bin = source_bins(i)
1564 : ELSE
1565 : ! Fill remainder with a cyclic distribution
1566 0 : bin = MOD(i, old_nbins)
1567 : END IF
1568 1977889 : virtual_bin = bin*multiplicity + bin_multiplier(bin)
1569 1977889 : new_bins(i) = virtual_bin/nimages
1570 1977889 : images(i) = 1 + MOD(virtual_bin, nimages)
1571 1977889 : bin_multiplier(bin) = bin_multiplier(bin) + 1
1572 2514177 : IF (bin_multiplier(bin) >= multiplicity) THEN
1573 841139 : bin_multiplier(bin) = 0
1574 : END IF
1575 : END DO
1576 536288 : 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 2177315 : 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 2177315 : INTEGER, DIMENSION(:), POINTER :: cd_data, rd_data
1606 :
1607 : ! Row sizes
1608 2177315 : IF (nrow_block == 0) THEN
1609 : nblkrows = 0
1610 : sz = 0
1611 : ELSE
1612 2177315 : nblkrows = nrows/nrow_block
1613 2177315 : sz = MOD(nrows, nrow_block)
1614 : END IF
1615 2177315 : IF (sz > 0) nblkrows = nblkrows + 1
1616 8706764 : ALLOCATE (row_blk_sizes(nblkrows), rd_data(nblkrows))
1617 8550210 : row_blk_sizes = nrow_block
1618 2177315 : IF (sz /= 0) row_blk_sizes(nblkrows) = sz
1619 :
1620 : ! Column sizes
1621 2177315 : IF (ncol_block == 0) THEN
1622 : nblkcols = 0
1623 : sz = 0
1624 : ELSE
1625 2177315 : nblkcols = ncolumns/ncol_block
1626 2177315 : sz = MOD(ncolumns, ncol_block)
1627 : END IF
1628 2177315 : IF (sz > 0) nblkcols = nblkcols + 1
1629 8702028 : ALLOCATE (col_blk_sizes(nblkcols), cd_data(nblkcols))
1630 6213096 : col_blk_sizes = ncol_block
1631 2177315 : 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 2177315 : nprows = SIZE(pgrid, 1)
1641 6447266 : DO pdim = 0, MIN(nprows - 1, nblkrows - 1)
1642 12820161 : rd_data(1 + pdim:nblkrows:nprows) = pdim
1643 : END DO
1644 : ! Calculate process column distribution
1645 2177315 : npcols = SIZE(pgrid, 2)
1646 4351014 : DO pdim = 0, MIN(npcols - 1, nblkcols - 1)
1647 8386795 : 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 2177315 : reuse_arrays=.TRUE.)
1662 :
1663 2177315 : 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 205215 : 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 205215 : IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1679 1421112 : ALLOCATE (matrix_set(nmatrix))
1680 1010682 : DO imatrix = 1, nmatrix
1681 1010682 : NULLIFY (matrix_set(imatrix)%matrix)
1682 : END DO
1683 205215 : 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 199954 : 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 199954 : IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1700 5714454 : ALLOCATE (matrix_set(nmatrix, mmatrix))
1701 2500336 : DO jmatrix = 1, mmatrix
1702 5114592 : DO imatrix = 1, nmatrix
1703 4914638 : NULLIFY (matrix_set(imatrix, jmatrix)%matrix)
1704 : END DO
1705 : END DO
1706 199954 : 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 328 : 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 328 : IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1724 7662 : ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix))
1725 668 : DO kmatrix = 1, pmatrix
1726 2562 : DO jmatrix = 1, mmatrix
1727 6022 : DO imatrix = 1, nmatrix
1728 5682 : NULLIFY (matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1729 : END DO
1730 : END DO
1731 : END DO
1732 328 : 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 0 : 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 0 : IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1751 0 : ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix))
1752 0 : DO lmatrix = 1, qmatrix
1753 0 : DO kmatrix = 1, pmatrix
1754 0 : DO jmatrix = 1, mmatrix
1755 0 : DO imatrix = 1, nmatrix
1756 0 : NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1757 : END DO
1758 : END DO
1759 : END DO
1760 : END DO
1761 0 : 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 0 : 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 0 : IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1784 0 : ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix, smatrix))
1785 0 : DO hmatrix = 1, smatrix
1786 0 : DO lmatrix = 1, qmatrix
1787 0 : DO kmatrix = 1, pmatrix
1788 0 : DO jmatrix = 1, mmatrix
1789 0 : DO imatrix = 1, nmatrix
1790 0 : 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 0 : 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 205364 : SUBROUTINE deallocate_dbcsr_matrix_set_1d(matrix_set)
1805 :
1806 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_set
1807 :
1808 : INTEGER :: imatrix
1809 :
1810 205364 : IF (ASSOCIATED(matrix_set)) THEN
1811 1414976 : DO imatrix = LBOUND(matrix_set, 1), UBOUND(matrix_set, 1)
1812 1009836 : CALL dbcsr_deallocate_matrix(matrix_set(imatrix)%matrix)
1813 : END DO
1814 202570 : DEALLOCATE (matrix_set)
1815 : END IF
1816 :
1817 205364 : 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 204645 : 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 204645 : IF (ASSOCIATED(matrix_set)) THEN
1832 2915390 : DO jmatrix = LBOUND(matrix_set, 2), UBOUND(matrix_set, 2)
1833 9763223 : DO imatrix = LBOUND(matrix_set, 1), UBOUND(matrix_set, 1)
1834 4940548 : CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix)%matrix)
1835 : END DO
1836 : END DO
1837 201621 : DEALLOCATE (matrix_set)
1838 : END IF
1839 204645 : 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 328 : 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 328 : IF (ASSOCIATED(matrix_set)) THEN
1854 1324 : DO kmatrix = LBOUND(matrix_set, 3), UBOUND(matrix_set, 3)
1855 3242 : DO jmatrix = LBOUND(matrix_set, 2), UBOUND(matrix_set, 2)
1856 9810 : DO imatrix = LBOUND(matrix_set, 1), UBOUND(matrix_set, 1)
1857 5682 : CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1858 : END DO
1859 : END DO
1860 : END DO
1861 328 : DEALLOCATE (matrix_set)
1862 : END IF
1863 328 : 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 0 : 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 0 : IF (ASSOCIATED(matrix_set)) THEN
1878 0 : DO lmatrix = LBOUND(matrix_set, 4), UBOUND(matrix_set, 4)
1879 0 : DO kmatrix = LBOUND(matrix_set, 3), UBOUND(matrix_set, 3)
1880 0 : DO jmatrix = LBOUND(matrix_set, 2), UBOUND(matrix_set, 2)
1881 0 : DO imatrix = LBOUND(matrix_set, 1), UBOUND(matrix_set, 1)
1882 0 : CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1883 : END DO
1884 : END DO
1885 : END DO
1886 : END DO
1887 0 : DEALLOCATE (matrix_set)
1888 : END IF
1889 0 : 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 0 : 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 0 : IF (ASSOCIATED(matrix_set)) THEN
1906 0 : DO hmatrix = LBOUND(matrix_set, 5), UBOUND(matrix_set, 5)
1907 0 : DO lmatrix = LBOUND(matrix_set, 4), UBOUND(matrix_set, 4)
1908 0 : DO kmatrix = LBOUND(matrix_set, 3), UBOUND(matrix_set, 3)
1909 0 : DO jmatrix = LBOUND(matrix_set, 2), UBOUND(matrix_set, 2)
1910 0 : DO imatrix = LBOUND(matrix_set, 1), UBOUND(matrix_set, 1)
1911 0 : 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 0 : DEALLOCATE (matrix_set)
1918 : END IF
1919 0 : END SUBROUTINE deallocate_dbcsr_matrix_set_5d
1920 :
1921 0 : END MODULE cp_dbcsr_operations
|
