Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE cp_dbcsr_contrib
9 : USE OMP_LIB, ONLY: omp_get_num_threads
10 : USE cp_dbcsr_api, ONLY: &
11 : dbcsr_clear, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_type, &
12 : dbcsr_finalize, dbcsr_get_data_size, dbcsr_get_info, dbcsr_get_num_blocks, &
13 : dbcsr_get_occupation, dbcsr_get_readonly_block_p, dbcsr_get_stored_coordinates, &
14 : dbcsr_has_symmetry, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
15 : dbcsr_iterator_readonly_start, dbcsr_iterator_start, dbcsr_iterator_stop, &
16 : dbcsr_iterator_type, dbcsr_p_type, dbcsr_put_block, dbcsr_reserve_blocks, dbcsr_type
17 : USE dbm_tests, ONLY: generate_larnv_seed
18 : USE kinds, ONLY: default_string_length,&
19 : dp
20 : USE machine, ONLY: default_output_unit
21 : USE message_passing, ONLY: mp_comm_type
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 :
27 : PUBLIC :: dbcsr_hadamard_product
28 : PUBLIC :: dbcsr_maxabs
29 : PUBLIC :: dbcsr_frobenius_norm
30 : PUBLIC :: dbcsr_gershgorin_norm
31 : PUBLIC :: dbcsr_init_random
32 : PUBLIC :: dbcsr_reserve_diag_blocks
33 : PUBLIC :: dbcsr_reserve_all_blocks
34 : PUBLIC :: dbcsr_add_on_diag
35 : PUBLIC :: dbcsr_dot, dbcsr_dot_batch
36 : PUBLIC :: dbcsr_trace
37 : PUBLIC :: dbcsr_get_block_diag
38 : PUBLIC :: dbcsr_scale_by_vector
39 : PUBLIC :: dbcsr_get_diag
40 : PUBLIC :: dbcsr_set_diag
41 : PUBLIC :: dbcsr_checksum
42 : PUBLIC :: dbcsr_print
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Hadamard product: C = A . B (C needs to be different from A and B)
48 : !> \param matrix_a ...
49 : !> \param matrix_b ...
50 : !> \param matrix_c ...
51 : ! **************************************************************************************************
52 397520 : SUBROUTINE dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c)
53 : TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b
54 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_c
55 :
56 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_hadamard_product'
57 :
58 : INTEGER :: col, handle, nblkrows_tot_a, &
59 : nblkrows_tot_b, nblkrows_tot_c, row
60 : LOGICAL :: found_b
61 79504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
62 : TYPE(dbcsr_iterator_type) :: iter
63 :
64 79504 : CALL timeset(routineN, handle)
65 79504 : CPASSERT(omp_get_num_threads() == 1)
66 :
67 79504 : CALL dbcsr_get_info(matrix_a, nblkrows_total=nblkrows_tot_a)
68 79504 : CALL dbcsr_get_info(matrix_b, nblkrows_total=nblkrows_tot_b)
69 79504 : CALL dbcsr_get_info(matrix_c, nblkrows_total=nblkrows_tot_c)
70 79504 : IF (nblkrows_tot_a /= nblkrows_tot_b .OR. nblkrows_tot_a /= nblkrows_tot_c) THEN
71 0 : CPABORT("matrices not consistent")
72 : END IF
73 :
74 79504 : CALL dbcsr_clear(matrix_c)
75 79504 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
76 179315 : DO WHILE (dbcsr_iterator_blocks_left(iter))
77 99811 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
78 99811 : CALL dbcsr_get_readonly_block_p(matrix_b, row, col, block_b, found_b)
79 179315 : IF (found_b) THEN
80 19326781 : CALL dbcsr_put_block(matrix_c, row, col, block_a*block_b)
81 : END IF
82 : END DO
83 79504 : CALL dbcsr_iterator_stop(iter)
84 79504 : CALL dbcsr_finalize(matrix_c)
85 79504 : CALL timestop(handle)
86 79504 : END SUBROUTINE dbcsr_hadamard_product
87 :
88 : ! **************************************************************************************************
89 : !> \brief Compute the maxabs norm of a dbcsr matrix
90 : !> \param matrix ...
91 : !> \return ...
92 : ! **************************************************************************************************
93 27044 : FUNCTION dbcsr_maxabs(matrix) RESULT(norm)
94 : TYPE(dbcsr_type), INTENT(IN) :: matrix
95 : REAL(dp) :: norm
96 :
97 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_maxabs'
98 :
99 : INTEGER :: handle
100 13522 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
101 : TYPE(dbcsr_iterator_type) :: iter
102 : TYPE(mp_comm_type) :: group
103 :
104 13522 : CALL timeset(routineN, handle)
105 13522 : CPASSERT(omp_get_num_threads() == 1)
106 :
107 : norm = 0.0_dp
108 13522 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix) PRIVATE(iter, block) REDUCTION(MAX:norm)
109 : CALL dbcsr_iterator_readonly_start(iter, matrix)
110 : DO WHILE (dbcsr_iterator_blocks_left(iter))
111 : CALL dbcsr_iterator_next_block(iter, block=block)
112 : norm = MAX(norm, MAXVAL(ABS(block)))
113 : END DO
114 : CALL dbcsr_iterator_stop(iter)
115 : !$OMP END PARALLEL
116 :
117 13522 : CALL dbcsr_get_info(matrix, group=group)
118 13522 : CALL group%max(norm)
119 :
120 13522 : CALL timestop(handle)
121 13522 : END FUNCTION dbcsr_maxabs
122 :
123 : ! **************************************************************************************************
124 : !> \brief Compute the frobenius norm of a dbcsr matrix
125 : !> \param matrix ...
126 : !> \return ...
127 : ! **************************************************************************************************
128 388756 : FUNCTION dbcsr_frobenius_norm(matrix) RESULT(norm)
129 : TYPE(dbcsr_type), INTENT(IN) :: matrix
130 : REAL(dp) :: norm
131 :
132 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_frobenius_norm'
133 :
134 : INTEGER :: col, handle, row
135 : LOGICAL :: has_symmetry
136 194378 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
137 : TYPE(dbcsr_iterator_type) :: iter
138 : TYPE(mp_comm_type) :: group
139 :
140 194378 : CALL timeset(routineN, handle)
141 194378 : CPASSERT(omp_get_num_threads() == 1)
142 :
143 194378 : has_symmetry = dbcsr_has_symmetry(matrix)
144 : norm = 0.0_dp
145 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix, has_symmetry) &
146 194378 : !$OMP PRIVATE(iter, row, col, block) REDUCTION(+:norm)
147 : CALL dbcsr_iterator_readonly_start(iter, matrix)
148 : DO WHILE (dbcsr_iterator_blocks_left(iter))
149 : CALL dbcsr_iterator_next_block(iter, row, col, block)
150 : IF (has_symmetry .AND. row /= col) THEN
151 : norm = norm + 2.0_dp*SUM(block**2)
152 : ELSE
153 : norm = norm + SUM(block**2)
154 : END IF
155 : END DO
156 : CALL dbcsr_iterator_stop(iter)
157 : !$OMP END PARALLEL
158 :
159 194378 : CALL dbcsr_get_info(matrix, group=group)
160 194378 : CALL group%sum(norm)
161 194378 : norm = SQRT(norm)
162 :
163 194378 : CALL timestop(handle)
164 194378 : END FUNCTION dbcsr_frobenius_norm
165 :
166 : ! **************************************************************************************************
167 : !> \brief Compute the gershgorin norm of a dbcsr matrix
168 : !> \param matrix ...
169 : !> \return ...
170 : ! **************************************************************************************************
171 9906 : FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm)
172 : TYPE(dbcsr_type), INTENT(IN) :: matrix
173 : REAL(dp) :: norm
174 :
175 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_gershgorin_norm'
176 :
177 : INTEGER :: col, col_offset, handle, i, j, ncol, &
178 : nrow, row, row_offset
179 : LOGICAL :: has_symmetry
180 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buffer
181 9906 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
182 : TYPE(dbcsr_iterator_type) :: iter
183 : TYPE(mp_comm_type) :: group
184 :
185 9906 : CALL timeset(routineN, handle)
186 9906 : CPASSERT(omp_get_num_threads() == 1)
187 :
188 9906 : has_symmetry = dbcsr_has_symmetry(matrix)
189 9906 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
190 9906 : CPASSERT(nrow == ncol)
191 29710 : ALLOCATE (buffer(nrow))
192 9906 : buffer = 0.0_dp
193 :
194 9906 : CALL dbcsr_iterator_readonly_start(iter, matrix)
195 204271 : DO WHILE (dbcsr_iterator_blocks_left(iter))
196 194365 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_offset=row_offset, col_offset=col_offset)
197 650608 : DO j = 1, SIZE(block, 2)
198 4951087 : DO i = 1, SIZE(block, 1)
199 4310385 : buffer(row_offset + i - 1) = buffer(row_offset + i - 1) + ABS(block(i, j))
200 4756722 : IF (has_symmetry .AND. row /= col) THEN
201 53764 : buffer(col_offset + j - 1) = buffer(col_offset + j - 1) + ABS(block(i, j))
202 : END IF
203 : END DO
204 : END DO
205 : END DO
206 9906 : CALL dbcsr_iterator_stop(iter)
207 :
208 9906 : CALL dbcsr_get_info(matrix, group=group)
209 9906 : CALL group%sum(buffer)
210 159604 : norm = MAXVAL(buffer)
211 9906 : DEALLOCATE (buffer)
212 :
213 9906 : CALL timestop(handle)
214 39624 : END FUNCTION dbcsr_gershgorin_norm
215 :
216 : ! **************************************************************************************************
217 : !> \brief Fills the given matrix with random numbers.
218 : !> \param matrix ...
219 : !> \param keep_sparsity ...
220 : ! **************************************************************************************************
221 1302 : SUBROUTINE dbcsr_init_random(matrix, keep_sparsity)
222 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
223 : LOGICAL, OPTIONAL :: keep_sparsity
224 :
225 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_init_random'
226 :
227 : INTEGER :: col, col_size, handle, ncol, nrow, row, &
228 : row_size
229 : INTEGER, DIMENSION(4) :: iseed
230 : LOGICAL :: my_keep_sparsity
231 434 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
232 : TYPE(dbcsr_iterator_type) :: iter
233 :
234 434 : CALL timeset(routineN, handle)
235 434 : CPASSERT(omp_get_num_threads() == 1)
236 :
237 434 : my_keep_sparsity = .FALSE.
238 434 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
239 434 : IF (.NOT. my_keep_sparsity) CALL dbcsr_reserve_all_blocks(matrix)
240 434 : CALL dbcsr_get_info(matrix, nblkrows_total=nrow, nblkcols_total=ncol)
241 :
242 434 : CALL dbcsr_iterator_start(iter, matrix)
243 982 : DO WHILE (dbcsr_iterator_blocks_left(iter))
244 548 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
245 : ! set the seed for dlarnv, is here to guarantee same value of the random numbers
246 : ! for all layouts (and block distributions)
247 548 : iseed = generate_larnv_seed(row, nrow, col, ncol, 1)
248 982 : CALL dlarnv(1, iseed, row_size*col_size, block(1, 1))
249 : END DO
250 434 : CALL dbcsr_iterator_stop(iter)
251 :
252 434 : CALL timestop(handle)
253 434 : END SUBROUTINE dbcsr_init_random
254 :
255 : ! **************************************************************************************************
256 : !> \brief Reserves all diagonal blocks.
257 : !> \param matrix ...
258 : ! **************************************************************************************************
259 356508 : SUBROUTINE dbcsr_reserve_diag_blocks(matrix)
260 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
261 :
262 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reserve_diag_blocks'
263 :
264 : INTEGER :: handle, i, k, mynode, nblkcols_total, &
265 : nblkrows_total, owner
266 356508 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_diag
267 356508 : INTEGER, DIMENSION(:), POINTER :: local_rows
268 : TYPE(dbcsr_distribution_type) :: dist
269 :
270 356508 : CALL timeset(routineN, handle)
271 356508 : CPASSERT(omp_get_num_threads() == 1)
272 :
273 356508 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total)
274 356508 : CPASSERT(nblkrows_total == nblkcols_total)
275 :
276 356508 : CALL dbcsr_get_info(matrix, local_rows=local_rows, distribution=dist)
277 356508 : CALL dbcsr_distribution_get(dist, mynode=mynode)
278 990645 : ALLOCATE (local_diag(SIZE(local_rows)))
279 :
280 1088315 : k = 0
281 1088315 : DO i = 1, SIZE(local_rows)
282 731807 : CALL dbcsr_get_stored_coordinates(matrix, row=local_rows(i), column=local_rows(i), processor=owner)
283 1088315 : IF (owner == mynode) THEN
284 731807 : k = k + 1
285 731807 : local_diag(k) = local_rows(i)
286 : END IF
287 : END DO
288 :
289 356508 : CALL dbcsr_reserve_blocks(matrix, rows=local_diag(1:k), cols=local_diag(1:k))
290 356508 : DEALLOCATE (local_diag)
291 :
292 356508 : CALL timestop(handle)
293 1069524 : END SUBROUTINE dbcsr_reserve_diag_blocks
294 :
295 : ! **************************************************************************************************
296 : !> \brief Reserves all blocks.
297 : !> \param matrix ...
298 : ! **************************************************************************************************
299 1971573 : SUBROUTINE dbcsr_reserve_all_blocks(matrix)
300 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
301 :
302 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reserve_all_blocks'
303 :
304 : INTEGER :: handle, i, j, k, n
305 1971573 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
306 1971573 : INTEGER, DIMENSION(:), POINTER :: local_cols, local_rows
307 :
308 1971573 : CALL timeset(routineN, handle)
309 1971573 : CPASSERT(omp_get_num_threads() == 1)
310 :
311 1971573 : CALL dbcsr_get_info(matrix, local_rows=local_rows, local_cols=local_cols)
312 1971573 : n = SIZE(local_rows)*SIZE(local_cols)
313 7530506 : ALLOCATE (rows(n), cols(n))
314 :
315 4550754 : k = 0
316 4550754 : DO i = 1, SIZE(local_rows)
317 9123586 : DO j = 1, SIZE(local_cols)
318 4572832 : k = k + 1
319 4572832 : rows(k) = local_rows(i)
320 7152013 : cols(k) = local_cols(j)
321 : END DO
322 : END DO
323 :
324 1971573 : CALL dbcsr_reserve_blocks(matrix, rows=rows(1:k), cols=cols(1:k))
325 1971573 : DEALLOCATE (rows, cols)
326 :
327 1971573 : CALL timestop(handle)
328 1971573 : END SUBROUTINE dbcsr_reserve_all_blocks
329 :
330 : ! **************************************************************************************************
331 : !> \brief Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
332 : !> \param matrix ...
333 : !> \param alpha ...
334 : ! **************************************************************************************************
335 654172 : SUBROUTINE dbcsr_add_on_diag(matrix, alpha)
336 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
337 : REAL(kind=dp), INTENT(IN) :: alpha
338 :
339 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag'
340 :
341 : INTEGER :: col, col_size, handle, i, row, row_size
342 327086 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
343 : TYPE(dbcsr_iterator_type) :: iter
344 :
345 327086 : CALL timeset(routineN, handle)
346 327086 : CPASSERT(omp_get_num_threads() == 1)
347 :
348 327086 : CALL dbcsr_reserve_diag_blocks(matrix)
349 :
350 327086 : CALL dbcsr_iterator_start(iter, matrix)
351 6024907 : DO WHILE (dbcsr_iterator_blocks_left(iter))
352 5697821 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
353 6024907 : IF (row == col) THEN
354 675373 : CPASSERT(row_size == col_size)
355 3103721 : DO i = 1, row_size
356 3103721 : block(i, i) = block(i, i) + alpha
357 : END DO
358 : END IF
359 : END DO
360 327086 : CALL dbcsr_iterator_stop(iter)
361 :
362 327086 : CALL timestop(handle)
363 327086 : END SUBROUTINE dbcsr_add_on_diag
364 :
365 : ! **************************************************************************************************
366 : !> \brief Computes the dot product of two matrices, also known as the trace of their matrix product.
367 : !> \param matrix_a ...
368 : !> \param matrix_b ...
369 : !> \param trace ...
370 : ! **************************************************************************************************
371 783456 : SUBROUTINE dbcsr_dot(matrix_a, matrix_b, trace)
372 : TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_a, matrix_b
373 : REAL(KIND=dp), INTENT(OUT) :: trace
374 :
375 : REAL(KIND=dp) :: values(1)
376 3133824 : TYPE(dbcsr_p_type) :: a(1), b(1)
377 :
378 783456 : a(1)%matrix => matrix_a
379 783456 : b(1)%matrix => matrix_b
380 783456 : CALL dbcsr_dot_batch(a, b, values)
381 783456 : trace = values(1)
382 783456 : END SUBROUTINE dbcsr_dot
383 :
384 : ! **************************************************************************************************
385 : !> \brief Evaluate several dot products in one OpenMP region; MPI reductions remain outside it.
386 : !> \param matrix_a first operands, in reduction order
387 : !> \param matrix_b second operands
388 : !> \param trace dot products
389 : !> \param comm optional common communicator; all operands must use this exact communicator
390 : ! **************************************************************************************************
391 1586205 : SUBROUTINE dbcsr_dot_batch(matrix_a, matrix_b, trace, comm)
392 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_a(:), matrix_b(:)
393 : REAL(KIND=dp), INTENT(OUT) :: trace(SIZE(matrix_a))
394 :
395 : CLASS(mp_comm_type), INTENT(IN), OPTIONAL :: comm
396 :
397 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_dot_batch'
398 :
399 : INTEGER :: col, handle, i, row
400 : LOGICAL :: found_b, has_symmetry
401 1586205 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
402 : TYPE(dbcsr_iterator_type) :: iter
403 : TYPE(mp_comm_type) :: group
404 :
405 1586205 : CALL timeset(routineN, handle)
406 1586205 : CPASSERT(omp_get_num_threads() == 1)
407 :
408 1586205 : CPASSERT(SIZE(matrix_a) == SIZE(matrix_b))
409 8407806 : DO i = 1, SIZE(matrix_a)
410 6821601 : has_symmetry = dbcsr_has_symmetry(matrix_a(i)%matrix)
411 6821601 : CPASSERT(has_symmetry .EQV. dbcsr_has_symmetry(matrix_b(i)%matrix))
412 8407806 : IF (PRESENT(comm)) THEN
413 6038145 : CALL dbcsr_get_info(matrix_a(i)%matrix, group=group)
414 6038145 : CPASSERT(group == comm)
415 : END IF
416 : END DO
417 :
418 8407806 : trace = 0.0_dp
419 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix_a, matrix_b) &
420 1586205 : !$OMP PRIVATE(i, iter, row, col, block_a, block_b, found_b, has_symmetry) REDUCTION(+:trace)
421 : DO i = 1, SIZE(matrix_a)
422 : has_symmetry = dbcsr_has_symmetry(matrix_a(i)%matrix)
423 : CALL dbcsr_iterator_readonly_start(iter, matrix_a(i)%matrix)
424 : DO WHILE (dbcsr_iterator_blocks_left(iter))
425 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
426 : IF (SIZE(block_a) == 0) CYCLE ! Skip zero-sized blocks.
427 : CALL dbcsr_get_readonly_block_p(matrix_b(i)%matrix, row, col, block_b, found_b)
428 : IF (found_b) THEN
429 : IF (has_symmetry .AND. row /= col) THEN
430 : trace(i) = trace(i) + 2.0_dp*SUM(block_a*block_b)
431 : ELSE
432 : trace(i) = trace(i) + SUM(block_a*block_b)
433 : END IF
434 : END IF
435 : END DO
436 : CALL dbcsr_iterator_stop(iter)
437 : END DO
438 : !$OMP END PARALLEL
439 :
440 1586205 : IF (PRESENT(comm)) THEN
441 1585075 : IF (SIZE(trace) > 0 .AND. comm%num_pe > 1) CALL comm%sum(trace)
442 : ELSE
443 1566912 : DO i = 1, SIZE(matrix_a)
444 783456 : CALL dbcsr_get_info(matrix_a(i)%matrix, group=group)
445 1566912 : IF (group%num_pe > 1) CALL group%sum(trace(i))
446 : END DO
447 : END IF
448 :
449 1586205 : CALL timestop(handle)
450 1586205 : END SUBROUTINE dbcsr_dot_batch
451 :
452 : ! **************************************************************************************************
453 : !> \brief Computes the trace of the given matrix, also known as the sum of its diagonal elements.
454 : !> \param matrix ...
455 : !> \param trace ...
456 : ! **************************************************************************************************
457 79977 : SUBROUTINE dbcsr_trace(matrix, trace)
458 : TYPE(dbcsr_type), INTENT(IN) :: matrix
459 : REAL(KIND=dp), INTENT(OUT) :: trace
460 :
461 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace'
462 :
463 : INTEGER :: col, col_size, handle, i, row, row_size
464 26659 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
465 : TYPE(dbcsr_iterator_type) :: iter
466 : TYPE(mp_comm_type) :: group
467 :
468 26659 : CALL timeset(routineN, handle)
469 26659 : CPASSERT(omp_get_num_threads() == 1)
470 :
471 26659 : trace = 0.0_dp
472 26659 : CALL dbcsr_iterator_readonly_start(iter, matrix)
473 1648909 : DO WHILE (dbcsr_iterator_blocks_left(iter))
474 1622250 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
475 1648909 : IF (row == col) THEN
476 144064 : CPASSERT(row_size == col_size)
477 617441 : DO i = 1, row_size
478 617441 : trace = trace + block(i, i)
479 : END DO
480 : END IF
481 : END DO
482 26659 : CALL dbcsr_iterator_stop(iter)
483 :
484 26659 : CALL dbcsr_get_info(matrix, group=group)
485 26659 : CALL group%sum(trace)
486 :
487 26659 : CALL timestop(handle)
488 26659 : END SUBROUTINE dbcsr_trace
489 :
490 : ! **************************************************************************************************
491 : !> \brief Copies the diagonal blocks of matrix into diag.
492 : !> \param matrix ...
493 : !> \param diag ...
494 : ! **************************************************************************************************
495 232940 : SUBROUTINE dbcsr_get_block_diag(matrix, diag)
496 : TYPE(dbcsr_type), INTENT(IN) :: matrix
497 : TYPE(dbcsr_type), INTENT(INOUT) :: diag
498 :
499 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_block_diag'
500 :
501 : CHARACTER(len=default_string_length) :: name
502 : INTEGER :: col, handle, row
503 116470 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
504 : TYPE(dbcsr_iterator_type) :: iter
505 :
506 116470 : CALL timeset(routineN, handle)
507 116470 : CPASSERT(omp_get_num_threads() == 1)
508 :
509 116470 : CALL dbcsr_get_info(matrix, name=name)
510 116470 : CALL dbcsr_create(diag, template=matrix, name='diag of '//name)
511 :
512 116470 : CALL dbcsr_iterator_readonly_start(iter, matrix)
513 8649083 : DO WHILE (dbcsr_iterator_blocks_left(iter))
514 8532613 : CALL dbcsr_iterator_next_block(iter, row, col, block)
515 8649083 : IF (row == col) THEN
516 474651 : CALL dbcsr_put_block(diag, row, col, block)
517 : END IF
518 : END DO
519 116470 : CALL dbcsr_iterator_stop(iter)
520 116470 : CALL dbcsr_finalize(diag)
521 :
522 116470 : CALL timestop(handle)
523 116470 : END SUBROUTINE dbcsr_get_block_diag
524 :
525 : ! **************************************************************************************************
526 : !> \brief Scales the rows/columns of given matrix.
527 : !> \param matrix Matrix to be scaled in-place.
528 : !> \param alpha Vector with scaling factors.
529 : !> \param side Side from which to apply the vector. Allowed values are 'right' and 'left'.
530 : ! **************************************************************************************************
531 227399 : SUBROUTINE dbcsr_scale_by_vector(matrix, alpha, side)
532 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
533 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha
534 : CHARACTER(LEN=*), INTENT(IN) :: side
535 :
536 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_by_vector'
537 :
538 : INTEGER :: col_offset, col_size, handle, i, &
539 : nfullcols_total, nfullrows_total, &
540 : row_offset, row_size
541 : LOGICAL :: right
542 227399 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
543 : TYPE(dbcsr_iterator_type) :: iter
544 :
545 227399 : CALL timeset(routineN, handle)
546 227399 : CPASSERT(omp_get_num_threads() == 1)
547 :
548 227399 : IF (side == 'right') THEN
549 : right = .TRUE.
550 0 : ELSE IF (side == 'left') THEN
551 : right = .FALSE.
552 : ELSE
553 0 : CPABORT("Unknown side: "//TRIM(side))
554 : END IF
555 :
556 : ! Check that alpha and matrix have matching sizes.
557 227399 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
558 227399 : IF (right) THEN
559 227399 : CPASSERT(nfullcols_total == SIZE(alpha))
560 : ELSE
561 0 : CPASSERT(nfullrows_total == SIZE(alpha))
562 : END IF
563 :
564 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix, alpha, right) &
565 227399 : !$OMP PRIVATE(iter, block, row_size, col_size, row_offset, col_offset, i)
566 : CALL dbcsr_iterator_start(iter, matrix)
567 : DO WHILE (dbcsr_iterator_blocks_left(iter))
568 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
569 : row_offset=row_offset, col_offset=col_offset)
570 : IF (SIZE(block) == 0) CYCLE ! Skip zero-sized blocks.
571 : IF (right) THEN
572 : DO i = 1, col_size
573 : block(:, i) = block(:, i)*alpha(col_offset + i - 1)
574 : END DO
575 : ELSE
576 : DO i = 1, row_size
577 : block(i, :) = block(i, :)*alpha(row_offset + i - 1)
578 : END DO
579 : END IF
580 : END DO
581 : CALL dbcsr_iterator_stop(iter)
582 : !$OMP END PARALLEL
583 :
584 227399 : CALL timestop(handle)
585 227399 : END SUBROUTINE dbcsr_scale_by_vector
586 :
587 : ! **************************************************************************************************
588 : !> \brief Copies the diagonal elements from the given matrix into the given array.
589 : !> \param matrix ...
590 : !> \param diag ...
591 : ! **************************************************************************************************
592 11356 : SUBROUTINE dbcsr_get_diag(matrix, diag)
593 : TYPE(dbcsr_type), INTENT(IN) :: matrix
594 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
595 :
596 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag'
597 :
598 : INTEGER :: col, col_size, handle, i, &
599 : nfullcols_total, nfullrows_total, row, &
600 : row_offset, row_size
601 11356 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
602 : TYPE(dbcsr_iterator_type) :: iter
603 :
604 11356 : CALL timeset(routineN, handle)
605 11356 : CPASSERT(omp_get_num_threads() == 1)
606 :
607 11356 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
608 11356 : CPASSERT(nfullrows_total == nfullcols_total)
609 11356 : CPASSERT(nfullrows_total == SIZE(diag))
610 :
611 1229026 : diag(:) = 0.0_dp
612 11356 : CALL dbcsr_iterator_readonly_start(iter, matrix)
613 279568 : DO WHILE (dbcsr_iterator_blocks_left(iter))
614 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
615 268212 : col_size=col_size, row_offset=row_offset)
616 279568 : IF (row == col) THEN
617 34615 : CPASSERT(row_size == col_size)
618 646412 : DO i = 1, row_size
619 646412 : diag(row_offset + i - 1) = block(i, i)
620 : END DO
621 : END IF
622 : END DO
623 11356 : CALL dbcsr_iterator_stop(iter)
624 :
625 11356 : CALL timestop(handle)
626 11356 : END SUBROUTINE dbcsr_get_diag
627 :
628 : ! **************************************************************************************************
629 : !> \brief Copies the diagonal elements from the given array into the given matrix.
630 : !> \param matrix ...
631 : !> \param diag ...
632 : ! **************************************************************************************************
633 6020 : SUBROUTINE dbcsr_set_diag(matrix, diag)
634 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
635 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag
636 :
637 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag'
638 :
639 : INTEGER :: col, col_size, handle, i, &
640 : nfullcols_total, nfullrows_total, row, &
641 : row_offset, row_size
642 6020 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
643 : TYPE(dbcsr_iterator_type) :: iter
644 :
645 6020 : CALL timeset(routineN, handle)
646 6020 : CPASSERT(omp_get_num_threads() == 1)
647 :
648 6020 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
649 6020 : CPASSERT(nfullrows_total == nfullcols_total)
650 6020 : CPASSERT(nfullrows_total == SIZE(diag))
651 :
652 6020 : CALL dbcsr_reserve_diag_blocks(matrix)
653 :
654 6020 : CALL dbcsr_iterator_start(iter, matrix)
655 153759 : DO WHILE (dbcsr_iterator_blocks_left(iter))
656 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
657 147739 : col_size=col_size, row_offset=row_offset)
658 153759 : IF (row == col) THEN
659 16572 : CPASSERT(row_size == col_size)
660 95784 : DO i = 1, row_size
661 95784 : block(i, i) = diag(row_offset + i - 1)
662 : END DO
663 : END IF
664 : END DO
665 6020 : CALL dbcsr_iterator_stop(iter)
666 :
667 6020 : CALL timestop(handle)
668 6020 : END SUBROUTINE dbcsr_set_diag
669 :
670 : ! **************************************************************************************************
671 : !> \brief Calculates the checksum of a DBCSR matrix.
672 : !> \param matrix ...
673 : !> \param pos Enable position-dependent checksum.
674 : !> \return ...
675 : ! **************************************************************************************************
676 55680 : FUNCTION dbcsr_checksum(matrix, pos) RESULT(checksum)
677 :
678 : TYPE(dbcsr_type), INTENT(IN) :: matrix
679 : LOGICAL, INTENT(IN), OPTIONAL :: pos
680 : REAL(KIND=dp) :: checksum
681 :
682 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_checksum'
683 :
684 : INTEGER :: col_offset, col_size, handle, i, j, &
685 : row_offset, row_size
686 : LOGICAL :: my_pos
687 : REAL(KIND=dp) :: position_factor
688 18560 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
689 : TYPE(dbcsr_iterator_type) :: iter
690 : TYPE(mp_comm_type) :: group
691 :
692 18560 : CALL timeset(routineN, handle)
693 18560 : CPASSERT(omp_get_num_threads() == 1)
694 :
695 18560 : my_pos = .FALSE.
696 18560 : IF (PRESENT(pos)) THEN
697 286 : my_pos = pos
698 : END IF
699 :
700 18560 : checksum = 0.0_dp
701 18560 : CALL dbcsr_iterator_readonly_start(iter, matrix)
702 115034 : DO WHILE (dbcsr_iterator_blocks_left(iter))
703 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
704 96474 : row_offset=row_offset, col_offset=col_offset)
705 115034 : IF (my_pos) THEN
706 21372 : DO i = 1, row_size
707 218025 : DO j = 1, col_size
708 196653 : position_factor = LOG(REAL((row_offset + i - 1)*(col_offset + j - 1), KIND=dp))
709 215930 : checksum = checksum + block(i, j)*position_factor
710 : END DO
711 : END DO
712 : ELSE
713 4419117 : checksum = checksum + SUM(block**2)
714 : END IF
715 : END DO
716 18560 : CALL dbcsr_iterator_stop(iter)
717 :
718 18560 : CALL dbcsr_get_info(matrix, group=group)
719 18560 : CALL group%sum(checksum)
720 :
721 18560 : CALL timestop(handle)
722 18560 : END FUNCTION dbcsr_checksum
723 :
724 : ! **************************************************************************************************
725 : !> \brief Prints given matrix in matlab format (only present blocks).
726 : !> \param matrix ...
727 : !> \param variable_name ...
728 : !> \param unit_nr ...
729 : ! **************************************************************************************************
730 0 : SUBROUTINE dbcsr_print(matrix, variable_name, unit_nr)
731 : TYPE(dbcsr_type), INTENT(IN) :: matrix
732 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: variable_name
733 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
734 :
735 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_print'
736 :
737 : CHARACTER(len=default_string_length) :: my_variable_name, name
738 : INTEGER :: col_offset, col_size, handle, i, iw, j, nblkcols_total, nblkrows_total, &
739 : nfullcols_total, nfullrows_total, row_offset, row_size
740 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
741 : TYPE(dbcsr_iterator_type) :: iter
742 :
743 0 : CALL timeset(routineN, handle)
744 0 : CPASSERT(omp_get_num_threads() == 1)
745 :
746 0 : iw = default_output_unit
747 0 : IF (PRESENT(unit_nr)) iw = unit_nr
748 :
749 0 : my_variable_name = 'a'
750 0 : IF (PRESENT(variable_name)) my_variable_name = variable_name
751 :
752 : ! Print matrix properties.
753 : CALL dbcsr_get_info(matrix, name=name, &
754 : nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total, &
755 0 : nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
756 0 : WRITE (iw, *) "===", routineN, "==="
757 0 : WRITE (iw, *) "Name:", name
758 0 : WRITE (iw, *) "Symmetry:", dbcsr_has_symmetry(matrix)
759 0 : WRITE (iw, *) "Number of blocks:", dbcsr_get_num_blocks(matrix)
760 0 : WRITE (iw, *) "Data size:", dbcsr_get_data_size(matrix)
761 0 : WRITE (iw, *) "Occupation:", dbcsr_get_occupation(matrix)
762 0 : WRITE (iw, *) "Full size:", nfullrows_total, "x", nfullcols_total
763 0 : WRITE (iw, *) "Blocked size:", nblkrows_total, "x", nblkcols_total
764 :
765 : ! Print matrix blocks.
766 0 : CALL dbcsr_iterator_readonly_start(iter, matrix)
767 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
768 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
769 0 : row_offset=row_offset, col_offset=col_offset)
770 0 : DO i = 1, row_size
771 0 : DO j = 1, col_size
772 0 : WRITE (iw, '(A,I4,A,I4,A,E23.16,A)') TRIM(my_variable_name)//'(', &
773 0 : row_offset + i - 1, ',', col_offset + j - 1, ')=', block(i, j), ';'
774 : END DO
775 : END DO
776 : END DO
777 0 : CALL dbcsr_iterator_stop(iter)
778 :
779 0 : CALL timestop(handle)
780 0 : END SUBROUTINE dbcsr_print
781 :
782 : END MODULE cp_dbcsr_contrib
|
