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 Full-matrix operations not provided by the CP2K FM packages.
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_struct, ONLY: cp_fm_struct_equivalent
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_get_info,&
25 : cp_fm_release,&
26 : cp_fm_to_fm,&
27 : cp_fm_type
28 : USE kinds, ONLY: dp
29 : USE mathconstants, ONLY: z_one,&
30 : z_zero
31 : USE mathlib, ONLY: diag_complex,&
32 : gemm_square
33 : USE parallel_gemm_api, ONLY: parallel_gemm
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 : PRIVATE
38 :
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils_fm'
40 :
41 : PUBLIC :: cfm_contract_ABA, fm_contract_ABA, fm_invert, fm_sqrt, &
42 : local_add_on_diag, local_complex_power
43 :
44 : INTERFACE fm_contract_ABA
45 : MODULE PROCEDURE fm_contract_ABA, fm_contract_ABA_set
46 : END INTERFACE fm_contract_ABA
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Computes A^H B A for complex full matrices.
52 : !> \param matrix_A left and right matrix A; H denotes its Hermitian transpose
53 : !> \param matrix_B input matrix B; overwritten by A^H B A if matrix_C is absent
54 : !> \param matrix_C optional output matrix C=A^H B A
55 : !> \param workspace optional reusable scratch matrix with the same layout as A;
56 : !> its contents are overwritten and it must not alias the input or output matrices
57 : ! **************************************************************************************************
58 44836 : SUBROUTINE cfm_contract_ABA(matrix_A, matrix_B, matrix_C, workspace)
59 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_A
60 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_B
61 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: matrix_C, workspace
62 :
63 : INTEGER :: ncol, nrow
64 : LOGICAL :: compatible_workspace
65 : TYPE(cp_cfm_type) :: work
66 :
67 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
68 22418 : CPASSERT(ASSOCIATED(matrix_B%matrix_struct))
69 22418 : CALL cp_cfm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
70 22418 : CPASSERT(nrow == ncol)
71 22418 : IF (PRESENT(workspace)) THEN
72 5458 : compatible_workspace = cp_fm_struct_equivalent(matrix_A%matrix_struct, workspace%matrix_struct)
73 : compatible_workspace = compatible_workspace .AND. &
74 54580 : ALL(matrix_A%matrix_struct%descriptor == workspace%matrix_struct%descriptor)
75 5458 : CPASSERT(compatible_workspace)
76 5458 : CPASSERT(ASSOCIATED(workspace%local_data))
77 5458 : CPASSERT(.NOT. ASSOCIATED(workspace%local_data, matrix_A%local_data))
78 5458 : CPASSERT(.NOT. ASSOCIATED(workspace%local_data, matrix_B%local_data))
79 5458 : IF (PRESENT(matrix_C)) THEN
80 0 : CPASSERT(.NOT. ASSOCIATED(workspace%local_data, matrix_C%local_data))
81 : END IF
82 5458 : work = workspace
83 : ELSE
84 16960 : CALL cp_cfm_create(work, matrix_A%matrix_struct)
85 : END IF
86 :
87 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, z_one, matrix_B, matrix_A, &
88 22418 : z_zero, work)
89 22418 : IF (PRESENT(matrix_C)) THEN
90 2496 : CPASSERT(ASSOCIATED(matrix_C%matrix_struct))
91 : CALL parallel_gemm("C", "N", nrow, nrow, nrow, z_one, matrix_A, work, &
92 2496 : z_zero, matrix_C)
93 : ELSE
94 : CALL parallel_gemm("C", "N", nrow, nrow, nrow, z_one, matrix_A, work, &
95 19922 : z_zero, matrix_B)
96 : END IF
97 :
98 22418 : IF (.NOT. PRESENT(workspace)) CALL cp_cfm_release(work)
99 :
100 22418 : END SUBROUTINE cfm_contract_ABA
101 :
102 : ! **************************************************************************************************
103 : !> \brief Add a scalar to the diagonal of a local complex matrix.
104 : !> \param matrix input/output matrix
105 : !> \param alpha scalar added to each diagonal element
106 : ! **************************************************************************************************
107 66560 : SUBROUTINE local_add_on_diag(matrix, alpha)
108 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
109 : COMPLEX(KIND=dp), INTENT(IN) :: alpha
110 :
111 : INTEGER :: i, n
112 :
113 66560 : n = SIZE(matrix, 1)
114 66560 : CPASSERT(n == SIZE(matrix, 2))
115 :
116 453120 : DO i = 1, n
117 453120 : matrix(i, i) = matrix(i, i) + alpha
118 : END DO
119 :
120 66560 : END SUBROUTINE local_add_on_diag
121 :
122 : ! **************************************************************************************************
123 : !> \brief Compute a spectral power of a local complex Hermitian matrix.
124 : !> Eigenvalues not larger than eps are discarded.
125 : !> \param matrix input/output matrix
126 : !> \param exponent matrix exponent
127 : !> \param eps eigenvalue threshold
128 : !> \param cond_nr optional condition number
129 : !> \param min_ev optional minimum absolute eigenvalue
130 : !> \param max_ev optional maximum absolute eigenvalue
131 : ! **************************************************************************************************
132 82880 : SUBROUTINE local_complex_power(matrix, exponent, eps, cond_nr, min_ev, max_ev)
133 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
134 : REAL(KIND=dp), INTENT(IN) :: exponent, eps
135 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: cond_nr, min_ev, max_ev
136 :
137 : CHARACTER(len=*), PARAMETER :: routineN = 'local_complex_power'
138 :
139 82880 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvectors
140 : INTEGER :: handle, i, n
141 : REAL(KIND=dp) :: sqrt_power_eval
142 82880 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
143 :
144 82880 : CALL timeset(routineN, handle)
145 :
146 82880 : n = SIZE(matrix, 1)
147 82880 : CPASSERT(n == SIZE(matrix, 2))
148 6687296 : matrix(:, :) = 0.5_dp*(matrix(:, :) + CONJG(TRANSPOSE(matrix(:, :))))
149 :
150 497280 : ALLOCATE (eigenvalues(n), eigenvectors(n, n))
151 82880 : CALL diag_complex(matrix, eigenvectors, eigenvalues)
152 :
153 143040 : IF (PRESENT(cond_nr)) cond_nr = MAXVAL(ABS(eigenvalues))/MINVAL(ABS(eigenvalues))
154 112960 : IF (PRESENT(min_ev)) min_ev = MINVAL(ABS(eigenvalues))
155 112960 : IF (PRESENT(max_ev)) max_ev = MAXVAL(ABS(eigenvalues))
156 :
157 : ! Reconstruct matrix**exponent as X*X^H, where
158 : ! X = U*diag(eigenvalues**(exponent/2)).
159 565264 : DO i = 1, n
160 482384 : IF (eps < eigenvalues(i)) THEN
161 482384 : sqrt_power_eval = eigenvalues(i)**(0.5_dp*exponent)
162 : ELSE
163 : sqrt_power_eval = 0.0_dp
164 : END IF
165 3385088 : eigenvectors(:, i) = eigenvectors(:, i)*sqrt_power_eval
166 : END DO
167 :
168 : ! X*X^H = U*diag(eigenvalues**exponent)*U^H.
169 82880 : CALL gemm_square(eigenvectors, 'N', eigenvectors, 'C', matrix)
170 :
171 82880 : DEALLOCATE (eigenvalues, eigenvectors)
172 :
173 82880 : CALL timestop(handle)
174 :
175 82880 : END SUBROUTINE local_complex_power
176 :
177 : ! **************************************************************************************************
178 : !> \brief Inverts a symmetric matrix. First, Cholesky decomposition is tried.
179 : !> If it fails, the matrix is diagonalized and inverted by taking
180 : !> inverse eigenvalues.
181 : !> \param matrix_A input A; output inverse or filtered pseudoinverse
182 : !> \param eigenvalue_threshold optional relative eigenvalue threshold
183 : !> \param unit_nr optional output unit for a failed Cholesky decomposition
184 : ! **************************************************************************************************
185 2004 : SUBROUTINE fm_invert(matrix_A, eigenvalue_threshold, unit_nr)
186 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_A
187 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eigenvalue_threshold
188 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
189 :
190 : INTEGER :: info, n_dependent
191 : LOGICAL :: cholesky_successful
192 : REAL(KIND=dp) :: threshold
193 : TYPE(cp_fm_type) :: work
194 :
195 668 : threshold = 0.0_dp
196 668 : IF (PRESENT(eigenvalue_threshold)) threshold = eigenvalue_threshold
197 668 : CALL cp_fm_create(work, matrix_A%matrix_struct)
198 :
199 : ! Preserve A because a failed Cholesky decomposition overwrites it.
200 668 : CALL cp_fm_to_fm(matrix_A, work)
201 668 : CALL cp_fm_cholesky_decompose(matrix_A, info_out=info)
202 668 : cholesky_successful = info == 0
203 : n_dependent = 0
204 668 : IF (cholesky_successful) THEN
205 668 : CALL cp_fm_cholesky_invert(matrix_A)
206 668 : CALL cp_fm_uplo_to_full(matrix_A, work)
207 : ELSE
208 0 : CALL cp_fm_to_fm(work, matrix_A)
209 0 : CALL cp_fm_power(matrix_A, work, -1.0_dp, threshold, n_dependent)
210 0 : IF (PRESENT(unit_nr)) THEN
211 0 : IF (unit_nr > 0) THEN
212 : WRITE (unit_nr, '(T2,A)') &
213 0 : 'Cholesky decomposition failed; matrix inverted by diagonalization.'
214 0 : WRITE (unit_nr, '(T2,A,T72,I9)') 'Discarded eigenmodes:', n_dependent
215 : END IF
216 : END IF
217 : END IF
218 668 : CALL cp_fm_release(work)
219 :
220 668 : END SUBROUTINE fm_invert
221 :
222 : ! **************************************************************************************************
223 : !> \brief For input A, computes B such that B^T B=A.
224 : !> First, Cholesky decomposition is tried. If it fails, B is computed
225 : !> by diagonalizing A.
226 : !> \param matrix_A symmetric input matrix A, retained unchanged
227 : !> \param matrix_B output matrix B
228 : !> \param eigenvalue_threshold optional relative eigenvalue threshold
229 : !> \param unit_nr optional output unit for a failed Cholesky decomposition
230 : ! **************************************************************************************************
231 120 : SUBROUTINE fm_sqrt(matrix_A, matrix_B, eigenvalue_threshold, unit_nr)
232 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
233 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_B
234 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eigenvalue_threshold
235 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
236 :
237 : INTEGER :: i_row, info, j_col, n_dependent, &
238 : ncol_local, nrow_local
239 40 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
240 : LOGICAL :: cholesky_successful
241 : REAL(KIND=dp) :: threshold
242 : TYPE(cp_fm_type) :: work
243 :
244 40 : threshold = 0.0_dp
245 40 : IF (PRESENT(eigenvalue_threshold)) threshold = eigenvalue_threshold
246 40 : CALL cp_fm_create(work, matrix_A%matrix_struct)
247 :
248 40 : CALL cp_fm_to_fm(matrix_A, matrix_B)
249 40 : CALL cp_fm_cholesky_decompose(matrix_B, info_out=info)
250 40 : cholesky_successful = info == 0
251 40 : IF (cholesky_successful) THEN
252 : n_dependent = 0
253 : CALL cp_fm_get_info(matrix=matrix_B, nrow_local=nrow_local, ncol_local=ncol_local, &
254 40 : row_indices=row_indices, col_indices=col_indices)
255 1208 : DO j_col = 1, ncol_local
256 21402 : DO i_row = 1, nrow_local
257 21362 : IF (row_indices(i_row) > col_indices(j_col)) THEN
258 9805 : matrix_B%local_data(i_row, j_col) = 0.0_dp
259 : END IF
260 : END DO
261 : END DO
262 : ELSE
263 0 : CALL cp_fm_to_fm(matrix_A, matrix_B)
264 0 : CALL cp_fm_power(matrix_B, work, 0.5_dp, threshold, n_dependent)
265 0 : IF (PRESENT(unit_nr)) THEN
266 0 : IF (unit_nr > 0) THEN
267 : WRITE (unit_nr, '(T2,A)') &
268 0 : 'Cholesky decomposition failed, using diagonalization.'
269 0 : WRITE (unit_nr, '(T2,A,T72,I9)') 'Discarded eigenmodes:', n_dependent
270 : END IF
271 : END IF
272 : END IF
273 40 : CALL cp_fm_release(work)
274 :
275 40 : END SUBROUTINE fm_sqrt
276 :
277 : ! **************************************************************************************************
278 : !> \brief Computes A^T B A.
279 : !> \param matrix_A left and right matrix A
280 : !> \param matrix_B input matrix B; overwritten by A^T B A if matrix_C is absent
281 : !> \param matrix_C optional output matrix C=A^T B A
282 : ! **************************************************************************************************
283 7086 : SUBROUTINE fm_contract_ABA(matrix_A, matrix_B, matrix_C)
284 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
285 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_B
286 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: matrix_C
287 :
288 : INTEGER :: ncol, nrow
289 : TYPE(cp_fm_type) :: work
290 :
291 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
292 2362 : CPASSERT(ASSOCIATED(matrix_B%matrix_struct))
293 2362 : CALL cp_fm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
294 2362 : CPASSERT(nrow == ncol)
295 2362 : CALL cp_fm_create(work, matrix_A%matrix_struct)
296 :
297 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, 1.0_dp, matrix_B, matrix_A, &
298 2362 : 0.0_dp, work)
299 2362 : IF (PRESENT(matrix_C)) THEN
300 2348 : CPASSERT(ASSOCIATED(matrix_C%matrix_struct))
301 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
302 2348 : 0.0_dp, matrix_C)
303 : ELSE
304 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
305 14 : 0.0_dp, matrix_B)
306 : END IF
307 :
308 2362 : CALL cp_fm_release(work)
309 :
310 2362 : END SUBROUTINE fm_contract_ABA
311 :
312 : ! **************************************************************************************************
313 : !> \brief Computes A^T B_i A for a set of matrices.
314 : !> \param matrix_A left and right matrix A
315 : !> \param matrix_B input matrices B_i; overwritten if matrix_C is absent
316 : !> \param matrix_C optional output matrices C_i=A^T B_i A
317 : ! **************************************************************************************************
318 274 : SUBROUTINE fm_contract_ABA_set(matrix_A, matrix_B, matrix_C)
319 : TYPE(cp_fm_type), INTENT(IN) :: matrix_A
320 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: matrix_B
321 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
322 : OPTIONAL :: matrix_C
323 :
324 : INTEGER :: i, ncol, nrow
325 : TYPE(cp_fm_type) :: work
326 :
327 0 : CPASSERT(ASSOCIATED(matrix_A%matrix_struct))
328 274 : CALL cp_fm_get_info(matrix_A, nrow_global=nrow, ncol_global=ncol)
329 274 : CPASSERT(nrow == ncol)
330 274 : IF (PRESENT(matrix_C)) THEN
331 0 : CPASSERT(SIZE(matrix_C) == SIZE(matrix_B))
332 : END IF
333 274 : CALL cp_fm_create(work, matrix_A%matrix_struct)
334 :
335 2348 : DO i = 1, SIZE(matrix_B)
336 2074 : CPASSERT(ASSOCIATED(matrix_B(i)%matrix_struct))
337 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, 1.0_dp, matrix_B(i), matrix_A, &
338 2074 : 0.0_dp, work)
339 2348 : IF (PRESENT(matrix_C)) THEN
340 0 : CPASSERT(ASSOCIATED(matrix_C(i)%matrix_struct))
341 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
342 0 : 0.0_dp, matrix_C(i))
343 : ELSE
344 : CALL parallel_gemm("T", "N", nrow, nrow, nrow, 1.0_dp, matrix_A, work, &
345 2074 : 0.0_dp, matrix_B(i))
346 : END IF
347 : END DO
348 :
349 274 : CALL cp_fm_release(work)
350 :
351 274 : END SUBROUTINE fm_contract_ABA_set
352 :
353 : END MODULE gw_utils_fm
|
