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