Branch data Line data Source code
1 : : !!****f* ABINIT/m_slicewf
2 : : !! NAME
3 : : !! m_slicewf
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains a routine updating the whole wave functions at a given k-point,
7 : : !! using the Spectrum Slicing method (2021 implementation using xG abstraction layer)
8 : : !! for a given spin-polarization, from a fixed Hamiltonian but might also simply compute
9 : : !! eigenvectors and eigenvalues at this k point. it will also update the matrix elements
10 : : !! of the Hamiltonian.
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 2018-2025 ABINIT group (BS, IML)
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 : : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
18 : : !!
19 : : !! SOURCE
20 : :
21 : : #if defined HAVE_CONFIG_H
22 : : #include "config.h"
23 : : #endif
24 : :
25 : : #include "abi_common.h"
26 : :
27 : : ! nvtx related macro definition
28 : : #include "nvtx_macros.h"
29 : :
30 : : module m_slicewf
31 : :
32 : : use defs_abitypes
33 : : use defs_basis
34 : : use m_abicore
35 : : use m_errors
36 : : use m_fstrings
37 : : use m_time
38 : :
39 : : use m_slice
40 : : !use m_chebfi
41 : : !use m_chebfi2
42 : : use m_invovl
43 : :
44 : : use m_cgtools, only : dotprod_g
45 : : use m_dtset, only : dataset_type
46 : :
47 : : use m_hamiltonian, only : gs_hamiltonian_type
48 : : use m_pawcprj, only : pawcprj_type
49 : : use m_nonlop, only : nonlop
50 : : use m_prep_kgb, only : prep_getghc, prep_nonlop
51 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
52 : : use m_getghc, only : multithreaded_getghc
53 : : use m_gemm_nonlop_projectors , only : gemm_nonlop_use_gemm
54 : : use m_mpinfo, only : copy_mpi_enreg
55 : :
56 : : use m_xg
57 : : use m_xgTransposer
58 : :
59 : : #if defined(HAVE_GPU_MARKERS)
60 : : use m_nvtx_data
61 : : #endif
62 : :
63 : : #if defined(HAVE_GPU)
64 : : use m_gpu_toolbox
65 : : #endif
66 : :
67 : : #if defined(HAVE_YAKL)
68 : : use gator_mod
69 : : #endif
70 : :
71 : : use, intrinsic :: iso_c_binding, only: c_associated,c_loc,c_ptr,c_f_pointer,c_double,c_size_t
72 : :
73 : : use m_xmpi
74 : : use m_xomp
75 : : #ifdef HAVE_OPENMP
76 : : use omp_lib
77 : : #endif
78 : :
79 : : implicit none
80 : :
81 : : private
82 : :
83 : : integer, parameter :: l_tim_getghc=7
84 : : real(dp), parameter :: inv_sqrt2 = 1/sqrt2
85 : :
86 : : ! For use in getghc_gsc1
87 : : integer, save :: l_cpopt
88 : : integer, save :: l_icplx
89 : : integer, save :: l_npw
90 : : integer, save :: l_nband_filter
91 : : integer, save :: l_nspinor
92 : : logical, save :: l_paw
93 : : integer, save :: l_prtvol
94 : : integer, save :: l_sij_opt
95 : : integer, save :: l_paral_kgb
96 : : integer, save :: l_useria
97 : : integer, save :: l_block_sliced
98 : : real(dp), save :: l_invovl_precision
99 : :
100 : : #ifdef HAVE_OPENMP_OFFLOAD
101 : : ! For use in getBm1X
102 : : integer, save, private :: l_blockdim_evil = 0
103 : : #endif
104 : :
105 : : type(mpi_type),pointer,save :: l_mpi_enreg
106 : : type(gs_hamiltonian_type),pointer,save :: l_gs_hamk
107 : :
108 : : integer, parameter :: DEBUG_ROWS = 5
109 : : integer, parameter :: DEBUG_COLUMNS = 5
110 : :
111 : : public :: slicewf
112 : :
113 : : CONTAINS !========================================================================================
114 : : !!***
115 : :
116 : : !!****f* m_slicewf/slicewf
117 : : !! NAME
118 : : !! slicewf
119 : : !!
120 : : !! FUNCTION
121 : : !! This routine updates the whole wave functions set at a given k-point,
122 : : !! using the Spectrum Slicing method (2025 version using xG abstraction layer)
123 : : !!
124 : : !! INPUTS
125 : : !! dtset= input variables for this dataset
126 : : !! mpi_enreg= MPI-parallelisation information
127 : : !! nband= number of bands at this k point
128 : : !! npw= number of plane waves at this k point
129 : : !! nspinor= number of spinorial components of the wavefunctions
130 : : !! prtvol= control print volume and debugging
131 : : !!
132 : : !! OUTPUT
133 : : !! eig(nband)= eigenvalues (hartree) for all bands
134 : : !! enl_out(nband)= contribution of each band to the nl part of energy
135 : : !! resid(nband)= residuals for each band
136 : : !!
137 : : !! SIDE EFFECTS
138 : : !! cg(2,npw*nspinor*nband)= planewave coefficients of wavefunctions
139 : : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
140 : : !!
141 : : !! SOURCE
142 : :
143 : 18 : subroutine slicewf(cg,dtset,eig,enl_out,gs_hamk,mpi_enreg,&
144 : 18 : & nband,npw,nspinor,prtvol,resid)
145 : :
146 : : implicit none
147 : :
148 : : ! Arguments ------------------------------------
149 : : integer,intent(in) :: nband,npw,prtvol,nspinor
150 : : type(mpi_type),target,intent(in) :: mpi_enreg
151 : : real(dp),target,intent(inout) :: cg(2,npw*nspinor*nband)
152 : : real(dp),target,intent(out) :: resid(nband)
153 : : real(dp),intent(out) :: enl_out(nband)
154 : : real(dp),target,intent(out) :: eig(nband)
155 : : !real(dp),target,intent(in) :: occ(nband) ! kept for reference
156 : : type(dataset_type),intent(in) :: dtset
157 : : type(gs_hamiltonian_type),target,intent(inout) :: gs_hamk
158 : :
159 : : ! Local variables-------------------------------
160 : : ! scalars
161 : : integer, parameter :: tim_slicewf = 2160
162 : : integer, parameter :: tim_nonlop = 1753
163 : : integer :: iband,shift,space,blockdim
164 : : integer :: spacedim,spacecom,gpu_option
165 : : integer :: me_g0,me_g0_fft
166 : : logical :: transfer_cg
167 : : !integer(kind=c_size_t) :: localMem
168 : 18 : type(slice_t) :: slice
169 : : type(xgBlock_t) :: xgx0,xgeigen,xgresidu
170 : : ! arrays
171 : : real(dp) :: tsec(2)
172 : : !integer(kind=c_size_t) :: sliceMem(2)
173 : 18 : real(dp), allocatable :: l_gvnlxc(:,:)
174 : :
175 : : ! Parameters for nonlop call in NC
176 : : integer,parameter :: choice=1, paw_opt=0, signs=1
177 : : real(dp) :: gsc_dummy(1,1)
178 [ + + + + : 2862 : type(pawcprj_type) :: cprj_dum(gs_hamk%natom,1)
+ + + + +
- + - ]
179 : :
180 : : ! *********************************************************************
181 : :
182 : : !################ INITIALIZATION #####################################
183 : : !######################################################################
184 : :
185 : 18 : call timab(tim_slicewf,1,tsec)
186 : :
187 : : !Set module variables
188 : 18 : l_paw = (gs_hamk%usepaw==1)
189 [ + - ]: 18 : l_cpopt=-1;l_sij_opt=0;if (l_paw) l_sij_opt=1
190 : 18 : l_npw = npw
191 : 18 : l_nspinor = nspinor
192 : 18 : l_prtvol = prtvol
193 : 18 : l_mpi_enreg => mpi_enreg
194 : 18 : l_gs_hamk => gs_hamk
195 : 18 : l_nband_filter = nband
196 : 18 : l_paral_kgb = dtset%paral_kgb
197 : 18 : l_block_sliced = dtset%invovl_blksliced
198 : 18 : l_invovl_precision = dtset%invovl_precision
199 : :
200 : : !Variables
201 : 18 : spacedim = l_npw*l_nspinor
202 : 18 : spacecom = l_mpi_enreg%comm_bandspinorfft
203 : 18 : gpu_option = dtset%gpu_option
204 : 18 : blockdim=l_mpi_enreg%nproc_band*l_mpi_enreg%bandpp
205 : : !for debug
206 : 18 : l_useria=dtset%useria
207 : :
208 : : !Depends on istwfk
209 [ - + ]: 18 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
210 : 0 : space = SPACE_CR
211 : 0 : l_icplx = 2
212 : : else ! complex
213 : 18 : space = SPACE_C
214 : 18 : l_icplx = 1
215 : : end if
216 : :
217 : 18 : me_g0 = -1
218 : 18 : me_g0_fft = -1
219 [ - + ]: 18 : if (space==SPACE_CR) then
220 : 0 : me_g0 = 0
221 : 0 : me_g0_fft = 0
222 [ # # ]: 0 : if (gs_hamk%istwf_k == 2) then
223 [ # # ]: 0 : if (l_mpi_enreg%me_g0 == 1) me_g0 = 1
224 [ # # ]: 0 : if (l_mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
225 : : end if
226 : : end if
227 : :
228 : 18 : transfer_cg = .false.
229 : : #ifdef HAVE_OPENMP_OFFLOAD
230 : : !$OMP TARGET ENTER DATA MAP(to:eig,resid) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
231 : : transfer_cg = .not. xomp_target_is_present(c_loc(cg))
232 : : !$OMP TARGET ENTER DATA MAP(to:cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. transfer_cg)
233 : : #endif
234 : :
235 : 18 : call xgBlock_map(xgx0,cg,space,spacedim,nband,comm=spacecom,me_g0=me_g0,gpu_option=gpu_option)
236 : :
237 : 18 : call xgBlock_map_1d(xgeigen,eig,SPACE_R,nband,gpu_option=gpu_option)
238 : :
239 : 18 : call xgBlock_map_1d(xgresidu,resid,SPACE_R,nband,gpu_option=gpu_option)
240 : :
241 : 18 : write(std_out,*) 'calling slice_init'
242 : 18 : write(std_out,*) 'dtset%paral_kgb=', dtset%paral_kgb
243 : 18 : write(std_out,*) 'spacecom=', spacecom, xmpi_comm_size(spacecom)
244 : 18 : write(std_out,*) 'spacedim=', spacedim
245 : :
246 : : call slice_init(slice,dtset%nslice,nband,spacedim,dtset%tolwfr_diago,dtset%paral_kgb,&
247 : : dtset%paral_slice,dtset%mdeg_filter,dtset%nbdbuf,0.d0,dtset%ecut,l_mpi_enreg%bandpp,&
248 : : space,spacecom,me_g0,me_g0_fft,l_paw,l_mpi_enreg%comm_spinorfft,l_mpi_enreg%comm_band,&
249 : : 0,l_gs_hamk%gpu_option,gpu_kokkos_nthrd=dtset%gpu_kokkos_nthrd,&
250 : 18 : gpu_thread_limit=dtset%gpu_thread_limit)
251 : :
252 : : ! todo will be cleaned and input variables will be deleted
253 : : ! for debug purposes
254 : : ! ==================================================
255 : : ! Manually setting slices from .abi
256 : : ! Slice 1
257 : : ! slice%neigenpairs_per_slice(1) = dtset%slice1_m ! number of TRUE eigenvalues in (poly_low, poly_upp)
258 : : ! slice%poly_degrees(1) = dtset%slice1_deg ! filter degree
259 : : ! slice%low_bounds(1) = dtset%slice1_part_a ! first slice is lowpass so interval to suppress
260 : : ! slice%upp_bounds(1) = dtset%slice1_part_b ! used for convergence <----
261 : : ! slice%low_bounds(1) = dtset%slice1_poly_a ! with overlap
262 : : ! slice%upp_bounds(1) = dtset%slice1_poly_b ! with overlap
263 : : !
264 : : ! ! Slice 2
265 : : ! slice%neigenpairs_per_slice(2) = dtset%slice2_m ! number of TRUE eigenvalues in (poly_low, poly_upp)
266 : : ! slice%poly_degrees(2) = dtset%slice2_deg ! slice%ndeg_filter
267 : : ! slice%low_bounds(2) = dtset%slice2_part_a ! used for convergence <-----
268 : : ! slice%upp_bounds(2) = dtset%slice2_part_b ! used for convergence <-----
269 : : ! slice%low_bounds(2) = dtset%slice2_poly_a ! with overlap
270 : : ! slice%upp_bounds(2) = dtset%slice2_poly_b ! with overlap
271 : : !
272 : : ! ! Slice 3
273 : : ! if (dtset%nslice==3) then
274 : : ! slice%neigenpairs_per_slice(3) = dtset%slice3_m ! number of TRUE eigenvalues in (poly_low, poly_upp)
275 : : ! slice%poly_degrees(3) = dtset%slice3_deg ! slice%ndeg_filter
276 : : ! slice%low_bounds(3) = dtset%slice3_part_a ! used for convergence <-----
277 : : ! slice%upp_bounds(3) = dtset%slice3_part_b ! used for convergence <-----
278 : : ! slice%low_bounds(3) = dtset%slice3_poly_a ! with overlap
279 : : ! slice%upp_bounds(3) = dtset%slice3_poly_b ! with overlap
280 : : !
281 : : ! end if
282 : : ! ==================================================
283 : :
284 : : !################ RUUUUUUUN #####################################
285 : : !######################################################################
286 : :
287 : 18 : call slice_run(slice, xgx0, getghc_gsc1, getBm1X, xgeigen, xgresidu, nspinor)
288 : :
289 [ - + ]: 18 : if ( .not. l_paw ) then
290 : 0 : call timab(tim_nonlop,1,tsec)
291 : : #ifdef FC_CRAY
292 : : ABI_MALLOC(l_gvnlxc,(1,1))
293 : : #else
294 [ # # # # ]: 0 : ABI_MALLOC(l_gvnlxc,(0,0))
295 : : #endif
296 : : !end if
297 : :
298 : : ABI_NVTX_START_RANGE(NVTX_SLICE_NONLOP)
299 : : !Call nonlop
300 [ # # ]: 0 : if (l_paral_kgb==0) then
301 : :
302 : : call nonlop(choice,l_cpopt,cprj_dum,enl_out,l_gs_hamk,0,eig,mpi_enreg,nband,1,paw_opt,&
303 : 0 : & signs,gsc_dummy,l_tim_getghc,cg,l_gvnlxc)
304 : :
305 : : else
306 : : #ifdef HAVE_OPENMP_OFFLOAD
307 : : !$OMP TARGET UPDATE FROM(cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
308 : : #endif
309 [ # # ]: 0 : do iband=1,nband/blockdim
310 : 0 : shift = (iband-1)*blockdim*l_npw*l_nspinor
311 : : call prep_nonlop(choice,l_cpopt,cprj_dum, &
312 : : & enl_out((iband-1)*blockdim+1:iband*blockdim),l_gs_hamk,0,&
313 : : & eig((iband-1)*blockdim+1:iband*blockdim),blockdim,mpi_enreg,1,paw_opt,signs,&
314 : : & gsc_dummy,l_tim_getghc, &
315 : : & cg(:,shift+1:shift+blockdim*l_npw*l_nspinor),&
316 : : !& l_gvnlxc(:,shift+1:shift+blockdim*l_npw*l_nspinor),&
317 : : & l_gvnlxc(:,:),&
318 : 0 : & already_transposed=.false.)
319 : : end do
320 : : end if
321 : : ABI_NVTX_END_RANGE()
322 [ # # ]: 0 : ABI_FREE(l_gvnlxc)
323 : 0 : call timab(tim_nonlop,2,tsec)
324 : : end if
325 : :
326 : : !Free slice
327 : 18 : call slice_free(slice)
328 : :
329 : : #ifdef HAVE_OPENMP_OFFLOAD
330 : : !$OMP TARGET UPDATE FROM(eig,resid) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
331 : : !$OMP TARGET EXIT DATA MAP(delete:eig,resid) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
332 : : !$OMP TARGET EXIT DATA MAP(from:cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. transfer_cg)
333 : : #endif
334 : :
335 : 18 : call timab(tim_slicewf,2,tsec)
336 : :
337 : : DBG_EXIT("COLL")
338 : :
339 [ - + - + : 36 : end subroutine slicewf
- + - + -
+ ]
340 : : !!***
341 : :
342 : : !----------------------------------------------------------------------
343 : :
344 : : !!****f* m_slicewf/getghc_gsc1
345 : : !! NAME
346 : : !! getghc_gsc1
347 : : !!
348 : : !! FUNCTION
349 : : !! This routine computes H|C> and possibly S|C> for a given wave function C.
350 : : !! It acts as a driver for getghc, taken into account parallelism, multithreading, etc.
351 : : !!
352 : : !! SIDE EFFECTS
353 : : !! X <type(xgBlock_t)>= memory block containing |C>
354 : : !! AX <type(xgBlock_t)>= memory block containing H|C>
355 : : !! BX <type(xgBlock_t)>= memory block containing S|C>
356 : : !!
357 : : !! SOURCE
358 : :
359 : 1278 : subroutine getghc_gsc1(X,AX,BX)
360 : :
361 : : implicit none
362 : :
363 : : !Arguments ------------------------------------
364 : : type(xgBlock_t), intent(inout) :: X
365 : : type(xgBlock_t), intent(inout) :: AX
366 : : type(xgBlock_t), intent(inout) :: BX
367 : : integer :: blockdim
368 : : integer :: spacedim
369 [ + + + + : 203202 : type(pawcprj_type) :: cprj_dum(l_gs_hamk%natom,1)
+ + + + +
- + - ]
370 : :
371 : : !Local variables-------------------------------
372 : : !scalars
373 : : real(dp) :: eval
374 : : !arrays
375 : 1278 : real(dp), pointer :: cg(:,:)
376 : 1278 : real(dp), pointer :: ghc(:,:)
377 : 1278 : real(dp), pointer :: gsc(:,:)
378 : : real(dp) :: l_gvnlxc(1,1)
379 : :
380 : : ! *********************************************************************
381 : :
382 : : ABI_NVTX_START_RANGE(NVTX_GETGHC)
383 : :
384 : 1278 : call xgBlock_getSize(X,spacedim,blockdim)
385 : 1278 : call xgBlock_check(X,AX)
386 : 1278 : call xgBlock_check(X,BX)
387 : :
388 : 1278 : call xgBlock_reverseMap(X,cg,rows=1,cols=spacedim*blockdim)
389 : 1278 : call xgBlock_reverseMap(AX,ghc,rows=1,cols=spacedim*blockdim)
390 : 1278 : call xgBlock_reverseMap(BX,gsc,rows=1,cols=spacedim*blockdim)
391 : :
392 : : call multithreaded_getghc(l_cpopt,cg,cprj_dum,ghc,gsc,&
393 : 1278 : l_gs_hamk,l_gvnlxc,eval,l_mpi_enreg,blockdim,l_prtvol,l_sij_opt,l_tim_getghc,0)
394 : :
395 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
396 : : call gpu_device_synchronize()
397 : : #endif
398 : :
399 [ - + ]: 1278 : if ( .not. l_paw ) call xgBlock_copy(X,BX)
400 : :
401 : : ABI_NVTX_END_RANGE()
402 : :
403 : 2556 : end subroutine getghc_gsc1
404 : : !!***
405 : :
406 : : !----------------------------------------------------------------------
407 : :
408 : : !!****f* m_slicewf/getBm1X
409 : : !! NAME
410 : : !! getBm1X
411 : : !!
412 : : !! FUNCTION
413 : : !! This routine computes S^-1|C> for a given wave function C.
414 : : !! It acts as a driver for apply_invovl.
415 : : !!
416 : : !! SIDE EFFECTS
417 : : !! X <type(xgBlock_t)>= memory block containing |C>
418 : : !! Bm1X <type(xgBlock_t)>= memory block containing S^-1|C>
419 : : !!
420 : : !! SOURCE
421 : :
422 : 720 : subroutine getBm1X(X,Bm1X)
423 : :
424 : : implicit none
425 : :
426 : : !Arguments ------------------------------------
427 : : type(xgBlock_t), intent(inout) :: X
428 : : type(xgBlock_t), intent(inout) :: Bm1X
429 : :
430 : : !Local variables-------------------------------
431 : : !scalars
432 : : integer :: blockdim
433 : : integer :: spacedim
434 : : #ifdef HAVE_OPENMP_OFFLOAD
435 : : character(len=500) :: msg
436 : : #endif
437 : : !arrays
438 : 720 : real(dp), pointer :: ghc_filter(:,:)
439 : 720 : real(dp), pointer :: gsm1hc_filter(:,:)
440 : 720 : type(pawcprj_type), allocatable :: cwaveprj_next(:,:) !dummy
441 : :
442 : : ! *********************************************************************
443 : :
444 : : ! working bandpp will be equal to blockdim
445 : 720 : call xgBlock_getSize(X,spacedim,blockdim)
446 : :
447 : : #ifdef HAVE_OPENMP_OFFLOAD
448 : : if (l_blockdim_evil > 0) then
449 : : if (blockdim /= l_blockdim_evil) then
450 : : write(msg,'(2(a,i0))')' you are about to apply_invovl with ndat=',blockdim,&
451 : : ' different from allocated ndat=',l_blockdim_evil
452 : : ABI_WARNING(msg)
453 : : end if
454 : : if (blockdim > l_blockdim_evil) then
455 : : ABI_ERROR('apply_invovl will not work with ndat larger than allocated')
456 : : end if
457 : : else
458 : : l_blockdim_evil = blockdim
459 : : end if
460 : : #endif
461 : :
462 [ + - ]: 720 : if(l_paw) then
463 : :
464 : 720 : call xgBlock_reverseMap(X,ghc_filter,rows=1,cols=spacedim*blockdim)
465 : 720 : call xgBlock_reverseMap(Bm1X,gsm1hc_filter,rows=1,cols=spacedim*blockdim)
466 : :
467 : : !cwaveprj_next is dummy
468 [ - + ]: 720 : if(gemm_nonlop_use_gemm) then
469 [ # # # # : 0 : ABI_MALLOC(cwaveprj_next, (1,1))
# # ]
470 : : else
471 [ - + + - : 48960 : ABI_MALLOC(cwaveprj_next, (l_gs_hamk%natom,l_nspinor*blockdim))
- + + - +
- + - - +
- + + + +
+ ]
472 : 720 : call pawcprj_alloc(cwaveprj_next,0,l_gs_hamk%dimcprj)
473 : : end if
474 : :
475 : : ABI_NVTX_START_RANGE(NVTX_INVOVL)
476 : : call apply_invovl(l_gs_hamk, ghc_filter(:,:), gsm1hc_filter(:,:), cwaveprj_next(:,:), &
477 [ - + - + ]: 720 : spacedim/l_nspinor, blockdim, l_mpi_enreg, l_nspinor, l_block_sliced, l_invovl_precision)
478 : : ABI_NVTX_END_RANGE()
479 : :
480 : 720 : call pawcprj_free(cwaveprj_next)
481 [ + - + + : 45360 : ABI_FREE(cwaveprj_next)
- + - + -
+ ]
482 : :
483 : : else
484 : :
485 : 0 : call xgBlock_copy(X,Bm1X)
486 : :
487 : : end if
488 : :
489 [ - + - - : 1440 : end subroutine getBm1X
- - - - ]
490 : : !!***
491 : :
492 : : end module m_slicewf
493 : : !!***
|
