Branch data 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 : : !!***
|
