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 : 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 286000 : 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 57200 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
62 : TYPE(dbcsr_iterator_type) :: iter
63 :
64 57200 : CALL timeset(routineN, handle)
65 57200 : CPASSERT(omp_get_num_threads() == 1)
66 :
67 57200 : CALL dbcsr_get_info(matrix_a, nblkrows_total=nblkrows_tot_a)
68 57200 : CALL dbcsr_get_info(matrix_b, nblkrows_total=nblkrows_tot_b)
69 57200 : CALL dbcsr_get_info(matrix_c, nblkrows_total=nblkrows_tot_c)
70 57200 : 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 57200 : CALL dbcsr_clear(matrix_c)
75 57200 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
76 145199 : DO WHILE (dbcsr_iterator_blocks_left(iter))
77 87999 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
78 87999 : CALL dbcsr_get_readonly_block_p(matrix_b, row, col, block_b, found_b)
79 145199 : IF (found_b) THEN
80 18003585 : CALL dbcsr_put_block(matrix_c, row, col, block_a*block_b)
81 : END IF
82 : END DO
83 57200 : CALL dbcsr_iterator_stop(iter)
84 57200 : CALL dbcsr_finalize(matrix_c)
85 57200 : CALL timestop(handle)
86 57200 : END SUBROUTINE dbcsr_hadamard_product
87 :
88 : ! **************************************************************************************************
89 : !> \brief Compute the maxabs norm of a dbcsr matrix
90 : !> \param matrix ...
91 : !> \return ...
92 : ! **************************************************************************************************
93 38502 : 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 12834 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
101 : TYPE(dbcsr_iterator_type) :: iter
102 : TYPE(mp_comm_type) :: group
103 :
104 12834 : CALL timeset(routineN, handle)
105 12834 : CPASSERT(omp_get_num_threads() == 1)
106 :
107 12834 : norm = 0.0_dp
108 12834 : CALL dbcsr_iterator_readonly_start(iter, matrix)
109 172372 : DO WHILE (dbcsr_iterator_blocks_left(iter))
110 159538 : CALL dbcsr_iterator_next_block(iter, block=block)
111 4059289 : norm = MAX(norm, MAXVAL(ABS(block)))
112 : END DO
113 12834 : CALL dbcsr_iterator_stop(iter)
114 :
115 12834 : CALL dbcsr_get_info(matrix, group=group)
116 12834 : CALL group%max(norm)
117 :
118 12834 : CALL timestop(handle)
119 12834 : END FUNCTION dbcsr_maxabs
120 :
121 : ! **************************************************************************************************
122 : !> \brief Compute the frobenius norm of a dbcsr matrix
123 : !> \param matrix ...
124 : !> \return ...
125 : ! **************************************************************************************************
126 1113738 : FUNCTION dbcsr_frobenius_norm(matrix) RESULT(norm)
127 : TYPE(dbcsr_type), INTENT(IN) :: matrix
128 : REAL(dp) :: norm
129 :
130 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_frobenius_norm'
131 :
132 : INTEGER :: col, handle, row
133 : LOGICAL :: has_symmetry
134 371246 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
135 : TYPE(dbcsr_iterator_type) :: iter
136 : TYPE(mp_comm_type) :: group
137 :
138 371246 : CALL timeset(routineN, handle)
139 371246 : CPASSERT(omp_get_num_threads() == 1)
140 :
141 371246 : has_symmetry = dbcsr_has_symmetry(matrix)
142 371246 : norm = 0.0_dp
143 371246 : CALL dbcsr_iterator_readonly_start(iter, matrix)
144 6371885 : DO WHILE (dbcsr_iterator_blocks_left(iter))
145 6000639 : CALL dbcsr_iterator_next_block(iter, row, col, block)
146 6371885 : IF (has_symmetry .AND. row /= col) THEN
147 5588714 : norm = norm + 2.0_dp*SUM(block**2)
148 : ELSE
149 73606885 : norm = norm + SUM(block**2)
150 : END IF
151 : END DO
152 371246 : CALL dbcsr_iterator_stop(iter)
153 :
154 371246 : CALL dbcsr_get_info(matrix, group=group)
155 371246 : CALL group%sum(norm)
156 371246 : norm = SQRT(norm)
157 :
158 371246 : CALL timestop(handle)
159 371246 : END FUNCTION dbcsr_frobenius_norm
160 :
161 : ! **************************************************************************************************
162 : !> \brief Compute the gershgorin norm of a dbcsr matrix
163 : !> \param matrix ...
164 : !> \return ...
165 : ! **************************************************************************************************
166 7642 : FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm)
167 : TYPE(dbcsr_type), INTENT(IN) :: matrix
168 : REAL(dp) :: norm
169 :
170 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_gershgorin_norm'
171 :
172 : INTEGER :: col, col_offset, handle, i, j, ncol, &
173 : nrow, row, row_offset
174 : LOGICAL :: has_symmetry
175 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buffer
176 7642 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
177 : TYPE(dbcsr_iterator_type) :: iter
178 : TYPE(mp_comm_type) :: group
179 :
180 7642 : CALL timeset(routineN, handle)
181 7642 : CPASSERT(omp_get_num_threads() == 1)
182 :
183 7642 : has_symmetry = dbcsr_has_symmetry(matrix)
184 7642 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
185 7642 : CPASSERT(nrow == ncol)
186 22918 : ALLOCATE (buffer(nrow))
187 136856 : buffer = 0.0_dp
188 :
189 7642 : CALL dbcsr_iterator_readonly_start(iter, matrix)
190 194283 : DO WHILE (dbcsr_iterator_blocks_left(iter))
191 186641 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_offset=row_offset, col_offset=col_offset)
192 609890 : DO j = 1, SIZE(block, 2)
193 4481467 : DO i = 1, SIZE(block, 1)
194 3879219 : buffer(row_offset + i - 1) = buffer(row_offset + i - 1) + ABS(block(i, j))
195 4294826 : IF (has_symmetry .AND. row /= col) THEN
196 40164 : buffer(col_offset + j - 1) = buffer(col_offset + j - 1) + ABS(block(i, j))
197 : END IF
198 : END DO
199 : END DO
200 : END DO
201 7642 : CALL dbcsr_iterator_stop(iter)
202 :
203 7642 : CALL dbcsr_get_info(matrix, group=group)
204 7642 : CALL group%sum(buffer)
205 144498 : norm = MAXVAL(buffer)
206 7642 : DEALLOCATE (buffer)
207 :
208 7642 : CALL timestop(handle)
209 30568 : END FUNCTION dbcsr_gershgorin_norm
210 :
211 : ! **************************************************************************************************
212 : !> \brief Fills the given matrix with random numbers.
213 : !> \param matrix ...
214 : !> \param keep_sparsity ...
215 : ! **************************************************************************************************
216 708 : SUBROUTINE dbcsr_init_random(matrix, keep_sparsity)
217 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
218 : LOGICAL, OPTIONAL :: keep_sparsity
219 :
220 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_init_random'
221 :
222 : INTEGER :: col, col_size, handle, ncol, nrow, row, &
223 : row_size
224 : INTEGER, DIMENSION(4) :: iseed
225 : LOGICAL :: my_keep_sparsity
226 236 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
227 : TYPE(dbcsr_iterator_type) :: iter
228 :
229 236 : CALL timeset(routineN, handle)
230 236 : CPASSERT(omp_get_num_threads() == 1)
231 :
232 236 : my_keep_sparsity = .FALSE.
233 236 : IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
234 236 : IF (.NOT. my_keep_sparsity) CALL dbcsr_reserve_all_blocks(matrix)
235 236 : CALL dbcsr_get_info(matrix, nblkrows_total=nrow, nblkcols_total=ncol)
236 :
237 236 : CALL dbcsr_iterator_start(iter, matrix)
238 695 : DO WHILE (dbcsr_iterator_blocks_left(iter))
239 459 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
240 : ! set the seed for dlarnv, is here to guarantee same value of the random numbers
241 : ! for all layouts (and block distributions)
242 459 : iseed = generate_larnv_seed(row, nrow, col, ncol, 1)
243 695 : CALL dlarnv(1, iseed, row_size*col_size, block(1, 1))
244 : END DO
245 236 : CALL dbcsr_iterator_stop(iter)
246 :
247 236 : CALL timestop(handle)
248 236 : END SUBROUTINE dbcsr_init_random
249 :
250 : ! **************************************************************************************************
251 : !> \brief Reserves all diagonal blocks.
252 : !> \param matrix ...
253 : ! **************************************************************************************************
254 481506 : SUBROUTINE dbcsr_reserve_diag_blocks(matrix)
255 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
256 :
257 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reserve_diag_blocks'
258 :
259 : INTEGER :: handle, i, k, mynode, nblkcols_total, &
260 : nblkrows_total, owner
261 481506 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_diag
262 481506 : INTEGER, DIMENSION(:), POINTER :: local_rows
263 : TYPE(dbcsr_distribution_type) :: dist
264 :
265 481506 : CALL timeset(routineN, handle)
266 481506 : CPASSERT(omp_get_num_threads() == 1)
267 :
268 481506 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total)
269 481506 : CPASSERT(nblkrows_total == nblkcols_total)
270 :
271 481506 : CALL dbcsr_get_info(matrix, local_rows=local_rows, distribution=dist)
272 481506 : CALL dbcsr_distribution_get(dist, mynode=mynode)
273 1377073 : ALLOCATE (local_diag(SIZE(local_rows)))
274 :
275 1409832 : k = 0
276 1409832 : DO i = 1, SIZE(local_rows)
277 928326 : CALL dbcsr_get_stored_coordinates(matrix, row=local_rows(i), column=local_rows(i), processor=owner)
278 1409832 : IF (owner == mynode) THEN
279 928326 : k = k + 1
280 928326 : local_diag(k) = local_rows(i)
281 : END IF
282 : END DO
283 :
284 481506 : CALL dbcsr_reserve_blocks(matrix, rows=local_diag(1:k), cols=local_diag(1:k))
285 481506 : DEALLOCATE (local_diag)
286 :
287 481506 : CALL timestop(handle)
288 1444518 : END SUBROUTINE dbcsr_reserve_diag_blocks
289 :
290 : ! **************************************************************************************************
291 : !> \brief Reserves all blocks.
292 : !> \param matrix ...
293 : ! **************************************************************************************************
294 1599001 : SUBROUTINE dbcsr_reserve_all_blocks(matrix)
295 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
296 :
297 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reserve_all_blocks'
298 :
299 : INTEGER :: handle, i, j, k, n
300 1599001 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
301 1599001 : INTEGER, DIMENSION(:), POINTER :: local_cols, local_rows
302 :
303 1599001 : CALL timeset(routineN, handle)
304 1599001 : CPASSERT(omp_get_num_threads() == 1)
305 :
306 1599001 : CALL dbcsr_get_info(matrix, local_rows=local_rows, local_cols=local_cols)
307 1599001 : n = SIZE(local_rows)*SIZE(local_cols)
308 6119916 : ALLOCATE (rows(n), cols(n))
309 :
310 3703961 : k = 0
311 3703961 : DO i = 1, SIZE(local_rows)
312 7320085 : DO j = 1, SIZE(local_cols)
313 3616124 : k = k + 1
314 3616124 : rows(k) = local_rows(i)
315 5721084 : cols(k) = local_cols(j)
316 : END DO
317 : END DO
318 :
319 1599001 : CALL dbcsr_reserve_blocks(matrix, rows=rows(1:k), cols=cols(1:k))
320 1599001 : DEALLOCATE (rows, cols)
321 :
322 1599001 : CALL timestop(handle)
323 1599001 : END SUBROUTINE dbcsr_reserve_all_blocks
324 :
325 : ! **************************************************************************************************
326 : !> \brief Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
327 : !> \param matrix ...
328 : !> \param alpha ...
329 : ! **************************************************************************************************
330 862892 : SUBROUTINE dbcsr_add_on_diag(matrix, alpha)
331 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
332 : REAL(kind=dp), INTENT(IN) :: alpha
333 :
334 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag'
335 :
336 : INTEGER :: col, col_size, handle, i, row, row_size
337 431446 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
338 : TYPE(dbcsr_iterator_type) :: iter
339 :
340 431446 : CALL timeset(routineN, handle)
341 431446 : CPASSERT(omp_get_num_threads() == 1)
342 :
343 431446 : CALL dbcsr_reserve_diag_blocks(matrix)
344 :
345 431446 : CALL dbcsr_iterator_start(iter, matrix)
346 6412214 : DO WHILE (dbcsr_iterator_blocks_left(iter))
347 5980768 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
348 6412214 : IF (row == col) THEN
349 841838 : CPASSERT(row_size == col_size)
350 4128318 : DO i = 1, row_size
351 4128318 : block(i, i) = block(i, i) + alpha
352 : END DO
353 : END IF
354 : END DO
355 431446 : CALL dbcsr_iterator_stop(iter)
356 :
357 431446 : CALL timestop(handle)
358 431446 : END SUBROUTINE dbcsr_add_on_diag
359 :
360 : ! **************************************************************************************************
361 : !> \brief Computes the dot product of two matrices, also known as the trace of their matrix product.
362 : !> \param matrix_a ...
363 : !> \param matrix_b ...
364 : !> \param trace ...
365 : ! **************************************************************************************************
366 5713107 : SUBROUTINE dbcsr_dot(matrix_a, matrix_b, trace)
367 : TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b
368 : REAL(KIND=dp), INTENT(OUT) :: trace
369 :
370 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_dot'
371 :
372 : INTEGER :: col, handle, row
373 : LOGICAL :: found_b, has_symmetry
374 1904369 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_a, block_b
375 : TYPE(dbcsr_iterator_type) :: iter
376 : TYPE(mp_comm_type) :: group
377 :
378 1904369 : CALL timeset(routineN, handle)
379 1904369 : CPASSERT(omp_get_num_threads() == 1)
380 :
381 1904369 : CPASSERT(dbcsr_has_symmetry(matrix_a) .EQV. dbcsr_has_symmetry(matrix_b))
382 1904369 : has_symmetry = dbcsr_has_symmetry(matrix_a)
383 :
384 1904369 : trace = 0.0_dp
385 1904369 : CALL dbcsr_iterator_readonly_start(iter, matrix_a)
386 35190884 : DO WHILE (dbcsr_iterator_blocks_left(iter))
387 33286515 : CALL dbcsr_iterator_next_block(iter, row, col, block_a)
388 99859545 : IF (SIZE(block_a) == 0) CYCLE ! Skip zero-sized blocks.
389 33285559 : CALL dbcsr_get_readonly_block_p(matrix_b, row, col, block_b, found_b)
390 35189928 : IF (found_b) THEN
391 32911852 : IF (has_symmetry .AND. row /= col) THEN
392 957524007 : trace = trace + 2.0_dp*SUM(block_a*block_b)
393 : ELSE
394 631094713 : trace = trace + SUM(block_a*block_b)
395 : END IF
396 : END IF
397 : END DO
398 1904369 : CALL dbcsr_iterator_stop(iter)
399 :
400 1904369 : CALL dbcsr_get_info(matrix_a, group=group)
401 1904369 : CALL group%sum(trace)
402 :
403 1904369 : CALL timestop(handle)
404 1904369 : END SUBROUTINE dbcsr_dot
405 :
406 : ! **************************************************************************************************
407 : !> \brief Computes the trace of the given matrix, also known as the sum of its diagonal elements.
408 : !> \param matrix ...
409 : !> \param trace ...
410 : ! **************************************************************************************************
411 48033 : SUBROUTINE dbcsr_trace(matrix, trace)
412 : TYPE(dbcsr_type), INTENT(IN) :: matrix
413 : REAL(KIND=dp), INTENT(OUT) :: trace
414 :
415 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace'
416 :
417 : INTEGER :: col, col_size, handle, i, row, row_size
418 16011 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
419 : TYPE(dbcsr_iterator_type) :: iter
420 : TYPE(mp_comm_type) :: group
421 :
422 16011 : CALL timeset(routineN, handle)
423 16011 : CPASSERT(omp_get_num_threads() == 1)
424 :
425 16011 : trace = 0.0_dp
426 16011 : CALL dbcsr_iterator_readonly_start(iter, matrix)
427 1431818 : DO WHILE (dbcsr_iterator_blocks_left(iter))
428 1415807 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, col_size=col_size)
429 1431818 : IF (row == col) THEN
430 122371 : CPASSERT(row_size == col_size)
431 509140 : DO i = 1, row_size
432 509140 : trace = trace + block(i, i)
433 : END DO
434 : END IF
435 : END DO
436 16011 : CALL dbcsr_iterator_stop(iter)
437 :
438 16011 : CALL dbcsr_get_info(matrix, group=group)
439 16011 : CALL group%sum(trace)
440 :
441 16011 : CALL timestop(handle)
442 16011 : END SUBROUTINE dbcsr_trace
443 :
444 : ! **************************************************************************************************
445 : !> \brief Copies the diagonal blocks of matrix into diag.
446 : !> \param matrix ...
447 : !> \param diag ...
448 : ! **************************************************************************************************
449 199616 : SUBROUTINE dbcsr_get_block_diag(matrix, diag)
450 : TYPE(dbcsr_type), INTENT(IN) :: matrix
451 : TYPE(dbcsr_type), INTENT(INOUT) :: diag
452 :
453 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_block_diag'
454 :
455 : CHARACTER(len=default_string_length) :: name
456 : INTEGER :: col, handle, row
457 99808 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
458 : TYPE(dbcsr_iterator_type) :: iter
459 :
460 99808 : CALL timeset(routineN, handle)
461 99808 : CPASSERT(omp_get_num_threads() == 1)
462 :
463 99808 : CALL dbcsr_get_info(matrix, name=name)
464 99808 : CALL dbcsr_create(diag, template=matrix, name='diag of '//name)
465 :
466 99808 : CALL dbcsr_iterator_readonly_start(iter, matrix)
467 8478720 : DO WHILE (dbcsr_iterator_blocks_left(iter))
468 8378912 : CALL dbcsr_iterator_next_block(iter, row, col, block)
469 8478720 : IF (row == col) THEN
470 445001 : CALL dbcsr_put_block(diag, row, col, block)
471 : END IF
472 : END DO
473 99808 : CALL dbcsr_iterator_stop(iter)
474 99808 : CALL dbcsr_finalize(diag)
475 :
476 99808 : CALL timestop(handle)
477 99808 : END SUBROUTINE dbcsr_get_block_diag
478 :
479 : ! **************************************************************************************************
480 : !> \brief Scales the rows/columns of given matrix.
481 : !> \param matrix Matrix to be scaled in-place.
482 : !> \param alpha Vector with scaling factors.
483 : !> \param side Side from which to apply the vector. Allowed values are 'right' and 'left'.
484 : ! **************************************************************************************************
485 166260 : SUBROUTINE dbcsr_scale_by_vector(matrix, alpha, side)
486 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
487 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha
488 : CHARACTER(LEN=*), INTENT(IN) :: side
489 :
490 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_by_vector'
491 :
492 : INTEGER :: col_offset, col_size, handle, i, &
493 : nfullcols_total, nfullrows_total, &
494 : row_offset, row_size
495 : LOGICAL :: right
496 166260 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
497 : TYPE(dbcsr_iterator_type) :: iter
498 :
499 166260 : CALL timeset(routineN, handle)
500 166260 : CPASSERT(omp_get_num_threads() == 1)
501 :
502 166260 : IF (side == 'right') THEN
503 : right = .TRUE.
504 0 : ELSE IF (side == 'left') THEN
505 : right = .FALSE.
506 : ELSE
507 0 : CPABORT("Unknown side: "//TRIM(side))
508 : END IF
509 :
510 : ! Check that alpha and matrix have matching sizes.
511 166260 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
512 166260 : IF (right) THEN
513 166260 : CPASSERT(nfullcols_total == SIZE(alpha))
514 : ELSE
515 0 : CPASSERT(nfullrows_total == SIZE(alpha))
516 : END IF
517 :
518 166260 : CALL dbcsr_iterator_start(iter, matrix)
519 739998 : DO WHILE (dbcsr_iterator_blocks_left(iter))
520 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
521 573738 : row_offset=row_offset, col_offset=col_offset)
522 1721214 : IF (SIZE(block) == 0) CYCLE ! Skip zero-sized blocks.
523 739998 : IF (right) THEN
524 12748749 : DO i = 1, col_size
525 68650030 : block(:, i) = block(:, i)*alpha(col_offset + i - 1)
526 : END DO
527 : ELSE
528 0 : DO i = 1, row_size
529 0 : block(i, :) = block(i, :)*alpha(row_offset + i - 1)
530 : END DO
531 : END IF
532 : END DO
533 166260 : CALL dbcsr_iterator_stop(iter)
534 :
535 166260 : CALL timestop(handle)
536 166260 : END SUBROUTINE dbcsr_scale_by_vector
537 :
538 : ! **************************************************************************************************
539 : !> \brief Copies the diagonal elements from the given matrix into the given array.
540 : !> \param matrix ...
541 : !> \param diag ...
542 : ! **************************************************************************************************
543 2696 : SUBROUTINE dbcsr_get_diag(matrix, diag)
544 : TYPE(dbcsr_type), INTENT(IN) :: matrix
545 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
546 :
547 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag'
548 :
549 : INTEGER :: col, col_size, handle, i, &
550 : nfullcols_total, nfullrows_total, row, &
551 : row_offset, row_size
552 2696 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
553 : TYPE(dbcsr_iterator_type) :: iter
554 :
555 2696 : CALL timeset(routineN, handle)
556 2696 : CPASSERT(omp_get_num_threads() == 1)
557 :
558 2696 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
559 2696 : CPASSERT(nfullrows_total == nfullcols_total)
560 2696 : CPASSERT(nfullrows_total == SIZE(diag))
561 :
562 75640 : diag(:) = 0.0_dp
563 2696 : CALL dbcsr_iterator_readonly_start(iter, matrix)
564 67153 : DO WHILE (dbcsr_iterator_blocks_left(iter))
565 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
566 64457 : col_size=col_size, row_offset=row_offset)
567 67153 : IF (row == col) THEN
568 8676 : CPASSERT(row_size == col_size)
569 45188 : DO i = 1, row_size
570 45188 : diag(row_offset + i - 1) = block(i, i)
571 : END DO
572 : END IF
573 : END DO
574 2696 : CALL dbcsr_iterator_stop(iter)
575 :
576 2696 : CALL timestop(handle)
577 2696 : END SUBROUTINE dbcsr_get_diag
578 :
579 : ! **************************************************************************************************
580 : !> \brief Copies the diagonal elements from the given array into the given matrix.
581 : !> \param matrix ...
582 : !> \param diag ...
583 : ! **************************************************************************************************
584 2954 : SUBROUTINE dbcsr_set_diag(matrix, diag)
585 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
586 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: diag
587 :
588 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag'
589 :
590 : INTEGER :: col, col_size, handle, i, &
591 : nfullcols_total, nfullrows_total, row, &
592 : row_offset, row_size
593 2954 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
594 : TYPE(dbcsr_iterator_type) :: iter
595 :
596 2954 : CALL timeset(routineN, handle)
597 2954 : CPASSERT(omp_get_num_threads() == 1)
598 :
599 2954 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
600 2954 : CPASSERT(nfullrows_total == nfullcols_total)
601 2954 : CPASSERT(nfullrows_total == SIZE(diag))
602 :
603 2954 : CALL dbcsr_reserve_diag_blocks(matrix)
604 :
605 2954 : CALL dbcsr_iterator_start(iter, matrix)
606 137056 : DO WHILE (dbcsr_iterator_blocks_left(iter))
607 : CALL dbcsr_iterator_next_block(iter, row, col, block, row_size=row_size, &
608 134102 : col_size=col_size, row_offset=row_offset)
609 137056 : IF (row == col) THEN
610 10966 : CPASSERT(row_size == col_size)
611 51978 : DO i = 1, row_size
612 51978 : block(i, i) = diag(row_offset + i - 1)
613 : END DO
614 : END IF
615 : END DO
616 2954 : CALL dbcsr_iterator_stop(iter)
617 :
618 2954 : CALL timestop(handle)
619 2954 : END SUBROUTINE dbcsr_set_diag
620 :
621 : ! **************************************************************************************************
622 : !> \brief Calculates the checksum of a DBCSR matrix.
623 : !> \param matrix ...
624 : !> \param pos Enable position-dependent checksum.
625 : !> \return ...
626 : ! **************************************************************************************************
627 47322 : FUNCTION dbcsr_checksum(matrix, pos) RESULT(checksum)
628 :
629 : TYPE(dbcsr_type), INTENT(IN) :: matrix
630 : LOGICAL, INTENT(IN), OPTIONAL :: pos
631 : REAL(KIND=dp) :: checksum
632 :
633 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_checksum'
634 :
635 : INTEGER :: col_offset, col_size, handle, i, j, &
636 : row_offset, row_size
637 : LOGICAL :: my_pos
638 : REAL(KIND=dp) :: position_factor
639 15774 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
640 : TYPE(dbcsr_iterator_type) :: iter
641 : TYPE(mp_comm_type) :: group
642 :
643 15774 : CALL timeset(routineN, handle)
644 15774 : CPASSERT(omp_get_num_threads() == 1)
645 :
646 15774 : my_pos = .FALSE.
647 15774 : IF (PRESENT(pos)) THEN
648 314 : my_pos = pos
649 : END IF
650 :
651 15774 : checksum = 0.0_dp
652 15774 : CALL dbcsr_iterator_readonly_start(iter, matrix)
653 102956 : DO WHILE (dbcsr_iterator_blocks_left(iter))
654 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
655 87182 : row_offset=row_offset, col_offset=col_offset)
656 102956 : IF (my_pos) THEN
657 21996 : DO i = 1, row_size
658 222881 : DO j = 1, col_size
659 200885 : position_factor = LOG(REAL((row_offset + i - 1)*(col_offset + j - 1), KIND=dp))
660 220714 : checksum = checksum + block(i, j)*position_factor
661 : END DO
662 : END DO
663 : ELSE
664 3909479 : checksum = checksum + SUM(block**2)
665 : END IF
666 : END DO
667 15774 : CALL dbcsr_iterator_stop(iter)
668 :
669 15774 : CALL dbcsr_get_info(matrix, group=group)
670 15774 : CALL group%sum(checksum)
671 :
672 15774 : CALL timestop(handle)
673 15774 : END FUNCTION dbcsr_checksum
674 :
675 : ! **************************************************************************************************
676 : !> \brief Prints given matrix in matlab format (only present blocks).
677 : !> \param matrix ...
678 : !> \param variable_name ...
679 : !> \param unit_nr ...
680 : ! **************************************************************************************************
681 0 : SUBROUTINE dbcsr_print(matrix, variable_name, unit_nr)
682 : TYPE(dbcsr_type), INTENT(IN) :: matrix
683 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: variable_name
684 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
685 :
686 : CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_print'
687 :
688 : CHARACTER(len=default_string_length) :: my_variable_name, name
689 : INTEGER :: col_offset, col_size, handle, i, iw, j, nblkcols_total, nblkrows_total, &
690 : nfullcols_total, nfullrows_total, row_offset, row_size
691 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
692 : TYPE(dbcsr_iterator_type) :: iter
693 :
694 0 : CALL timeset(routineN, handle)
695 0 : CPASSERT(omp_get_num_threads() == 1)
696 :
697 0 : iw = default_output_unit
698 0 : IF (PRESENT(unit_nr)) iw = unit_nr
699 :
700 0 : my_variable_name = 'a'
701 0 : IF (PRESENT(variable_name)) my_variable_name = variable_name
702 :
703 : ! Print matrix properties.
704 : CALL dbcsr_get_info(matrix, name=name, &
705 : nblkrows_total=nblkrows_total, nblkcols_total=nblkcols_total, &
706 0 : nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
707 0 : WRITE (iw, *) "===", routineN, "==="
708 0 : WRITE (iw, *) "Name:", name
709 0 : WRITE (iw, *) "Symmetry:", dbcsr_has_symmetry(matrix)
710 0 : WRITE (iw, *) "Number of blocks:", dbcsr_get_num_blocks(matrix)
711 0 : WRITE (iw, *) "Data size:", dbcsr_get_data_size(matrix)
712 0 : WRITE (iw, *) "Occupation:", dbcsr_get_occupation(matrix)
713 0 : WRITE (iw, *) "Full size:", nfullrows_total, "x", nfullcols_total
714 0 : WRITE (iw, *) "Blocked size:", nblkrows_total, "x", nblkcols_total
715 :
716 : ! Print matrix blocks.
717 0 : CALL dbcsr_iterator_readonly_start(iter, matrix)
718 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
719 : CALL dbcsr_iterator_next_block(iter, block=block, row_size=row_size, col_size=col_size, &
720 0 : row_offset=row_offset, col_offset=col_offset)
721 0 : DO i = 1, row_size
722 0 : DO j = 1, col_size
723 0 : WRITE (iw, '(A,I4,A,I4,A,E23.16,A)') TRIM(my_variable_name)//'(', &
724 0 : row_offset + i - 1, ',', col_offset + j - 1, ')=', block(i, j), ';'
725 : END DO
726 : END DO
727 : END DO
728 0 : CALL dbcsr_iterator_stop(iter)
729 :
730 0 : CALL timestop(handle)
731 0 : END SUBROUTINE dbcsr_print
732 :
733 : END MODULE cp_dbcsr_contrib
|
