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 full-matrix operations used by GW modules.
10 : !> \par History
11 : !> 09.2026 created Jan Wilhelm
12 : ! **************************************************************************************************
13 : MODULE gw_utils_fm
14 : USE cp_cfm_types, ONLY: cp_cfm_create,&
15 : cp_cfm_get_info,&
16 : cp_cfm_release,&
17 : cp_cfm_type
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
19 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
20 : cp_fm_cholesky_invert
21 : USE cp_fm_diag, ONLY: cp_fm_power
22 : USE cp_fm_types, ONLY: cp_fm_create,&
23 : cp_fm_get_info,&
24 : cp_fm_release,&
25 : cp_fm_to_fm,&
26 : cp_fm_type
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: z_one,&
29 : z_zero
30 : USE parallel_gemm_api, ONLY: parallel_gemm
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 : PRIVATE
35 :
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils_fm'
37 :
38 : PUBLIC :: cfm_contract_ABA, fm_contract_ABA, fm_invert, fm_sqrt
39 :
40 : INTERFACE fm_contract_ABA
41 : MODULE PROCEDURE fm_contract_ABA, fm_contract_ABA_set
42 : END INTERFACE fm_contract_ABA
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Computes A^H B A for complex full matrices.
48 : !> \param matrix_A left and right matrix A; H denotes its Hermitian transpose
49 : !> \param matrix_B input matrix B; overwritten by A^H B A if matrix_C is absent
50 : !> \param matrix_C optional output matrix C=A^H B A
51 : ! **************************************************************************************************
52 39720 : SUBROUTINE cfm_contract_ABA(matrix_A, matrix_B, matrix_C)
53 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_A
54 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_B
55 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: matrix_C
56 :
57 : INTEGER :: ncol, nrow
58 : TYPE(cp_cfm_type) :: work
59 :
60 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
61 13240 : CPASSERT(ASSOCIATED(matrix_B%matrix_struct))
62 13240 : CALL cp_cfm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
63 13240 : CPASSERT(nrow == ncol)
64 13240 : CALL cp_cfm_create(work, matrix_A%matrix_struct)
65 :
66 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, z_one, matrix_B, matrix_A, &
67 13240 : z_zero, work)
68 13240 : IF (PRESENT(matrix_C)) THEN
69 3196 : CPASSERT(ASSOCIATED(matrix_C%matrix_struct))
70 : CALL parallel_gemm("C", "N", nrow, nrow, nrow, z_one, matrix_A, work, &
71 3196 : z_zero, matrix_C)
72 : ELSE
73 : CALL parallel_gemm("C", "N", nrow, nrow, nrow, z_one, matrix_A, work, &
74 10044 : z_zero, matrix_B)
75 : END IF
76 :
77 13240 : CALL cp_cfm_release(work)
78 :
79 13240 : END SUBROUTINE cfm_contract_ABA
80 :
81 : ! **************************************************************************************************
82 : !> \brief Inverts a symmetric matrix. First, Cholesky decomposition is tried.
83 : !> If it fails, the matrix is diagonalized and inverted by taking
84 : !> inverse eigenvalues.
85 : !> \param matrix_A input A; output inverse or filtered pseudoinverse
86 : !> \param eigenvalue_threshold optional relative eigenvalue threshold
87 : !> \param unit_nr optional output unit for a failed Cholesky decomposition
88 : ! **************************************************************************************************
89 1962 : SUBROUTINE fm_invert(matrix_A, eigenvalue_threshold, unit_nr)
90 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_A
91 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eigenvalue_threshold
92 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
93 :
94 : INTEGER :: info, n_dependent
95 : LOGICAL :: cholesky_successful
96 : REAL(KIND=dp) :: threshold
97 : TYPE(cp_fm_type) :: work
98 :
99 654 : threshold = 0.0_dp
100 654 : IF (PRESENT(eigenvalue_threshold)) threshold = eigenvalue_threshold
101 654 : CALL cp_fm_create(work, matrix_A%matrix_struct)
102 :
103 : ! Preserve A because a failed Cholesky decomposition overwrites it.
104 654 : CALL cp_fm_to_fm(matrix_A, work)
105 654 : CALL cp_fm_cholesky_decompose(matrix_A, info_out=info)
106 654 : cholesky_successful = info == 0
107 : n_dependent = 0
108 654 : IF (cholesky_successful) THEN
109 654 : CALL cp_fm_cholesky_invert(matrix_A)
110 654 : CALL cp_fm_uplo_to_full(matrix_A, work)
111 : ELSE
112 0 : CALL cp_fm_to_fm(work, matrix_A)
113 0 : CALL cp_fm_power(matrix_A, work, -1.0_dp, threshold, n_dependent)
114 0 : IF (PRESENT(unit_nr)) THEN
115 0 : IF (unit_nr > 0) THEN
116 : WRITE (unit_nr, '(T2,A)') &
117 0 : 'Cholesky decomposition failed; matrix inverted by diagonalization.'
118 0 : WRITE (unit_nr, '(T2,A,T72,I9)') 'Discarded eigenmodes:', n_dependent
119 : END IF
120 : END IF
121 : END IF
122 654 : CALL cp_fm_release(work)
123 :
124 654 : END SUBROUTINE fm_invert
125 :
126 : ! **************************************************************************************************
127 : !> \brief For input A, computes B such that B^T B=A.
128 : !> First, Cholesky decomposition is tried. If it fails, B is computed
129 : !> by diagonalizing A.
130 : !> \param matrix_A symmetric input matrix A, retained unchanged
131 : !> \param matrix_B output matrix B
132 : !> \param eigenvalue_threshold optional relative eigenvalue threshold
133 : !> \param unit_nr optional output unit for a failed Cholesky decomposition
134 : ! **************************************************************************************************
135 120 : SUBROUTINE fm_sqrt(matrix_A, matrix_B, eigenvalue_threshold, unit_nr)
136 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
137 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_B
138 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eigenvalue_threshold
139 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
140 :
141 : INTEGER :: i_row, info, j_col, n_dependent, &
142 : ncol_local, nrow_local
143 40 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
144 : LOGICAL :: cholesky_successful
145 : REAL(KIND=dp) :: threshold
146 : TYPE(cp_fm_type) :: work
147 :
148 40 : threshold = 0.0_dp
149 40 : IF (PRESENT(eigenvalue_threshold)) threshold = eigenvalue_threshold
150 40 : CALL cp_fm_create(work, matrix_A%matrix_struct)
151 :
152 40 : CALL cp_fm_to_fm(matrix_A, matrix_B)
153 40 : CALL cp_fm_cholesky_decompose(matrix_B, info_out=info)
154 40 : cholesky_successful = info == 0
155 40 : IF (cholesky_successful) THEN
156 : n_dependent = 0
157 : CALL cp_fm_get_info(matrix=matrix_B, nrow_local=nrow_local, ncol_local=ncol_local, &
158 40 : row_indices=row_indices, col_indices=col_indices)
159 1208 : DO j_col = 1, ncol_local
160 21402 : DO i_row = 1, nrow_local
161 21362 : IF (row_indices(i_row) > col_indices(j_col)) THEN
162 9805 : matrix_B%local_data(i_row, j_col) = 0.0_dp
163 : END IF
164 : END DO
165 : END DO
166 : ELSE
167 0 : CALL cp_fm_to_fm(matrix_A, matrix_B)
168 0 : CALL cp_fm_power(matrix_B, work, 0.5_dp, threshold, n_dependent)
169 0 : IF (PRESENT(unit_nr)) THEN
170 0 : IF (unit_nr > 0) THEN
171 : WRITE (unit_nr, '(T2,A)') &
172 0 : 'Cholesky decomposition failed, using diagonalization.'
173 0 : WRITE (unit_nr, '(T2,A,T72,I9)') 'Discarded eigenmodes:', n_dependent
174 : END IF
175 : END IF
176 : END IF
177 40 : CALL cp_fm_release(work)
178 :
179 40 : END SUBROUTINE fm_sqrt
180 :
181 : ! **************************************************************************************************
182 : !> \brief Computes A^T B A.
183 : !> \param matrix_A left and right matrix A
184 : !> \param matrix_B input matrix B; overwritten by A^T B A if matrix_C is absent
185 : !> \param matrix_C optional output matrix C=A^T B A
186 : ! **************************************************************************************************
187 7026 : SUBROUTINE fm_contract_ABA(matrix_A, matrix_B, matrix_C)
188 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
189 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_B
190 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: matrix_C
191 :
192 : INTEGER :: ncol, nrow
193 : TYPE(cp_fm_type) :: work
194 :
195 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
196 2342 : CPASSERT(ASSOCIATED(matrix_B%matrix_struct))
197 2342 : CALL cp_fm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
198 2342 : CPASSERT(nrow == ncol)
199 2342 : CALL cp_fm_create(work, matrix_A%matrix_struct)
200 :
201 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, 1.0_dp, matrix_B, matrix_A, &
202 2342 : 0.0_dp, work)
203 2342 : IF (PRESENT(matrix_C)) THEN
204 2342 : CPASSERT(ASSOCIATED(matrix_C%matrix_struct))
205 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
206 2342 : 0.0_dp, matrix_C)
207 : ELSE
208 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
209 0 : 0.0_dp, matrix_B)
210 : END IF
211 :
212 2342 : CALL cp_fm_release(work)
213 :
214 2342 : END SUBROUTINE fm_contract_ABA
215 :
216 : ! **************************************************************************************************
217 : !> \brief Computes A^T B_i A for a set of matrices.
218 : !> \param matrix_A left and right matrix A
219 : !> \param matrix_B input matrices B_i; overwritten if matrix_C is absent
220 : !> \param matrix_C optional output matrices C_i=A^T B_i A
221 : ! **************************************************************************************************
222 248 : SUBROUTINE fm_contract_ABA_set(matrix_A, matrix_B, matrix_C)
223 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
224 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: matrix_B
225 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
226 : OPTIONAL :: matrix_C
227 :
228 : INTEGER :: i, ncol, nrow
229 : TYPE(cp_fm_type) :: work
230 :
231 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
232 248 : CALL cp_fm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
233 248 : CPASSERT(nrow == ncol)
234 248 : IF (PRESENT(matrix_C)) THEN
235 0 : CPASSERT(SIZE(matrix_C) == SIZE(matrix_B))
236 : END IF
237 248 : CALL cp_fm_create(work, matrix_A%matrix_struct)
238 :
239 2088 : DO i = 1, SIZE(matrix_B)
240 1840 : CPASSERT(ASSOCIATED(matrix_B(i)%matrix_struct))
241 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, 1.0_dp, matrix_B(i), matrix_A, &
242 1840 : 0.0_dp, work)
243 2088 : IF (PRESENT(matrix_C)) THEN
244 0 : CPASSERT(ASSOCIATED(matrix_C(i)%matrix_struct))
245 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
246 0 : 0.0_dp, matrix_C(i))
247 : ELSE
248 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
249 1840 : 0.0_dp, matrix_B(i))
250 : END IF
251 : END DO
252 :
253 248 : CALL cp_fm_release(work)
254 :
255 248 : END SUBROUTINE fm_contract_ABA_set
256 :
257 : END MODULE gw_utils_fm
|
