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 : ! **************************************************************************************************
9 : !> \brief Common DBCSR matrix operations used by GW modules.
10 : !> \par History
11 : !> 09.2026 created Jan Wilhelm
12 : ! **************************************************************************************************
13 : MODULE gw_utils_dbcsr
14 : USE cp_dbcsr_api, ONLY: &
15 : dbcsr_copy, dbcsr_create, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
16 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
17 : dbcsr_multiply, dbcsr_release, dbcsr_type
18 : USE kinds, ONLY: dp
19 : #include "./base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 :
24 : PUBLIC :: dbcsr_contract_ABA, hadamard_product, hadamard_product_inplace
25 :
26 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils_dbcsr'
27 :
28 : CONTAINS
29 :
30 : ! **************************************************************************************************
31 : !> \brief Computes the scaled element-wise product C = factor (A ◦ B) while preserving the
32 : !> block structure of A. Blocks absent from B are retained in C with zero values.
33 : !> \param matrix_A First factor and source of the block structure
34 : !> \param matrix_B Second factor
35 : !> \param matrix_C Scaled element-wise product
36 : !> \param factor Scaling factor
37 : ! **************************************************************************************************
38 4796 : SUBROUTINE hadamard_product(matrix_A, matrix_B, matrix_C, factor)
39 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_C
40 : REAL(KIND=dp), INTENT(IN) :: factor
41 :
42 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product'
43 :
44 : INTEGER :: handle, icol, irow
45 : LOGICAL :: found
46 2398 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_B, block_C
47 : TYPE(dbcsr_iterator_type) :: iterator
48 :
49 2398 : CALL timeset(routineN, handle)
50 :
51 2398 : CALL dbcsr_copy(matrix_C, matrix_A)
52 2398 : CALL dbcsr_iterator_start(iterator, matrix_C)
53 301060 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
54 298662 : CALL dbcsr_iterator_next_block(iterator, irow, icol, block_C)
55 298662 : CALL dbcsr_get_block_p(matrix_B, irow, icol, block_B, found)
56 301060 : IF (found) THEN
57 530996716 : block_C(:, :) = factor*block_C(:, :)*block_B(:, :)
58 : ELSE
59 0 : block_C(:, :) = 0.0_dp
60 : END IF
61 : END DO
62 2398 : CALL dbcsr_iterator_stop(iterator)
63 :
64 2398 : CALL timestop(handle)
65 :
66 2398 : END SUBROUTINE hadamard_product
67 :
68 : ! **************************************************************************************************
69 : !> \brief Form A = factor * (A element-wise B) without changing A's block structure.
70 : !> \param matrix_A First factor, overwritten by the product.
71 : !> \param matrix_B Second factor; a missing block represents zero.
72 : !> \param factor Product scale factor.
73 : ! **************************************************************************************************
74 9200 : SUBROUTINE hadamard_product_inplace(matrix_A, matrix_B, factor)
75 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B
76 : REAL(KIND=dp), INTENT(IN) :: factor
77 :
78 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product_inplace'
79 :
80 : INTEGER :: handle, icol, irow
81 : LOGICAL :: found
82 4600 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_A, block_B
83 : TYPE(dbcsr_iterator_type) :: iterator
84 :
85 4600 : CALL timeset(routineN, handle)
86 :
87 4600 : CALL dbcsr_iterator_start(iterator, matrix_A)
88 246348 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
89 241748 : CALL dbcsr_iterator_next_block(iterator, irow, icol, block_A)
90 241748 : CALL dbcsr_get_block_p(matrix_B, irow, icol, block_B, found)
91 246348 : IF (found) THEN
92 385435240 : block_A(:, :) = factor*block_A(:, :)*block_B(:, :)
93 : ELSE
94 0 : block_A(:, :) = 0.0_dp
95 : END IF
96 : END DO
97 4600 : CALL dbcsr_iterator_stop(iterator)
98 :
99 4600 : CALL timestop(handle)
100 4600 : END SUBROUTINE hadamard_product_inplace
101 :
102 : ! **************************************************************************************************
103 : !> \brief Computes C=A B A^T or C=A^T B A for DBCSR matrices.
104 : !> \param trans_A_left transposition applied to the left occurrence of A
105 : !> \param trans_A_right transposition applied to the right occurrence of A
106 : !> \param matrix_A left and right matrix A
107 : !> \param matrix_B input matrix B
108 : !> \param matrix_C output matrix C
109 : !> \param eps_filter filtering threshold for both matrix multiplications
110 : !> \param retain_sparsity if true, only existing blocks of C are filled
111 : ! **************************************************************************************************
112 9800 : SUBROUTINE dbcsr_contract_ABA(trans_A_left, trans_A_right, matrix_A, matrix_B, matrix_C, &
113 : eps_filter, retain_sparsity)
114 : CHARACTER(LEN=1), INTENT(IN) :: trans_A_left, trans_A_right
115 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_C
116 : REAL(KIND=dp), INTENT(IN) :: eps_filter
117 : LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity
118 :
119 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbcsr_contract_ABA'
120 :
121 : INTEGER :: handle
122 : LOGICAL :: my_retain_sparsity
123 : TYPE(dbcsr_type) :: work
124 :
125 9800 : CALL timeset(routineN, handle)
126 :
127 9800 : my_retain_sparsity = .FALSE.
128 9800 : IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
129 :
130 9800 : CALL dbcsr_create(work, template=matrix_A)
131 :
132 9800 : IF (trans_A_left == "N" .AND. trans_A_right == "T") THEN
133 : ! C = A B A^T
134 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_A, matrix_B, &
135 5200 : 0.0_dp, work, filter_eps=eps_filter)
136 : CALL dbcsr_multiply("N", "T", 1.0_dp, work, matrix_A, &
137 : 0.0_dp, matrix_C, filter_eps=eps_filter, &
138 5200 : retain_sparsity=my_retain_sparsity)
139 4600 : ELSE IF (trans_A_left == "T" .AND. trans_A_right == "N") THEN
140 : ! C = A^T B A
141 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_B, matrix_A, &
142 4600 : 0.0_dp, work, filter_eps=eps_filter)
143 : CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_A, work, &
144 : 0.0_dp, matrix_C, filter_eps=eps_filter, &
145 4600 : retain_sparsity=my_retain_sparsity)
146 : ELSE
147 0 : CPABORT("Unsupported transposition pair in dbcsr_contract_ABA")
148 : END IF
149 9800 : CALL dbcsr_release(work)
150 :
151 9800 : CALL timestop(handle)
152 :
153 9800 : END SUBROUTINE dbcsr_contract_ABA
154 :
155 : END MODULE gw_utils_dbcsr
|
