(git:3f0d10b)
Loading...
Searching...
No Matches
preconditioner_apply.F
Go to the documentation of this file.
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 computes preconditioners, and implements methods to apply them
10!> currently used in qs_ot
11!> \par History
12!> - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
13!> \author Joost VandeVondele (09.2002)
14! **************************************************************************************************
15MODULE preconditioner_apply
16 USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm
17 USE cp_cfm_types, ONLY: cp_cfm_create,&
18 cp_cfm_get_info,&
19 cp_cfm_release,&
20 cp_cfm_to_fm,&
21 cp_cfm_type,&
22 cp_fm_to_cfm
23 USE cp_dbcsr_api, ONLY: &
24 dbcsr_add, dbcsr_copy, dbcsr_get_info, dbcsr_init_p, dbcsr_iterator_blocks_left, &
25 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
26 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 copy_dbcsr_to_local,&
29 copy_fm_to_dbcsr,&
30 copy_local_to_dbcsr,&
31 cp_dbcsr_m_by_n_from_template
32 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
33 cp_fm_scale_and_add
34 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_restore
35 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 cp_fm_struct_release,&
37 cp_fm_struct_type
38 USE cp_fm_types, ONLY: cp_fm_create,&
39 cp_fm_get_info,&
40 cp_fm_release,&
41 cp_fm_to_fm,&
42 cp_fm_type
43 USE input_constants, ONLY: &
44 ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, ot_precond_fermi_low_rank, &
45 ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_full_kinetic, &
46 ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, &
47 ot_precond_solver_chebyshev, ot_precond_solver_direct, ot_precond_solver_inv_chol, &
48 ot_precond_solver_update
49 USE kinds, ONLY: dp
50 USE lattice_preconditioner_operator, ONLY: apply_lattice_inverse_dense,&
51 apply_lattice_state_inverse_dense
52 USE local_gemm_api, ONLY: local_gemm_pu_gpu,&
53 local_gemm_ctxt_type
54 USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight
55 USE mathconstants, ONLY: z_one,&
56 z_zero
57 USE parallel_gemm_api, ONLY: parallel_gemm
58 USE preconditioner_types, ONLY: lattice_rhs_distribution_type,&
59 preconditioner_p_type,&
60 preconditioner_type,&
61 release_lattice_rhs_distribution
62
63!$ USE OMP_LIB, ONLY: omp_get_max_threads
64#include "./base/base_uses.f90"
65
66 IMPLICIT NONE
67 PRIVATE
68
69 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
70
71 PUBLIC :: apply_preconditioner_cfm_complex, apply_preconditioner_dbcsr_complex, &
72 apply_preconditioner_fm, apply_preconditioner_dbcsr, apply_preconditioner_complex_batch
73
74 ! Reused across channels in each worker's batch; no FM/CFM lifecycle in the local path.
75 TYPE complex_preconditioner_work_type
76 COMPLEX(dp), ALLOCATABLE :: buffer(:, :, :)
77 REAL(dp), ALLOCATABLE :: packed(:, :, :)
78 TYPE(local_gemm_ctxt_type) :: gemm_ctx
79 END TYPE complex_preconditioner_work_type
80
81CONTAINS
82
83! **************************************************************************************************
84!> \brief Apply single-rank complex preconditioners in bounded OpenMP batches.
85!> Packing and DBCSR writeback stay outside OpenMP; each worker owns its dense scratch.
86!> \param preconditioners channel preconditioners
87!> \param matrix_in real/imaginary input pointers, indexed by component and channel
88!> \param matrix_out real/imaginary output pointers, indexed by component and channel
89!> \param handled false when a channel requires the distributed implementation
90!> \param output_scale optional real scale per channel, applied before writeback
91! **************************************************************************************************
92 SUBROUTINE apply_preconditioner_complex_batch(preconditioners, matrix_in, matrix_out, handled, output_scale)
93 TYPE(preconditioner_p_type), INTENT(IN) :: preconditioners(:)
94 TYPE(dbcsr_p_type), INTENT(IN) :: matrix_in(:, :), matrix_out(:, :)
95 LOGICAL, INTENT(OUT) :: handled
96 REAL(dp), INTENT(IN), OPTIONAL :: output_scale(:)
97
98 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_preconditioner_complex_batch'
99
100 INTEGER :: first, handle, i, j, k, n, nbatch, &
101 nchannels, nworkers
102 TYPE(complex_preconditioner_work_type), &
103 ALLOCATABLE :: work(:)
104 TYPE(preconditioner_type), POINTER :: p
105
106 handled = .false.
107 nchannels = SIZE(preconditioners)
108 IF (nchannels == 0) RETURN
109 DO i = 1, nchannels
110 p => preconditioners(i)%preconditioner
111 IF (.NOT. ASSOCIATED(p)) RETURN
112 IF (.NOT. ASSOCIATED(p%complex_fm)) RETURN
113 IF (p%complex_fm%matrix_struct%para_env%num_pe /= 1) RETURN
114 SELECT CASE (p%in_use)
115 CASE (ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_fermi_low_rank, &
116 ot_precond_full_single, ot_precond_full_single_inverse, &
117 ot_precond_full_kinetic, ot_precond_s_inverse)
118 CASE DEFAULT
119 RETURN
120 END SELECT
121 END DO
122
123 CALL timeset(routinen, handle)
124 nworkers = 1
125!$ nworkers = MIN(nchannels, omp_get_max_threads())
126 ALLOCATE (work(nworkers))
127 DO j = 1, nworkers
128 CALL work(j)%gemm_ctx%create(local_gemm_pu_gpu, timing=.false.)
129 END DO
130 DO first = 1, nchannels, nworkers
131 nbatch = min(nworkers, nchannels - first + 1)
132 DO j = 1, nbatch
133 i = first + j - 1
134 p => preconditioners(i)%preconditioner
135 CALL dbcsr_get_info(matrix_in(1, i)%matrix, nfullrows_total=n, nfullcols_total=k)
136 cpassert(n == p%complex_fm%matrix_struct%nrow_global)
137 IF (ALLOCATED(work(j)%buffer)) THEN
138 IF (SIZE(work(j)%buffer, 1) /= n .OR. SIZE(work(j)%buffer, 2) /= k) THEN
139 DEALLOCATE (work(j)%buffer, work(j)%packed)
140 END IF
141 END IF
142 IF (.NOT. ALLOCATED(work(j)%buffer)) THEN
143 ALLOCATE (work(j)%buffer(n, k, 2), work(j)%packed(n, k, 2))
144 END IF
145 CALL copy_dbcsr_to_local(matrix_in(1, i)%matrix, work(j)%packed(:, :, 1))
146 CALL copy_dbcsr_to_local(matrix_in(2, i)%matrix, work(j)%packed(:, :, 2))
147 END DO
148!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) &
149!$OMP SHARED(work,preconditioners,first,nbatch,output_scale) PRIVATE(i)
150 DO j = 1, nbatch
151 i = first + j - 1
152 CALL apply_preconditioner_complex_local(preconditioners(i)%preconditioner, work(j))
153 IF (PRESENT(output_scale)) THEN
154 IF (output_scale(i) /= 1.0_dp) work(j)%packed(:, :, :) = output_scale(i)*work(j)%packed(:, :, :)
155 END IF
156 END DO
157!$OMP END PARALLEL DO
158 DO j = 1, nbatch
159 i = first + j - 1
160 CALL copy_local_to_dbcsr(matrix_out(1, i)%matrix, work(j)%packed(:, :, 1))
161 CALL copy_local_to_dbcsr(matrix_out(2, i)%matrix, work(j)%packed(:, :, 2))
162 END DO
163 END DO
164 DO j = 1, nworkers
165 CALL work(j)%gemm_ctx%destroy()
166 END DO
167 handled = .true.
168 CALL timestop(handle)
169 END SUBROUTINE apply_preconditioner_complex_batch
170
171! **************************************************************************************************
172!> \brief Apply the existing complex preconditioner with private local GEMM storage.
173!> \param p channel preconditioner; its matrices are read only
174!> \param w worker storage
175! **************************************************************************************************
176 SUBROUTINE apply_preconditioner_complex_local(p, w)
177 TYPE(preconditioner_type), INTENT(IN) :: p
178 TYPE(complex_preconditioner_work_type), &
179 INTENT(INOUT), TARGET :: w
180
181 COMPLEX(dp), CONTIGUOUS, POINTER :: basis(:, :), result(:, :), rhs(:, :), &
182 rotation(:, :)
183 INTEGER :: k, n
184
185 NULLIFY (rotation)
186 n = SIZE(w%buffer, 1)
187 k = SIZE(w%buffer, 2)
188 IF (n == 0 .OR. k == 0) RETURN
189 basis => p%complex_fm%local_data
190 rhs => w%buffer(:, :, 1)
191 result => w%buffer(:, :, 2)
192 rhs(:, :) = cmplx(w%packed(:, :, 1), w%packed(:, :, 2), kind=dp)
193 SELECT CASE (p%in_use)
194 CASE (ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_fermi_low_rank)
195 cpassert(ASSOCIATED(p%full_evals))
196 cpassert(ASSOCIATED(p%occ_evals))
197 cpassert(n == SIZE(p%full_evals))
198 IF (p%in_use /= ot_precond_fermi_low_rank) THEN
199 cpassert(k == SIZE(p%occ_evals))
200 END IF
201 IF (p%in_use == ot_precond_full_all_covariant) THEN
202 cpassert(ASSOCIATED(p%occ_rotation_complex))
203 rotation => p%occ_rotation_complex%local_data
204 CALL w%gemm_ctx%gemm('N', 'N', n, k, k, z_one, rhs, n, &
205 rotation, SIZE(rotation, 1), z_zero, result, n)
206 rhs => w%buffer(:, :, 2)
207 result => w%buffer(:, :, 1)
208 END IF
209 CALL w%gemm_ctx%gemm('C', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
210 rhs, n, z_zero, result, n)
211 CALL scale_complex_spectral(p, result, p%complex_fm%matrix_struct%row_indices, &
212 p%complex_fm%matrix_struct%col_indices)
213 CALL w%gemm_ctx%gemm('N', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
214 result, n, z_zero, rhs, n)
215 IF (p%in_use == ot_precond_full_all_covariant) THEN
216 CALL w%gemm_ctx%gemm('N', 'C', n, k, k, z_one, rhs, n, &
217 rotation, SIZE(rotation, 1), z_zero, result, n)
218 ELSE
219 result => rhs
220 END IF
221 CASE DEFAULT
222 CALL w%gemm_ctx%gemm('N', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
223 rhs, n, z_zero, result, n)
224 END SELECT
225 w%packed(:, :, 1) = real(result, kind=dp)
226 w%packed(:, :, 2) = aimag(result)
227 END SUBROUTINE apply_preconditioner_complex_local
228
229! **************************************************************************************************
230!> \brief Scale spectral coefficients identically for local and distributed applications.
231!> \param p spectral preconditioner
232!> \param DATA local spectral coefficients
233!> \param rows global row indices
234!> \param cols global column indices
235! **************************************************************************************************
236 SUBROUTINE scale_complex_spectral(p, DATA, rows, cols)
237 TYPE(preconditioner_type), INTENT(IN) :: p
238 COMPLEX(dp), INTENT(INOUT) :: data(:, :)
239 INTEGER, INTENT(IN) :: rows(:), cols(:)
240
241 INTEGER :: i, j, nocc, rank_used
242 REAL(dp) :: scale
243
244 IF (p%in_use == ot_precond_fermi_low_rank) THEN
245 nocc = SIZE(p%occ_evals)
246 rank_used = min(p%spectral_rank, SIZE(p%full_evals) - nocc)
247 DO i = 1, SIZE(DATA, 1)
248 scale = p%spectral_base_scale
249 IF (rows(i) > nocc .AND. rows(i) <= nocc + rank_used) THEN
250 scale = low_rank_inverse_weight(p%full_evals(rows(i)), p%spectral_reference, p%energy_gap)
251 END IF
252 DATA(i, :) = scale*DATA(i, :)
253 END DO
254 ELSE
255 DO j = 1, SIZE(DATA, 2)
256 DO i = 1, SIZE(DATA, 1)
257 scale = 1.0_dp/max(p%energy_gap, p%full_evals(rows(i)) - p%occ_evals(cols(j)))
258 DATA(i, j) = scale*DATA(i, j)
259 END DO
260 END DO
261 END IF
262 END SUBROUTINE scale_complex_spectral
263
264! **************************************************************************************************
265!> \brief applies a previously created preconditioner to a full matrix
266!> \param preconditioner_env ...
267!> \param matrix_in ...
268!> \param matrix_out ...
269! **************************************************************************************************
270 SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
271
272 TYPE(preconditioner_type) :: preconditioner_env
273 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
274
275 CHARACTER(len=*), PARAMETER :: routinen = 'apply_preconditioner_fm'
276
277 INTEGER :: handle
278
279 CALL timeset(routinen, handle)
280
281 SELECT CASE (preconditioner_env%in_use)
282 CASE (0)
283 cpabort("No preconditioner in use")
284 CASE (ot_precond_full_single)
285 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
286 CASE (ot_precond_full_all)
287 CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
288 CASE (ot_precond_fermi_low_rank)
289 CALL apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
290 CASE (ot_precond_full_all_covariant)
291 CALL apply_full_all_covariant(preconditioner_env, matrix_in, matrix_out)
292 CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
293 IF (preconditioner_env%lattice_fft_active) THEN
294 CALL apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out)
295 ELSE
296 SELECT CASE (preconditioner_env%solver)
297 CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
298 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
299 CASE (ot_precond_solver_direct)
300 CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
301 CASE (ot_precond_solver_chebyshev)
302 CALL apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
303 CASE DEFAULT
304 cpabort("Solver not implemented")
305 END SELECT
306 END IF
307 CASE DEFAULT
308 cpabort("Unknown preconditioner")
309 END SELECT
310
311 CALL timestop(handle)
312
313 END SUBROUTINE apply_preconditioner_fm
314
315! **************************************************************************************************
316!> \brief ...
317!> \param preconditioner_env ...
318!> \param matrix_in ...
319!> \param matrix_out ...
320! **************************************************************************************************
321 SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
322
323 TYPE(preconditioner_type) :: preconditioner_env
324 TYPE(dbcsr_type) :: matrix_in, matrix_out
325
326 CHARACTER(len=*), PARAMETER :: routinen = 'apply_preconditioner_dbcsr'
327
328 INTEGER :: handle
329
330 CALL timeset(routinen, handle)
331
332 SELECT CASE (preconditioner_env%in_use)
333 CASE (0)
334 cpabort("No preconditioner in use")
335 CASE (ot_precond_full_single)
336 CALL apply_single(preconditioner_env, matrix_in, matrix_out)
337 CASE (ot_precond_full_all)
338 CALL apply_all(preconditioner_env, matrix_in, matrix_out)
339 CASE (ot_precond_fermi_low_rank)
340 CALL apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
341 CASE (ot_precond_full_all_covariant)
342 CALL apply_all_covariant(preconditioner_env, matrix_in, matrix_out)
343 CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
344 IF (preconditioner_env%lattice_fft_active) THEN
345 CALL apply_lattice_operator_dbcsr(preconditioner_env, matrix_in, matrix_out)
346 ELSE
347 SELECT CASE (preconditioner_env%solver)
348 CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
349 CALL apply_single(preconditioner_env, matrix_in, matrix_out)
350 CASE (ot_precond_solver_direct)
351 cpabort("Apply_full_direct not supported with ot")
352 !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
353 CASE (ot_precond_solver_chebyshev)
354 CALL apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
355 CASE DEFAULT
356 cpabort("Wrong solver")
357 END SELECT
358 END IF
359 CASE DEFAULT
360 cpabort("Wrong preconditioner")
361 END SELECT
362
363 CALL timestop(handle)
364
365 END SUBROUTINE apply_preconditioner_dbcsr
366
367! **************************************************************************************************
368!> \brief Apply a complex k-point orbital preconditioner.
369!> \param preconditioner_env complex preconditioner storage
370!> \param matrix_in complex input channel
371!> \param matrix_out complex output channel
372! **************************************************************************************************
373 SUBROUTINE apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
374
375 TYPE(preconditioner_type) :: preconditioner_env
376 TYPE(cp_cfm_type), INTENT(IN) :: matrix_in
377 TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_out
378
379 CHARACTER(len=*), PARAMETER :: routinen = 'apply_preconditioner_cfm_complex'
380
381 INTEGER :: handle, k, n, npre
382 TYPE(cp_cfm_type) :: canonical_in, canonical_out
383
384 CALL timeset(routinen, handle)
385
386 SELECT CASE (preconditioner_env%in_use)
387 CASE (ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_all_covariant, &
388 ot_precond_full_single, ot_precond_full_single_inverse, &
389 ot_precond_full_kinetic, ot_precond_s_inverse)
390 CASE DEFAULT
391 cpabort("Unsupported complex K-point OT preconditioner")
392 END SELECT
393 cpassert(ASSOCIATED(preconditioner_env%complex_fm))
394 CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
395 CALL cp_cfm_get_info(preconditioner_env%complex_fm, nrow_global=npre)
396 cpassert(n == npre)
397
398 SELECT CASE (preconditioner_env%in_use)
399 CASE (ot_precond_full_all, ot_precond_fermi_low_rank)
400 CALL apply_complex_spectral(preconditioner_env, matrix_in, matrix_out)
401 CASE (ot_precond_full_all_covariant)
402 cpassert(ASSOCIATED(preconditioner_env%occ_rotation_complex))
403 CALL cp_cfm_create(canonical_in, matrix_in%matrix_struct, &
404 name='complex covariant FULL_ALL input')
405 CALL cp_cfm_create(canonical_out, matrix_out%matrix_struct, &
406 name='complex covariant FULL_ALL output')
407 CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_in, &
408 preconditioner_env%occ_rotation_complex, z_zero, canonical_in)
409 CALL apply_complex_spectral(preconditioner_env, canonical_in, canonical_out)
410 CALL cp_cfm_gemm('N', 'C', n, k, k, z_one, canonical_out, &
411 preconditioner_env%occ_rotation_complex, z_zero, matrix_out)
412 CALL cp_cfm_release(canonical_out)
413 CALL cp_cfm_release(canonical_in)
414 CASE DEFAULT
415 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
416 matrix_in, z_zero, matrix_out)
417 END SELECT
418
419 CALL timestop(handle)
420
421 END SUBROUTINE apply_preconditioner_cfm_complex
422
423! **************************************************************************************************
424!> \brief Apply the spectral inverse in its canonical complex basis.
425!> \param preconditioner_env complex spectral preconditioner
426!> \param matrix_in input orbital-gradient matrix
427!> \param matrix_out preconditioned orbital-gradient matrix
428! **************************************************************************************************
429 SUBROUTINE apply_complex_spectral(preconditioner_env, matrix_in, matrix_out)
430
431 TYPE(preconditioner_type) :: preconditioner_env
432 TYPE(cp_cfm_type), INTENT(IN) :: matrix_in
433 TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_out
434
435 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
436 POINTER :: local_data
437 INTEGER :: k, n, ncol_local, nrow_local
438 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
439 TYPE(cp_cfm_type) :: matrix_spectral
440
441 cpassert(ASSOCIATED(preconditioner_env%full_evals))
442 cpassert(ASSOCIATED(preconditioner_env%occ_evals))
443 CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
444 cpassert(n == SIZE(preconditioner_env%full_evals))
445 IF (preconditioner_env%in_use /= ot_precond_fermi_low_rank) THEN
446 cpassert(k == SIZE(preconditioner_env%occ_evals))
447 END IF
448 CALL cp_cfm_create(matrix_spectral, matrix_in%matrix_struct, &
449 name='complex spectral input')
450 CALL cp_cfm_gemm('C', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
451 matrix_in, z_zero, matrix_spectral)
452 CALL cp_cfm_get_info(matrix_spectral, nrow_local=nrow_local, ncol_local=ncol_local, &
453 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
454 CALL scale_complex_spectral(preconditioner_env, local_data(1:nrow_local, 1:ncol_local), &
455 row_indices, col_indices)
456 CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
457 matrix_spectral, z_zero, matrix_out)
458 CALL cp_cfm_release(matrix_spectral)
459
460 END SUBROUTINE apply_complex_spectral
461
462! **************************************************************************************************
463!> \brief Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
464!> \param preconditioner_env complex k-point preconditioner
465!> \param matrix_in_re real input channel
466!> \param matrix_in_im imaginary input channel
467!> \param matrix_out_re real output channel
468!> \param matrix_out_im imaginary output channel
469! **************************************************************************************************
470 SUBROUTINE apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, &
471 matrix_out_re, matrix_out_im)
472
473 TYPE(preconditioner_type), TARGET :: preconditioner_env
474 TYPE(dbcsr_type), TARGET :: matrix_in_re, matrix_in_im, &
475 matrix_out_re, matrix_out_im
476
477 CHARACTER(len=*), PARAMETER :: routinen = 'apply_preconditioner_dbcsr_complex'
478
479 INTEGER :: handle, k, n
480 LOGICAL :: handled
481 TYPE(cp_cfm_type) :: matrix_in, matrix_out
482 TYPE(cp_fm_type) :: matrix_in_im_fm, matrix_in_re_fm, &
483 matrix_out_im_fm, matrix_out_re_fm
484 TYPE(dbcsr_p_type) :: input(2, 1), output(2, 1)
485 TYPE(preconditioner_p_type) :: preconditioners(1)
486
487 preconditioners(1)%preconditioner => preconditioner_env
488 input(1, 1)%matrix => matrix_in_re
489 input(2, 1)%matrix => matrix_in_im
490 output(1, 1)%matrix => matrix_out_re
491 output(2, 1)%matrix => matrix_out_im
492 CALL apply_preconditioner_complex_batch(preconditioners, input, output, handled)
493 IF (handled) RETURN
494
495 CALL timeset(routinen, handle)
496
497 CALL dbcsr_get_info(matrix_in_re, nfullrows_total=n, nfullcols_total=k)
498 cpassert(ASSOCIATED(preconditioner_env%complex_fm))
499
500 CALL cp_fm_create(matrix_in_re_fm, preconditioner_env%complex_fm%matrix_struct, &
501 nrow=n, ncol=k, name='complex preconditioner input real')
502 CALL cp_fm_create(matrix_in_im_fm, preconditioner_env%complex_fm%matrix_struct, &
503 nrow=n, ncol=k, name='complex preconditioner input imaginary')
504 CALL cp_fm_create(matrix_out_re_fm, preconditioner_env%complex_fm%matrix_struct, &
505 nrow=n, ncol=k, name='complex preconditioner output real')
506 CALL cp_fm_create(matrix_out_im_fm, preconditioner_env%complex_fm%matrix_struct, &
507 nrow=n, ncol=k, name='complex preconditioner output imaginary')
508 CALL copy_dbcsr_to_fm(matrix_in_re, matrix_in_re_fm)
509 CALL copy_dbcsr_to_fm(matrix_in_im, matrix_in_im_fm)
510
511 CALL cp_cfm_create(matrix_in, matrix_in_re_fm%matrix_struct, &
512 name='complex preconditioner input')
513 CALL cp_cfm_create(matrix_out, matrix_in_re_fm%matrix_struct, &
514 name='complex preconditioner output')
515 CALL cp_fm_to_cfm(matrix_in_re_fm, matrix_in_im_fm, matrix_in)
516 CALL apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
517 CALL cp_cfm_to_fm(matrix_out, matrix_out_re_fm, matrix_out_im_fm)
518 CALL dbcsr_set(matrix_out_re, 0.0_dp)
519 CALL dbcsr_set(matrix_out_im, 0.0_dp)
520 CALL copy_fm_to_dbcsr(matrix_out_re_fm, matrix_out_re)
521 CALL copy_fm_to_dbcsr(matrix_out_im_fm, matrix_out_im)
522
523 CALL cp_cfm_release(matrix_out)
524 CALL cp_cfm_release(matrix_in)
525 CALL cp_fm_release(matrix_out_im_fm)
526 CALL cp_fm_release(matrix_out_re_fm)
527 CALL cp_fm_release(matrix_in_im_fm)
528 CALL cp_fm_release(matrix_in_re_fm)
529
530 CALL timestop(handle)
531
532 END SUBROUTINE apply_preconditioner_dbcsr_complex
533
534! **************************************************************************************************
535!> \brief apply to full matrix, complete inversion has already been done
536!> \param preconditioner_env ...
537!> \param matrix_in ...
538!> \param matrix_out ...
539! **************************************************************************************************
540 SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
541
542 TYPE(preconditioner_type) :: preconditioner_env
543 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
544
545 CHARACTER(len=*), PARAMETER :: routinen = 'apply_full_single'
546
547 INTEGER :: handle, k, n
548
549 CALL timeset(routinen, handle)
550
551 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
552 CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
553 matrix_in, 0.0_dp, matrix_out)
554 CALL timestop(handle)
555
556 END SUBROUTINE apply_full_single
557
558! **************************************************************************************************
559!> \brief apply to dbcsr matrix, complete inversion has already been done
560!> \param preconditioner_env ...
561!> \param matrix_in ...
562!> \param matrix_out ...
563! **************************************************************************************************
564 SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
565
566 TYPE(preconditioner_type) :: preconditioner_env
567 TYPE(dbcsr_type) :: matrix_in, matrix_out
568
569 CHARACTER(len=*), PARAMETER :: routinen = 'apply_single'
570
571 INTEGER :: handle
572
573 CALL timeset(routinen, handle)
574
575 IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
576 cpabort("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
577 END IF
578 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
579 0.0_dp, matrix_out)
580
581 CALL timestop(handle)
582
583 END SUBROUTINE apply_single
584
585! **************************************************************************************************
586!> \brief preconditioner contains the factorization, application done by
587!> solving the linear system
588!> \param preconditioner_env ...
589!> \param matrix_in ...
590!> \param matrix_out ...
591! **************************************************************************************************
592 SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
593
594 TYPE(preconditioner_type) :: preconditioner_env
595 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
596
597 CHARACTER(len=*), PARAMETER :: routinen = 'apply_full_direct'
598
599 INTEGER :: handle, k, n
600 TYPE(cp_fm_type) :: work
601
602 CALL timeset(routinen, handle)
603
604 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
605 CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single")
606 CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
607 & "SOLVE", transa="T")
608 CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
609 & "SOLVE", transa="N")
610 CALL cp_fm_release(work)
611
612 CALL timestop(handle)
613
614 END SUBROUTINE apply_full_direct
615
616! **************************************************************************************************
617!> \brief Apply a Chebyshev approximation to the inverse of the stored dense SPD operator.
618!> \param preconditioner_env ...
619!> \param matrix_in ...
620!> \param matrix_out ...
621! **************************************************************************************************
622 SUBROUTINE apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
623
624 TYPE(preconditioner_type) :: preconditioner_env
625 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
626
627 INTEGER :: degree, iteration, k, n
628 REAL(kind=dp) :: delta, rho, rho_previous, sigma, theta
629 TYPE(cp_fm_type) :: direction, residual
630
631 cpassert(ASSOCIATED(preconditioner_env%fm))
632 degree = preconditioner_env%polynomial_degree
633 cpassert(degree >= 1)
634 theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
635 delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
636 sigma = theta/delta
637 rho_previous = 1.0_dp/sigma
638 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
639 CALL cp_fm_create(direction, matrix_in%matrix_struct, name='Chebyshev direction')
640 CALL cp_fm_create(residual, matrix_in%matrix_struct, name='Chebyshev residual')
641 CALL cp_fm_to_fm(matrix_in, direction)
642 CALL cp_fm_scale(1.0_dp/theta, direction)
643 CALL cp_fm_to_fm(direction, matrix_out)
644 DO iteration = 2, degree
645 CALL cp_fm_to_fm(matrix_in, residual)
646 CALL parallel_gemm('N', 'N', n, k, n, -1.0_dp, preconditioner_env%fm, &
647 matrix_out, 1.0_dp, residual)
648 rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
649 CALL cp_fm_scale(rho*rho_previous, direction)
650 CALL cp_fm_scale_and_add(1.0_dp, direction, 2.0_dp*rho/delta, residual)
651 CALL cp_fm_scale_and_add(1.0_dp, matrix_out, 1.0_dp, direction)
652 rho_previous = rho
653 END DO
654 CALL cp_fm_release(residual)
655 CALL cp_fm_release(direction)
656
657 END SUBROUTINE apply_chebyshev_fm
658
659! **************************************************************************************************
660!> \brief Apply a Chebyshev approximation to the inverse of the stored sparse SPD operator.
661!> \param preconditioner_env ...
662!> \param matrix_in ...
663!> \param matrix_out ...
664! **************************************************************************************************
665 SUBROUTINE apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
666
667 TYPE(preconditioner_type) :: preconditioner_env
668 TYPE(dbcsr_type) :: matrix_in, matrix_out
669
670 INTEGER :: degree, iteration
671 REAL(kind=dp) :: delta, rho, rho_previous, sigma, theta
672 TYPE(dbcsr_type) :: direction, residual
673
674 cpassert(ASSOCIATED(preconditioner_env%dbcsr_matrix))
675 degree = preconditioner_env%polynomial_degree
676 cpassert(degree >= 1)
677 theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
678 delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
679 sigma = theta/delta
680 rho_previous = 1.0_dp/sigma
681 CALL dbcsr_copy(direction, matrix_in, name='Chebyshev direction')
682 CALL dbcsr_copy(residual, matrix_in, name='Chebyshev residual')
683 CALL dbcsr_set(matrix_out, 0.0_dp)
684 CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp/theta)
685 CALL dbcsr_set(direction, 0.0_dp)
686 CALL dbcsr_add(direction, matrix_in, 1.0_dp, 1.0_dp/theta)
687 DO iteration = 2, degree
688 CALL dbcsr_set(residual, 0.0_dp)
689 CALL dbcsr_add(residual, matrix_in, 1.0_dp, 1.0_dp)
690 CALL dbcsr_multiply('N', 'N', -1.0_dp, preconditioner_env%dbcsr_matrix, &
691 matrix_out, 1.0_dp, residual)
692 rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
693 CALL dbcsr_add(direction, residual, rho*rho_previous, 2.0_dp*rho/delta)
694 CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp)
695 rho_previous = rho
696 END DO
697 CALL dbcsr_release(residual)
698 CALL dbcsr_release(direction)
699
700 END SUBROUTINE apply_chebyshev_dbcsr
701
702! **************************************************************************************************
703!> \brief Applies the overlap-inverse base and a bounded low-rank spectral correction.
704!> \param preconditioner_env ...
705!> \param matrix_in ...
706!> \param matrix_out ...
707! **************************************************************************************************
708 SUBROUTINE apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
709
710 TYPE(preconditioner_type) :: preconditioner_env
711 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
712
713 CHARACTER(len=*), PARAMETER :: routinen = 'apply_low_rank_fm'
714
715 INTEGER :: handle, i, j, k, n, ncol_local, nocc, &
716 nrow_local, rank_used
717 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
718 REAL(kind=dp) :: correction
719 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
720 POINTER :: local_data
721 TYPE(cp_fm_type) :: matrix_tmp
722
723 CALL timeset(routinen, handle)
724
725 cpassert(ASSOCIATED(preconditioner_env%fm))
726 cpassert(ASSOCIATED(preconditioner_env%full_evals))
727 cpassert(ASSOCIATED(preconditioner_env%occ_evals))
728
729 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
730 nocc = SIZE(preconditioner_env%occ_evals)
731 cpassert(n == SIZE(preconditioner_env%full_evals))
732 rank_used = min(preconditioner_env%spectral_rank, n - nocc)
733
734 ! With the exact base V V^T = S^-1. The optional lattice base replaces
735 ! S^-1 by the inverse of its block-circulant translational projection:
736 ! P = B^-1/window + V_r diag(1/max(gap, epsilon_a-mu) - 1/window) V_r^T.
737 ! The common reference mu makes the operator right-covariant for any number
738 ! of orbital columns. The lattice path retains only non-negative corrections.
739 SELECT CASE (preconditioner_env%low_rank_base)
740 CASE (ot_low_rank_base_overlap)
741 cpassert(ASSOCIATED(preconditioner_env%base_fm))
742 CALL parallel_gemm('N', 'N', n, k, n, preconditioner_env%spectral_base_scale, &
743 preconditioner_env%base_fm, matrix_in, 0.0_dp, matrix_out)
744 CASE (ot_low_rank_base_lattice_fft)
745 CALL apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out, &
746 preconditioner_env%spectral_base_scale)
747 CASE DEFAULT
748 cpabort("Unknown FERMI_LOW_RANK base operator during application")
749 END SELECT
750
751 IF (rank_used > 0) THEN
752 CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name=routinen, &
753 nrow=rank_used, ncol=k)
754 CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
755 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
756
757 CALL parallel_gemm('T', 'N', rank_used, k, n, 1.0_dp, preconditioner_env%fm, &
758 matrix_in, 0.0_dp, matrix_tmp, a_first_col=nocc + 1)
759
760 DO j = 1, ncol_local
761 DO i = 1, nrow_local
762 correction = low_rank_inverse_weight( &
763 preconditioner_env%full_evals(nocc + row_indices(i)), &
764 preconditioner_env%spectral_reference, &
765 preconditioner_env%energy_gap) &
766 - preconditioner_env%spectral_base_scale
767 local_data(i, j) = local_data(i, j)*correction
768 END DO
769 END DO
770
771 CALL parallel_gemm('N', 'N', n, k, rank_used, 1.0_dp, preconditioner_env%fm, &
772 matrix_tmp, 1.0_dp, matrix_out, a_first_col=nocc + 1)
773 CALL cp_fm_release(matrix_tmp)
774 END IF
775
776 CALL timestop(handle)
777
778 END SUBROUTINE apply_low_rank_fm
779
780! **************************************************************************************************
781!> \brief Apply a replicated-cell inverse operator to a distributed full matrix.
782!> \param preconditioner_env ...
783!> \param matrix_in ...
784!> \param matrix_out ...
785!> \param scale_factor ...
786! **************************************************************************************************
787 SUBROUTINE apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out, scale_factor)
788
789 TYPE(preconditioner_type) :: preconditioner_env
790 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
791 REAL(kind=dp), INTENT(IN), OPTIONAL :: scale_factor
792
793 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_lattice_operator_fm'
794
795 INTEGER :: fft_ok, handle, i
796 LOGICAL :: used_fft
797 REAL(kind=dp) :: imaginary_residual, output_max, scale
798 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
799 coarse_residual, coarse_rhs, &
800 local_input, local_output
801
802 CALL timeset(routinen, handle)
803 cpassert(ALLOCATED(preconditioner_env%lattice_inverse_k))
804 scale = 1.0_dp
805 IF (PRESENT(scale_factor)) scale = scale_factor
806 CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
807 preconditioner_env%lattice_rhs_distribution, local_input)
808 ALLOCATE (local_output(SIZE(local_input, 1), SIZE(local_input, 2)))
809 IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
810 cpassert(ALLOCATED(preconditioner_env%lattice_local_operator))
811 cpassert(ALLOCATED(preconditioner_env%lattice_local_coarse_inverse))
812 ALLOCATE (base_input(SIZE(local_input, 1), SIZE(local_input, 2)), &
813 coarse_projected(SIZE(preconditioner_env%lattice_local_indices), SIZE(local_input, 2)), &
814 coarse_rhs(SIZE(preconditioner_env%lattice_local_indices), SIZE(local_input, 2)))
815 DO i = 1, SIZE(preconditioner_env%lattice_local_indices)
816 coarse_projected(i, :) = local_input(preconditioner_env%lattice_local_indices(i), :)
817 END DO
818 coarse_rhs(:, :) = matmul(preconditioner_env%lattice_local_coarse_inverse, &
819 coarse_projected)
820 base_input(:, :) = local_input - &
821 matmul(preconditioner_env%lattice_local_operator, coarse_rhs)
822 ELSE
823 ALLOCATE (base_input(SIZE(local_input, 1), SIZE(local_input, 2)), source=local_input)
824 END IF
825 used_fft = .true.
826 imaginary_residual = 0.0_dp
827 IF (SIZE(local_input, 2) > 0) THEN
828 CALL apply_lattice_inverse_dense(preconditioner_env%lattice_inverse_k, &
829 preconditioner_env%lattice_dims, base_input, local_output, &
830 used_fft, imaginary_residual)
831 END IF
832 IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
833 ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
834 coarse_projected(:, :) = matmul(transpose(preconditioner_env%lattice_local_operator), &
835 local_output)
836 coarse_residual(:, :) = matmul(preconditioner_env%lattice_local_coarse_inverse, &
837 coarse_projected)
838 local_output(preconditioner_env%lattice_local_indices, :) = &
839 local_output(preconditioner_env%lattice_local_indices, :) + coarse_rhs - coarse_residual
840 DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
841 END IF
842 fft_ok = merge(1, 0, used_fft)
843 output_max = 0.0_dp
844 IF (SIZE(local_output) > 0) output_max = maxval(abs(local_output))
845 CALL preconditioner_env%para_env%min(fft_ok)
846 CALL preconditioner_env%para_env%max(imaginary_residual)
847 CALL preconditioner_env%para_env%max(output_max)
848 IF (fft_ok == 0) cpabort("Lattice preconditioner inverse failed to use FFTs")
849 IF (imaginary_residual > 1.0e-10_dp*max(1.0_dp, output_max)) THEN
850 cpabort("Lattice preconditioner inverse produced a non-real Gamma-point result")
851 END IF
852 CALL collect_lattice_rhs(preconditioner_env, preconditioner_env%lattice_rhs_distribution, &
853 local_output, matrix_out, scale)
854 DEALLOCATE (base_input, local_input, local_output)
855 CALL timestop(handle)
856
857 END SUBROUTINE apply_lattice_operator_fm
858
859! **************************************************************************************************
860!> \brief Apply state-dependent replicated-cell inverse operators to a distributed full matrix.
861!> \param preconditioner_env ...
862!> \param matrix_in ...
863!> \param matrix_out ...
864! **************************************************************************************************
865 SUBROUTINE apply_lattice_state_operator_fm(preconditioner_env, matrix_in, matrix_out)
866
867 TYPE(preconditioner_type) :: preconditioner_env
868 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
869
870 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_lattice_state_operator_fm'
871
872 INTEGER :: fft_ok, global_column, handle, i, j, k, &
873 n, state
874 LOGICAL :: used_fft
875 REAL(kind=dp) :: imaginary_residual, output_max
876 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
877 coarse_residual, coarse_rhs, complement_input, local_input, local_output, &
878 occupied_component, occupied_update
879
880 CALL timeset(routinen, handle)
881 cpassert(ALLOCATED(preconditioner_env%lattice_state_inverse_k))
882 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
883 cpassert(SIZE(preconditioner_env%lattice_state_inverse_k, 4) == k)
884 CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
885 preconditioner_env%lattice_rhs_distribution, local_input)
886 cpassert(ALLOCATED(preconditioner_env%lattice_occ_vectors))
887 cpassert(ALLOCATED(preconditioner_env%lattice_occ_dual))
888 ALLOCATE (complement_input(n, SIZE(local_input, 2)), &
889 local_output(n, SIZE(local_input, 2)), &
890 occupied_component(k, SIZE(local_input, 2)), &
891 occupied_update(n, SIZE(local_input, 2)))
892 occupied_component(:, :) = matmul(transpose(preconditioner_env%lattice_occ_vectors), &
893 local_input)
894 complement_input(:, :) = local_input - &
895 matmul(preconditioner_env%lattice_occ_dual, occupied_component)
896 ALLOCATE (base_input(n, SIZE(local_input, 2)), source=complement_input)
897 IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
898 cpassert(ALLOCATED(preconditioner_env%lattice_state_local_operator))
899 cpassert(ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse))
900 ALLOCATE (coarse_rhs(SIZE(preconditioner_env%lattice_state_local_indices, 1), &
901 SIZE(local_input, 2)), &
902 coarse_projected(SIZE(preconditioner_env%lattice_state_local_indices, 1), &
903 SIZE(local_input, 2)))
904 DO j = 1, SIZE(local_input, 2)
905 state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
906 DO i = 1, SIZE(preconditioner_env%lattice_state_local_indices, 1)
907 coarse_projected(i, j) = &
908 complement_input(preconditioner_env%lattice_state_local_indices(i, state), j)
909 END DO
910 coarse_rhs(:, j) = matmul( &
911 preconditioner_env%lattice_state_local_coarse_inverse(:, :, state), &
912 coarse_projected(:, j))
913 base_input(:, j) = complement_input(:, j) - &
914 matmul(preconditioner_env%lattice_state_local_operator(:, :, state), &
915 coarse_rhs(:, j))
916 END DO
917 END IF
918 used_fft = .true.
919 imaginary_residual = 0.0_dp
920 IF (SIZE(local_input, 2) > 0) THEN
921 CALL apply_lattice_state_inverse_dense(preconditioner_env%lattice_state_inverse_k, &
922 preconditioner_env%lattice_dims, &
923 base_input, local_output, used_fft, imaginary_residual, &
924 preconditioner_env%lattice_rhs_distribution%global_columns)
925 END IF
926 IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
927 ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
928 DO j = 1, SIZE(local_input, 2)
929 state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
930 coarse_projected(:, j) = &
931 matmul(transpose(preconditioner_env%lattice_state_local_operator(:, :, state)), &
932 local_output(:, j))
933 coarse_residual(:, j) = matmul( &
934 preconditioner_env%lattice_state_local_coarse_inverse(:, :, state), &
935 coarse_projected(:, j))
936 local_output(preconditioner_env%lattice_state_local_indices(:, state), j) = &
937 local_output(preconditioner_env%lattice_state_local_indices(:, state), j) + &
938 coarse_rhs(:, j) - coarse_residual(:, j)
939 END DO
940 DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
941 END IF
942 fft_ok = merge(1, 0, used_fft)
943 output_max = 0.0_dp
944 IF (SIZE(local_output) > 0) output_max = maxval(abs(local_output))
945 CALL preconditioner_env%para_env%min(fft_ok)
946 CALL preconditioner_env%para_env%max(imaginary_residual)
947 CALL preconditioner_env%para_env%max(output_max)
948 IF (fft_ok == 0) THEN
949 cpabort("State-dependent lattice preconditioner inverse failed to use FFTs")
950 END IF
951 IF (imaginary_residual > 1.0e-10_dp*max(1.0_dp, output_max)) THEN
952 cpabort("State-dependent lattice preconditioner produced a non-real Gamma-point result")
953 END IF
954 occupied_component(:, :) = matmul(transpose(preconditioner_env%lattice_occ_dual), &
955 local_output)
956 occupied_update(:, :) = matmul(preconditioner_env%lattice_occ_vectors, occupied_component)
957 local_output(:, :) = local_output - occupied_update
958 occupied_component(:, :) = matmul(transpose(preconditioner_env%lattice_occ_vectors), &
959 local_input)
960 DO j = 1, SIZE(local_input, 2)
961 global_column = preconditioner_env%lattice_rhs_distribution%global_columns(j)
962 DO i = 1, k
963 occupied_component(i, j) = occupied_component(i, j)/max( &
964 preconditioner_env%energy_gap, &
965 preconditioner_env%full_evals(i) - &
966 preconditioner_env%occ_evals(global_column))
967 END DO
968 END DO
969 occupied_update(:, :) = matmul(preconditioner_env%lattice_occ_vectors, occupied_component)
970 local_output(:, :) = local_output + occupied_update
971 CALL collect_lattice_rhs(preconditioner_env, preconditioner_env%lattice_rhs_distribution, &
972 local_output, matrix_out)
973 DEALLOCATE (base_input, complement_input, local_input, local_output, occupied_component, occupied_update)
974 CALL timestop(handle)
975
976 END SUBROUTINE apply_lattice_state_operator_fm
977
978! **************************************************************************************************
979!> \brief Route complete orbital columns to distinct ranks for lattice-FFT application.
980!> \param preconditioner_env ...
981!> \param matrix_in distributed input matrix
982!> \param distribution communication map retained for the reverse route
983!> \param local_input complete columns owned by this rank
984! **************************************************************************************************
985 SUBROUTINE distribute_lattice_rhs(preconditioner_env, matrix_in, distribution, local_input)
986
987 TYPE(preconditioner_type) :: preconditioner_env
988 TYPE(cp_fm_type), INTENT(IN) :: matrix_in
989 TYPE(lattice_rhs_distribution_type), INTENT(INOUT) :: distribution
990 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
991 INTENT(OUT) :: local_input
992
993 INTEGER :: destination, i, ientry, j, k, n, &
994 ncol_local, nowned, nrow_local, &
995 num_pe, slot
996 INTEGER, ALLOCATABLE, DIMENSION(:) :: cursor, send_global_columns, &
997 send_global_rows
998 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
999 LOGICAL :: rebuild
1000 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_values, send_values
1001 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1002 POINTER :: input_data
1003
1004 num_pe = preconditioner_env%para_env%num_pe
1005 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k, &
1006 nrow_local=nrow_local, ncol_local=ncol_local, &
1007 row_indices=row_indices, col_indices=col_indices, local_data=input_data)
1008 rebuild = distribution%nrow_global /= n .OR. distribution%ncol_global /= k .OR. &
1009 distribution%nrow_local /= nrow_local .OR. distribution%ncol_local /= ncol_local .OR. &
1010 distribution%num_pe /= num_pe
1011 IF (.NOT. rebuild) THEN
1012 rebuild = .NOT. ALLOCATED(distribution%local_rows) .OR. &
1013 .NOT. ALLOCATED(distribution%local_columns)
1014 END IF
1015 IF (.NOT. rebuild) THEN
1016 rebuild = any(distribution%local_rows /= row_indices) .OR. &
1017 any(distribution%local_columns /= col_indices)
1018 END IF
1019 IF (rebuild) THEN
1020 CALL release_lattice_rhs_distribution(distribution)
1021 distribution%nrow_global = n
1022 distribution%ncol_global = k
1023 distribution%nrow_local = nrow_local
1024 distribution%ncol_local = ncol_local
1025 distribution%num_pe = num_pe
1026 ALLOCATE (distribution%local_rows(nrow_local), distribution%local_columns(ncol_local))
1027 distribution%local_rows(:) = row_indices
1028 distribution%local_columns(:) = col_indices
1029 ALLOCATE (distribution%send_counts(num_pe), distribution%send_displacements(num_pe), &
1030 distribution%recv_counts(num_pe), distribution%recv_displacements(num_pe), &
1031 cursor(num_pe))
1032 distribution%send_counts = 0
1033 DO j = 1, ncol_local
1034 destination = modulo(col_indices(j) - 1, num_pe) + 1
1035 distribution%send_counts(destination) = &
1036 distribution%send_counts(destination) + nrow_local
1037 END DO
1038 CALL preconditioner_env%para_env%alltoall(distribution%send_counts, &
1039 distribution%recv_counts, 1)
1040 distribution%send_displacements(1) = 0
1041 distribution%recv_displacements(1) = 0
1042 DO i = 2, num_pe
1043 distribution%send_displacements(i) = distribution%send_displacements(i - 1) + &
1044 distribution%send_counts(i - 1)
1045 distribution%recv_displacements(i) = distribution%recv_displacements(i - 1) + &
1046 distribution%recv_counts(i - 1)
1047 END DO
1048
1049 ALLOCATE (send_global_rows(sum(distribution%send_counts)), &
1050 send_global_columns(sum(distribution%send_counts)), &
1051 distribution%send_local_rows(sum(distribution%send_counts)), &
1052 distribution%send_local_columns(sum(distribution%send_counts)), &
1053 distribution%recv_global_rows(sum(distribution%recv_counts)), &
1054 distribution%recv_global_columns(sum(distribution%recv_counts)))
1055 cursor(:) = distribution%send_displacements
1056 DO j = 1, ncol_local
1057 destination = modulo(col_indices(j) - 1, num_pe) + 1
1058 DO i = 1, nrow_local
1059 cursor(destination) = cursor(destination) + 1
1060 ientry = cursor(destination)
1061 send_global_rows(ientry) = row_indices(i)
1062 send_global_columns(ientry) = col_indices(j)
1063 distribution%send_local_rows(ientry) = i
1064 distribution%send_local_columns(ientry) = j
1065 END DO
1066 END DO
1067 CALL preconditioner_env%para_env%alltoall( &
1068 send_global_rows, distribution%send_counts, distribution%send_displacements, &
1069 distribution%recv_global_rows, distribution%recv_counts, distribution%recv_displacements)
1070 CALL preconditioner_env%para_env%alltoall( &
1071 send_global_columns, distribution%send_counts, distribution%send_displacements, &
1072 distribution%recv_global_columns, distribution%recv_counts, distribution%recv_displacements)
1073
1074 nowned = 0
1075 DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
1076 nowned = nowned + 1
1077 END DO
1078 ALLOCATE (distribution%global_columns(nowned))
1079 nowned = 0
1080 DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
1081 nowned = nowned + 1
1082 distribution%global_columns(nowned) = j
1083 END DO
1084 DEALLOCATE (cursor, send_global_columns, send_global_rows)
1085 END IF
1086
1087 nowned = SIZE(distribution%global_columns)
1088 ALLOCATE (local_input(n, nowned), source=0.0_dp)
1089 ALLOCATE (send_values(sum(distribution%send_counts)), &
1090 recv_values(sum(distribution%recv_counts)))
1091 DO ientry = 1, SIZE(send_values)
1092 send_values(ientry) = input_data(distribution%send_local_rows(ientry), &
1093 distribution%send_local_columns(ientry))
1094 END DO
1095 CALL preconditioner_env%para_env%alltoall( &
1096 send_values, distribution%send_counts, distribution%send_displacements, &
1097 recv_values, distribution%recv_counts, distribution%recv_displacements)
1098 DO ientry = 1, SIZE(recv_values)
1099 slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
1100 local_input(distribution%recv_global_rows(ientry), slot) = recv_values(ientry)
1101 END DO
1102
1103 DEALLOCATE (recv_values, send_values)
1104
1105 END SUBROUTINE distribute_lattice_rhs
1106
1107! **************************************************************************************************
1108!> \brief Return rank-owned lattice-FFT columns to their original full-matrix layout.
1109!> \param preconditioner_env ...
1110!> \param distribution communication map from distribute_lattice_rhs
1111!> \param local_output complete preconditioned columns owned by this rank
1112!> \param matrix_out distributed output matrix
1113!> \param scale_factor optional final scaling
1114! **************************************************************************************************
1115 SUBROUTINE collect_lattice_rhs(preconditioner_env, distribution, local_output, matrix_out, scale_factor)
1116
1117 TYPE(preconditioner_type) :: preconditioner_env
1118 TYPE(lattice_rhs_distribution_type), INTENT(IN) :: distribution
1119 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: local_output
1120 TYPE(cp_fm_type), INTENT(IN) :: matrix_out
1121 REAL(kind=dp), INTENT(IN), OPTIONAL :: scale_factor
1122
1123 INTEGER :: ientry, num_pe, slot
1124 REAL(kind=dp) :: scale
1125 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_values, send_values
1126 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1127 POINTER :: output_data
1128
1129 num_pe = preconditioner_env%para_env%num_pe
1130 cpassert(SIZE(local_output, 1) == distribution%nrow_global)
1131 cpassert(SIZE(local_output, 2) == SIZE(distribution%global_columns))
1132 scale = 1.0_dp
1133 IF (PRESENT(scale_factor)) scale = scale_factor
1134 ALLOCATE (send_values(sum(distribution%recv_counts)), &
1135 recv_values(sum(distribution%send_counts)))
1136 DO ientry = 1, SIZE(send_values)
1137 slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
1138 send_values(ientry) = local_output(distribution%recv_global_rows(ientry), slot)
1139 END DO
1140 CALL preconditioner_env%para_env%alltoall( &
1141 send_values, distribution%recv_counts, distribution%recv_displacements, &
1142 recv_values, distribution%send_counts, distribution%send_displacements)
1143
1144 CALL cp_fm_get_info(matrix_out, local_data=output_data)
1145 output_data = 0.0_dp
1146 DO ientry = 1, SIZE(recv_values)
1147 output_data(distribution%send_local_rows(ientry), &
1148 distribution%send_local_columns(ientry)) = scale*recv_values(ientry)
1149 END DO
1150 DEALLOCATE (recv_values, send_values)
1151
1152 END SUBROUTINE collect_lattice_rhs
1153
1154! **************************************************************************************************
1155!> \brief Bridge a DBCSR OT gradient to the full-matrix lattice-FFT operator.
1156!> \param preconditioner_env ...
1157!> \param matrix_in ...
1158!> \param matrix_out ...
1159! **************************************************************************************************
1160 SUBROUTINE apply_lattice_operator_dbcsr(preconditioner_env, matrix_in, matrix_out)
1161
1162 TYPE(preconditioner_type) :: preconditioner_env
1163 TYPE(dbcsr_type) :: matrix_in, matrix_out
1164
1165 INTEGER :: k, n
1166 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1167 TYPE(cp_fm_type) :: input_fm, output_fm
1168
1169 NULLIFY (fm_struct)
1170 CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1171 CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1172 context=preconditioner_env%ctxt, &
1173 para_env=preconditioner_env%para_env)
1174 CALL cp_fm_create(input_fm, fm_struct, name="lattice preconditioner input")
1175 CALL cp_fm_create(output_fm, fm_struct, name="lattice preconditioner output")
1176 CALL cp_fm_struct_release(fm_struct)
1177 CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1178 CALL apply_lattice_operator_fm(preconditioner_env, input_fm, output_fm)
1179 CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1180 CALL cp_fm_release(input_fm)
1181 CALL cp_fm_release(output_fm)
1182
1183 END SUBROUTINE apply_lattice_operator_dbcsr
1184
1185! **************************************************************************************************
1186!> \brief DBCSR variant of the bounded low-rank spectral correction.
1187!> \param preconditioner_env ...
1188!> \param matrix_in ...
1189!> \param matrix_out ...
1190! **************************************************************************************************
1191 SUBROUTINE apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
1192
1193 TYPE(preconditioner_type) :: preconditioner_env
1194 TYPE(dbcsr_type) :: matrix_in, matrix_out
1195
1196 CHARACTER(len=*), PARAMETER :: routinen = 'apply_low_rank_dbcsr'
1197
1198 INTEGER :: col, col_offset, col_size, handle, i, j, &
1199 k, n, nocc, rank_used, row, &
1200 row_offset, row_size
1201 REAL(kind=dp) :: correction
1202 REAL(kind=dp), DIMENSION(:, :), POINTER :: DATA
1203 TYPE(dbcsr_iterator_type) :: iter
1204 TYPE(dbcsr_type) :: matrix_tmp
1205
1206 CALL timeset(routinen, handle)
1207
1208 IF (preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
1209 CALL apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
1210 CALL timestop(handle)
1211 RETURN
1212 END IF
1213
1214 cpassert(ASSOCIATED(preconditioner_env%base_dbcsr_matrix))
1215 cpassert(ASSOCIATED(preconditioner_env%dbcsr_matrix))
1216 cpassert(ASSOCIATED(preconditioner_env%full_evals))
1217 cpassert(ASSOCIATED(preconditioner_env%occ_evals))
1218
1219 CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1220 nocc = SIZE(preconditioner_env%occ_evals)
1221 cpassert(n == SIZE(preconditioner_env%full_evals))
1222 rank_used = min(preconditioner_env%spectral_rank, n - nocc)
1223
1224 ! DBCSR form of the same overlap-inverse base plus positive low-rank update.
1225 CALL dbcsr_multiply('N', 'N', preconditioner_env%spectral_base_scale, &
1226 preconditioner_env%base_dbcsr_matrix, matrix_in, 0.0_dp, matrix_out)
1227
1228 IF (rank_used > 0) THEN
1229 CALL dbcsr_copy(matrix_tmp, matrix_in, name=routinen)
1230 CALL dbcsr_set(matrix_tmp, 0.0_dp)
1231 CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1232 matrix_in, 0.0_dp, matrix_tmp, &
1233 first_row=nocc + 1, last_row=nocc + rank_used)
1234
1235 CALL dbcsr_iterator_start(iter, matrix_tmp)
1236 DO WHILE (dbcsr_iterator_blocks_left(iter))
1237 CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
1238 row_size=row_size, col_size=col_size, &
1239 row_offset=row_offset, col_offset=col_offset)
1240 DO j = 1, col_size
1241 DO i = 1, row_size
1242 IF (row_offset + i - 1 >= nocc + 1 .AND. &
1243 row_offset + i - 1 <= nocc + rank_used) THEN
1244 correction = low_rank_inverse_weight( &
1245 preconditioner_env%full_evals(row_offset + i - 1), &
1246 preconditioner_env%spectral_reference, &
1247 preconditioner_env%energy_gap) &
1248 - preconditioner_env%spectral_base_scale
1249 DATA(i, j) = DATA(i, j)*correction
1250 END IF
1251 END DO
1252 END DO
1253 END DO
1254 CALL dbcsr_iterator_stop(iter)
1255
1256 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1257 matrix_tmp, 1.0_dp, matrix_out, &
1258 first_k=nocc + 1, last_k=nocc + rank_used)
1259 CALL dbcsr_release(matrix_tmp)
1260 END IF
1261
1262 CALL timestop(handle)
1263
1264 END SUBROUTINE apply_low_rank_dbcsr
1265
1266! **************************************************************************************************
1267!> \brief Bridge DBCSR OT gradients to the full-matrix lattice-FFT implementation.
1268!> \param preconditioner_env ...
1269!> \param matrix_in ...
1270!> \param matrix_out ...
1271! **************************************************************************************************
1272 SUBROUTINE apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
1273
1274 TYPE(preconditioner_type) :: preconditioner_env
1275 TYPE(dbcsr_type) :: matrix_in, matrix_out
1276
1277 INTEGER :: k, n
1278 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1279 TYPE(cp_fm_type) :: input_fm, output_fm
1280
1281 NULLIFY (fm_struct)
1282 CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1283 CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1284 context=preconditioner_env%ctxt, &
1285 para_env=preconditioner_env%para_env)
1286 CALL cp_fm_create(input_fm, fm_struct, name="lattice low-rank input")
1287 CALL cp_fm_create(output_fm, fm_struct, name="lattice low-rank output")
1288 CALL cp_fm_struct_release(fm_struct)
1289 CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1290 CALL apply_low_rank_fm(preconditioner_env, input_fm, output_fm)
1291 CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1292 CALL cp_fm_release(input_fm)
1293 CALL cp_fm_release(output_fm)
1294
1295 END SUBROUTINE apply_low_rank_dbcsr_via_fm
1296
1297! **************************************************************************************************
1298!> \brief full all to a full matrix
1299!> \param preconditioner_env ...
1300!> \param matrix_in ...
1301!> \param matrix_out ...
1302! **************************************************************************************************
1303 SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
1304
1305 TYPE(preconditioner_type) :: preconditioner_env
1306 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
1307
1308 CHARACTER(len=*), PARAMETER :: routinen = 'apply_full_all'
1309
1310 INTEGER :: handle, i, j, k, n, ncol_local, &
1311 nrow_local
1312 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1313 REAL(kind=dp) :: dum
1314 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1315 POINTER :: local_data
1316 TYPE(cp_fm_type) :: matrix_tmp
1317
1318 CALL timeset(routinen, handle)
1319
1320 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1321
1322 CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
1323 CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
1324 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
1325
1326 !
1327 CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
1328 matrix_in, 0.0_dp, matrix_tmp)
1329
1330 ! do the right scaling
1331 DO j = 1, ncol_local
1332 DO i = 1, nrow_local
1333 dum = 1.0_dp/max(preconditioner_env%energy_gap, &
1334 preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
1335 local_data(i, j) = local_data(i, j)*dum
1336 END DO
1337 END DO
1338
1339 ! mult back
1340 CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
1341 matrix_tmp, 0.0_dp, matrix_out)
1342
1343 CALL cp_fm_release(matrix_tmp)
1344
1345 CALL timestop(handle)
1346
1347 END SUBROUTINE apply_full_all
1348
1349! **************************************************************************************************
1350!> \brief full all to a dbcsr matrix
1351!> \param preconditioner_env ...
1352!> \param matrix_in ...
1353!> \param matrix_out ...
1354! **************************************************************************************************
1355 SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
1356
1357 TYPE(preconditioner_type) :: preconditioner_env
1358 TYPE(dbcsr_type) :: matrix_in, matrix_out
1359
1360 CHARACTER(len=*), PARAMETER :: routinen = 'apply_all'
1361
1362 INTEGER :: col, col_offset, col_size, handle, i, j, &
1363 row, row_offset, row_size
1364 REAL(kind=dp) :: dum
1365 REAL(kind=dp), DIMENSION(:, :), POINTER :: DATA
1366 TYPE(dbcsr_iterator_type) :: iter
1367 TYPE(dbcsr_type) :: matrix_tmp
1368
1369 CALL timeset(routinen, handle)
1370
1371 CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
1372 CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1373 matrix_in, 0.0_dp, matrix_tmp)
1374 ! do the right scaling
1375 CALL dbcsr_iterator_start(iter, matrix_tmp)
1376 DO WHILE (dbcsr_iterator_blocks_left(iter))
1377 CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
1378 row_size=row_size, col_size=col_size, &
1379 row_offset=row_offset, col_offset=col_offset)
1380 DO j = 1, col_size
1381 DO i = 1, row_size
1382 dum = 1.0_dp/max(preconditioner_env%energy_gap, &
1383 preconditioner_env%full_evals(row_offset + i - 1) &
1384 - preconditioner_env%occ_evals(col_offset + j - 1))
1385 DATA(i, j) = DATA(i, j)*dum
1386 END DO
1387 END DO
1388 END DO
1389 CALL dbcsr_iterator_stop(iter)
1390
1391 ! mult back
1392 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1393 matrix_tmp, 0.0_dp, matrix_out)
1394 CALL dbcsr_release(matrix_tmp)
1395 CALL timestop(handle)
1396
1397 END SUBROUTINE apply_all
1398
1399! **************************************************************************************************
1400!> \brief Apply the state-selective inverse in a private canonical occupied gauge.
1401!> \param preconditioner_env ...
1402!> \param matrix_in ...
1403!> \param matrix_out ...
1404! **************************************************************************************************
1405 SUBROUTINE apply_full_all_covariant(preconditioner_env, matrix_in, matrix_out)
1406
1407 TYPE(preconditioner_type) :: preconditioner_env
1408 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
1409
1410 CHARACTER(len=*), PARAMETER :: routinen = 'apply_full_all_covariant'
1411
1412 INTEGER :: handle, k, n
1413 TYPE(cp_fm_type) :: canonical_in, canonical_out
1414
1415 CALL timeset(routinen, handle)
1416 cpassert(ASSOCIATED(preconditioner_env%occ_rotation))
1417 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1418
1419 CALL cp_fm_create(canonical_in, matrix_in%matrix_struct, name=routinen)
1420 CALL cp_fm_create(canonical_out, matrix_out%matrix_struct, name=routinen)
1421 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, matrix_in, &
1422 preconditioner_env%occ_rotation, 0.0_dp, canonical_in)
1423 IF (preconditioner_env%lattice_fft_active) THEN
1424 CALL apply_lattice_state_operator_fm(preconditioner_env, canonical_in, canonical_out)
1425 ELSE
1426 CALL apply_full_all(preconditioner_env, canonical_in, canonical_out)
1427 END IF
1428 CALL parallel_gemm('N', 'T', n, k, k, 1.0_dp, canonical_out, &
1429 preconditioner_env%occ_rotation, 0.0_dp, matrix_out)
1430 CALL cp_fm_release(canonical_in)
1431 CALL cp_fm_release(canonical_out)
1432 CALL timestop(handle)
1433
1434 END SUBROUTINE apply_full_all_covariant
1435
1436! **************************************************************************************************
1437!> \brief DBCSR form of the rotation-covariant state-selective inverse.
1438!> \param preconditioner_env ...
1439!> \param matrix_in ...
1440!> \param matrix_out ...
1441! **************************************************************************************************
1442 SUBROUTINE apply_all_covariant(preconditioner_env, matrix_in, matrix_out)
1443
1444 TYPE(preconditioner_type) :: preconditioner_env
1445 TYPE(dbcsr_type) :: matrix_in, matrix_out
1446
1447 CHARACTER(len=*), PARAMETER :: routinen = 'apply_all_covariant'
1448
1449 INTEGER :: handle, k, n
1450 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1451 TYPE(cp_fm_type) :: input_fm, output_fm
1452 TYPE(dbcsr_type) :: canonical_in, canonical_out
1453
1454 CALL timeset(routinen, handle)
1455 cpassert(ASSOCIATED(preconditioner_env%occ_rotation))
1456 CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1457
1458 IF (preconditioner_env%lattice_fft_active) THEN
1459 NULLIFY (fm_struct)
1460 CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1461 context=preconditioner_env%ctxt, &
1462 para_env=preconditioner_env%para_env)
1463 CALL cp_fm_create(input_fm, fm_struct, name=routinen)
1464 CALL cp_fm_create(output_fm, fm_struct, name=routinen)
1465 CALL cp_fm_struct_release(fm_struct)
1466 CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1467 CALL apply_full_all_covariant(preconditioner_env, input_fm, output_fm)
1468 CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1469 CALL cp_fm_release(input_fm)
1470 CALL cp_fm_release(output_fm)
1471 CALL timestop(handle)
1472 RETURN
1473 END IF
1474
1475 IF (.NOT. ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1476 CALL dbcsr_init_p(preconditioner_env%occ_rotation_dbcsr)
1477 CALL cp_dbcsr_m_by_n_from_template(preconditioner_env%occ_rotation_dbcsr, &
1478 template=matrix_in, m=k, n=k, &
1479 sym=dbcsr_type_no_symmetry)
1480 CALL copy_fm_to_dbcsr(preconditioner_env%occ_rotation, &
1481 preconditioner_env%occ_rotation_dbcsr)
1482 END IF
1483
1484 CALL dbcsr_copy(canonical_in, matrix_in, name=routinen)
1485 CALL dbcsr_set(canonical_in, 0.0_dp)
1486 CALL dbcsr_copy(canonical_out, matrix_out, name=routinen)
1487 CALL dbcsr_set(canonical_out, 0.0_dp)
1488 CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_in, &
1489 preconditioner_env%occ_rotation_dbcsr, 0.0_dp, canonical_in)
1490 CALL apply_all(preconditioner_env, canonical_in, canonical_out)
1491 CALL dbcsr_multiply('N', 'T', 1.0_dp, canonical_out, &
1492 preconditioner_env%occ_rotation_dbcsr, 0.0_dp, matrix_out)
1493 CALL dbcsr_release(canonical_in)
1494 CALL dbcsr_release(canonical_out)
1495 CALL timestop(handle)
1496
1497 END SUBROUTINE apply_all_covariant
1498
1499END MODULE preconditioner_apply
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:563
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:845
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:793
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
cp_dbcsr_operations::copy_local_to_dbcsr
subroutine, public copy_local_to_dbcsr(matrix, dense)
Store a single-rank, nonsymmetric matrix without redistribution.
Definition cp_dbcsr_operations.F:201
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:1284
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:240
cp_dbcsr_operations::copy_dbcsr_to_local
subroutine, public copy_dbcsr_to_local(matrix, dense)
Load a single-rank, nonsymmetric matrix without redistribution.
Definition cp_dbcsr_operations.F:172
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1503
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:319
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1046
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ot_precond_solver_chebyshev
integer, parameter, public ot_precond_solver_chebyshev
Definition input_constants.F:565
input_constants::ot_low_rank_base_overlap
integer, parameter, public ot_low_rank_base_overlap
Definition input_constants.F:560
input_constants::ot_precond_full_all_covariant
integer, parameter, public ot_precond_full_all_covariant
Definition input_constants.F:551
input_constants::ot_low_rank_base_lattice_fft
integer, parameter, public ot_low_rank_base_lattice_fft
Definition input_constants.F:560
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:551
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:551
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:565
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:551
input_constants::ot_precond_fermi_low_rank
integer, parameter, public ot_precond_fermi_low_rank
Definition input_constants.F:551
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:551
input_constants::ot_precond_solver_update
integer, parameter, public ot_precond_solver_update
Definition input_constants.F:565
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:551
input_constants::ot_precond_solver_direct
integer, parameter, public ot_precond_solver_direct
Definition input_constants.F:565
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lattice_preconditioner_operator
Block-circulant inverse operators for replicated-cell OT preconditioners.
Definition lattice_preconditioner_operator.F:11
lattice_preconditioner_operator::apply_lattice_inverse_dense
subroutine, public apply_lattice_inverse_dense(inverse_k, lattice_dims, matrix_in, matrix_out, used_fft, imaginary_residual)
Apply inverse Fourier blocks to a dense cell-major real matrix.
Definition lattice_preconditioner_operator.F:506
lattice_preconditioner_operator::apply_lattice_state_inverse_dense
subroutine, public apply_lattice_state_inverse_dense(inverse_k, lattice_dims, matrix_in, matrix_out, used_fft, imaginary_residual, state_indices)
Apply one inverse Fourier operator to every column of a dense real matrix.
Definition lattice_preconditioner_operator.F:566
local_gemm_api
Definition local_gemm_api.F:8
local_gemm_api::local_gemm_pu_gpu
integer, parameter, public local_gemm_pu_gpu
Definition local_gemm_api.F:41
low_rank_preconditioner_model
Small algebraic helpers for the rotationally covariant low-rank OT preconditioner.
Definition low_rank_preconditioner_model.F:11
low_rank_preconditioner_model::low_rank_inverse_weight
pure elemental real(kind=dp) function, public low_rank_inverse_weight(eigenvalue, reference, energy_gap)
Spectral inverse weight relative to a common occupied reference level.
Definition low_rank_preconditioner_model.F:30
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::degree
real(kind=dp), parameter, public degree
Definition mathconstants.F:140
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_apply
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_apply.F:15
preconditioner_apply::apply_preconditioner_dbcsr_complex
subroutine, public apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, matrix_out_re, matrix_out_im)
Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
Definition preconditioner_apply.F:472
preconditioner_apply::apply_preconditioner_fm
subroutine, public apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
applies a previously created preconditioner to a full matrix
Definition preconditioner_apply.F:271
preconditioner_apply::apply_preconditioner_cfm_complex
subroutine, public apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
Apply a complex k-point orbital preconditioner.
Definition preconditioner_apply.F:374
preconditioner_apply::apply_preconditioner_dbcsr
subroutine, public apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
...
Definition preconditioner_apply.F:322
preconditioner_apply::apply_preconditioner_complex_batch
subroutine, public apply_preconditioner_complex_batch(preconditioners, matrix_in, matrix_out, handled, output_scale)
Apply single-rank complex preconditioners in bounded OpenMP batches. Packing and DBCSR writeback stay...
Definition preconditioner_apply.F:93
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner_types::release_lattice_rhs_distribution
subroutine, public release_lattice_rhs_distribution(distribution)
Release a cached lattice right-hand-side communication map.
Definition preconditioner_types.F:308
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
local_gemm_api::local_gemm_ctxt_type
Definition local_gemm_api.F:47
preconditioner_types::lattice_rhs_distribution_type
Definition preconditioner_types.F:48
preconditioner_types::preconditioner_p_type
Definition preconditioner_types.F:110
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:60