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_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
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 426340 : 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 85268 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
62 : TYPE(dbcsr_iterator_type) :: iter
63 :
64 85268 : CALL timeset(routineN, handle)
65 85268 : CPASSERT(omp_get_num_threads() == 1)
66 :
67 85268 : CALL dbcsr_get_info(matrix_a, nblkrows_total=nblkrows_tot_a)
68 85268 : CALL dbcsr_get_info(matrix_b, nblkrows_total=nblkrows_tot_b)
69 85268 : CALL dbcsr_get_info(matrix_c, nblkrows_total=nblkrows_tot_c)
70 85268 : 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 85268 : CALL dbcsr_clear(matrix_c)
75 85268 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
76 191107 : DO WHILE (dbcsr_iterator_blocks_left(iter))
77 105839 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
78 105839 : CALL dbcsr_get_readonly_block_p(matrix_b, row, col, block_b, found_b)
79 191107 : IF (found_b) THEN
80 19520289 : CALL dbcsr_put_block(matrix_c, row, col, block_a*block_b)
81 : END IF
82 : END DO
83 85268 : CALL dbcsr_iterator_stop(iter)
84 85268 : CALL dbcsr_finalize(matrix_c)
85 85268 : CALL timestop(handle)
86 85268 : 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 396194 : 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 198097 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
137 : TYPE(dbcsr_iterator_type) :: iter
138 : TYPE(mp_comm_type) :: group
139 :
140 198097 : CALL timeset(routineN, handle)
141 198097 : CPASSERT(omp_get_num_threads() == 1)
142 :
143 198097 : has_symmetry = dbcsr_has_symmetry(matrix)
144 : norm = 0.0_dp
145 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix, has_symmetry) &
146 198097 : !$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 198097 : CALL dbcsr_get_info(matrix, group=group)
160 198097 : CALL group%sum(norm)
161 198097 : norm = SQRT(norm)
162 :
163 198097 : CALL timestop(handle)
164 198097 : END FUNCTION dbcsr_frobenius_norm
165 :
166 : ! **************************************************************************************************
167 : !> \brief Compute the gershgorin norm of a dbcsr matrix
168 : !> \param matrix ...
169 : !> \return ...
170 : ! **************************************************************************************************
171 9321 : 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 9321 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
182 : TYPE(dbcsr_iterator_type) :: iter
183 : TYPE(mp_comm_type) :: group
184 :
185 9321 : CALL timeset(routineN, handle)
186 9321 : CPASSERT(omp_get_num_threads() == 1)
187 :
188 9321 : has_symmetry = dbcsr_has_symmetry(matrix)
189 9321 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
190 9321 : CPASSERT(nrow == ncol)
191 27955 : ALLOCATE (buffer(nrow))
192 9321 : buffer = 0.0_dp
193 :
194 9321 : CALL dbcsr_iterator_readonly_start(iter, matrix)
195 203472 : DO WHILE (dbcsr_iterator_blocks_left(iter))
196 194151 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_offset=row_offset, col_offset=col_offset)
197 647599 : DO j = 1, SIZE(block, 2)
198 4897901 : DO i = 1, SIZE(block, 1)
199 4259623 : buffer(row_offset + i - 1) = buffer(row_offset + i - 1) + ABS(block(i, j))
200 4703750 : 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 9321 : CALL dbcsr_iterator_stop(iter)
207 :
208 9321 : CALL dbcsr_get_info(matrix, group=group)
209 9321 : CALL group%sum(buffer)
210 157394 : norm = MAXVAL(buffer)
211 9321 : DEALLOCATE (buffer)
212 :
213 9321 : CALL timestop(handle)
214 37284 : 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 1350 : 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 450 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
232 : TYPE(dbcsr_iterator_type) :: iter
233 :
234 450 : CALL timeset(routineN, handle)
235 450 : CPASSERT(omp_get_num_threads() == 1)
236 :
237 450 : my_keep_sparsity = .FALSE.
238 450 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
239 450 : IF (.NOT. my_keep_sparsity) CALL dbcsr_reserve_all_blocks(matrix)
240 450 : CALL dbcsr_get_info(matrix, nblkrows_total=nrow, nblkcols_total=ncol)
241 :
242 450 : CALL dbcsr_iterator_start(iter, matrix)
243 1006 : DO WHILE (dbcsr_iterator_blocks_left(iter))
244 556 : 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 556 : iseed = generate_larnv_seed(row, nrow, col, ncol, 1)
248 1006 : CALL dlarnv(1, iseed, row_size*col_size, block(1, 1))
249 : END DO
250 450 : CALL dbcsr_iterator_stop(iter)
251 :
252 450 : CALL timestop(handle)
253 450 : END SUBROUTINE dbcsr_init_random
254 :
255 : ! **************************************************************************************************
256 : !> \brief Reserves all diagonal blocks.
257 : !> \param matrix ...
258 : ! **************************************************************************************************
259 358191 : 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 358191 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_diag
267 358191 : INTEGER, DIMENSION(:), POINTER :: local_rows
268 : TYPE(dbcsr_distribution_type) :: dist
269 :
270 358191 : CALL timeset(routineN, handle)
271 358191 : CPASSERT(omp_get_num_threads() == 1)
272 :
273 358191 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total)
274 358191 : CPASSERT(nblkrows_total == nblkcols_total)
275 :
276 358191 : CALL dbcsr_get_info(matrix, local_rows=local_rows, distribution=dist)
277 358191 : CALL dbcsr_distribution_get(dist, mynode=mynode)
278 996719 : ALLOCATE (local_diag(SIZE(local_rows)))
279 :
280 1092850 : k = 0
281 1092850 : DO i = 1, SIZE(local_rows)
282 734659 : CALL dbcsr_get_stored_coordinates(matrix, row=local_rows(i), column=local_rows(i), processor=owner)
283 1092850 : IF (owner == mynode) THEN
284 734659 : k = k + 1
285 734659 : local_diag(k) = local_rows(i)
286 : END IF
287 : END DO
288 :
289 358191 : CALL dbcsr_reserve_blocks(matrix, rows=local_diag(1:k), cols=local_diag(1:k))
290 358191 : DEALLOCATE (local_diag)
291 :
292 358191 : CALL timestop(handle)
293 1074573 : END SUBROUTINE dbcsr_reserve_diag_blocks
294 :
295 : ! **************************************************************************************************
296 : !> \brief Reserves all blocks.
297 : !> \param matrix ...
298 : ! **************************************************************************************************
299 1968914 : 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 1968914 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
306 1968914 : INTEGER, DIMENSION(:), POINTER :: local_cols, local_rows
307 :
308 1968914 : CALL timeset(routineN, handle)
309 1968914 : CPASSERT(omp_get_num_threads() == 1)
310 :
311 1968914 : CALL dbcsr_get_info(matrix, local_rows=local_rows, local_cols=local_cols)
312 1968914 : n = SIZE(local_rows)*SIZE(local_cols)
313 7522490 : ALLOCATE (rows(n), cols(n))
314 :
315 4575676 : k = 0
316 4575676 : DO i = 1, SIZE(local_rows)
317 9236137 : DO j = 1, SIZE(local_cols)
318 4660461 : k = k + 1
319 4660461 : rows(k) = local_rows(i)
320 7267223 : cols(k) = local_cols(j)
321 : END DO
322 : END DO
323 :
324 1968914 : CALL dbcsr_reserve_blocks(matrix, rows=rows(1:k), cols=cols(1:k))
325 1968914 : DEALLOCATE (rows, cols)
326 :
327 1968914 : CALL timestop(handle)
328 1968914 : 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 657554 : 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 328777 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
343 : TYPE(dbcsr_iterator_type) :: iter
344 :
345 328777 : CALL timeset(routineN, handle)
346 328777 : CPASSERT(omp_get_num_threads() == 1)
347 :
348 328777 : CALL dbcsr_reserve_diag_blocks(matrix)
349 :
350 328777 : CALL dbcsr_iterator_start(iter, matrix)
351 6030328 : DO WHILE (dbcsr_iterator_blocks_left(iter))
352 5701551 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
353 6030328 : IF (row == col) THEN
354 678237 : CPASSERT(row_size == col_size)
355 3144545 : DO i = 1, row_size
356 3144545 : block(i, i) = block(i, i) + alpha
357 : END DO
358 : END IF
359 : END DO
360 328777 : CALL dbcsr_iterator_stop(iter)
361 :
362 328777 : CALL timestop(handle)
363 328777 : 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 13299802 : SUBROUTINE dbcsr_dot(matrix_a, matrix_b, trace)
372 : TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b
373 : REAL(KIND=dp), INTENT(OUT) :: trace
374 :
375 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_dot'
376 :
377 : INTEGER :: col, handle, row
378 : LOGICAL :: found_b, has_symmetry
379 6649901 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
380 : TYPE(dbcsr_iterator_type) :: iter
381 : TYPE(mp_comm_type) :: group
382 :
383 6649901 : CALL timeset(routineN, handle)
384 6649901 : CPASSERT(omp_get_num_threads() == 1)
385 :
386 6649901 : CPASSERT(dbcsr_has_symmetry(matrix_a) .EQV. dbcsr_has_symmetry(matrix_b))
387 6649901 : has_symmetry = dbcsr_has_symmetry(matrix_a)
388 :
389 6649901 : trace = 0.0_dp
390 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix_a, matrix_b, has_symmetry) &
391 6649901 : !$OMP PRIVATE(iter, row, col, block_a, block_b, found_b) REDUCTION(+:trace)
392 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
393 : DO WHILE (dbcsr_iterator_blocks_left(iter))
394 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
395 : IF (SIZE(block_a) == 0) CYCLE ! Skip zero-sized blocks.
396 : CALL dbcsr_get_readonly_block_p(matrix_b, row, col, block_b, found_b)
397 : IF (found_b) THEN
398 : IF (has_symmetry .AND. row /= col) THEN
399 : trace = trace + 2.0_dp*SUM(block_a*block_b)
400 : ELSE
401 : trace = trace + SUM(block_a*block_b)
402 : END IF
403 : END IF
404 : END DO
405 : CALL dbcsr_iterator_stop(iter)
406 : !$OMP END PARALLEL
407 :
408 6649901 : CALL dbcsr_get_info(matrix_a, group=group)
409 6649901 : CALL group%sum(trace)
410 :
411 6649901 : CALL timestop(handle)
412 6649901 : END SUBROUTINE dbcsr_dot
413 :
414 : ! **************************************************************************************************
415 : !> \brief Computes the trace of the given matrix, also known as the sum of its diagonal elements.
416 : !> \param matrix ...
417 : !> \param trace ...
418 : ! **************************************************************************************************
419 79977 : SUBROUTINE dbcsr_trace(matrix, trace)
420 : TYPE(dbcsr_type), INTENT(IN) :: matrix
421 : REAL(KIND=dp), INTENT(OUT) :: trace
422 :
423 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace'
424 :
425 : INTEGER :: col, col_size, handle, i, row, row_size
426 26659 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
427 : TYPE(dbcsr_iterator_type) :: iter
428 : TYPE(mp_comm_type) :: group
429 :
430 26659 : CALL timeset(routineN, handle)
431 26659 : CPASSERT(omp_get_num_threads() == 1)
432 :
433 26659 : trace = 0.0_dp
434 26659 : CALL dbcsr_iterator_readonly_start(iter, matrix)
435 1648909 : DO WHILE (dbcsr_iterator_blocks_left(iter))
436 1622250 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
437 1648909 : IF (row == col) THEN
438 144064 : CPASSERT(row_size == col_size)
439 617441 : DO i = 1, row_size
440 617441 : trace = trace + block(i, i)
441 : END DO
442 : END IF
443 : END DO
444 26659 : CALL dbcsr_iterator_stop(iter)
445 :
446 26659 : CALL dbcsr_get_info(matrix, group=group)
447 26659 : CALL group%sum(trace)
448 :
449 26659 : CALL timestop(handle)
450 26659 : END SUBROUTINE dbcsr_trace
451 :
452 : ! **************************************************************************************************
453 : !> \brief Copies the diagonal blocks of matrix into diag.
454 : !> \param matrix ...
455 : !> \param diag ...
456 : ! **************************************************************************************************
457 232812 : SUBROUTINE dbcsr_get_block_diag(matrix, diag)
458 : TYPE(dbcsr_type), INTENT(IN) :: matrix
459 : TYPE(dbcsr_type), INTENT(INOUT) :: diag
460 :
461 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_block_diag'
462 :
463 : CHARACTER(len=default_string_length) :: name
464 : INTEGER :: col, handle, row
465 116406 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
466 : TYPE(dbcsr_iterator_type) :: iter
467 :
468 116406 : CALL timeset(routineN, handle)
469 116406 : CPASSERT(omp_get_num_threads() == 1)
470 :
471 116406 : CALL dbcsr_get_info(matrix, name=name)
472 116406 : CALL dbcsr_create(diag, template=matrix, name='diag of '//name)
473 :
474 116406 : CALL dbcsr_iterator_readonly_start(iter, matrix)
475 8649187 : DO WHILE (dbcsr_iterator_blocks_left(iter))
476 8532781 : CALL dbcsr_iterator_next_block(iter, row, col, block)
477 8649187 : IF (row == col) THEN
478 474399 : CALL dbcsr_put_block(diag, row, col, block)
479 : END IF
480 : END DO
481 116406 : CALL dbcsr_iterator_stop(iter)
482 116406 : CALL dbcsr_finalize(diag)
483 :
484 116406 : CALL timestop(handle)
485 116406 : END SUBROUTINE dbcsr_get_block_diag
486 :
487 : ! **************************************************************************************************
488 : !> \brief Scales the rows/columns of given matrix.
489 : !> \param matrix Matrix to be scaled in-place.
490 : !> \param alpha Vector with scaling factors.
491 : !> \param side Side from which to apply the vector. Allowed values are 'right' and 'left'.
492 : ! **************************************************************************************************
493 267997 : SUBROUTINE dbcsr_scale_by_vector(matrix, alpha, side)
494 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
495 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha
496 : CHARACTER(LEN=*), INTENT(IN) :: side
497 :
498 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_by_vector'
499 :
500 : INTEGER :: col_offset, col_size, handle, i, &
501 : nfullcols_total, nfullrows_total, &
502 : row_offset, row_size
503 : LOGICAL :: right
504 267997 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
505 : TYPE(dbcsr_iterator_type) :: iter
506 :
507 267997 : CALL timeset(routineN, handle)
508 267997 : CPASSERT(omp_get_num_threads() == 1)
509 :
510 267997 : IF (side == 'right') THEN
511 : right = .TRUE.
512 0 : ELSE IF (side == 'left') THEN
513 : right = .FALSE.
514 : ELSE
515 0 : CPABORT("Unknown side: "//TRIM(side))
516 : END IF
517 :
518 : ! Check that alpha and matrix have matching sizes.
519 267997 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
520 267997 : IF (right) THEN
521 267997 : CPASSERT(nfullcols_total == SIZE(alpha))
522 : ELSE
523 0 : CPASSERT(nfullrows_total == SIZE(alpha))
524 : END IF
525 :
526 : !$OMP PARALLEL DEFAULT(NONE) SHARED(matrix, alpha, right) &
527 267997 : !$OMP PRIVATE(iter, block, row_size, col_size, row_offset, col_offset, i)
528 : CALL dbcsr_iterator_start(iter, matrix)
529 : DO WHILE (dbcsr_iterator_blocks_left(iter))
530 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
531 : row_offset=row_offset, col_offset=col_offset)
532 : IF (SIZE(block) == 0) CYCLE ! Skip zero-sized blocks.
533 : IF (right) THEN
534 : DO i = 1, col_size
535 : block(:, i) = block(:, i)*alpha(col_offset + i - 1)
536 : END DO
537 : ELSE
538 : DO i = 1, row_size
539 : block(i, :) = block(i, :)*alpha(row_offset + i - 1)
540 : END DO
541 : END IF
542 : END DO
543 : CALL dbcsr_iterator_stop(iter)
544 : !$OMP END PARALLEL
545 :
546 267997 : CALL timestop(handle)
547 267997 : END SUBROUTINE dbcsr_scale_by_vector
548 :
549 : ! **************************************************************************************************
550 : !> \brief Copies the diagonal elements from the given matrix into the given array.
551 : !> \param matrix ...
552 : !> \param diag ...
553 : ! **************************************************************************************************
554 11270 : SUBROUTINE dbcsr_get_diag(matrix, diag)
555 : TYPE(dbcsr_type), INTENT(IN) :: matrix
556 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
557 :
558 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag'
559 :
560 : INTEGER :: col, col_size, handle, i, &
561 : nfullcols_total, nfullrows_total, row, &
562 : row_offset, row_size
563 11270 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
564 : TYPE(dbcsr_iterator_type) :: iter
565 :
566 11270 : CALL timeset(routineN, handle)
567 11270 : CPASSERT(omp_get_num_threads() == 1)
568 :
569 11270 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
570 11270 : CPASSERT(nfullrows_total == nfullcols_total)
571 11270 : CPASSERT(nfullrows_total == SIZE(diag))
572 :
573 1228468 : diag(:) = 0.0_dp
574 11270 : CALL dbcsr_iterator_readonly_start(iter, matrix)
575 279385 : DO WHILE (dbcsr_iterator_blocks_left(iter))
576 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
577 268115 : col_size=col_size, row_offset=row_offset)
578 279385 : IF (row == col) THEN
579 34526 : CPASSERT(row_size == col_size)
580 646099 : DO i = 1, row_size
581 646099 : diag(row_offset + i - 1) = block(i, i)
582 : END DO
583 : END IF
584 : END DO
585 11270 : CALL dbcsr_iterator_stop(iter)
586 :
587 11270 : CALL timestop(handle)
588 11270 : END SUBROUTINE dbcsr_get_diag
589 :
590 : ! **************************************************************************************************
591 : !> \brief Copies the diagonal elements from the given array into the given matrix.
592 : !> \param matrix ...
593 : !> \param diag ...
594 : ! **************************************************************************************************
595 6012 : SUBROUTINE dbcsr_set_diag(matrix, diag)
596 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
597 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag
598 :
599 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag'
600 :
601 : INTEGER :: col, col_size, handle, i, &
602 : nfullcols_total, nfullrows_total, row, &
603 : row_offset, row_size
604 6012 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
605 : TYPE(dbcsr_iterator_type) :: iter
606 :
607 6012 : CALL timeset(routineN, handle)
608 6012 : CPASSERT(omp_get_num_threads() == 1)
609 :
610 6012 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
611 6012 : CPASSERT(nfullrows_total == nfullcols_total)
612 6012 : CPASSERT(nfullrows_total == SIZE(diag))
613 :
614 6012 : CALL dbcsr_reserve_diag_blocks(matrix)
615 :
616 6012 : CALL dbcsr_iterator_start(iter, matrix)
617 153739 : DO WHILE (dbcsr_iterator_blocks_left(iter))
618 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
619 147727 : col_size=col_size, row_offset=row_offset)
620 153739 : IF (row == col) THEN
621 16560 : CPASSERT(row_size == col_size)
622 95680 : DO i = 1, row_size
623 95680 : block(i, i) = diag(row_offset + i - 1)
624 : END DO
625 : END IF
626 : END DO
627 6012 : CALL dbcsr_iterator_stop(iter)
628 :
629 6012 : CALL timestop(handle)
630 6012 : END SUBROUTINE dbcsr_set_diag
631 :
632 : ! **************************************************************************************************
633 : !> \brief Calculates the checksum of a DBCSR matrix.
634 : !> \param matrix ...
635 : !> \param pos Enable position-dependent checksum.
636 : !> \return ...
637 : ! **************************************************************************************************
638 55728 : FUNCTION dbcsr_checksum(matrix, pos) RESULT(checksum)
639 :
640 : TYPE(dbcsr_type), INTENT(IN) :: matrix
641 : LOGICAL, INTENT(IN), OPTIONAL :: pos
642 : REAL(KIND=dp) :: checksum
643 :
644 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_checksum'
645 :
646 : INTEGER :: col_offset, col_size, handle, i, j, &
647 : row_offset, row_size
648 : LOGICAL :: my_pos
649 : REAL(KIND=dp) :: position_factor
650 18576 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
651 : TYPE(dbcsr_iterator_type) :: iter
652 : TYPE(mp_comm_type) :: group
653 :
654 18576 : CALL timeset(routineN, handle)
655 18576 : CPASSERT(omp_get_num_threads() == 1)
656 :
657 18576 : my_pos = .FALSE.
658 18576 : IF (PRESENT(pos)) THEN
659 286 : my_pos = pos
660 : END IF
661 :
662 18576 : checksum = 0.0_dp
663 18576 : CALL dbcsr_iterator_readonly_start(iter, matrix)
664 115098 : DO WHILE (dbcsr_iterator_blocks_left(iter))
665 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
666 96522 : row_offset=row_offset, col_offset=col_offset)
667 115098 : IF (my_pos) THEN
668 21372 : DO i = 1, row_size
669 218025 : DO j = 1, col_size
670 196653 : position_factor = LOG(REAL((row_offset + i - 1)*(col_offset + j - 1), KIND=dp))
671 215930 : checksum = checksum + block(i, j)*position_factor
672 : END DO
673 : END DO
674 : ELSE
675 4422423 : checksum = checksum + SUM(block**2)
676 : END IF
677 : END DO
678 18576 : CALL dbcsr_iterator_stop(iter)
679 :
680 18576 : CALL dbcsr_get_info(matrix, group=group)
681 18576 : CALL group%sum(checksum)
682 :
683 18576 : CALL timestop(handle)
684 18576 : END FUNCTION dbcsr_checksum
685 :
686 : ! **************************************************************************************************
687 : !> \brief Prints given matrix in matlab format (only present blocks).
688 : !> \param matrix ...
689 : !> \param variable_name ...
690 : !> \param unit_nr ...
691 : ! **************************************************************************************************
692 0 : SUBROUTINE dbcsr_print(matrix, variable_name, unit_nr)
693 : TYPE(dbcsr_type), INTENT(IN) :: matrix
694 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: variable_name
695 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
696 :
697 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_print'
698 :
699 : CHARACTER(len=default_string_length) :: my_variable_name, name
700 : INTEGER :: col_offset, col_size, handle, i, iw, j, nblkcols_total, nblkrows_total, &
701 : nfullcols_total, nfullrows_total, row_offset, row_size
702 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
703 : TYPE(dbcsr_iterator_type) :: iter
704 :
705 0 : CALL timeset(routineN, handle)
706 0 : CPASSERT(omp_get_num_threads() == 1)
707 :
708 0 : iw = default_output_unit
709 0 : IF (PRESENT(unit_nr)) iw = unit_nr
710 :
711 0 : my_variable_name = 'a'
712 0 : IF (PRESENT(variable_name)) my_variable_name = variable_name
713 :
714 : ! Print matrix properties.
715 : CALL dbcsr_get_info(matrix, name=name, &
716 : nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total, &
717 0 : nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
718 0 : WRITE (iw, *) "===", routineN, "==="
719 0 : WRITE (iw, *) "Name:", name
720 0 : WRITE (iw, *) "Symmetry:", dbcsr_has_symmetry(matrix)
721 0 : WRITE (iw, *) "Number of blocks:", dbcsr_get_num_blocks(matrix)
722 0 : WRITE (iw, *) "Data size:", dbcsr_get_data_size(matrix)
723 0 : WRITE (iw, *) "Occupation:", dbcsr_get_occupation(matrix)
724 0 : WRITE (iw, *) "Full size:", nfullrows_total, "x", nfullcols_total
725 0 : WRITE (iw, *) "Blocked size:", nblkrows_total, "x", nblkcols_total
726 :
727 : ! Print matrix blocks.
728 0 : CALL dbcsr_iterator_readonly_start(iter, matrix)
729 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
730 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
731 0 : row_offset=row_offset, col_offset=col_offset)
732 0 : DO i = 1, row_size
733 0 : DO j = 1, col_size
734 0 : WRITE (iw, '(A,I4,A,I4,A,E23.16,A)') TRIM(my_variable_name)//'(', &
735 0 : row_offset + i - 1, ',', col_offset + j - 1, ')=', block(i, j), ';'
736 : END DO
737 : END DO
738 : END DO
739 0 : CALL dbcsr_iterator_stop(iter)
740 :
741 0 : CALL timestop(handle)
742 0 : END SUBROUTINE dbcsr_print
743 :
744 : END MODULE cp_dbcsr_contrib
|
