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