Branch data Line data Source code
1 : : !!****m* ABINIT/m_invovl
2 : : !! NAME
3 : : !! m_invovl
4 : : !!
5 : : !! FUNCTION
6 : : !! Provides functions to invert the overlap matrix S. Used by Chebyshev in PAW
7 : : !! See paper by A. Levitt and M. Torrent for details
8 : : !! S = 1 + projs * s_projs * projs'
9 : : !! S^-1 = 1 + projs * inv_s_projs * projs', with
10 : : !! inv_s_projs = - (s_projs^-1 + projs'*projs)^-1
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 2013-2026 ABINIT group (AL)
14 : : !! This file is distributed under the terms of the
15 : : !! GNU General Public License, see ~abinit/COPYING
16 : : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : ! nvtx related macro definition
27 : : #include "nvtx_macros.h"
28 : :
29 : : MODULE m_invovl
30 : :
31 : : use defs_basis
32 : : use m_errors
33 : : use m_xmpi
34 : : use m_xomp
35 : : use m_abicore
36 : : use m_abi_linalg
37 : :
38 : : use defs_abitypes, only : mpi_type
39 : : use m_time, only : timab
40 : : use m_hamiltonian, only : gs_hamiltonian_type
41 : : use m_bandfft_kpt, only : bandfft_kpt_get_ikpt
42 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_axpby
43 : : use m_gemm_nonlop_projectors, only : gemm_nonlop_use_gemm
44 : : use m_nonlop, only : nonlop
45 : : use m_prep_kgb, only : prep_nonlop
46 : :
47 : : #ifdef HAVE_FC_ISO_C_BINDING
48 : : ! FIXME Don't know what's wrong with GCC when OpenMP GPU Offload is enabled here...
49 : : #ifdef FC_GNU
50 : : use, intrinsic :: iso_c_binding, only : c_int32_t, c_int64_t, c_float, c_double, c_size_t, c_loc
51 : : #else
52 : : use, intrinsic :: iso_c_binding, only : c_ptr, c_int32_t, c_int64_t, c_float, c_double, c_size_t, c_loc
53 : : #endif
54 : : #endif
55 : :
56 : : #if defined(HAVE_GPU_MARKERS)
57 : : use m_nvtx_data
58 : : #endif
59 : :
60 : : #ifdef HAVE_GPU
61 : : use m_gputk
62 : : #endif
63 : :
64 : : #ifdef HAVE_KOKKOS
65 : : use m_manage_kokkos, only : add_array_kokkos
66 : : #endif
67 : :
68 : : implicit none
69 : :
70 : : private
71 : :
72 : : !public procedures.
73 : : public :: init_invovl
74 : : public :: make_invovl
75 : : public :: apply_invovl
76 : : public :: destroy_invovl
77 : :
78 : : ! Those routines are here to assess memory requirements
79 : : public :: invovl_ompgpu_work_mem
80 : : public :: invovl_ompgpu_static_mem
81 : : !!***
82 : :
83 : : !!****t* m_invovl/invovl_kpt_type
84 : : !! NAME
85 : : !! invovl_kpt_type
86 : : !!
87 : : !! FUNCTION
88 : : !! Contains information needed to invert the overlap matrix S
89 : : !!
90 : : !! SOURCE
91 : :
92 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
93 : :
94 : : type, public :: invovl_kpt_type
95 : :
96 : : integer(kind=c_int32_t) :: nprojs
97 : : ! total number of projectors
98 : : ! nlmn for a specific atom = count(indlmn(3,:,itypat)>0)
99 : : ! A value of -1 means that the following arrays are not allocated
100 : :
101 : : real(kind=c_double), allocatable :: gram_projs(:,:,:)
102 : : ! gram_projs(cplx, nprojs, nprojs)
103 : : ! projs' * projs
104 : :
105 : : real(kind=c_double), allocatable :: inv_sij(:,:,:,:)
106 : : ! inv_sij(cplx, lmnmax, lmnmax, ntypat)
107 : : ! inverse of ham%sij
108 : :
109 : : real(kind=c_double), allocatable :: inv_s_approx(:,:,:,:)
110 : : ! inv_s_approx(cplx, lmnmax, lmnmax, ntypat)
111 : : ! preconditionner
112 : :
113 : : end type invovl_kpt_type
114 : :
115 : : !> companion type to invovl_kpt_type to pass data to gpu/cuda
116 : : type, bind(c), public :: invovl_kpt_gpu_type
117 : :
118 : : integer(kind=c_int32_t) :: nprojs
119 : : ! total number of projectors
120 : : ! nlmn for a specific atom = count(indlmn(3,:,itypat)>0)
121 : : ! A value of -1 means that the following arrays are not allocated
122 : :
123 : : type(c_ptr) :: gram_projs
124 : : ! gram_projs(cplx, nprojs, nprojs)
125 : : ! projs' * projs
126 : :
127 : : integer(kind=c_int32_t) :: gram_projs_dim(3)
128 : :
129 : : type(c_ptr) :: inv_sij
130 : : ! inv_sij(cplx, lmnmax, lmnmax, ntypat)
131 : : ! inverse of ham%sij
132 : :
133 : : integer(kind=c_int32_t) :: inv_sij_dim(4)
134 : :
135 : : type(c_ptr) :: inv_s_approx
136 : : ! inv_s_approx(cplx, lmnmax, lmnmax, ntypat)
137 : : ! preconditionner
138 : :
139 : : integer(kind=c_int32_t) :: inv_s_approx_dim(4)
140 : :
141 : : end type invovl_kpt_gpu_type
142 : :
143 : : #else
144 : :
145 : : type, public :: invovl_kpt_type
146 : :
147 : : integer :: nprojs
148 : : ! total number of projectors
149 : : ! nlmn for a specific atom = count(indlmn(3,:,itypat)>0)
150 : : ! A value of -1 means that the following arrays are not allocated
151 : :
152 : : real(dp), allocatable :: gram_projs(:,:,:)
153 : : ! gram_projs(cplx, nprojs, nprojs)
154 : : ! projs' * projs
155 : :
156 : : real(dp), allocatable :: inv_sij(:,:,:,:)
157 : : ! inv_sij(cplx, lmnmax, lmnmax, ntypat)
158 : : ! inverse of ham%sij
159 : :
160 : : real(dp), allocatable :: inv_s_approx(:,:,:,:)
161 : : ! inv_s_approx(cplx, lmnmax, lmnmax, ntypat)
162 : : ! preconditionner
163 : :
164 : : end type invovl_kpt_type
165 : :
166 : : #endif
167 : :
168 : : !!***
169 : :
170 : : integer, parameter :: tim_nonlop = 13
171 : :
172 : : integer, parameter :: &
173 : : & timer_apply_inv_ovl_opernla = 1630, &
174 : : & timer_apply_inv_ovl_opernlb = 1631, &
175 : : & timer_apply_inv_ovl_inv_s = 1632
176 : :
177 : : type(invovl_kpt_type), public,save,allocatable, target :: invovl_kpt(:)
178 : : #ifdef HAVE_OPENMP_OFFLOAD
179 : : !Module variable keeping track of which K-point data is so=tored on GPU
180 : : integer, save :: current_ikpt_in_gpu=-1
181 : : #endif
182 : :
183 : : #if defined(HAVE_GPU_CUDA)
184 : :
185 : : !> this interface is only useful when gpu is enabled
186 : : !! these functions are defined in 46_manage_gpu/gpu_apply_invovl_inner.cu
187 : : !! these functions are entry point for calling cuda implemented functions
188 : : interface
189 : :
190 : : !> allocate GPU workspace for apply_invovl
191 : : subroutine f_gpu_apply_invovl_inner_alloc(proj_dim, ntypat, realloc) bind(c, name='gpu_apply_invovl_inner_alloc')
192 : : use, intrinsic :: iso_c_binding
193 : : integer(kind=c_int32_t), intent(in) :: proj_dim(3)
194 : : integer(kind=c_int32_t), value, intent(in) :: ntypat
195 : : integer(kind=c_int32_t), value, intent(in) :: realloc
196 : : end subroutine f_gpu_apply_invovl_inner_alloc
197 : :
198 : : !> deallocate GPU workspace for apply_invovl
199 : : subroutine f_gpu_apply_invovl_inner_dealloc() bind(c, name='gpu_apply_invovl_inner_dealloc')
200 : : end subroutine f_gpu_apply_invovl_inner_dealloc
201 : :
202 : : !> allocate GPU workspace for make_invovl (sij and s_approx)
203 : : subroutine f_gpu_apply_invovl_matrix_alloc(cplx, nprojs, lmnmax, ntypat, realloc) bind(c, name='gpu_apply_invovl_matrix_alloc')
204 : : use, intrinsic :: iso_c_binding
205 : : integer(kind=c_int32_t), value, intent(in) :: cplx
206 : : integer(kind=c_int32_t), value, intent(in) :: nprojs
207 : : integer(kind=c_int32_t), value, intent(in) :: ntypat
208 : : integer(kind=c_int32_t), value, intent(in) :: lmnmax
209 : : integer(kind=c_int32_t), value, intent(in) :: realloc
210 : : end subroutine f_gpu_apply_invovl_matrix_alloc
211 : :
212 : : !> deallocate GPU workspace for make_invovl (sij and s_approx)
213 : : subroutine f_gpu_apply_invovl_matrix_dealloc() bind(c, name='gpu_apply_invovl_matrix_dealloc')
214 : : end subroutine f_gpu_apply_invovl_matrix_dealloc
215 : :
216 : : !> init data for GPU
217 : : subroutine f_gpu_init_invovl_data(indlmn_dim, indlmn_ptr) bind(c, name='init_invovl_data')
218 : : use, intrinsic :: iso_c_binding
219 : : integer(kind=c_int32_t), intent(in) :: indlmn_dim(3)
220 : : type(c_ptr) , value :: indlmn_ptr
221 : : end subroutine f_gpu_init_invovl_data
222 : :
223 : : !> upload inverse overlap matrices
224 : : subroutine f_upload_inverse_overlap(invovl_gpu, cplx, nprojs, lmnmax, ntypat) bind(c, name='upload_inverse_overlap')
225 : : use, intrinsic :: iso_c_binding
226 : : import invovl_kpt_gpu_type
227 : : type(invovl_kpt_gpu_type), value, intent(in) :: invovl_gpu
228 : : integer(kind=c_int32_t), value, intent(in) :: cplx
229 : : integer(kind=c_int32_t), value, intent(in) :: nprojs
230 : : integer(kind=c_int32_t), value, intent(in) :: lmnmax
231 : : integer(kind=c_int32_t), value, intent(in) :: ntypat
232 : : end subroutine f_upload_inverse_overlap
233 : :
234 : : !> solve_inner on GPU
235 : : subroutine f_solve_inner_gpu(proj_dim, proj_ptr, sm1proj_ptr, ptp_sm1proj_ptr, &
236 : : & nattyp_dim, nattyp_ptr, ntypat, lmnmax, cplx, block_sliced) bind(c, name='solve_inner_gpu')
237 : :
238 : : use, intrinsic :: iso_c_binding
239 : : integer(kind=c_int32_t), intent(in) :: proj_dim(3)
240 : : type(c_ptr) , value :: proj_ptr
241 : : type(c_ptr) , value :: sm1proj_ptr
242 : : type(c_ptr) , value :: ptp_sm1proj_ptr
243 : : integer(kind=c_int32_t), value, intent(in) :: nattyp_dim
244 : : type(c_ptr) , value :: nattyp_ptr
245 : : integer(kind=c_int32_t), value, intent(in) :: ntypat
246 : : integer(kind=c_int32_t), value, intent(in) :: lmnmax
247 : : integer(kind=c_int32_t), value, intent(in) :: cplx
248 : : integer(kind=c_int32_t), value, intent(in) :: block_sliced
249 : : end subroutine f_solve_inner_gpu
250 : :
251 : : end interface
252 : :
253 : : #endif
254 : :
255 : : CONTAINS
256 : :
257 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
258 : : !!****f* m_invovl/make_invovl_kpt_gpu
259 : : !! NAME
260 : : !! make_invovl_kpt
261 : : !!
262 : : !! FUNCTION
263 : : !! Create a invovl_pkt_gpu_type from a cpu counter part for cuda interoperability
264 : : !! SOURCE
265 : : function make_invovl_kpt_gpu(invovl) result(invovl_gpu)
266 : : type(invovl_kpt_type), intent(inout),target :: invovl
267 : : type(invovl_kpt_gpu_type) :: invovl_gpu
268 : :
269 : : invovl_gpu%nprojs = invovl%nprojs
270 : :
271 : : invovl_gpu%gram_projs = c_loc(invovl%gram_projs(1,1,1))
272 : : invovl_gpu%gram_projs_dim = (/ &
273 : : & size(invovl%gram_projs,1), &
274 : : & size(invovl%gram_projs,2), &
275 : : & size(invovl%gram_projs,3) &
276 : : & /)
277 : :
278 : : invovl_gpu%inv_sij = c_loc(invovl%inv_sij(1,1,1,1))
279 : : invovl_gpu%inv_sij_dim = (/ &
280 : : & size(invovl%inv_sij,1), &
281 : : & size(invovl%inv_sij,2), &
282 : : & size(invovl%inv_sij,3), &
283 : : & size(invovl%inv_sij,4) &
284 : : & /)
285 : :
286 : : invovl_gpu%inv_s_approx = c_loc(invovl%inv_s_approx(1,1,1,1))
287 : : invovl_gpu%inv_s_approx_dim = (/ &
288 : : & size(invovl%inv_s_approx,1), &
289 : : & size(invovl%inv_s_approx,2), &
290 : : & size(invovl%inv_s_approx,3), &
291 : : & size(invovl%inv_s_approx,4) &
292 : : & /)
293 : :
294 : : end function make_invovl_kpt_gpu
295 : : !!***
296 : : #endif
297 : :
298 : : #ifdef HAVE_OPENMP_OFFLOAD
299 : :
300 : : ! NOTE cumbersome implementation to avoid using derived-type syntax in OMP directives (loose support in compilers)
301 : : subroutine alloc_invovl_arrays_ompgpu(gram_projs, inv_s_approx, inv_sij)
302 : : real(dp), intent(inout) :: gram_projs(:,:,:)
303 : : real(dp), intent(inout) :: inv_sij(:,:,:,:)
304 : : real(dp), intent(inout) :: inv_s_approx(:,:,:,:)
305 : : !$OMP TARGET ENTER DATA MAP(to:gram_projs,inv_sij,inv_s_approx)
306 : : end subroutine alloc_invovl_arrays_ompgpu
307 : :
308 : : subroutine delete_invovl_arrays_ompgpu(gram_projs, inv_s_approx, inv_sij)
309 : : real(dp), intent(inout) :: gram_projs(:,:,:)
310 : : real(dp), intent(inout) :: inv_sij(:,:,:,:)
311 : : real(dp), intent(inout) :: inv_s_approx(:,:,:,:)
312 : : !$OMP TARGET EXIT DATA MAP(delete:gram_projs,inv_sij,inv_s_approx)
313 : : end subroutine delete_invovl_arrays_ompgpu
314 : :
315 : : subroutine refresh_invovl_ompgpu_kpt(ikpt)
316 : : integer,intent(in) :: ikpt
317 : :
318 : : ! *************************************************************************
319 : :
320 : : if(ikpt < 0) then
321 : : ABI_BUG("Requested GPU upload of a negative K-point index !")
322 : : end if
323 : :
324 : : if(current_ikpt_in_gpu /= -1) then
325 : : call delete_invovl_arrays_ompgpu(invovl_kpt(current_ikpt_in_gpu)%gram_projs, &
326 : : & invovl_kpt(current_ikpt_in_gpu)%inv_s_approx, &
327 : : & invovl_kpt(current_ikpt_in_gpu)%inv_sij)
328 : : end if
329 : :
330 : : call alloc_invovl_arrays_ompgpu(invovl_kpt(ikpt)%gram_projs, &
331 : : & invovl_kpt(ikpt)%inv_s_approx, &
332 : : & invovl_kpt(ikpt)%inv_sij)
333 : : current_ikpt_in_gpu=ikpt
334 : :
335 : : end subroutine refresh_invovl_ompgpu_kpt
336 : : #endif
337 : :
338 : : !!****f* m_invovl/init_invovl
339 : : !! NAME
340 : : !! init_invovl
341 : : !!
342 : : !! FUNCTION
343 : : !! Initalization of the invovl_kpt array
344 : : !!
345 : : !! INPUTS
346 : : !! nkpt= number of k-points
347 : : !!
348 : : !! SOURCE
349 : :
350 : 126 : subroutine init_invovl(nkpt)
351 : :
352 : : integer, intent(in) :: nkpt
353 : : integer :: ikpt
354 : :
355 : : ! *************************************************************************
356 : :
357 [ - + + - : 1820 : ABI_MALLOC(invovl_kpt, (nkpt))
+ - + - -
+ - + - +
+ + ]
358 : : ! TODO add cycling if kpt parallelism
359 [ + + ]: 1568 : do ikpt=1,nkpt
360 : 1568 : invovl_kpt(ikpt)%nprojs = -1
361 : : end do
362 : :
363 : 126 : end subroutine init_invovl
364 : : !!***
365 : :
366 : : !!****f* m_invovl/destroy_invovl_ikpt
367 : : !! NAME
368 : : !! destroy_invovl_ikpt
369 : : !!
370 : : !! FUNCTION
371 : : !! Destruction of the ikpt-th member of invovl array
372 : : !!
373 : : !! INPUTS
374 : : !! ikpt= index of k-point
375 : : !!
376 : : !! SOURCE
377 : 16747 : subroutine destroy_invovl_ikpt(ikpt, gpu_option)
378 : :
379 : : integer, intent(in) :: ikpt
380 : : integer, intent(in) :: gpu_option
381 : :
382 : : ! *************************************************************************
383 : :
384 : : if(gpu_option==ABI_GPU_OPENMP) then
385 : : #ifdef HAVE_OPENMP_OFFLOAD
386 : : if(current_ikpt_in_gpu == ikpt) then
387 : : call delete_invovl_arrays_ompgpu(invovl_kpt(ikpt)%gram_projs, &
388 : : & invovl_kpt(ikpt)%inv_s_approx, invovl_kpt(ikpt)%inv_sij)
389 : : current_ikpt_in_gpu = -1
390 : : end if
391 : : #endif
392 : : end if
393 : :
394 [ - + ]: 16747 : ABI_FREE(invovl_kpt(ikpt)%gram_projs)
395 [ - + ]: 16747 : ABI_FREE(invovl_kpt(ikpt)%inv_sij)
396 [ - + ]: 16747 : ABI_FREE(invovl_kpt(ikpt)%inv_s_approx)
397 : 16747 : invovl_kpt(ikpt)%nprojs = -1
398 : :
399 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_FC_ISO_C_BINDING)
400 : : if (gpu_option == ABI_GPU_LEGACY .or. gpu_option == ABI_GPU_KOKKOS) then
401 : : call f_gpu_apply_invovl_inner_dealloc()
402 : : call f_gpu_apply_invovl_matrix_dealloc()
403 : : end if
404 : : #endif
405 : :
406 : 16747 : end subroutine destroy_invovl_ikpt
407 : : !!***
408 : :
409 : : !!****f* m_invovl/destroy_invovl
410 : : !! NAME
411 : : !! destroy_invovl
412 : : !!
413 : : !! FUNCTION
414 : : !! Destruction of the invovl_kpt array
415 : : !!
416 : : !! INPUTS
417 : : !! nkpt= number of k-points
418 : : !!
419 : : !! SOURCE
420 : 125 : subroutine destroy_invovl(nkpt, gpu_option)
421 : :
422 : : integer, intent(in) :: nkpt
423 : : integer, intent(in) :: gpu_option
424 : : integer :: ikpt
425 : :
426 : : ! *************************************************************************
427 : :
428 : : ! TODO add cycling if kpt parallelism
429 [ + + ]: 1539 : do ikpt=1,nkpt
430 [ + + ]: 1414 : if(invovl_kpt(ikpt)%nprojs == -1) then
431 : : ! write(0, *) 'ERROR invovl_kpt is unallocated'
432 : : cycle
433 : : end if
434 : 1539 : call destroy_invovl_ikpt(ikpt, gpu_option)
435 : : end do
436 : :
437 [ + - + + : 1539 : ABI_FREE(invovl_kpt)
- + - + -
+ - + ]
438 : :
439 : 125 : end subroutine destroy_invovl
440 : : !!***
441 : :
442 : : !!****f* m_invovl/make_invovl
443 : : !! NAME
444 : : !! make_invovl
445 : : !!
446 : : !! FUNCTION
447 : : !! Builds of the invovl structure
448 : : !!
449 : : !! INPUTS
450 : : !!
451 : : !! SOURCE
452 : :
453 : 16775 : subroutine make_invovl(ham, dimffnl, ffnl, ph3d, mpi_enreg)
454 : :
455 : : use m_abi_linalg
456 : :
457 : : type(gs_hamiltonian_type),intent(in), target :: ham
458 : : integer, intent(in) :: dimffnl
459 : : real(dp),intent(in) :: ffnl(ham%npw_k,dimffnl,ham%lmnmax,ham%ntypat)
460 : : real(dp),intent(in) :: ph3d(2,ham%npw_k,ham%matblk)
461 : : type(mpi_type) :: mpi_enreg
462 : :
463 : 33550 : real(dp) :: atom_projs(2, ham%npw_k, ham%lmnmax)
464 : 33550 : real(dp) :: temp(ham%npw_k)
465 : 16775 : complex(dp), allocatable :: work(:)
466 : 16775 : real(dp), allocatable,target :: projs(:,:,:)
467 : 16775 : real(dp), allocatable :: gram_proj(:,:,:)
468 : 16775 : integer, allocatable :: ipiv(:)
469 : :
470 : : integer :: itypat, ilmn, nlmn, jlmn, ia, iaph3d, shift
471 : : integer :: il, ilm, jlm, ipw, info, ierr, cplx
472 : : integer :: ikpt_this_proc,cplex_dij
473 : : logical :: parity
474 : : real(dp) :: tsec(2)
475 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
476 : : character(len=500) :: message
477 : : #endif
478 : : character :: blas_transpose
479 : :
480 : : type(invovl_kpt_type), pointer :: invovl
481 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
482 : : type(invovl_kpt_gpu_type) :: invovl_gpu
483 : : #endif
484 : 33550 : integer :: array_nprojs_pp(mpi_enreg%nproc_fft)
485 : : integer :: iproc, slice_size
486 : 16775 : real(dp), allocatable :: gramwork(:,:,:)
487 : :
488 : : integer, parameter :: timer_mkinvovl = 1620, timer_mkinvovl_build_d = 1621, timer_mkinvovl_build_ptp = 1622
489 : :
490 : : ! *************************************************************************
491 : :
492 : : !! S = 1 + PDP', so
493 : : !! S^-1 = 1 + P inv_s_projs P', with
494 : : !! inv_s_projs = - (D^-1 + P'*P)^-1
495 : :
496 : : ABI_NVTX_START_RANGE(NVTX_MAKE_INVOVL)
497 : :
498 [ + + ]: 16775 : if(ham%istwf_k == 1) then
499 : 16639 : cplx = 2
500 : 16639 : blas_transpose = 'c'
501 : : else
502 : 136 : cplx = 1
503 : 136 : blas_transpose = 't'
504 : : end if
505 : :
506 : 16775 : ikpt_this_proc=bandfft_kpt_get_ikpt()
507 : 16775 : invovl => invovl_kpt(ikpt_this_proc)
508 : :
509 [ + + ]: 16775 : if(invovl%nprojs /= -1) then
510 : : ! We have been here before, cleanup before remaking
511 : 15736 : call destroy_invovl_ikpt(ikpt_this_proc, ham%gpu_option)
512 : : end if
513 : :
514 : 16775 : iaph3d = 1
515 : :
516 : 16775 : call timab(timer_mkinvovl,1,tsec)
517 : 16775 : call timab(timer_mkinvovl_build_d,1,tsec)
518 : :
519 : : ! build nprojs
520 : 16775 : invovl%nprojs = 0
521 [ + + ]: 34144 : do itypat=1,ham%ntypat
522 [ + + + + ]: 173466 : invovl%nprojs = invovl%nprojs + count(ham%indlmn(3,:,itypat)>0)*ham%nattyp(itypat)
523 : : end do
524 : :
525 [ - + + - : 67100 : ABI_MALLOC(projs, (2, ham%npw_k, invovl%nprojs))
- + + - +
- + - - +
- + ]
526 [ - + + - : 100650 : ABI_MALLOC(invovl%inv_sij, (cplx, ham%lmnmax, ham%lmnmax, ham%ntypat))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
527 [ - + + - : 100650 : ABI_MALLOC(invovl%inv_s_approx, (cplx, ham%lmnmax, ham%lmnmax, ham%ntypat))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
528 : : ! workspace for inversion
529 [ - + + - : 50325 : ABI_MALLOC(ipiv, (ham%lmnmax))
+ - + - -
+ - + ]
530 [ + - - + : 50325 : ABI_MALLOC(work, (64*ham%lmnmax))
+ - + - +
- - + -
+ ]
531 : :
532 [ + + + + : 133788887 : projs = zero
+ + ]
533 [ + + + + : 3519766 : invovl%inv_sij = zero
+ + + + ]
534 [ + + + + : 3519766 : invovl%inv_s_approx = zero
+ + + + ]
535 : :
536 : : shift = 0
537 [ + + ]: 34144 : do itypat = 1, ham%ntypat
538 [ + + + + ]: 156691 : nlmn = count(ham%indlmn(3,:,itypat)>0)
539 [ + + ]: 17369 : if (size(ham%sij(:,itypat))==ham%lmnmax*(ham%lmnmax+1)/2) then
540 : : cplex_dij = 1
541 [ - + ]: 192 : else if (size(ham%sij(:,itypat))==ham%lmnmax*(ham%lmnmax+1)) then
542 : : cplex_dij = 2
543 : : else
544 : 0 : ABI_ERROR('sij size not recognize')
545 : : end if
546 : : !! unpack ham%sij into inv_sij
547 [ + + ]: 156671 : do jlmn = 1, nlmn
548 [ + + ]: 139302 : if (cplex_dij==1) then
549 : 137766 : invovl%inv_sij(1, jlmn, jlmn, itypat) = ham%sij(jlmn*(jlmn-1)/2 + jlmn, itypat)
550 : : else
551 : 1536 : invovl%inv_sij(1, jlmn, jlmn, itypat) = ham%sij(jlmn*(jlmn-1) + 2*jlmn-1, itypat)
552 : : end if
553 [ + + ]: 647378 : do ilmn = 1, jlmn-1
554 [ + + ]: 630009 : if (cplex_dij==1) then
555 : 485331 : invovl%inv_sij(1, ilmn, jlmn, itypat) = ham%sij(jlmn*(jlmn-1)/2 + ilmn, itypat)
556 : 485331 : invovl%inv_sij(1, jlmn, ilmn, itypat) = ham%sij(jlmn*(jlmn-1)/2 + ilmn, itypat)
557 : : else
558 : 5376 : invovl%inv_sij(1, ilmn, jlmn, itypat) = ham%sij(jlmn*(jlmn-1) + 2*ilmn-1, itypat)
559 : 5376 : invovl%inv_sij(1, jlmn, ilmn, itypat) = ham%sij(jlmn*(jlmn-1) + 2*ilmn-1, itypat)
560 : :
561 : : end if
562 : : end do
563 : : end do
564 : :
565 : : ! Invert sij
566 [ + + ]: 17369 : if(cplx == 2) then
567 : 17097 : call ZHETRF('U', nlmn, invovl%inv_sij(:,:,:,itypat), ham%lmnmax, ipiv, work, (64*nlmn), info)
568 : 17097 : call ZHETRI('U', nlmn, invovl%inv_sij(:,:,:,itypat), ham%lmnmax, ipiv, work, info)
569 : : else
570 : 272 : call DSYTRF('U', nlmn, invovl%inv_sij(:,:,:,itypat), ham%lmnmax, ipiv, work, (64*nlmn), info)
571 : 272 : call DSYTRI('U', nlmn, invovl%inv_sij(:,:,:,itypat), ham%lmnmax, ipiv, work, info)
572 : : end if
573 : : ! complete the matrix
574 [ + + ]: 156671 : do ilm=1, nlmn
575 [ + + ]: 647378 : do jlm=1, ilm-1
576 : 630009 : invovl%inv_sij(1,ilm,jlm,itypat) = invovl%inv_sij(1,jlm,ilm,itypat)
577 : : end do
578 : : end do
579 : :
580 : : !! loop on atoms to build atom_projs and fill projs, s_projs
581 [ + + ]: 68019 : do ia = 1, ham%nattyp(itypat)
582 : :
583 : : !! build atom_projs, from opernlb
584 : : !! P = 4pi/sqrt(ucvol)* conj(diag(ph3d)) * ffnl * diag(parity), with parity = (-i)^l
585 [ + + + + : 133871987 : atom_projs(:,:,:) = zero
+ + ]
586 : :
587 : : ! start from 4pi/sqrt(ucvol)*ffnl
588 : : ! atom_projs(1, :, 1:nlmn) = four_pi/sqrt(ham%ucvol) * ffnl(:, 1, 1:nlmn)
589 : : ! TODO vectorize (DCOPY with stride)
590 [ + + ]: 5109778 : do ipw=1, ham%npw_k
591 [ + + ]: 49608872 : atom_projs(1,ipw, 1:nlmn) = four_pi/sqrt(ham%ucvol) * ffnl(ipw, 1, 1:nlmn, itypat)
592 : : end do
593 : :
594 : : ! multiply by (-i)^l
595 [ + + ]: 308705 : do ilmn=1,nlmn
596 : 274830 : il=mod(ham%indlmn(1,ilmn, itypat),4);
597 : 274830 : parity=(mod(il,2)==0)
598 [ + + ]: 274830 : if (il>1) then
599 : : ! multiply by -1
600 [ + + + + ]: 11679440 : atom_projs(:,:,ilmn) = -atom_projs(:,:,ilmn)
601 : : end if
602 [ + + ]: 308705 : if(.not. parity) then
603 : : ! multiply by -i
604 [ + + ]: 30658668 : temp = atom_projs(2,:,ilmn)
605 [ + + ]: 30658668 : atom_projs(2,:,ilmn) = -atom_projs(1,:,ilmn)
606 [ + + ]: 30658668 : atom_projs(1,:,ilmn) = temp
607 : : end if
608 : : end do
609 : :
610 : : ! multiply by conj(ph3d)
611 [ + + ]: 308705 : do ilmn=1,nlmn
612 [ + + ]: 44773924 : temp = atom_projs(1, :, ilmn)
613 [ + + ]: 44773924 : atom_projs(1, :, ilmn) = atom_projs(1, :, ilmn) * ph3d(1, :, iaph3d) + atom_projs(2, :, ilmn) * ph3d(2, :, iaph3d)
614 [ + + ]: 44807799 : atom_projs(2, :, ilmn) = atom_projs(2, :, ilmn) * ph3d(1, :, iaph3d) - temp * ph3d(2, :, iaph3d)
615 : : end do
616 : :
617 : : ! me_g0 trick
618 [ + + + + ]: 33875 : if(ham%istwf_k == 2 .and. mpi_enreg%me_g0 == 1) then
619 [ + + ]: 270 : atom_projs(1,1,:) = atom_projs(1,1,:) / sqrt2
620 [ + + ]: 270 : atom_projs(2,1,:) = zero
621 : : end if
622 [ + + ]: 33875 : if(ham%istwf_k > 1) then
623 [ + + + + : 224880 : atom_projs(:,:,:) = atom_projs(:,:,:) * sqrt2
+ + ]
624 : : end if
625 : :
626 : :
627 : : !! atom_projs and typat_s_projs are built, copy them to projs and inv_s_projs
628 [ + + + + : 133805987 : projs(:, :, shift+1:shift+nlmn) = atom_projs(:, :, 1:nlmn)
+ + ]
629 : 33875 : shift = shift + nlmn
630 : :
631 : : ! build inv_s_approx = (D^-1+PtP)^-1 restricted to a single atom block
632 : : ! can be optimized (real, build directly from ffnl)
633 [ + + ]: 33875 : if(ia == 1) then
634 : : ! D^-1
635 [ + + + + ]: 1277343 : invovl%inv_s_approx(1, :, :, itypat) = invovl%inv_sij(1, :, :, itypat)
636 : : ! + PtP
637 [ - + + - : 86805 : ABI_MALLOC(gram_proj, (cplx, nlmn, nlmn))
- + + - -
+ + - + -
+ - - + -
+ - + ]
638 : : call abi_xgemm(blas_transpose,'N', nlmn, nlmn, (3-cplx)*ham%npw_k, cone, atom_projs(:,:,1), (3-cplx)*ham%npw_k, &
639 : 17361 : & atom_projs(:,:,1), (3-cplx)*ham%npw_k, czero, gram_proj(:,:,1), nlmn,x_cplx=cplx)
640 : 17361 : call xmpi_sum(gram_proj,mpi_enreg%comm_bandspinorfft,ierr)
641 [ + + + + : 3499803 : invovl%inv_s_approx(:,1:nlmn,1:nlmn,itypat) = invovl%inv_s_approx(:,1:nlmn,1:nlmn,itypat) + gram_proj(:,:,:)
+ + ]
642 [ - + ]: 17361 : ABI_FREE(gram_proj)
643 : : ! ^-1
644 [ + + ]: 17361 : if(cplx == 2) then
645 : 17089 : call ZHETRF('U', nlmn, invovl%inv_s_approx(:,:,:,itypat), ham%lmnmax, ipiv, work, (64*nlmn), info)
646 : 17089 : call ZHETRI('U', nlmn, invovl%inv_s_approx(:,:,:,itypat), ham%lmnmax, ipiv, work, info)
647 : : else
648 : 272 : call DSYTRF('U', nlmn, invovl%inv_s_approx(:,:,:,itypat), ham%lmnmax, ipiv, work, (64*nlmn), info)
649 : 272 : call DSYTRI('U', nlmn, invovl%inv_s_approx(:,:,:,itypat), ham%lmnmax, ipiv, work, info)
650 : : end if
651 : : ! complete lower triangle of matrix
652 [ + + ]: 173960 : do ilm=1, nlmn
653 [ + + ]: 647082 : do jlm=1, ilm-1
654 : 490483 : invovl%inv_s_approx(1, ilm, jlm, itypat) = invovl%inv_s_approx(1, jlm, ilm, itypat)
655 [ + + ]: 629721 : if(cplx == 2) then
656 : 482867 : invovl%inv_s_approx(2, ilm, jlm, itypat) = -invovl%inv_s_approx(2, jlm, ilm, itypat)
657 : : end if
658 : : end do
659 : : end do
660 : : end if
661 : :
662 : 51244 : iaph3d = iaph3d + 1
663 : : end do
664 : : end do
665 [ - + ]: 16775 : ABI_FREE(ipiv)
666 [ - + ]: 16775 : ABI_FREE(work)
667 : :
668 : 16775 : call timab(timer_mkinvovl_build_d, 2, tsec)
669 : 16775 : call timab(timer_mkinvovl_build_ptp, 1, tsec)
670 : :
671 : : ! Compute P'P one column slice at a time (might be too big to fit in one proc)
672 [ + + ]: 16775 : if(mpi_enreg%paral_kgb == 1) then
673 : : ! Split the work evenly the fft processors
674 [ + + ]: 13464 : array_nprojs_pp(:) = invovl%nprojs / mpi_enreg%nproc_fft
675 : : ! not enough work, there's MOD(nprojs,mpi_enreg%nproc_fft) tasks left
676 : : ! assign them to the first ones
677 [ - + ]: 6728 : array_nprojs_pp(1:MOD(invovl%nprojs,mpi_enreg%nproc_fft)) = array_nprojs_pp(1:MOD(invovl%nprojs,mpi_enreg%nproc_fft)) + 1
678 : : else
679 [ + + ]: 20094 : array_nprojs_pp = invovl%nprojs
680 : : end if
681 [ - + + - : 83875 : ABI_MALLOC(invovl%gram_projs, (cplx,invovl%nprojs,array_nprojs_pp(mpi_enreg%me_fft+1)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
682 : 16775 : shift = 0
683 [ - + - - ]: 16775 : if (ham%gpu_option==ABI_GPU_OPENMP .and. mpi_enreg%nproc_fft==1) then
684 : : #ifdef HAVE_OPENMP_OFFLOAD
685 : : ! compute gram_projs in one GEMM, only one FFT proc expected in GPU mode
686 : : slice_size = array_nprojs_pp(1)
687 : : !$OMP TARGET ENTER DATA MAP(alloc:invovl%gram_projs)
688 : : !$OMP TARGET ENTER DATA MAP(to:projs)
689 : :
690 : : call abi_xgemm(blas_transpose,'N', invovl%nprojs, slice_size, (3-cplx)*ham%npw_k, cone, &
691 : : & projs, (3-cplx)*ham%npw_k, &
692 : : & projs, (3-cplx)*ham%npw_k, czero, &
693 : : & invovl%gram_projs, invovl%nprojs, &
694 : : & x_cplx=cplx, gpu_option=ham%gpu_option)
695 : : call xmpi_sum(invovl%gram_projs,mpi_enreg%comm_band,ierr,use_omp_map=.true.)
696 : : !$OMP TARGET EXIT DATA MAP(from:invovl%gram_projs)
697 : : !$OMP TARGET EXIT DATA MAP(delete:projs)
698 : : #endif
699 : : else
700 [ + + ]: 33558 : do iproc = 1, mpi_enreg%nproc_fft
701 : : ! compute local contribution to slice iproc of gram_projs
702 : 16783 : slice_size = array_nprojs_pp(iproc)
703 [ - + + - : 83915 : ABI_MALLOC(gramwork, (cplx,invovl%nprojs,slice_size))
- + + - -
+ + - + -
+ - - + -
+ - + ]
704 : : call abi_xgemm(blas_transpose,'N', invovl%nprojs, slice_size, (3-cplx)*ham%npw_k, cone,&
705 : : & projs(:,:,1), (3-cplx)*ham%npw_k, &
706 : : & projs(:, :, shift+1), (3-cplx)*ham%npw_k, czero, &
707 : : & gramwork(:,:,1), invovl%nprojs,&
708 : 16783 : & x_cplx=cplx)
709 : 16783 : shift = shift + slice_size
710 : : ! reduce on proc i
711 : 16783 : call xmpi_sum_master(gramwork, iproc-1, mpi_enreg%comm_fft, ierr)
712 [ + + ]: 16783 : if(iproc == mpi_enreg%me_fft+1) then
713 [ + - + - : 21993224 : invovl%gram_projs = gramwork
+ - - + -
- - - - -
- - + + +
+ + + ]
714 : : end if
715 [ - + ]: 33558 : ABI_FREE(gramwork)
716 : : end do
717 : 16775 : call xmpi_sum(invovl%gram_projs,mpi_enreg%comm_band,ierr)
718 : : end if
719 : :
720 : 16775 : call timab(timer_mkinvovl_build_ptp, 2, tsec)
721 : 16775 : call timab(timer_mkinvovl,2,tsec)
722 : :
723 [ - + ]: 16775 : ABI_FREE(projs)
724 : :
725 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
726 : :
727 : : ! upload inverse overlap matrices (sij and s_approx) to GPU memory
728 : : if (ham%gpu_option==ABI_GPU_LEGACY .or. ham%gpu_option==ABI_GPU_KOKKOS) then
729 : : ! allocate memory for sij and s_approx on GPU
730 : : write(message,'(a,a,i12,a,a,i6,a,a,i6,a,a,es12.4,a)') &
731 : : & 'Allocate GPU memory for inverse overlap computations (sij and s_approx) : ',&
732 : : & 'nprojs=',invovl%nprojs,ch10,&
733 : : & 'nlmnax=',ham%lmnmax,ch10,&
734 : : & 'ntypat=',ham%ntypat,ch10,&
735 : : & 'gram_projs_gpu_size (GBytes)=',1e-9*cplx*invovl%nprojs*invovl%nprojs*dp,ch10
736 : : call wrtout(std_out,message,'COLL')
737 : : call f_gpu_apply_invovl_matrix_alloc(cplx, invovl%nprojs, ham%lmnmax, ham%ntypat, 0)
738 : :
739 : : invovl_gpu = make_invovl_kpt_gpu(invovl)
740 : : call f_upload_inverse_overlap(invovl_gpu, cplx, invovl%nprojs, ham%lmnmax, ham%ntypat)
741 : : write(message,*) 'Invovl uploaded to GPU memory'
742 : : call wrtout(std_out,message,'COLL')
743 : : end if
744 : :
745 : : #endif
746 : :
747 : : #ifdef HAVE_OPENMP_OFFLOAD
748 : : if (ham%gpu_option==ABI_GPU_OPENMP) then
749 : : call refresh_invovl_ompgpu_kpt(ikpt_this_proc)
750 : : end if
751 : : #endif
752 : :
753 : : ! LB-10/06/24: This message is too verbose on some cases (for example many k-points)
754 : : ! write(message,*) 'Invovl built'
755 : : ! call wrtout(std_out,message,'COLL')
756 : :
757 : : ABI_NVTX_END_RANGE()
758 : :
759 [ - + - + ]: 16775 : end subroutine make_invovl
760 : : !!***
761 : :
762 : : !!****f* m_invovl/apply_invovl
763 : : !! NAME
764 : : !! apply_invovl
765 : : !!
766 : : !! FUNCTION
767 : : !! Applies the inverse of the overlap matrix to cwavef
768 : : !!
769 : : !! INPUTS
770 : : !!
771 : : !! SOURCE
772 : :
773 [ - + ]: 178662 : subroutine apply_invovl(ham, cwavef, sm1cwavef, cwaveprj, npw, ndat, mpi_enreg, nspinor, block_sliced, invovl_precision, algo)
774 : :
775 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU)
776 : : use, intrinsic :: iso_c_binding
777 : : #endif
778 : :
779 : : ! args
780 : : type(gs_hamiltonian_type), intent(in), target :: ham
781 : : integer, intent(in) :: npw, ndat
782 : : integer, intent(in) :: nspinor
783 : : integer, intent(in) :: block_sliced
784 : : integer, intent(in), optional :: algo
785 : : real(dp), intent(in) :: invovl_precision
786 : : real(dp), intent(inout), target :: cwavef(2, npw*nspinor*ndat) ! TODO should be in, fix nonlop
787 : : type(mpi_type) :: mpi_enreg
788 : : real(dp), intent(inout), target :: sm1cwavef(2, npw*nspinor*ndat)
789 : : type(pawcprj_type), intent(inout) :: cwaveprj(:,:)
790 : : #ifdef HAVE_OPENMP_OFFLOAD
791 : : logical :: transfer_omp_args
792 : : #endif
793 : :
794 : 178662 : real(dp),allocatable, target :: proj(:,:,:), sm1proj(:,:,:), PtPsm1proj(:,:,:)
795 : :
796 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_FC_ISO_C_BINDING)
797 : : ! used to pass proj dimensions to cuda
798 : : integer(kind=c_int32_t) :: proj_dim(3)
799 : : integer(kind=c_int32_t) :: nattyp_dim
800 : : integer(kind=c_int32_t) :: indlmn_dim(3)
801 : : #endif
802 : :
803 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS) && defined(HAVE_YAKL)
804 : : integer(kind=c_int32_t) :: cwavef_size
805 : : #endif
806 : :
807 : : integer :: idat, iatom, nprojs, nlmn, shift, l_algo
808 : : #ifdef HAVE_OPENMP_OFFLOAD
809 : : integer :: icplx, iproj
810 : : #endif
811 : : real(dp) :: tsec(2)
812 : :
813 : : integer :: choice, cpopt, paw_opt , cplx, old_me_g0
814 : 178662 : type(pawcprj_type),allocatable :: cwaveprj_in(:,:)
815 : :
816 : : integer :: ikpt_this_proc
817 : : ! dummies
818 : 357324 : real(dp) :: enlout(ndat), lambda_block(1), gvnlxc(1,1)
819 : : integer, parameter :: nnlout = 0, idir = 0, signs = 2
820 : :
821 : : type(invovl_kpt_type), pointer :: invovl
822 : :
823 : : ! *************************************************************************
824 : :
825 : :
826 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_PREP)
827 [ + + ]: 178662 : l_algo=0; if(present(algo)) l_algo=algo
828 : 178662 : ikpt_this_proc=bandfft_kpt_get_ikpt()
829 : 178662 : invovl => invovl_kpt(ikpt_this_proc)
830 : 178662 : nprojs=invovl%nprojs
831 : :
832 [ + + ]: 178662 : if(ham%istwf_k == 1) then
833 : 160803 : cplx = 2
834 : : else
835 : 17859 : cplx = 1
836 : : end if
837 : :
838 [ - + + - : 893310 : ABI_MALLOC(proj, (cplx,invovl%nprojs,nspinor*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
839 [ - + + - : 714648 : ABI_MALLOC(sm1proj, (cplx,invovl%nprojs,nspinor*ndat))
- + + - -
+ + - - +
- + ]
840 [ - + + - : 714648 : ABI_MALLOC(PtPsm1proj, (cplx,invovl%nprojs,nspinor*ndat))
- + + - -
+ + - - +
- + ]
841 [ + - ]: 178662 : if(ham%gpu_option/=ABI_GPU_OPENMP) then
842 [ + + + + : 73471254 : proj = zero
+ + ]
843 [ + + + + : 73471254 : sm1proj = zero
+ + ]
844 [ + + + + : 73471254 : PtPsm1proj = zero
+ + ]
845 : : end if
846 : :
847 : : #ifdef HAVE_OPENMP_OFFLOAD
848 : : if(ham%gpu_option == ABI_GPU_OPENMP) then
849 : : if(ikpt_this_proc /= current_ikpt_in_gpu) call refresh_invovl_ompgpu_kpt(ikpt_this_proc)
850 : : !$OMP TARGET ENTER DATA MAP(alloc:proj,sm1proj,PtPsm1proj)
851 : :
852 : : transfer_omp_args = .not. ( xomp_target_is_present(c_loc(sm1cwavef)) &
853 : : .and. xomp_target_is_present(c_loc(cwavef)))
854 : : !$OMP TARGET ENTER DATA MAP(alloc:gvnlxc)
855 : : if(transfer_omp_args) then
856 : : !$OMP TARGET ENTER DATA MAP(alloc:sm1cwavef,cwavef)
857 : : !$OMP TARGET UPDATE TO(sm1cwavef,cwavef)
858 : : end if
859 : : end if
860 : : #endif
861 : :
862 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_FC_ISO_C_BINDING)
863 : :
864 : : proj_dim = (/ size(proj,1), size(proj,2), size(proj,3) /)
865 : : nattyp_dim = size(ham%nattyp)
866 : : indlmn_dim = (/ size(ham%indlmn,1), size(ham%indlmn,2), size(ham%indlmn,3) /)
867 : :
868 : : !! memory allocation of data used in solve_inner_gpu
869 : : !! note : this is actually done only once
870 : : if (ham%gpu_option==ABI_GPU_LEGACY .or. ham%gpu_option==ABI_GPU_KOKKOS) then
871 : :
872 : : #ifdef DEBUG_VERBOSE_GPU
873 : : if(xmpi_comm_rank(xmpi_world) == 0) then
874 : : call check_gpu_mem("gpu_apply_invovl_inner_alloc begin")
875 : : end if
876 : : #endif
877 : :
878 : : ! make sure to use sizes from apply_invovl
879 : : call f_gpu_apply_invovl_inner_alloc(proj_dim, ham%ntypat, 0)
880 : :
881 : : #ifdef DEBUG_VERBOSE_GPU
882 : : if(xmpi_comm_rank(xmpi_world) == 0) then
883 : : call check_gpu_mem("gpu_apply_invovl_inner_alloc end")
884 : : end if
885 : : #endif
886 : :
887 : : ! TODO find a better place to put that initialization
888 : : call f_gpu_init_invovl_data(indlmn_dim, c_loc(ham%indlmn(1,1,1)))
889 : :
890 : : end if
891 : :
892 : : #endif
893 : :
894 : 178662 : call timab(timer_apply_inv_ovl_opernla, 1, tsec)
895 : :
896 : : ! cwaveprj may be dummy or unused if gemm nonlop is turned on
897 [ + + + + : 535986 : if((.not. (gemm_nonlop_use_gemm .or. ham%gpu_option==ABI_GPU_OPENMP)) .or. size(cwaveprj) > 1) then
- + + + ]
898 [ - + + - : 4778856 : ABI_MALLOC(cwaveprj_in, (ham%natom,nspinor*ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
899 : 178374 : call pawcprj_alloc(cwaveprj_in,0,ham%dimcprj)
900 : : else
901 [ - + - + : 864 : ABI_MALLOC(cwaveprj_in, (1,1))
+ + + + ]
902 : 288 : call pawcprj_alloc(cwaveprj_in,0,(/1/))
903 : : end if
904 : : ABI_NVTX_END_RANGE()
905 : :
906 : : ! get the cprj
907 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_NONLOP1)
908 : 178662 : choice = 0 ! only compute cprj, nothing else
909 : 178662 : cpopt = 0 ! compute and save cprj
910 : 178662 : paw_opt = 3 ! S nonlocal operator
911 [ + + + - ]: 178662 : if (mpi_enreg%paral_kgb==1 .and. ham%gpu_option/=ABI_GPU_OPENMP) then
912 : : call prep_nonlop(choice,cpopt,cwaveprj_in,enlout,ham,idir,lambda_block,ndat,mpi_enreg,&
913 : : & nnlout,paw_opt,signs,sm1cwavef,tim_nonlop,cwavef,gvnlxc,&
914 : : & already_transposed=.true.,&
915 : : & gpu_option=ham%gpu_option,&
916 : 62460 : & vectproj=proj)
917 : : else
918 [ + + - + ]: 116202 : if(ham%istwf_k==2 .and. ham%gpu_option==ABI_GPU_OPENMP) then
919 : 0 : old_me_g0=mpi_enreg%me_g0
920 [ # # ]: 0 : if (mpi_enreg%me_fft==0) then
921 : 0 : mpi_enreg%me_g0=1
922 : : else
923 : 0 : mpi_enreg%me_g0=0
924 : : end if
925 : : end if
926 : : call nonlop(choice,cpopt,cwaveprj_in,enlout,ham,idir,lambda_block,mpi_enreg,ndat,&
927 : 116202 : & nnlout,paw_opt,signs,sm1cwavef,tim_nonlop,cwavef,gvnlxc,vectproj=proj)
928 : : end if
929 : : ABI_NVTX_END_RANGE()
930 : :
931 : 178662 : call timab(timer_apply_inv_ovl_opernla, 2, tsec)
932 : 178662 : call timab(timer_apply_inv_ovl_inv_s, 1, tsec)
933 : :
934 : : ! If using GEMM nonlop, proj array is used directly instead of writing in cwaveprj_in content, so skip this copy
935 [ + + + - ]: 178662 : if(.not. (gemm_nonlop_use_gemm .or. ham%gpu_option==ABI_GPU_OPENMP)) then
936 : : ! copy cwaveprj_in to proj(:,:)
937 [ + + ]: 1475318 : do idat=1, ndat*nspinor
938 : 1297232 : shift = 0
939 [ + + ]: 4238262 : do iatom = 1, ham%natom
940 : 2762944 : nlmn = cwaveprj_in(iatom, idat)%nlmn
941 [ + + + + ]: 74631296 : proj(1:cplx, shift+1:shift+nlmn, idat) = cwaveprj_in(iatom, idat)%cp(1:cplx, 1:nlmn)
942 : 4060176 : shift = shift + nlmn
943 : : end do
944 : : end do
945 : : end if
946 : :
947 : : !multiply by S^1
948 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_INNER)
949 : : ! TODO : when solve_inner_gpu is ready, update the following to activate GPU computation
950 [ + - ]: 178662 : if (ham%gpu_option == ABI_GPU_LEGACY .or. ham%gpu_option==ABI_GPU_KOKKOS) then
951 : :
952 : : #if defined(HAVE_FC_ISO_C_BINDING) && defined(HAVE_GPU_CUDA)
953 : :
954 : : if (mpi_enreg%nproc_fft /= 1) then
955 : : ABI_ERROR("[66_wfs/m_invovl.F90:apply_invovl] nproc_fft must be 1, when GPU/CUDA is activated")
956 : : end if
957 : :
958 : : call f_solve_inner_gpu(proj_dim, c_loc(proj(1,1,1)), &
959 : : & c_loc(sm1proj(1,1,1)), c_loc(PtPsm1proj(1,1,1)), &
960 : : & nattyp_dim, c_loc(ham%nattyp(1)), ham%ntypat, &
961 : : & ham%lmnmax, cplx, block_sliced)
962 : :
963 : : #endif
964 : :
965 : : else
966 : :
967 [ + + ]: 178662 : if(l_algo == 0) then
968 : : call solve_inner(invovl, ham, cplx, mpi_enreg, proj, ndat*nspinor, sm1proj, PtPsm1proj, &
969 : 171686 : & block_sliced, ham%gpu_option, invovl_precision)
970 [ + - ]: 6976 : else if(l_algo == 1) then
971 : : call solve_inner_sqmr(invovl, ham, cplx, mpi_enreg, proj, ndat*nspinor, sm1proj, PtPsm1proj, &
972 : 6976 : & block_sliced, ham%gpu_option, invovl_precision)
973 : : else
974 : 0 : ABI_BUG("Incorrect algo value")
975 : : end if
976 : :
977 [ + - ]: 178662 : if(ham%gpu_option == ABI_GPU_DISABLED) then
978 [ + + + + : 73471254 : sm1proj = - sm1proj
+ + ]
979 [ + + + + : 73471254 : PtPsm1proj = - PtPsm1proj
+ + ]
980 : : else if(ham%gpu_option == ABI_GPU_OPENMP) then
981 : : #ifdef HAVE_OPENMP_OFFLOAD
982 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:sm1proj,PtPsm1proj)
983 : : do idat =1, ndat*nspinor
984 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(iproj,icplx)
985 : : do iproj = 1, nprojs
986 : : do icplx = 1, cplx
987 : : sm1proj(icplx,iproj,idat) = - sm1proj(icplx,iproj,idat)
988 : : PtPsm1proj(icplx,iproj,idat) = - PtPsm1proj(icplx,iproj,idat)
989 : : end do
990 : : end do
991 : : end do
992 : : #endif
993 : : end if
994 : : end if
995 : :
996 : : ABI_NVTX_END_RANGE()
997 : :
998 : : ! If using GEMM nonlop, sm1proj array is used directly instead of reading cwaveprj content, so skip this copy
999 [ + + + - ]: 178662 : if(.not. (gemm_nonlop_use_gemm .or. ham%gpu_option==ABI_GPU_OPENMP)) then
1000 : : ! copy sm1proj to cwaveprj(:,:)
1001 [ + + ]: 1475318 : do idat=1, ndat*nspinor
1002 : 1297232 : shift = 0
1003 [ + + ]: 4238262 : do iatom = 1, ham%natom
1004 : 2762944 : nlmn = cwaveprj(iatom, idat)%nlmn
1005 [ + + + + ]: 74631296 : cwaveprj(iatom, idat)%cp(1:cplx, 1:nlmn) = sm1proj(1:cplx, shift+1:shift+nlmn, idat)
1006 : 4060176 : shift = shift + nlmn
1007 : : end do
1008 : : end do
1009 : : end if
1010 : 178662 : call timab(timer_apply_inv_ovl_inv_s, 2, tsec)
1011 : 178662 : call timab(timer_apply_inv_ovl_opernlb, 1, tsec)
1012 : :
1013 : : ! get the corresponding wf
1014 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_NONLOP2)
1015 : 178662 : cpopt = 2 ! reuse cprj
1016 : 178662 : choice = 7 ! get wf from cprj, without the application of S
1017 : : paw_opt = 3
1018 [ + + + - ]: 178662 : if (mpi_enreg%paral_kgb==1 .and. ham%gpu_option/=ABI_GPU_OPENMP) then
1019 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,ham,idir,lambda_block,ndat,mpi_enreg,nnlout,&
1020 : : & paw_opt,signs,sm1cwavef,tim_nonlop,cwavef,gvnlxc,already_transposed=.true.,&
1021 : 62460 : & gpu_option=ham%gpu_option,vectproj=sm1proj)
1022 : : else
1023 : : call nonlop(choice,cpopt,cwaveprj,enlout,ham,idir,lambda_block,mpi_enreg,ndat,nnlout,paw_opt,&
1024 : 116202 : & signs,sm1cwavef,tim_nonlop,cwavef,gvnlxc,vectproj=sm1proj)
1025 : : end if
1026 : : ABI_NVTX_END_RANGE()
1027 : :
1028 : 178662 : call timab(timer_apply_inv_ovl_opernlb, 2, tsec)
1029 [ + + - + ]: 178662 : if (ham%istwf_k==2 .and. ham%gpu_option==ABI_GPU_OPENMP) mpi_enreg%me_g0=old_me_g0
1030 : :
1031 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_POST2)
1032 [ + + + + ]: 535986 : if(size(cwaveprj) > 1) then
1033 : : #ifdef HAVE_OPENMP_OFFLOAD
1034 : : !$OMP TARGET UPDATE FROM(PtPsm1proj,proj) IF(ham%gpu_option==ABI_GPU_OPENMP)
1035 : : #endif
1036 : : ! copy PtPsm1proj to cwaveprj(:,:)
1037 [ + + ]: 1477334 : do idat=1, ndat*nspinor
1038 : 1298960 : shift = 0
1039 [ + + ]: 4243734 : do iatom = 1, ham%natom
1040 : 2766400 : nlmn = cwaveprj(iatom, idat)%nlmn
1041 [ + + + + ]: 74717696 : cwaveprj(iatom, idat)%cp(1:cplx, 1:nlmn) = PtPsm1proj(1:cplx, shift+1:shift+nlmn, idat)
1042 : 4065360 : shift = shift + nlmn
1043 : : end do
1044 : : end do
1045 : : !cwaveprj_in is empty if GEMM nonlop is being used, so populate it here
1046 [ + + - + ]: 178374 : if(gemm_nonlop_use_gemm .or. ham%gpu_option == ABI_GPU_OPENMP) then
1047 [ + + ]: 2016 : do idat=1, ndat*nspinor
1048 : 1728 : shift = 0
1049 [ + + ]: 183558 : do iatom = 1, ham%natom
1050 : 3456 : nlmn = cwaveprj_in(iatom, idat)%nlmn
1051 [ + + + + ]: 86400 : cwaveprj_in(iatom, idat)%cp(1:cplx, 1:nlmn) = proj(1:cplx, shift+1:shift+nlmn, idat)
1052 : 5184 : shift = shift + nlmn
1053 : : end do
1054 : : end do
1055 : : end if
1056 : 178374 : call pawcprj_axpby(one, one, cwaveprj_in, cwaveprj)
1057 : : end if
1058 : 178662 : call pawcprj_free(cwaveprj_in)
1059 [ + - + + : 2945350 : ABI_FREE(cwaveprj_in)
- + - + -
+ ]
1060 : :
1061 [ + - ]: 178662 : if (ham%gpu_option == ABI_GPU_LEGACY .or. ham%gpu_option==ABI_GPU_KOKKOS) then
1062 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS) && defined(HAVE_YAKL)
1063 : : cwavef_size = 2*npw*nspinor*ndat
1064 : : call add_array_kokkos(c_loc(sm1cwavef), c_loc(cwavef), cwavef_size)
1065 : : #endif
1066 [ - + ]: 178662 : else if(ham%gpu_option == ABI_GPU_OPENMP) then
1067 : 0 : call abi_xaxpy(2*npw*nspinor*ndat, cone, cwavef, 1, sm1cwavef, 1, gpu_option=ABI_GPU_OPENMP)
1068 : : else
1069 [ + + + + ]: 735307302 : sm1cwavef = cwavef + sm1cwavef
1070 : : end if
1071 : :
1072 : : ABI_NVTX_END_RANGE()
1073 : :
1074 : : if(ham%gpu_option==ABI_GPU_OPENMP) then
1075 : : #ifdef HAVE_OPENMP_OFFLOAD
1076 : : if(transfer_omp_args) then
1077 : : !$OMP TARGET UPDATE FROM(sm1cwavef,cwavef)
1078 : : !$OMP TARGET EXIT DATA MAP(delete:sm1cwavef,cwavef)
1079 : : end if
1080 : : !$OMP TARGET EXIT DATA MAP(delete:gvnlxc)
1081 : : !$OMP TARGET EXIT DATA MAP(delete:proj,sm1proj,PtPsm1proj)
1082 : : #endif
1083 : : end if
1084 : :
1085 [ - + ]: 178662 : ABI_FREE(proj)
1086 [ - + ]: 178662 : ABI_FREE(sm1proj)
1087 [ - + ]: 178662 : ABI_FREE(PtPsm1proj)
1088 : :
1089 [ - + - - : 178662 : end subroutine apply_invovl
- - - - ]
1090 : : !!***
1091 : :
1092 : : !!****f* m_invovl/solve_inner
1093 : : !! NAME
1094 : : !! solve_inner
1095 : : !!
1096 : : !! FUNCTION
1097 : : !! Helper function: iteratively solves the inner system
1098 : : !!
1099 : : !! INPUTS
1100 : : !!
1101 : : !! SOURCE
1102 : 171686 : subroutine solve_inner(invovl, ham, cplx, mpi_enreg, proj, ndat, sm1proj, PtPsm1proj, block_sliced, gpu_option, invovl_precision)
1103 : :
1104 : : use m_abi_linalg
1105 : :
1106 : : integer,intent(in) :: ndat,cplx
1107 : : type(invovl_kpt_type), intent(in), target :: invovl
1108 : : real(dp), intent(in) :: invovl_precision
1109 : : real(dp), intent(inout) :: proj(cplx, invovl%nprojs,ndat)
1110 : : real(dp), intent(inout), target :: sm1proj(cplx, invovl%nprojs, ndat)
1111 : : real(dp), intent(inout), target :: PtPsm1proj(cplx, invovl%nprojs, ndat)
1112 : 171686 : real(dp), allocatable :: temp_proj(:,:,:)
1113 : : type(mpi_type), intent(in) :: mpi_enreg
1114 : : type(gs_hamiltonian_type),intent(in) :: ham
1115 : : integer, intent(in) :: block_sliced, gpu_option
1116 : :
1117 : 343372 : integer :: array_nlmntot_pp(mpi_enreg%nproc_fft)
1118 : : integer :: nlmntot_this_proc, ibeg, iend, ierr, i, nprojs
1119 : : integer :: nstagnated
1120 : 343372 : real(dp) :: resid(cplx, invovl%nprojs,ndat), precondresid(cplx, invovl%nprojs,ndat)
1121 : 343372 : real(dp) :: normprojs(ndat), errs(ndat), maxerr, previous_maxerr
1122 : : ! Stagnation check: if the (global, over ndat) residual fails to shrink by
1123 : : ! at least a factor stagnation_tol over max_stagnated_it consecutive
1124 : : ! iterations, further iterations are extremely unlikely to help (this
1125 : : ! typically happens once invovl_precision is set close to, or below, the
1126 : : ! relative residual level attainable in double precision for the given
1127 : : ! conditioning of A) : give up early instead of iterating all the way to
1128 : : ! maxit.
1129 : : real(dp), parameter :: stagnation_tol = 0.99_dp
1130 : : integer, parameter :: max_stagnated_it = 3
1131 : : character(len=500) :: message
1132 : :
1133 : : real(dp) :: convergence_rate
1134 : : integer :: additional_steps_to_take
1135 : : #ifdef HAVE_OPENMP_OFFLOAD
1136 : : real(dp) :: sum_tmp
1137 : : integer :: idat,iproj,icplx
1138 : : #endif
1139 : :
1140 : : ! *************************************************************************
1141 : :
1142 : 171686 : nprojs = invovl%nprojs
1143 [ + - ]: 171686 : if(gpu_option == ABI_GPU_DISABLED) then
1144 [ + + + + : 69497302 : normprojs = SUM(SUM(proj**2, 1),1)
+ + ]
1145 : 171686 : call xmpi_sum(normprojs, mpi_enreg%comm_fft, ierr)
1146 : : else if(gpu_option == ABI_GPU_OPENMP) then
1147 : : #ifdef HAVE_OPENMP_OFFLOAD
1148 : : !$OMP TARGET ENTER DATA MAP(alloc:errs,precondresid,resid,normprojs)
1149 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:normprojs,proj) PRIVATE(idat,sum_tmp)
1150 : : do idat = 1,ndat
1151 : : sum_tmp=0
1152 : : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:sum_tmp) PRIVATE(iproj,icplx)
1153 : : do iproj = 1,nprojs
1154 : : do icplx = 1,cplx
1155 : : sum_tmp = sum_tmp + proj(icplx,iproj,idat)**2
1156 : : end do
1157 : : end do
1158 : : normprojs(idat)=sum_tmp
1159 : : end do
1160 : : !$OMP TARGET UPDATE FROM(normprojs)
1161 : : #endif
1162 : : end if
1163 : :
1164 : : ! Compute work distribution : split nprojs evenly between the fft processors
1165 : : ! On GPU, work isn't shared on npfft so path is the same wether paral_kgb={0,1}
1166 [ + + + - ]: 171686 : if(mpi_enreg%paral_kgb == 1 .and. gpu_option == ABI_GPU_DISABLED) then
1167 [ + + ]: 125208 : array_nlmntot_pp(:) = nprojs / mpi_enreg%nproc_fft
1168 : : ! not enough work, there's MOD(nprojs,mpi_enreg%nproc_fft) tasks left
1169 : : ! assign them to the first ones
1170 [ - + ]: 62460 : array_nlmntot_pp(1:MOD(nprojs,mpi_enreg%nproc_fft)) = array_nlmntot_pp(1:MOD(nprojs,mpi_enreg%nproc_fft)) + 1
1171 [ + + ]: 62604 : ibeg = SUM(array_nlmntot_pp(1:mpi_enreg%me_fft)) + 1
1172 : 62460 : iend = ibeg + array_nlmntot_pp(1+mpi_enreg%me_fft) - 1
1173 : 62460 : nlmntot_this_proc = iend - ibeg + 1
1174 : : else
1175 : 109226 : ibeg = 1
1176 : 109226 : iend = nprojs
1177 : 109226 : nlmntot_this_proc = nprojs
1178 : : end if
1179 : :
1180 [ + - ]: 171686 : if(gpu_option == ABI_GPU_DISABLED) then
1181 [ - + + - : 858430 : ABI_MALLOC(temp_proj, (cplx, nlmntot_this_proc, ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
1182 : : end if
1183 : :
1184 : : ! first guess for sm1proj
1185 : 171686 : call apply_block(ham, cplx, invovl%inv_s_approx, nprojs, ndat, proj, sm1proj, block_sliced, gpu_option)
1186 : :
1187 : 171686 : previous_maxerr = 1.0 ! Avoid uninit value
1188 : :
1189 : : ! Iterative refinement
1190 : : ! TODO use a more efficient iterative algorithm than iterative refinement, use locking
1191 : 171686 : additional_steps_to_take = -1
1192 [ + - ]: 1550561 : do i=1, 30
1193 : : #ifdef FC_NVHPC
1194 : : ! Silly fix for NVHPC 25.1
1195 : : if(ndat == -42) write(100,*) ndat
1196 : : #endif
1197 : : ! compute resid = proj - (D^-1 + PtP)sm1proj
1198 : 1550561 : call apply_block(ham, cplx, invovl%inv_sij, nprojs, ndat, sm1proj, resid, block_sliced, gpu_option)
1199 : :
1200 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_INNER_GEMM)
1201 [ + - ]: 1550561 : if(gpu_option == ABI_GPU_DISABLED) then
1202 [ + - + - : 753397000 : temp_proj = sm1proj(:,ibeg:iend,:)
+ - - + -
- + + + +
+ + ]
1203 : :
1204 : : ! compute matrix multiplication : PtPsm1proj(:,:,1) = invovl%gram * temp_proj(:,:,1)
1205 : : call abi_xgemm('N', 'N', nprojs, ndat, nlmntot_this_proc, cone, &
1206 : : & invovl%gram_projs(:,:,1), nprojs, &
1207 : : & temp_proj(:,:,1), nlmntot_this_proc, czero, &
1208 : : & PtPsm1proj(:,:,1), nprojs, &
1209 : 1550561 : & x_cplx=cplx)
1210 : 1550561 : call xmpi_sum(PtPsm1proj, mpi_enreg%comm_fft, ierr)
1211 [ + + + + : 752330567 : resid = proj - resid - PtPsm1proj
+ + ]
1212 [ # # ]: 0 : else if(gpu_option == ABI_GPU_OPENMP) then
1213 : : ! compute matrix multiplication : PtPsm1proj(:,:,1) = invovl%gram * sm1proj(:,:,1)
1214 : : call abi_xgemm('N', 'N', nprojs, ndat, nlmntot_this_proc, cone, &
1215 : : & invovl%gram_projs, nprojs,&
1216 : : & sm1proj, nlmntot_this_proc, czero, &
1217 : : & PtPsm1proj, nprojs, &
1218 : 0 : & x_cplx=cplx, gpu_option=ABI_GPU_OPENMP)
1219 : :
1220 : : #ifdef HAVE_OPENMP_OFFLOAD
1221 : : !$OMP TARGET TEAMS DISTRIBUTE &
1222 : : !$OMP& PRIVATE(idat) MAP(to:proj,resid,PtPsm1proj)
1223 : : do idat =1, ndat
1224 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(iproj,icplx)
1225 : : do iproj =1, nprojs
1226 : : do icplx = 1,cplx
1227 : : resid(icplx, iproj, idat) = proj(icplx, iproj, idat) - resid(icplx, iproj, idat) - PtPsm1proj(icplx, iproj, idat)
1228 : : end do
1229 : : end do
1230 : : end do
1231 : : #endif
1232 : : end if
1233 : :
1234 : : ! exit check
1235 [ + - ]: 1550561 : if(gpu_option == ABI_GPU_DISABLED) then
1236 [ + + + + : 752330567 : errs = SUM(SUM(resid**2, 1),1)
+ + ]
1237 : 1550561 : call xmpi_sum(errs, mpi_enreg%comm_fft, ierr)
1238 : : else if (gpu_option == ABI_GPU_OPENMP) then
1239 : : #ifdef HAVE_OPENMP_OFFLOAD
1240 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:errs,resid) PRIVATE(idat,sum_tmp)
1241 : : do idat = 1,ndat
1242 : : sum_tmp=0
1243 : : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:sum_tmp) PRIVATE(iproj,icplx)
1244 : : do iproj = 1,nprojs
1245 : : do icplx = 1,cplx
1246 : : sum_tmp = sum_tmp + resid(icplx,iproj,idat)**2
1247 : : end do
1248 : : end do
1249 : : errs(idat)=sum_tmp
1250 : : end do
1251 : : !$OMP TARGET UPDATE FROM(errs)
1252 : : #endif
1253 : : end if
1254 : : ABI_NVTX_END_RANGE()
1255 : :
1256 [ + - - + : 14298218 : maxerr = sqrt(MAXVAL(errs/normprojs))
- - + + ]
1257 [ + + + + ]: 1550561 : if(maxerr < invovl_precision .or. additional_steps_to_take == 1) then
1258 : : exit
1259 : : ! We might stall and never get to the specified precision because of machine errors.
1260 : : ! If we got to 1e-10, extrapolate convergence rate and determine the number of additional
1261 : : ! steps to take to reach precision
1262 [ + + + + ]: 1462593 : else if(maxerr < 1e-10 .and. additional_steps_to_take == -1) then
1263 : 165918 : convergence_rate = -LOG(1e-10) / i
1264 [ - + ]: 165918 : additional_steps_to_take = CEILING(-LOG(invovl_precision/1e-10)/convergence_rate) + 1
1265 [ + + ]: 1296675 : else if(additional_steps_to_take > 0) then
1266 [ + + ]: 541092 : if(previous_maxerr<maxerr)exit
1267 : 457374 : additional_steps_to_take = additional_steps_to_take - 1
1268 : : end if
1269 : :
1270 : : ! Stagnation check (see comment on stagnation_tol/max_stagnated_it above).
1271 [ + + ]: 1378875 : if(maxerr > stagnation_tol*previous_maxerr) then
1272 : 2437 : nstagnated = nstagnated + 1
1273 [ - + ]: 2437 : if(nstagnated >= max_stagnated_it) exit
1274 : : else
1275 : : nstagnated = 0
1276 : : end if
1277 : :
1278 : 1378875 : previous_maxerr=maxerr
1279 : :
1280 : : ! add preconditionned residual
1281 : 1378875 : call apply_block(ham, cplx, invovl%inv_s_approx, nprojs, ndat, resid, precondresid, block_sliced, gpu_option)
1282 [ + - ]: 1550561 : if(gpu_option == ABI_GPU_DISABLED) then
1283 [ + + + + : 682833265 : sm1proj = sm1proj + precondresid
+ + ]
1284 [ # # ]: 0 : else if(gpu_option == ABI_GPU_OPENMP) then
1285 : 0 : call abi_xaxpy(cplx*nprojs*ndat, cone, precondresid, 1, sm1proj, 1, gpu_option=ABI_GPU_OPENMP)
1286 : : end if
1287 : : end do
1288 : :
1289 [ + + - + ]: 171686 : if(maxerr >= invovl_precision .and. maxerr >= 1e-10) then
1290 : 0 : write(message, *) 'In invovl, max error was', maxerr, ' after 30 iterations'
1291 : 0 : ABI_WARNING(message)
1292 : : else
1293 : 171686 : write(message,'(a,i2,a,es13.5)') 'Iterative solver in invovl finished in ', i, ' iterations, error', maxerr
1294 : 171686 : call wrtout(std_out,message,'COLL')
1295 : : end if
1296 : :
1297 [ + - ]: 171686 : if(gpu_option == ABI_GPU_DISABLED) then
1298 [ - + ]: 171686 : ABI_FREE(temp_proj)
1299 : : else if(gpu_option == ABI_GPU_OPENMP) then
1300 : : #ifdef HAVE_OPENMP_OFFLOAD
1301 : : !$OMP TARGET EXIT DATA MAP(delete:errs,resid,precondresid,normprojs)
1302 : : #endif
1303 : : end if
1304 : :
1305 [ - + ]: 171686 : end subroutine solve_inner
1306 : : !!***
1307 : :
1308 : : !!****f* m_invovl/solve_inner_sqmr
1309 : : !! NAME
1310 : : !! solve_inner_sqmr
1311 : : !!
1312 : : !! FUNCTION
1313 : : !! Helper function: iteratively solves the inner system using SQMR algorithm
1314 : : !!
1315 : : !! INPUTS
1316 : : !!
1317 : : !! SOURCE
1318 : 6976 : subroutine solve_inner_sqmr(invovl, ham, cplx, mpi_enreg, proj, ndat, sm1proj, PtPsm1proj, block_sliced, gpu_option, invovl_precision)
1319 : :
1320 : : use m_abi_linalg
1321 : :
1322 : : integer,intent(in) :: ndat,cplx
1323 : : type(invovl_kpt_type), intent(in), target :: invovl
1324 : : real(dp), intent(in) :: invovl_precision
1325 : : real(dp), intent(inout) :: proj(cplx, invovl%nprojs,ndat)
1326 : : real(dp), intent(inout), target :: sm1proj(cplx, invovl%nprojs, ndat)
1327 : : real(dp), intent(inout), target :: PtPsm1proj(cplx, invovl%nprojs, ndat)
1328 : 6976 : real(dp), allocatable :: temp_proj(:,:,:)
1329 : : type(mpi_type), intent(in) :: mpi_enreg
1330 : : type(gs_hamiltonian_type),intent(in) :: ham
1331 : : integer, intent(in) :: block_sliced, gpu_option
1332 : :
1333 : 13952 : integer :: array_nlmntot_pp(mpi_enreg%nproc_fft)
1334 : : integer :: nlmntot_this_proc, ibeg, iend, ierr, i, nprojs
1335 : : integer :: nstagnated
1336 : : integer, parameter :: maxit = 30
1337 : : real(dp), parameter :: breakdown_tol = 1.0e-300_dp
1338 : : ! Stagnation check: if the (global, over ndat) residual fails to shrink by
1339 : : ! at least a factor stagnation_tol over max_stagnated_it consecutive
1340 : : ! iterations, further iterations are extremely unlikely to help (this
1341 : : ! typically happens once invovl_precision is set close to, or below, the
1342 : : ! relative residual level attainable in double precision for the given
1343 : : ! conditioning of A) : give up early instead of iterating all the way to
1344 : : ! maxit. Requiring several consecutive non-improving iterations (rather
1345 : : ! than a single one) avoids bailing out on transient slow steps, which are
1346 : : ! common with Krylov-type methods and do not by themselves indicate a
1347 : : ! genuine stall.
1348 : : real(dp), parameter :: stagnation_tol = 0.99_dp
1349 : : integer, parameter :: max_stagnated_it = 3
1350 : : real(dp) :: previous_maxerr
1351 : : ! SQMR work vectors: r is the underlying (BiCG-type) Lanczos residual, t is
1352 : : ! the preconditioned residual (M^-1 r, also reused as the "u" vector of the
1353 : : ! algorithm), q is the search direction, Aq = A q. d is the quasi-minimized
1354 : : ! smoothing correction added to sm1proj at each step, and Ad/rx track the
1355 : : ! *exact* residual (rx = proj - A*sm1proj) of the returned iterate at no
1356 : : ! extra matrix-vector product cost (see comments below).
1357 : 13952 : real(dp) :: r(cplx, invovl%nprojs,ndat), t(cplx, invovl%nprojs,ndat)
1358 : 13952 : real(dp) :: q(cplx, invovl%nprojs,ndat), Aq(cplx, invovl%nprojs,ndat)
1359 : 13952 : real(dp) :: d(cplx, invovl%nprojs,ndat), Ad(cplx, invovl%nprojs,ndat)
1360 : 13952 : real(dp) :: rx(cplx, invovl%nprojs,ndat)
1361 : 13952 : real(dp) :: normprojs(ndat), errs(ndat), maxerr
1362 : 13952 : real(dp) :: tau(ndat), tnorm2(ndat), rho(ndat), theta(ndat), sigma(ndat)
1363 : 13952 : real(dp) :: alpha(ndat), beta(ndat), psi(ndat), theta_old(ndat), rho_old(ndat)
1364 : : character(len=500) :: message
1365 : :
1366 : : #ifdef HAVE_OPENMP_OFFLOAD
1367 : : real(dp) :: sum_tmp
1368 : : integer :: idat,iproj,icplx
1369 : : #endif
1370 : :
1371 : : ! *************************************************************************
1372 : :
1373 : 6976 : nprojs = invovl%nprojs
1374 [ + - ]: 6976 : if(gpu_option == ABI_GPU_DISABLED) then
1375 [ + + + + : 3973952 : normprojs = SUM(SUM(proj**2, 1),1)
+ + ]
1376 : 6976 : call xmpi_sum(normprojs, mpi_enreg%comm_fft, ierr)
1377 : : else if(gpu_option == ABI_GPU_OPENMP) then
1378 : : #ifdef HAVE_OPENMP_OFFLOAD
1379 : : !$OMP TARGET ENTER DATA MAP(alloc:errs,normprojs,r,t,q,Aq,d,Ad,rx, &
1380 : : !$OMP& tau,tnorm2,rho,theta,sigma,alpha,beta,psi,theta_old,rho_old)
1381 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:normprojs,proj) PRIVATE(idat,sum_tmp)
1382 : : do idat = 1,ndat
1383 : : sum_tmp=0
1384 : : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:sum_tmp) PRIVATE(iproj,icplx)
1385 : : do iproj = 1,nprojs
1386 : : do icplx = 1,cplx
1387 : : sum_tmp = sum_tmp + proj(icplx,iproj,idat)**2
1388 : : end do
1389 : : end do
1390 : : normprojs(idat)=sum_tmp
1391 : : end do
1392 : : !$OMP TARGET UPDATE FROM(normprojs)
1393 : : #endif
1394 : : end if
1395 : :
1396 : : ! Compute work distribution : split nprojs evenly between the fft processors
1397 : : ! On GPU, work isn't shared on npfft so path is the same wether paral_kgb={0,1}
1398 [ - + - - ]: 6976 : if(mpi_enreg%paral_kgb == 1 .and. gpu_option == ABI_GPU_DISABLED) then
1399 [ # # ]: 0 : array_nlmntot_pp(:) = nprojs / mpi_enreg%nproc_fft
1400 : : ! not enough work, there's MOD(nprojs,mpi_enreg%nproc_fft) tasks left
1401 : : ! assign them to the first ones
1402 [ # # ]: 0 : array_nlmntot_pp(1:MOD(nprojs,mpi_enreg%nproc_fft)) = array_nlmntot_pp(1:MOD(nprojs,mpi_enreg%nproc_fft)) + 1
1403 [ # # ]: 0 : ibeg = SUM(array_nlmntot_pp(1:mpi_enreg%me_fft)) + 1
1404 : 0 : iend = ibeg + array_nlmntot_pp(1+mpi_enreg%me_fft) - 1
1405 : 0 : nlmntot_this_proc = iend - ibeg + 1
1406 : : else
1407 : 6976 : ibeg = 1
1408 : 6976 : iend = nprojs
1409 : 6976 : nlmntot_this_proc = nprojs
1410 : : end if
1411 : :
1412 [ + - ]: 6976 : if(gpu_option == ABI_GPU_DISABLED) then
1413 [ - + + - : 34880 : ABI_MALLOC(temp_proj, (cplx, nlmntot_this_proc, ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
1414 : : end if
1415 : :
1416 : : ! first guess for sm1proj
1417 : 6976 : call apply_block(ham, cplx, invovl%inv_s_approx, nprojs, ndat, proj, sm1proj, block_sliced, gpu_option)
1418 : :
1419 : : ! Solve (D^-1 + PtP) sm1proj = proj with SQMR (symmetric quasi-minimal
1420 : : ! residual, Freund & Nachtigal, "A new Krylov-subspace method for symmetric
1421 : : ! indefinite linear systems", 1994), left-preconditioned by inv_s_approx.
1422 : : ! A = D^-1+PtP is Hermitian but indefinite, and so is the block-diagonal
1423 : : ! preconditioner (built with ZHETRF/ZHETRI in mkinvovl): plain (P)CG/MINRES
1424 : : ! are not guaranteed to work in that setting, whereas SQMR tolerates an
1425 : : ! indefinite preconditioner and converges in markedly fewer iterations than
1426 : : ! the previous Richardson-type iterative refinement.
1427 : : ! See A. Levitt & M. Torrent, https://arxiv.org/abs/1406.4350, section 5.1.
1428 : : !
1429 : : ! All the scalar quantities below (tau, rho, sigma, theta...) are, in exact
1430 : : ! arithmetic, real numbers even though cplx=2 (complex Hermitian case),
1431 : : ! because A and the preconditioner are Hermitian: they are computed with
1432 : : ! dot_per_idat(), which sums v(icplx,iproj,idat)*w(icplx,iproj,idat) over
1433 : : ! icplx without conjugation, which is exactly Re(v^H w) for such quantities.
1434 : : !
1435 : : ! d is the quasi-minimized smoothing correction added to sm1proj at each
1436 : : ! step. Ad and rx are kept consistent with d and sm1proj (rx = proj -
1437 : : ! (D^-1+PtP)*sm1proj exactly) through the same linear recurrence as d,
1438 : : ! reusing Aq that is already computed for the algorithm: this gives the
1439 : : ! *exact* residual of the returned solution at every step, at no extra
1440 : : ! matrix-vector product cost, so that the convergence check below has
1441 : : ! exactly the same meaning as with the previous algorithm.
1442 : :
1443 : 6976 : call compute_Amult(sm1proj, rx)
1444 : 6976 : call vec_sub(proj, rx, r)
1445 : 6976 : call vec_copy(r, rx)
1446 : :
1447 : 6976 : call apply_block(ham, cplx, invovl%inv_s_approx, nprojs, ndat, r, t, block_sliced, gpu_option)
1448 : 6976 : call vec_copy(t, q)
1449 : :
1450 : 6976 : call zero_vec(d)
1451 : 6976 : call zero_vec(Ad)
1452 : 6976 : call zero_scalar(theta)
1453 : :
1454 : 6976 : call dot_per_idat(t, t, tau, need_host=.false.)
1455 : 6976 : call dot_per_idat(t, r, rho, need_host=.false.)
1456 : :
1457 : 6976 : call dot_per_idat(rx, rx, errs)
1458 [ + - - + : 95872 : maxerr = sqrt(MAXVAL(errs/normprojs))
- - + + ]
1459 : :
1460 : 6976 : i = 0
1461 : 6976 : previous_maxerr = maxerr
1462 : 6976 : nstagnated = 0
1463 [ + - ]: 6976 : if(maxerr >= invovl_precision) then
1464 [ + - ]: 40050 : do i=1, maxit
1465 : : #ifdef FC_NVHPC
1466 : : ! Silly fix for NVHPC 25.1
1467 : : if(ndat == -42) write(100,*) ndat
1468 : : #endif
1469 : : ! Aq = A q ; sigma = q . Aq ; alpha = rho/sigma
1470 : 40050 : call compute_Amult(q, Aq)
1471 : 40050 : call dot_per_idat(q, Aq, sigma, need_host=.false.)
1472 : 40050 : call compute_alpha()
1473 : :
1474 : : ! r = r - alpha*Aq
1475 : 40050 : call update_r()
1476 : :
1477 : 40050 : call copy_scalar(theta_old, theta)
1478 : : ! t = M^-1 r (preconditioned residual, doubles as "u" below since we use
1479 : : ! left preconditioning)
1480 : 40050 : call apply_block(ham, cplx, invovl%inv_s_approx, nprojs, ndat, r, t, block_sliced, gpu_option)
1481 : 40050 : call dot_per_idat(t, t, tnorm2, need_host=.false.)
1482 : 40050 : call compute_theta_psi_tau()
1483 : :
1484 : : ! d = psi*theta_old*d + psi*alpha*q ; sm1proj = sm1proj + d
1485 : 40050 : call update_d_x()
1486 : : ! Ad = psi*theta_old*Ad + psi*alpha*Aq ; rx = rx - Ad (keeps rx exact)
1487 : 40050 : call update_Ad_rx()
1488 : :
1489 : 40050 : call dot_per_idat(rx, rx, errs)
1490 [ + - - + : 570204 : maxerr = sqrt(MAXVAL(errs/normprojs))
- - + + ]
1491 [ + + ]: 40050 : if(maxerr < invovl_precision) exit
1492 : :
1493 : : ! Stagnation check (see comment on stagnation_tol/max_stagnated_it above).
1494 [ + + ]: 33074 : if(maxerr > stagnation_tol*previous_maxerr) then
1495 : 8 : nstagnated = nstagnated + 1
1496 [ - + ]: 8 : if(nstagnated >= max_stagnated_it) exit
1497 : : else
1498 : : nstagnated = 0
1499 : : end if
1500 : 33074 : previous_maxerr = maxerr
1501 : :
1502 : 33074 : call copy_scalar(rho_old, rho)
1503 : : ! rho = t . r ; beta = rho/rho_old ; q = t + beta*q
1504 : 33074 : call dot_per_idat(t, r, rho, need_host=.false.)
1505 : 33074 : call compute_beta()
1506 : 40050 : call update_q()
1507 : : end do
1508 : : end if
1509 : :
1510 : : ! Recover PtPsm1proj = PtP*sm1proj from bookkeeping instead of an extra GEMM:
1511 : : ! proj - rx = (D^-1+PtP)*sm1proj, and D^-1*sm1proj is cheap (block-diagonal).
1512 : 6976 : call apply_block(ham, cplx, invovl%inv_sij, nprojs, ndat, sm1proj, r, block_sliced, gpu_option)
1513 : 6976 : call vec_sub(proj, rx, t)
1514 : 6976 : call vec_sub(t, r, PtPsm1proj)
1515 : :
1516 [ - + - - ]: 6976 : if(maxerr >= invovl_precision .and. maxerr >= 1e-10) then
1517 : 0 : write(message, *) 'In invovl, max error was', maxerr, ' after ', i, ' iterations'
1518 : 0 : ABI_WARNING(message)
1519 : : else
1520 : 6976 : write(message,'(a,i2,a,es13.5)') 'Iterative solver in invovl finished in ', i, ' iterations, error', maxerr
1521 : : end if
1522 : 6976 : call wrtout(std_out,message,'COLL')
1523 : 6976 : flush(std_out)
1524 : :
1525 [ + - - + ]: 13952 : if(gpu_option == ABI_GPU_DISABLED) then
1526 [ - + ]: 6976 : ABI_FREE(temp_proj)
1527 : : else if(gpu_option == ABI_GPU_OPENMP) then
1528 : : #ifdef HAVE_OPENMP_OFFLOAD
1529 : : !$OMP TARGET EXIT DATA MAP(delete:errs,normprojs,r,t,q,Aq,d,Ad,rx, &
1530 : : !$OMP& tau,tnorm2,rho,theta,sigma,alpha,beta,psi,theta_old,rho_old)
1531 : : #endif
1532 : : end if
1533 : :
1534 : : contains
1535 : :
1536 : 47026 : subroutine compute_Amult(x, y)
1537 : : ! y := (D^-1 + PtP) x. Uses PtPsm1proj as the GEMM scratch buffer, exactly
1538 : : ! as the previous iterative refinement did.
1539 : : real(dp), intent(inout) :: x(cplx,nprojs,ndat)
1540 : : real(dp), intent(inout) :: y(cplx,nprojs,ndat)
1541 : : integer :: ierr_l
1542 : :
1543 : 47026 : call apply_block(ham, cplx, invovl%inv_sij, nprojs, ndat, x, y, block_sliced, gpu_option)
1544 : :
1545 : : ABI_NVTX_START_RANGE(NVTX_INVOVL_INNER_GEMM)
1546 [ + - ]: 47026 : if(gpu_option == ABI_GPU_DISABLED) then
1547 [ + - + - : 27819868 : temp_proj = x(:,ibeg:iend,:)
+ - - + -
- + + + +
+ + ]
1548 : :
1549 : : ! compute matrix multiplication : PtPsm1proj(:,:,1) = invovl%gram * temp_proj(:,:,1)
1550 : : call abi_xgemm('N', 'N', nprojs, ndat, nlmntot_this_proc, cone, &
1551 : : & invovl%gram_projs(:,:,1), nprojs, &
1552 : : & temp_proj(:,:,1), nlmntot_this_proc, czero, &
1553 : : & PtPsm1proj(:,:,1), nprojs, &
1554 : 47026 : & x_cplx=cplx)
1555 : 47026 : call xmpi_sum(PtPsm1proj, mpi_enreg%comm_fft, ierr_l)
1556 [ + + + + : 27772842 : y = y + PtPsm1proj
+ + ]
1557 [ # # ]: 0 : else if(gpu_option == ABI_GPU_OPENMP) then
1558 : : ! compute matrix multiplication : PtPsm1proj(:,:,1) = invovl%gram * x(:,:,1)
1559 : : call abi_xgemm('N', 'N', nprojs, ndat, nlmntot_this_proc, cone, &
1560 : : & invovl%gram_projs, nprojs,&
1561 : : & x, nlmntot_this_proc, czero, &
1562 : : & PtPsm1proj, nprojs, &
1563 : 0 : & x_cplx=cplx, gpu_option=ABI_GPU_OPENMP)
1564 : :
1565 : 0 : call abi_xaxpy(cplx*nprojs*ndat, cone, PtPsm1proj, 1, y, 1, gpu_option=gpu_option)
1566 : : end if
1567 : : ABI_NVTX_END_RANGE()
1568 : 47026 : end subroutine compute_Amult
1569 : :
1570 : 174152 : subroutine dot_per_idat(v, w, out_arr, need_host)
1571 : : ! out_arr(idat) = sum_{iproj,icplx} v(icplx,iproj,idat)*w(icplx,iproj,idat)
1572 : : ! need_host (GPU only, default .true.): whether out_arr must be made
1573 : : ! available on the host after this call. Most call sites only feed the
1574 : : ! result into further device-resident scalar arithmetic (see
1575 : : ! compute_alpha/compute_theta_psi_tau/compute_beta below) and can skip the
1576 : : ! blocking host<->device round trip entirely; only the calls computing
1577 : : ! errs (used for the host-side convergence/stagnation checks) need it.
1578 : : real(dp), intent(in) :: v(cplx,nprojs,ndat), w(cplx,nprojs,ndat)
1579 : : real(dp), intent(out) :: out_arr(ndat)
1580 : : logical, intent(in), optional :: need_host
1581 : : integer :: ierr_l
1582 : : logical :: need_host_
1583 : : #ifdef HAVE_OPENMP_OFFLOAD
1584 : : integer :: idat_l, iproj_l, icplx_l
1585 : : real(dp) :: sum_tmp_l
1586 : : #endif
1587 : :
1588 : 174152 : need_host_ = .true. ; if(present(need_host)) need_host_ = need_host
1589 [ + - ]: 174152 : if(gpu_option == ABI_GPU_DISABLED) then
1590 [ + + + + : 103143464 : out_arr = SUM(SUM(v*w, 1), 1)
+ + ]
1591 : 174152 : call xmpi_sum(out_arr, mpi_enreg%comm_fft, ierr_l)
1592 : : else if(gpu_option == ABI_GPU_OPENMP) then
1593 : : #ifdef HAVE_OPENMP_OFFLOAD
1594 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:out_arr,v,w) PRIVATE(idat_l,sum_tmp_l)
1595 : : do idat_l = 1,ndat
1596 : : sum_tmp_l=0
1597 : : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:sum_tmp_l) PRIVATE(iproj_l,icplx_l)
1598 : : do iproj_l = 1,nprojs
1599 : : do icplx_l = 1,cplx
1600 : : sum_tmp_l = sum_tmp_l + v(icplx_l,iproj_l,idat_l)*w(icplx_l,iproj_l,idat_l)
1601 : : end do
1602 : : end do
1603 : : out_arr(idat_l)=sum_tmp_l
1604 : : end do
1605 : : if(need_host_) then
1606 : : !$OMP TARGET UPDATE FROM(out_arr)
1607 : : end if
1608 : : #endif
1609 : : end if
1610 : 174152 : end subroutine dot_per_idat
1611 : :
1612 : 6976 : subroutine zero_scalar(y)
1613 : : ! y := 0 (ndat-sized scalar work array, e.g. theta)
1614 : : real(dp), intent(out) :: y(ndat)
1615 : : #ifdef HAVE_OPENMP_OFFLOAD
1616 : : integer :: idat_l
1617 : : #endif
1618 : :
1619 [ + - ]: 6976 : if(gpu_option == ABI_GPU_DISABLED) then
1620 [ + + ]: 88896 : y = zero
1621 : : else if(gpu_option == ABI_GPU_OPENMP) then
1622 : : #ifdef HAVE_OPENMP_OFFLOAD
1623 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(idat_l) MAP(to:y)
1624 : : do idat_l = 1,ndat
1625 : : y(idat_l) = zero
1626 : : end do
1627 : : #endif
1628 : : end if
1629 : 6976 : end subroutine zero_scalar
1630 : :
1631 : 73124 : subroutine copy_scalar(y, x)
1632 : : ! y := x (ndat-sized scalar work arrays, e.g. theta_old := theta). Done on
1633 : : ! device (instead of a plain host assignment) so the device-resident copy
1634 : : ! stays in sync without an extra host<->device round trip.
1635 : : real(dp), intent(out) :: y(ndat)
1636 : : real(dp), intent(in) :: x(ndat)
1637 : : #ifdef HAVE_OPENMP_OFFLOAD
1638 : : integer :: idat_l
1639 : : #endif
1640 : :
1641 [ + - ]: 73124 : if(gpu_option == ABI_GPU_DISABLED) then
1642 [ + + ]: 971412 : y = x
1643 : : else if(gpu_option == ABI_GPU_OPENMP) then
1644 : : #ifdef HAVE_OPENMP_OFFLOAD
1645 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(idat_l) MAP(to:y,x)
1646 : : do idat_l = 1,ndat
1647 : : y(idat_l) = x(idat_l)
1648 : : end do
1649 : : #endif
1650 : : end if
1651 : 73124 : end subroutine copy_scalar
1652 : :
1653 : 40050 : subroutine compute_alpha()
1654 : : ! sigma := sign(breakdown_tol,sigma) if too small ; alpha := rho/sigma.
1655 : : ! Kept entirely device-resident: sigma/rho/alpha are all already present
1656 : : ! on the device (never synced to host), so no round trip is needed here.
1657 : : #ifdef HAVE_OPENMP_OFFLOAD
1658 : : integer :: idat_l
1659 : : #endif
1660 : :
1661 [ + - ]: 40050 : if(gpu_option == ABI_GPU_DISABLED) then
1662 [ + + - + ]: 530154 : where(abs(sigma) < breakdown_tol) sigma = sign(breakdown_tol, sigma)
1663 [ + + ]: 530154 : alpha = rho / sigma
1664 : : else if(gpu_option == ABI_GPU_OPENMP) then
1665 : : #ifdef HAVE_OPENMP_OFFLOAD
1666 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(idat_l) MAP(to:sigma,rho,alpha)
1667 : : do idat_l = 1,ndat
1668 : : if(abs(sigma(idat_l)) < breakdown_tol) sigma(idat_l) = sign(breakdown_tol, sigma(idat_l))
1669 : : alpha(idat_l) = rho(idat_l) / sigma(idat_l)
1670 : : end do
1671 : : #endif
1672 : : end if
1673 : 40050 : end subroutine compute_alpha
1674 : :
1675 : 40050 : subroutine compute_theta_psi_tau()
1676 : : ! tau := sign(breakdown_tol,tau) if too small ; theta := tnorm2/tau ;
1677 : : ! psi := 1/(1+theta) ; tau := tnorm2*psi. Kept entirely device-resident.
1678 : : #ifdef HAVE_OPENMP_OFFLOAD
1679 : : integer :: idat_l
1680 : : #endif
1681 : :
1682 [ + - ]: 40050 : if(gpu_option == ABI_GPU_DISABLED) then
1683 [ + + - + ]: 530154 : where(abs(tau) < breakdown_tol) tau = sign(breakdown_tol, tau)
1684 [ + + ]: 530154 : theta = tnorm2 / tau
1685 [ + + ]: 530154 : psi = one / (one + theta)
1686 [ + + ]: 530154 : tau = tnorm2 * psi
1687 : : else if(gpu_option == ABI_GPU_OPENMP) then
1688 : : #ifdef HAVE_OPENMP_OFFLOAD
1689 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(idat_l) &
1690 : : !$OMP& MAP(to:tau,theta,psi,tnorm2)
1691 : : do idat_l = 1,ndat
1692 : : if(abs(tau(idat_l)) < breakdown_tol) tau(idat_l) = sign(breakdown_tol, tau(idat_l))
1693 : : theta(idat_l) = tnorm2(idat_l) / tau(idat_l)
1694 : : psi(idat_l) = one / (one + theta(idat_l))
1695 : : tau(idat_l) = tnorm2(idat_l) * psi(idat_l)
1696 : : end do
1697 : : #endif
1698 : : end if
1699 : 40050 : end subroutine compute_theta_psi_tau
1700 : :
1701 : 33074 : subroutine compute_beta()
1702 : : ! rho_old := sign(breakdown_tol,rho_old) if too small ; beta := rho/rho_old.
1703 : : ! Kept entirely device-resident.
1704 : : #ifdef HAVE_OPENMP_OFFLOAD
1705 : : integer :: idat_l
1706 : : #endif
1707 : :
1708 [ + - ]: 33074 : if(gpu_option == ABI_GPU_DISABLED) then
1709 [ + + - + ]: 441258 : where(abs(rho_old) < breakdown_tol) rho_old = sign(breakdown_tol, rho_old)
1710 [ + + ]: 441258 : beta = rho / rho_old
1711 : : else if(gpu_option == ABI_GPU_OPENMP) then
1712 : : #ifdef HAVE_OPENMP_OFFLOAD
1713 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(idat_l) MAP(to:rho_old,rho,beta)
1714 : : do idat_l = 1,ndat
1715 : : if(abs(rho_old(idat_l)) < breakdown_tol) rho_old(idat_l) = sign(breakdown_tol, rho_old(idat_l))
1716 : : beta(idat_l) = rho(idat_l) / rho_old(idat_l)
1717 : : end do
1718 : : #endif
1719 : : end if
1720 : 33074 : end subroutine compute_beta
1721 : :
1722 : 20928 : subroutine vec_sub(a, b, y)
1723 : : ! y := a - b (full array, no per-idat scaling)
1724 : : real(dp), intent(in) :: a(cplx,nprojs,ndat), b(cplx,nprojs,ndat)
1725 : : real(dp), intent(out) :: y(cplx,nprojs,ndat)
1726 : : #ifdef HAVE_OPENMP_OFFLOAD
1727 : : integer :: idat_l, iproj_l, icplx_l
1728 : : #endif
1729 : :
1730 [ + - ]: 20928 : if(gpu_option == ABI_GPU_DISABLED) then
1731 [ + + + + : 11921856 : y = a - b
+ + ]
1732 : : else if(gpu_option == ABI_GPU_OPENMP) then
1733 : : #ifdef HAVE_OPENMP_OFFLOAD
1734 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1735 : : !$OMP& PRIVATE(idat_l,iproj_l,icplx_l) MAP(to:y,a,b)
1736 : : do idat_l = 1,ndat
1737 : : do iproj_l = 1,nprojs
1738 : : do icplx_l = 1,cplx
1739 : : y(icplx_l,iproj_l,idat_l) = a(icplx_l,iproj_l,idat_l) - b(icplx_l,iproj_l,idat_l)
1740 : : end do
1741 : : end do
1742 : : end do
1743 : : #endif
1744 : : end if
1745 : 20928 : end subroutine vec_sub
1746 : :
1747 : 13952 : subroutine vec_copy(a, y)
1748 : : ! y := a
1749 : : real(dp), intent(in) :: a(cplx,nprojs,ndat)
1750 : : real(dp), intent(out) :: y(cplx,nprojs,ndat)
1751 : :
1752 [ + - ]: 13952 : if(gpu_option == ABI_GPU_DISABLED) then
1753 [ + + + + : 7947904 : y = a
+ + ]
1754 : : else if(gpu_option == ABI_GPU_OPENMP) then
1755 : : #ifdef HAVE_OPENMP_OFFLOAD
1756 : : call gpu_copy(y, a, int(cplx,c_size_t)*nprojs*ndat)
1757 : : #endif
1758 : : end if
1759 : 13952 : end subroutine vec_copy
1760 : :
1761 : 13952 : subroutine zero_vec(y)
1762 : : ! y := 0
1763 : : real(dp), intent(out) :: y(cplx,nprojs,ndat)
1764 : :
1765 [ + - ]: 13952 : if(gpu_option == ABI_GPU_DISABLED) then
1766 [ + + + + : 7947904 : y = zero
+ + ]
1767 : : else if(gpu_option == ABI_GPU_OPENMP) then
1768 : : #ifdef HAVE_OPENMP_OFFLOAD
1769 : : call gpu_set_to_zero(y,int(cplx,c_size_t)*nprojs*ndat)
1770 : : #endif
1771 : : end if
1772 : 13952 : end subroutine zero_vec
1773 : :
1774 : 40050 : subroutine update_r()
1775 : : ! r := r - alpha(idat)*Aq. alpha/r/Aq are all host-associated, already
1776 : : ! device-resident arrays: no coefficient argument is passed (an expression
1777 : : ! like "-alpha" would force a host-side compiler temporary, defeating the
1778 : : ! whole point of keeping alpha device-resident), and there is no
1779 : : ! TARGET UPDATE TO left to do.
1780 : : integer :: idat_l
1781 : : #ifdef HAVE_OPENMP_OFFLOAD
1782 : : integer :: iproj_l, icplx_l
1783 : : #endif
1784 : :
1785 [ + - ]: 40050 : if(gpu_option == ABI_GPU_DISABLED) then
1786 [ + + ]: 530154 : do idat_l = 1,ndat
1787 [ + + + + ]: 23798890 : r(:,:,idat_l) = r(:,:,idat_l) - alpha(idat_l)*Aq(:,:,idat_l)
1788 : : end do
1789 : : else if(gpu_option == ABI_GPU_OPENMP) then
1790 : : #ifdef HAVE_OPENMP_OFFLOAD
1791 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) PRIVATE(idat_l,iproj_l,icplx_l) MAP(to:r,Aq,alpha)
1792 : : do idat_l = 1,ndat
1793 : : do iproj_l = 1,nprojs
1794 : : do icplx_l = 1,cplx
1795 : : r(icplx_l,iproj_l,idat_l) = r(icplx_l,iproj_l,idat_l) - alpha(idat_l)*Aq(icplx_l,iproj_l,idat_l)
1796 : : end do
1797 : : end do
1798 : : end do
1799 : : #endif
1800 : : end if
1801 : 40050 : end subroutine update_r
1802 : :
1803 : 40050 : subroutine update_d_x()
1804 : : ! d := psi(idat)*theta_old(idat)*d + psi(idat)*alpha(idat)*q ;
1805 : : ! sm1proj := sm1proj + d. No arguments: psi/theta_old/alpha are
1806 : : ! host-associated device-resident arrays, and the products are formed
1807 : : ! inside the kernel instead of being passed in (which would force
1808 : : ! host-side temporaries for "psi*theta_old"/"psi*alpha" and lose device
1809 : : ! residency).
1810 : : integer :: idat_l
1811 : : #ifdef HAVE_OPENMP_OFFLOAD
1812 : : integer :: iproj_l, icplx_l
1813 : : #endif
1814 : :
1815 [ + - ]: 40050 : if(gpu_option == ABI_GPU_DISABLED) then
1816 [ + + ]: 530154 : do idat_l = 1,ndat
1817 [ + + + + ]: 23798890 : d(:,:,idat_l) = psi(idat_l)*theta_old(idat_l)*d(:,:,idat_l) + psi(idat_l)*alpha(idat_l)*q(:,:,idat_l)
1818 : : end do
1819 [ + + + + : 23798890 : sm1proj = sm1proj + d
+ + ]
1820 : : else if(gpu_option == ABI_GPU_OPENMP) then
1821 : : #ifdef HAVE_OPENMP_OFFLOAD
1822 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1823 : : !$OMP& PRIVATE(idat_l,iproj_l,icplx_l) MAP(to:d,q,sm1proj,psi,theta_old,alpha)
1824 : : do idat_l = 1,ndat
1825 : : do iproj_l = 1,nprojs
1826 : : do icplx_l = 1,cplx
1827 : : d(icplx_l,iproj_l,idat_l) = psi(idat_l)*theta_old(idat_l)*d(icplx_l,iproj_l,idat_l) + &
1828 : : & psi(idat_l)*alpha(idat_l)*q(icplx_l,iproj_l,idat_l)
1829 : : sm1proj(icplx_l,iproj_l,idat_l) = sm1proj(icplx_l,iproj_l,idat_l) + d(icplx_l,iproj_l,idat_l)
1830 : : end do
1831 : : end do
1832 : : end do
1833 : : #endif
1834 : : end if
1835 : 40050 : end subroutine update_d_x
1836 : :
1837 : 40050 : subroutine update_Ad_rx()
1838 : : ! Ad := psi(idat)*theta_old(idat)*Ad + psi(idat)*alpha(idat)*Aq ;
1839 : : ! rx := rx - Ad. Same no-argument rationale as update_d_x.
1840 : : integer :: idat_l
1841 : : #ifdef HAVE_OPENMP_OFFLOAD
1842 : : integer :: iproj_l, icplx_l
1843 : : #endif
1844 : :
1845 [ + - ]: 40050 : if(gpu_option == ABI_GPU_DISABLED) then
1846 [ + + ]: 530154 : do idat_l = 1,ndat
1847 [ + + + + ]: 23798890 : Ad(:,:,idat_l) = psi(idat_l)*theta_old(idat_l)*Ad(:,:,idat_l) + psi(idat_l)*alpha(idat_l)*Aq(:,:,idat_l)
1848 : : end do
1849 [ + + + + : 23798890 : rx = rx - Ad
+ + ]
1850 : : else if(gpu_option == ABI_GPU_OPENMP) then
1851 : : #ifdef HAVE_OPENMP_OFFLOAD
1852 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1853 : : !$OMP& PRIVATE(idat_l,iproj_l,icplx_l) MAP(to:Ad,Aq,rx,psi,theta_old,alpha)
1854 : : do idat_l = 1,ndat
1855 : : do iproj_l = 1,nprojs
1856 : : do icplx_l = 1,cplx
1857 : : Ad(icplx_l,iproj_l,idat_l) = psi(idat_l)*theta_old(idat_l)*Ad(icplx_l,iproj_l,idat_l) + &
1858 : : & psi(idat_l)*alpha(idat_l)*Aq(icplx_l,iproj_l,idat_l)
1859 : : rx(icplx_l,iproj_l,idat_l) = rx(icplx_l,iproj_l,idat_l) - Ad(icplx_l,iproj_l,idat_l)
1860 : : end do
1861 : : end do
1862 : : end do
1863 : : #endif
1864 : : end if
1865 : 40050 : end subroutine update_Ad_rx
1866 : :
1867 : 33074 : subroutine update_q()
1868 : : ! q := t + beta(idat)*q. beta is host-associated and already
1869 : : ! device-resident (computed in place by compute_beta).
1870 : : integer :: idat_l
1871 : : #ifdef HAVE_OPENMP_OFFLOAD
1872 : : integer :: iproj_l, icplx_l
1873 : : #endif
1874 : :
1875 [ + - ]: 33074 : if(gpu_option == ABI_GPU_DISABLED) then
1876 [ + + ]: 441258 : do idat_l = 1,ndat
1877 [ + + + + ]: 19824938 : q(:,:,idat_l) = t(:,:,idat_l) + beta(idat_l)*q(:,:,idat_l)
1878 : : end do
1879 : : else if(gpu_option == ABI_GPU_OPENMP) then
1880 : : #ifdef HAVE_OPENMP_OFFLOAD
1881 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1882 : : !$OMP& PRIVATE(idat_l,iproj_l,icplx_l) MAP(to:q,t,beta)
1883 : : do idat_l = 1,ndat
1884 : : do iproj_l = 1,nprojs
1885 : : do icplx_l = 1,cplx
1886 : : q(icplx_l,iproj_l,idat_l) = t(icplx_l,iproj_l,idat_l) + beta(idat_l)*q(icplx_l,iproj_l,idat_l)
1887 : : end do
1888 : : end do
1889 : : end do
1890 : : #endif
1891 : : end if
1892 : 33074 : end subroutine update_q
1893 : :
1894 : : end subroutine solve_inner_sqmr
1895 : : !!***
1896 : :
1897 : : !!****f* m_invovl/apply_block
1898 : : !! NAME
1899 : : !! apply_block
1900 : : !!
1901 : : !! FUNCTION
1902 : : !! Helper function: applies a block-diagonal matrix mat(lmnmax, lmnmax, ntypat)
1903 : : !!
1904 : : !! INPUTS
1905 : : !!
1906 : : !! SOURCE
1907 : 3209126 : subroutine apply_block(ham, cplx, mat, nprojs, ndat, x, y, block_sliced, gpu_option)
1908 : :
1909 : : use m_abi_linalg
1910 : :
1911 : : integer,intent(in) :: ndat, nprojs, cplx
1912 : : real(dp), intent(inout), target :: x(cplx, nprojs, ndat), y(cplx, nprojs, ndat)
1913 : : type(gs_hamiltonian_type),intent(in) :: ham
1914 : : real(dp), intent(in), target :: mat(cplx, ham%lmnmax, ham%lmnmax, ham%ntypat)
1915 : : integer, intent(in) :: block_sliced, gpu_option
1916 : :
1917 : : integer :: nlmn, shift, itypat, idat
1918 : 3209126 : real(dp),pointer :: work_x(:,:),work_y(:,:)
1919 : 3209126 : real(dp), ABI_CONTIGUOUS pointer :: mat_ptr(:,:,:)
1920 : :
1921 : : ! *************************************************************************
1922 : :
1923 : 3209126 : if (block_sliced == 1) then
1924 : :
1925 [ + + ]: 26208684 : do idat = 1, ndat
1926 : 23077616 : shift = 1
1927 [ + + ]: 57915764 : do itypat=1, ham%ntypat
1928 [ + + + + ]: 305616460 : nlmn = count(ham%indlmn(3,:,itypat)>0)
1929 : : !! apply mat to all atoms at once
1930 : : ! perform natom multiplications of size nlmn
1931 : : ! compute y = mat*x
1932 : 31707080 : work_x => x(:, shift:shift+nlmn*ham%nattyp(itypat)-1, idat)
1933 : 31707080 : work_y => y(:, shift:shift+nlmn*ham%nattyp(itypat)-1, idat)
1934 : 31707080 : mat_ptr => mat(:, :, :, itypat)
1935 [ + + ]: 31707080 : if(cplx == 2) then
1936 [ + - ]: 28701952 : if(gpu_option /= ABI_GPU_OPENMP) then
1937 : : call ZHEMM('L','U', nlmn, ham%nattyp(itypat), cone, &
1938 : : & mat_ptr, ham%lmnmax, &
1939 : : & work_x, nlmn, czero, &
1940 [ - + - + ]: 28701952 : & work_y, nlmn)
1941 : : else
1942 : : #ifdef HAVE_OPENMP_OFFLOAD
1943 : : !$OMP TARGET DATA USE_DEVICE_ADDR(mat_ptr,work_x,work_y)
1944 : : call abi_gpu_zhemm('L','U', nlmn, ham%nattyp(itypat), cone, &
1945 : : c_loc(mat_ptr), ham%lmnmax, &
1946 : : c_loc(work_x), nlmn, czero, &
1947 : : c_loc(work_y), nlmn)
1948 : : !$OMP END TARGET DATA
1949 : : #endif
1950 : : end if
1951 : : else
1952 [ + - ]: 3005128 : if(gpu_option /= ABI_GPU_OPENMP) then
1953 : : call DSYMM('L','U', nlmn, ham%nattyp(itypat), one, &
1954 : : & mat_ptr, ham%lmnmax, &
1955 : : & work_x, nlmn, zero, &
1956 [ - + - + ]: 3005128 : & work_y, nlmn)
1957 : : else
1958 : : #ifdef HAVE_OPENMP_OFFLOAD
1959 : : !$OMP TARGET DATA USE_DEVICE_ADDR(mat_ptr,work_x,work_y)
1960 : : call abi_gpu_xsymm(cplx, 'L','U', nlmn, ham%nattyp(itypat), cone, &
1961 : : c_loc(mat_ptr), ham%lmnmax, &
1962 : : c_loc(work_x), nlmn, czero, &
1963 : : c_loc(work_y), nlmn)
1964 : : !$OMP END TARGET DATA
1965 : : #endif
1966 : : end if
1967 : : end if
1968 : 54784696 : shift = shift + nlmn*ham%nattyp(itypat)
1969 : : end do
1970 : : end do
1971 : :
1972 : : else ! block_sliced = 0
1973 : :
1974 : : shift = 1
1975 [ + + ]: 234174 : do itypat=1, ham%ntypat
1976 [ + + + - ]: 1405044 : nlmn = count(ham%indlmn(3,:,itypat)>0)
1977 : 156116 : mat_ptr => mat(:, :, :, itypat)
1978 : : !! apply mat to all atoms at once, all idat at once
1979 : : ! perform natom multiplications of size nlmn
1980 : : ! be careful here matrix extracted from x and y are not memory contiguous
1981 : : ! the GPU version use batched GEMM with strides to workaround this
1982 [ + - ]: 156116 : if(gpu_option == ABI_GPU_DISABLED) then
1983 [ + + ]: 156116 : if(cplx == 2) then
1984 : : call ZHEMM('L','U', nlmn, ham%nattyp(itypat)*ndat, cone, &
1985 : : & mat(:, :, :, itypat), ham%lmnmax, &
1986 : : & x(:, shift:shift+nlmn*ham%nattyp(itypat)-1, 1:ndat), nlmn, czero, &
1987 [ + - + - : 117843120 : & y(:, shift:shift+nlmn*ham%nattyp(itypat)-1, 1:ndat), nlmn)
+ + + + +
+ + - + -
+ + + + +
+ + - + +
+ + + + +
- + + + +
+ + ]
1988 : : else
1989 : : call DSYMM('L','U', nlmn, ham%nattyp(itypat)*ndat, one, &
1990 : : & mat(:, :, :, itypat), ham%lmnmax, &
1991 : : & x(:, shift:shift+nlmn*ham%nattyp(itypat)-1, 1:ndat), nlmn, zero, &
1992 [ + - + - : 4902820 : & y(:, shift:shift+nlmn*ham%nattyp(itypat)-1, 1:ndat), nlmn)
+ + + + +
+ + - + -
+ + + + +
+ + - + +
+ + + + +
- + + + +
+ + ]
1993 : : end if
1994 : : else if(gpu_option == ABI_GPU_OPENMP) then
1995 : : #ifdef HAVE_OPENMP_OFFLOAD
1996 : : !$OMP TARGET DATA USE_DEVICE_ADDR(mat_ptr,x,y)
1997 : : call abi_gpu_xgemm_strided(cplx, 'N','N', &
1998 : : nlmn, ham%nattyp(itypat), nlmn, cone, &
1999 : : c_loc(mat_ptr), ham%lmnmax, 0, &
2000 : : c_loc(x(1,shift,1)), nlmn, nprojs, &
2001 : : czero, &
2002 : : c_loc(y(1,shift,1)), nlmn, nprojs, ndat)
2003 : : !$OMP END TARGET DATA
2004 : : #endif
2005 : : end if
2006 : 234174 : shift = shift + nlmn*ham%nattyp(itypat)
2007 : : end do
2008 : :
2009 : : end if
2010 : :
2011 [ + + ]: 3209126 : end subroutine apply_block
2012 : : !!***
2013 : :
2014 : 0 : function invovl_ompgpu_work_mem(ham, ndat) result(req_mem)
2015 : :
2016 : : type(gs_hamiltonian_type), intent(in) :: ham
2017 : : integer, intent(in) :: ndat
2018 : : integer :: nprojs, cplx, itypat
2019 : : integer(kind=c_size_t) :: req_mem
2020 : :
2021 : 0 : nprojs = 0
2022 [ # # ]: 0 : do itypat=1,ham%ntypat
2023 [ # # # # ]: 0 : nprojs = nprojs + count(ham%indlmn(3,:,itypat)>0)*ham%nattyp(itypat)
2024 : : end do
2025 [ # # ]: 0 : cplx = 2; if(ham%istwf_k == 2) cplx = 1
2026 : :
2027 : 0 : req_mem = 0
2028 : 0 : req_mem = req_mem + dp * cplx * int(nprojs, c_size_t) * int(ndat, c_size_t) ! proj
2029 : 0 : req_mem = req_mem + dp * cplx * int(nprojs, c_size_t) * int(ndat, c_size_t) ! sm1proj
2030 : 0 : req_mem = req_mem + dp * cplx * int(nprojs, c_size_t) * int(ndat, c_size_t) ! PtPsm1proj
2031 : : ! solve_inner (SQMR): r, t, q, Aq, d, Ad, rx
2032 : 0 : req_mem = req_mem + 7 * dp * cplx * int(nprojs, c_size_t) * int(ndat, c_size_t)
2033 : : ! solve_inner (SQMR scalars): tau, tnorm2, rho, theta, sigma, alpha, beta, psi, theta_old, rho_old, errs, normprojs
2034 : 0 : req_mem = req_mem + 12 * dp * int(ndat, c_size_t)
2035 : :
2036 : 0 : end function invovl_ompgpu_work_mem
2037 : :
2038 : 0 : function invovl_ompgpu_static_mem(ham) result(req_mem)
2039 : :
2040 : : type(gs_hamiltonian_type), intent(in) :: ham
2041 : : integer :: nprojs, cplx, itypat
2042 : : integer(kind=c_size_t) :: req_mem
2043 : :
2044 : 0 : nprojs = 0
2045 [ # # ]: 0 : do itypat=1,ham%ntypat
2046 [ # # # # ]: 0 : nprojs = nprojs + count(ham%indlmn(3,:,itypat)>0)*ham%nattyp(itypat)
2047 : : end do
2048 [ # # ]: 0 : cplx = 2; if(ham%istwf_k > 1) cplx = 1
2049 : :
2050 : 0 : req_mem = 0
2051 : 0 : req_mem = req_mem + dp * cplx * int(nprojs, c_size_t) * int(nprojs, c_size_t) ! gram_projs
2052 : : req_mem = req_mem + dp * cplx * int(ham%lmnmax, c_size_t) * &
2053 : 0 : & int(ham%lmnmax, c_size_t) * int(ham%ntypat, c_size_t) ! inv_sij
2054 : : req_mem = req_mem + dp * cplx * int(ham%lmnmax, c_size_t) * &
2055 : 0 : & int(ham%lmnmax, c_size_t) * int(ham%ntypat, c_size_t) ! inv_s_approx
2056 : :
2057 : 0 : end function invovl_ompgpu_static_mem
2058 : :
2059 [ # # # # : 0 : end MODULE m_invovl
# # # # #
# # # # #
# # # # ]
2060 : : !!***
|
