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