Line data Source code
1 : !!****f* ABINIT/m_chebfiwf
2 : !! NAME
3 : !! m_chebfiwf
4 : !!
5 : !! FUNCTION
6 : !! This module contains a routine updating the whole wave functions at a given k-point,
7 : !! using the Chebyshev filtering method (2021 implementation using xG abstraction layer)
8 : !! for a given spin-polarization, from a fixed hamiltonian
9 : !! but might also simply compute eigenvectors and eigenvalues at this k point.
10 : !! it will also update the matrix elements of the hamiltonian.
11 : !!
12 : !! COPYRIGHT
13 : !! Copyright (C) 2018-2026 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_chebfiwf
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_chebfi
40 : use m_chebfi2
41 : use m_invovl
42 :
43 : use m_cgtools, only : dotprod_g
44 : use m_dtset, only : dataset_type
45 :
46 : use m_hamiltonian, only : gs_hamiltonian_type
47 : use m_pawcprj, only : pawcprj_type
48 : use m_nonlop, only : nonlop
49 : use m_prep_kgb, only : prep_getghc, prep_nonlop
50 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
51 : use m_getghc, only : multithreaded_getghc
52 : use m_gemm_nonlop_projectors , only : gemm_nonlop_use_gemm
53 :
54 : use m_xg
55 : use m_xgTransposer
56 :
57 : #if defined(HAVE_GPU_MARKERS)
58 : use m_nvtx_data
59 : #endif
60 :
61 : #if defined(HAVE_GPU)
62 : use m_gpu_toolbox
63 : #endif
64 :
65 : #if defined(HAVE_YAKL)
66 : use gator_mod
67 : #endif
68 :
69 : use, intrinsic :: iso_c_binding, only: c_associated,c_loc,c_ptr,c_f_pointer,c_double,c_size_t
70 :
71 : use m_xmpi
72 : use m_xomp
73 : #ifdef HAVE_OPENMP
74 : use omp_lib
75 : #endif
76 :
77 : implicit none
78 :
79 : private
80 :
81 : integer, parameter :: l_tim_getghc=7
82 : real(dp), parameter :: inv_sqrt2 = 1/sqrt2
83 :
84 : ! For use in getghc_gsc1
85 : integer, save :: l_cpopt
86 : integer, save :: l_icplx
87 : integer, save :: l_nspinor
88 : logical, save :: l_paw
89 : integer, save :: l_prtvol
90 : integer, save :: l_sij_opt
91 : integer, save :: l_paral_kgb
92 : integer, save :: l_useria
93 : integer, save :: l_block_sliced
94 :
95 : type(mpi_type),pointer,save :: l_mpi_enreg
96 : type(gs_hamiltonian_type),pointer,save :: l_gs_hamk
97 :
98 : integer, parameter :: DEBUG_ROWS = 5
99 : integer, parameter :: DEBUG_COLUMNS = 5
100 :
101 : public :: chebfiwf2
102 :
103 : CONTAINS !========================================================================================
104 : !!***
105 :
106 : !!****f* m_chebfiwf/chebfiwf2
107 : !! NAME
108 : !! chebfiwf2
109 : !!
110 : !! FUNCTION
111 : !! This routine updates the whole wave functions set at a given k-point,
112 : !! using the Chebfi method (2021 version using xG abstraction layer)
113 : !!
114 : !! INPUTS
115 : !! dtset= input variables for this dataset
116 : !! mpi_enreg= MPI-parallelisation information
117 : !! nband= number of bands at this k point
118 : !! npw= number of plane waves at this k point
119 : !! nspinor= number of spinorial components of the wavefunctions
120 : !! prtvol= control print volume and debugging
121 : !!
122 : !! OUTPUT
123 : !! eig(nband)= eigenvalues (hartree) for all bands
124 : !! enl_out(nband)= contribution of each band to the nl part of energy
125 : !! resid(nband)= residuals for each band
126 : !!
127 : !! SIDE EFFECTS
128 : !! cg(2,npw*nspinor*nband)= planewave coefficients of wavefunctions
129 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
130 : !!
131 : !! SOURCE
132 :
133 9024 : subroutine chebfiwf2(cg,dtset,eig,occ,enl_out,gs_hamk,mpi_enreg,&
134 9024 : & nband,npw,nspinor,prtvol,resid)
135 :
136 : implicit none
137 :
138 : ! Arguments ------------------------------------
139 : integer,intent(in) :: nband,npw,prtvol,nspinor
140 : type(mpi_type),target,intent(in) :: mpi_enreg
141 : real(dp),target,intent(inout) :: cg(2,npw*nspinor*nband)
142 : real(dp),target,intent(out) :: resid(nband)
143 : real(dp),intent(out) :: enl_out(nband)
144 : real(dp),target,intent(out) :: eig(nband)
145 : real(dp),target,intent(in) :: occ(nband)
146 : type(dataset_type),intent(in) :: dtset
147 : type(gs_hamiltonian_type),target,intent(inout) :: gs_hamk
148 :
149 : ! Local variables-------------------------------
150 : ! scalars
151 : integer, parameter :: tim_chebfiwf2 = 1750
152 : integer, parameter :: tim_nonlop = 1753
153 : integer :: iband,shift,space,blockdim,total_spacedim,ierr
154 : integer :: me_g0,me_g0_fft
155 : logical :: transfer_cg
156 : integer(kind=c_size_t) :: localMem
157 9024 : type(chebfi_t) :: chebfi
158 : type(xgBlock_t) :: xgx0,xgeigen,xgocc,xgresidu
159 : ! arrays
160 : real(dp) :: tsec(2)
161 : integer(kind=c_size_t) :: chebfiMem(2)
162 9024 : real(dp), allocatable :: gvnlxc(:,:),occ_tmp(:)
163 :
164 : ! Parameters for nonlop call in NC
165 : integer,parameter :: choice=1, paw_opt=0, signs=1
166 : real(dp) :: gsc_dummy(1,1)
167 126976 : type(pawcprj_type) :: cprj_dum(gs_hamk%natom,1)
168 :
169 : ! *********************************************************************
170 :
171 : !################ INITIALIZATION #####################################
172 : !######################################################################
173 :
174 9024 : call timab(tim_chebfiwf2,1,tsec)
175 :
176 : !Set module variables
177 9024 : l_paw = (gs_hamk%usepaw==1)
178 9024 : l_cpopt=-1;l_sij_opt=0;if (l_paw) l_sij_opt=1
179 9024 : l_nspinor = nspinor
180 9024 : l_prtvol = prtvol
181 9024 : l_mpi_enreg => mpi_enreg
182 9024 : l_gs_hamk => gs_hamk
183 9024 : l_paral_kgb = dtset%paral_kgb
184 9024 : l_block_sliced = dtset%invovl_blksliced
185 :
186 : !Variables
187 9024 : blockdim=l_mpi_enreg%nproc_band*l_mpi_enreg%bandpp
188 : !for debug
189 9024 : l_useria=dtset%useria
190 :
191 : !Depends on istwfk
192 9024 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
193 : ! SPACE_CR mean that we have complex numbers but no re*im terms only re*re
194 : ! and im*im so that a vector of complex is consider as a long vector of real
195 : ! therefore the number of data is (2*npw*nspinor)*nband
196 : ! This space is completely equivalent to SPACE_R but will correctly set and
197 : ! get the array data into the xgBlock
198 3536 : space = SPACE_CR
199 3536 : l_icplx = 2
200 : else ! complex
201 5488 : space = SPACE_C
202 5488 : l_icplx = 1
203 : end if
204 :
205 9024 : me_g0 = -1
206 9024 : me_g0_fft = -1
207 9024 : if (space==SPACE_CR) then
208 3536 : me_g0 = 0
209 3536 : me_g0_fft = 0
210 3536 : if (gs_hamk%istwf_k == 2) then
211 442 : if (l_mpi_enreg%me_g0 == 1) me_g0 = 1
212 442 : if (l_mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
213 : end if
214 : end if
215 :
216 : !Memory info
217 9024 : if ( prtvol >= 3 ) then
218 0 : if (l_mpi_enreg%paral_kgb == 1) then
219 0 : total_spacedim = l_icplx*npw*nspinor
220 0 : call xmpi_sum(total_spacedim,l_mpi_enreg%comm_bandspinorfft,ierr)
221 : else
222 0 : total_spacedim = 0
223 : end if
224 : chebfiMem = chebfi_memInfo(nband,l_icplx*npw*nspinor,space,l_mpi_enreg%paral_kgb, &
225 0 : & total_spacedim,l_mpi_enreg%bandpp) !blockdim
226 0 : localMem = (int(2,c_size_t)*npw*nspinor*nband+3*nband)*kind(1.d0) !blockdim
227 0 : write(std_out,'(1x,A,F10.6,1x,A)') "Each MPI process calling chebfi should need around ", &
228 0 : (localMem+sum(chebfiMem))/1e9,"GB of peak memory as follows :"
229 0 : write(std_out,'(4x,A,F10.6,1x,A)') "Permanent memory in chebfiwf : ",real(localMem)/1e9,"GB"
230 0 : write(std_out,'(4x,A,F10.6,1x,A)') "Permanent memory in m_chebfi : ",real(chebfiMem(1))/1e9,"GB"
231 0 : write(std_out,'(4x,A,F10.6,1x,A)') "Temporary memory in m_chebfi : ",real(chebfiMem(2))/1e9,"GB"
232 : end if
233 :
234 9024 : transfer_cg = .false.
235 : #ifdef HAVE_OPENMP_OFFLOAD
236 : !$OMP TARGET ENTER DATA MAP(to:eig,resid,occ) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
237 : transfer_cg = .not. xomp_target_is_present(c_loc(cg))
238 : !$OMP TARGET ENTER DATA MAP(to:cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. transfer_cg)
239 : #endif
240 :
241 : call xgBlock_map(xgx0,cg,space,npw*nspinor,nband,comm=l_mpi_enreg%comm_bandspinorfft,me_g0=me_g0,&
242 9024 : & gpu_option=dtset%gpu_option)
243 :
244 9024 : call xgBlock_map_1d(xgeigen,eig,SPACE_R,nband,gpu_option=dtset%gpu_option)
245 :
246 9024 : call xgBlock_map_1d(xgresidu,resid,SPACE_R,nband,gpu_option=dtset%gpu_option)
247 :
248 : ! Occupancies in chebyshev are used for convergence criteria only
249 9024 : if (dtset%nbdbuf==-101.and.nspinor==1.and.dtset%nsppol==1) then
250 2448 : ABI_MALLOC(occ_tmp,(nband))
251 7344 : occ_tmp(:) = half*occ(:)
252 816 : call xgBlock_map_1d(xgocc,occ_tmp,SPACE_R,nband,gpu_option=dtset%gpu_option)
253 : else
254 8208 : call xgBlock_map_1d(xgocc,occ,SPACE_R,nband,gpu_option=dtset%gpu_option)
255 : end if
256 :
257 9024 : call timab(tim_chebfiwf2,2,tsec) ! IL TODO this should be deactivated
258 :
259 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_INIT)
260 : call chebfi_init(chebfi,nband,npw*nspinor,dtset%tolwfr_diago,dtset%ecut, &
261 : & dtset%paral_kgb,l_mpi_enreg%bandpp, &
262 : & dtset%mdeg_filter, dtset%nbdbuf, space,1, &
263 : & l_mpi_enreg%comm_bandspinorfft,me_g0,me_g0_fft,l_paw,&
264 : & l_mpi_enreg%comm_spinorfft,l_mpi_enreg%comm_band,&
265 : & dtset%chebfi_oracle,dtset%oracle_factor,dtset%oracle_min_occ,&
266 : & l_gs_hamk%gpu_option,gpu_kokkos_nthrd=dtset%gpu_kokkos_nthrd,&
267 9024 : & gpu_thread_limit=dtset%gpu_thread_limit)
268 : ABI_NVTX_END_RANGE()
269 :
270 : !################ RUUUUUUUN #####################################
271 : !######################################################################
272 :
273 9024 : call chebfi_run(chebfi,xgx0,getghc_gsc1,getBm1X,xgeigen,xgocc,xgresidu,nspinor)
274 :
275 9024 : if (allocated(occ_tmp)) then
276 816 : ABI_FREE(occ_tmp)
277 : end if
278 :
279 9024 : if ( .not. l_paw ) then
280 1496 : call timab(tim_nonlop,1,tsec)
281 : #ifdef FC_CRAY
282 : ABI_MALLOC(gvnlxc,(1,1))
283 : #else
284 1496 : ABI_MALLOC(gvnlxc,(0,0))
285 : #endif
286 : !end if
287 :
288 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NONLOP)
289 : !Call nonlop
290 1496 : if (l_paral_kgb==0) then
291 :
292 : call nonlop(choice,l_cpopt,cprj_dum,enl_out,l_gs_hamk,0,eig,mpi_enreg,nband,1,paw_opt,&
293 408 : & signs,gsc_dummy,l_tim_getghc,cg,gvnlxc)
294 :
295 : else
296 : #ifdef HAVE_OPENMP_OFFLOAD
297 : !$OMP TARGET UPDATE FROM(cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
298 : #endif
299 2176 : do iband=1,nband/blockdim
300 1088 : shift = (iband-1)*blockdim*npw*nspinor
301 : call prep_nonlop(choice,l_cpopt,cprj_dum, &
302 : & enl_out((iband-1)*blockdim+1:iband*blockdim),l_gs_hamk,0,&
303 : & eig((iband-1)*blockdim+1:iband*blockdim),blockdim,mpi_enreg,1,paw_opt,signs,&
304 : & gsc_dummy,l_tim_getghc,cg(:,shift+1:shift+blockdim*npw*nspinor),gvnlxc(:,:),&
305 2176 : & already_transposed=.false.)
306 : end do
307 : end if
308 : ABI_NVTX_END_RANGE()
309 1496 : ABI_FREE(gvnlxc)
310 1496 : call timab(tim_nonlop,2,tsec)
311 : end if
312 :
313 : !Free chebfi
314 9024 : call chebfi_free(chebfi)
315 :
316 : #ifdef HAVE_OPENMP_OFFLOAD
317 : !$OMP TARGET UPDATE FROM(eig,resid) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
318 : !$OMP TARGET EXIT DATA MAP(delete:eig,resid,occ) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
319 : !$OMP TARGET EXIT DATA MAP(from:cg) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. transfer_cg)
320 : #endif
321 :
322 9024 : call timab(tim_chebfiwf2,2,tsec)
323 :
324 : DBG_EXIT("COLL")
325 :
326 27072 : end subroutine chebfiwf2
327 : !!***
328 :
329 : !----------------------------------------------------------------------
330 :
331 : !!****f* m_chebfiwf/getghc_gsc1
332 : !! NAME
333 : !! getghc_gsc1
334 : !!
335 : !! FUNCTION
336 : !! This routine computes H|C> and possibly S|C> for a given wave function C.
337 : !! It acts as a driver for getghc, taken into account parallelism, multithreading, etc.
338 : !!
339 : !! SIDE EFFECTS
340 : !! X <type(xgBlock_t)>= memory block containing |C>
341 : !! AX <type(xgBlock_t)>= memory block containing H|C>
342 : !! BX <type(xgBlock_t)>= memory block containing S|C>
343 : !!
344 : !! SOURCE
345 :
346 72586 : subroutine getghc_gsc1(X,AX,BX)
347 :
348 : implicit none
349 :
350 : !Arguments ------------------------------------
351 : type(xgBlock_t), intent(inout) :: X
352 : type(xgBlock_t), intent(inout) :: AX
353 : type(xgBlock_t), intent(inout) :: BX
354 : integer :: blockdim
355 : integer :: spacedim
356 1020684 : type(pawcprj_type) :: cprj_dum(l_gs_hamk%natom,1)
357 :
358 : !Local variables-------------------------------
359 : !scalars
360 : real(dp) :: eval
361 : !arrays
362 72586 : real(dp), pointer :: cg(:,:)
363 72586 : real(dp), pointer :: ghc(:,:)
364 72586 : real(dp), pointer :: gsc(:,:)
365 : real(dp) :: gvnlxc(1,1)
366 :
367 : ! *********************************************************************
368 :
369 : ABI_NVTX_START_RANGE(NVTX_GETGHC)
370 :
371 72586 : call xgBlock_getSize(X,spacedim,blockdim)
372 72586 : call xgBlock_check(X,AX)
373 72586 : call xgBlock_check(X,BX)
374 :
375 72586 : call xgBlock_reverseMap(X,cg,rows=1,cols=spacedim*blockdim)
376 72586 : call xgBlock_reverseMap(AX,ghc,rows=1,cols=spacedim*blockdim)
377 72586 : call xgBlock_reverseMap(BX,gsc,rows=1,cols=spacedim*blockdim)
378 :
379 : call multithreaded_getghc(l_cpopt,cg,cprj_dum,ghc,gsc,&
380 72586 : l_gs_hamk,gvnlxc,eval,l_mpi_enreg,blockdim,l_prtvol,l_sij_opt,l_tim_getghc,0)
381 :
382 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
383 : call gpu_device_synchronize()
384 : #endif
385 :
386 72586 : if ( .not. l_paw ) call xgBlock_copy(X,BX)
387 :
388 : ABI_NVTX_END_RANGE()
389 :
390 145172 : end subroutine getghc_gsc1
391 : !!***
392 :
393 : !----------------------------------------------------------------------
394 :
395 : !!****f* m_chebfiwf/getBm1X
396 : !! NAME
397 : !! getBm1X
398 : !!
399 : !! FUNCTION
400 : !! This routine computes S^-1|C> for a given wave function C.
401 : !! It acts as a driver for apply_invovl.
402 : !!
403 : !! SIDE EFFECTS
404 : !! X <type(xgBlock_t)>= memory block containing |C>
405 : !! Bm1X <type(xgBlock_t)>= memory block containing S^-1|C>
406 : !!
407 : !! SOURCE
408 :
409 51718 : subroutine getBm1X(X,Bm1X)
410 :
411 : implicit none
412 :
413 : !Arguments ------------------------------------
414 : type(xgBlock_t), intent(inout) :: X
415 : type(xgBlock_t), intent(inout) :: Bm1X
416 :
417 : !Local variables-------------------------------
418 : !scalars
419 : integer :: blockdim
420 : integer :: spacedim
421 : !arrays
422 51718 : real(dp), pointer :: ghc_filter(:,:)
423 51718 : real(dp), pointer :: gsm1hc_filter(:,:)
424 51718 : type(pawcprj_type), allocatable :: cwaveprj_next(:,:) !dummy
425 :
426 : ! *********************************************************************
427 :
428 51718 : call xgBlock_getSize(X,spacedim,blockdim)
429 :
430 51718 : if(l_paw) then
431 :
432 51718 : call xgBlock_reverseMap(X,ghc_filter,rows=1,cols=spacedim*blockdim)
433 51718 : call xgBlock_reverseMap(Bm1X,gsm1hc_filter,rows=1,cols=spacedim*blockdim)
434 :
435 : !cwaveprj_next is dummy
436 51718 : if(gemm_nonlop_use_gemm) then
437 1152 : ABI_MALLOC(cwaveprj_next, (1,1))
438 : else
439 1247656 : ABI_MALLOC(cwaveprj_next, (l_gs_hamk%natom,l_nspinor*blockdim))
440 51430 : call pawcprj_alloc(cwaveprj_next,0,l_gs_hamk%dimcprj)
441 : end if
442 :
443 : ABI_NVTX_START_RANGE(NVTX_INVOVL)
444 : call apply_invovl(l_gs_hamk, ghc_filter(:,:), gsm1hc_filter(:,:), cwaveprj_next(:,:), &
445 51718 : spacedim/l_nspinor, blockdim, l_mpi_enreg, l_nspinor, l_block_sliced)
446 : ABI_NVTX_END_RANGE()
447 :
448 51718 : call pawcprj_free(cwaveprj_next)
449 761990 : ABI_FREE(cwaveprj_next)
450 :
451 : else
452 :
453 0 : call xgBlock_copy(X,Bm1X)
454 :
455 : end if
456 :
457 103436 : end subroutine getBm1X
458 : !!***
459 :
460 : end module m_chebfiwf
461 : !!***
|
