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