Line data Source code
1 : !!****m* ABINIT/m_mpi_setup
2 : !! NAME
3 : !! m_mpi_setup
4 : !!
5 : !! FUNCTION
6 : !! Initialize MPI parameters and datastructures for parallel execution
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (FJ, MT, FD)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_mpi_setup
23 :
24 : use defs_basis
25 : use m_distribfft
26 : use m_xmpi
27 : use m_xomp
28 : use m_hdr
29 : use m_sort
30 : use m_errors
31 : use m_abicore
32 :
33 : use defs_abitypes, only : MPI_type
34 : use m_fstrings, only : sjoin, itoa
35 : use m_time, only : abi_wtime
36 : use m_parser, only : intagm
37 : use m_geometry, only : mkrdim, metric
38 : use m_fftcore, only : fftalg_for_npfft, getng, kpgcount
39 : use m_mpinfo, only : init_mpi_enreg, mpi_distrib_is_ok, initmpi_atom, proc_distrb_cycle, &
40 : initmpi_grid, initmpi_pert, initmpi_img, distrb2, distrb2_hf, initmpi_world
41 : use m_libpaw_tools, only : libpaw_write_comm_set
42 : use m_dtset, only : dataset_type
43 : use m_kg, only : getmpw
44 : use m_dtfil, only : mkfilename
45 : use m_mep, only : NEB_CELL_ALGO_NONE
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : public :: mpi_setup
53 : !!***
54 :
55 : contains
56 : !!***
57 :
58 : !!****f* ABINIT/mpi_setup
59 : !! NAME
60 : !! mpi_setup
61 : !!
62 : !! FUNCTION
63 : !! Big loop on the datasets:
64 : !! - compute mgfft,mpw,nfft,... for this data set;
65 : !! - fill mpi_enreg
66 : !! *** At the output of this routine, all the dtsets input variables are known ***
67 : !! The content of dtsets should not be modified anymore afterwards.
68 : !!
69 : !! INPUTS
70 : !! filnam(5)=character strings giving file names
71 : !! ndtset= number of datasets to be read; if 0, no multi-dataset mode
72 : !! ndtset_alloc=number of datasets, corrected for allocation of at least
73 : !! one data set.
74 : !!
75 : !! OUTPUT
76 : !! dtsets(0:ndtset_alloc)=<type datafiles_type>contains all input variables,
77 : !! some of which are initialized here, while other were already
78 : !! initialized previously.
79 : !!
80 : !! SIDE EFFECTS
81 : !! mpi_enregs=information about MPI parallelization
82 : !!
83 : !! SOURCE
84 :
85 1428 : subroutine mpi_setup(dtsets,filnam,lenstr,mpi_enregs,ndtset,ndtset_alloc,string)
86 :
87 : !Arguments ------------------------------------
88 : !scalars
89 : integer,intent(in) :: lenstr,ndtset,ndtset_alloc
90 : type(MPI_type),intent(inout) :: mpi_enregs(0:ndtset_alloc)
91 : character(len=*),intent(in) :: string
92 : !arrays
93 : character(len=fnlen),intent(in) :: filnam(5)
94 : type(dataset_type),intent(inout) :: dtsets(0:ndtset_alloc)
95 :
96 : !Local variables -------------------------------
97 : !scalars
98 : integer :: blocksize,exchn2n3d,iband,idtset,iexit,ii,iikpt,iikpt_modulo, prtvol
99 : integer :: isppol,jdtset,marr,mband_lower,mband_upper
100 : integer :: me_fft,mgfft,mgfftdg,mkmem,mpw,mpw_k,max_mpw,optdriver
101 : integer :: mband_mem
102 : integer :: nfft,nfftdg,nkpt,nkpt_me,npert,nproc,nproc_fft,nqpt
103 : integer :: nspink,nsppol,nsym,nthreads,paral_fft,response,tnband,tread0,usepaw,vectsize
104 : integer :: fftalg,fftalga,fftalgc
105 : #ifdef HAVE_LINALG_ELPA
106 : integer :: icol,irow,np
107 : #endif
108 : logical :: fftalg_read,forbid_threads,ortalg_read,paral_kgb_read,wfoptalg_read,do_check
109 : real(dp) :: dilatmx,ecut,ecut_eff,ecutdg_eff,ucvol
110 : character(len=500) :: msg
111 : !arrays
112 : integer :: ngfft(18),ngfftdg(18),ngfftc(3),tread(12)
113 1428 : integer,allocatable :: intarr(:),istwfk(:),symrel(:,:,:)
114 1428 : integer,allocatable :: mybands(:)
115 1428 : integer,pointer :: nkpt_rbz(:)
116 : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
117 : real(dp) :: gmet(3,3),gprimd(3,3),kpt(3),qphon(3),rmet(3,3),rprimd(3,3)
118 1428 : real(dp),allocatable :: dprarr(:),kpt_with_shift(:,:)
119 1428 : real(dp),pointer :: nband_rbz(:,:)
120 : character(len=6) :: nm_mkmem(3)
121 :
122 : !*************************************************************************
123 :
124 : DBG_ENTER("COLL")
125 :
126 1428 : iexit=0;mpw_k=0
127 :
128 1428 : call init_mpi_enreg(mpi_enregs(0))
129 1428 : call initmpi_img(dtsets(0),mpi_enregs(0),-1)
130 1428 : nthreads=xomp_get_num_threads(open_parallel=.True.)
131 :
132 8387 : do idtset=1,ndtset_alloc
133 6959 : call init_mpi_enreg(mpi_enregs(idtset))
134 :
135 : ! Handy read-only variables.
136 6959 : optdriver = dtsets(idtset)%optdriver
137 6959 : prtvol = dtsets(idtset)%prtvol
138 :
139 : ! Read parallel input parameters
140 6959 : marr=max(12,dtsets(idtset)%npsp,dtsets(idtset)%nimage)
141 20877 : ABI_MALLOC(intarr,(marr))
142 20877 : ABI_MALLOC(dprarr,(marr))
143 6959 : nkpt =dtsets(idtset)%nkpt
144 6959 : nsppol=dtsets(idtset)%nsppol
145 6959 : jdtset=dtsets(idtset)%jdtset ; if(ndtset==0)jdtset=0
146 6959 : usepaw=dtsets(idtset)%usepaw
147 124359 : mband_upper=maxval(dtsets(idtset)%nband(1:nkpt*nsppol))
148 : mband_lower=minval(dtsets(idtset)%nband(1:nkpt*nsppol))
149 :
150 : ! Compute metric for this dataset
151 6959 : call mkrdim(dtsets(idtset)%acell_orig(1:3,1),dtsets(idtset)%rprim_orig(1:3,1:3,1),rprimd)
152 6959 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
153 :
154 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'max_ncpus',tread0,'INT')
155 6959 : if (tread0==1) dtsets(idtset)%max_ncpus=intarr(1)
156 :
157 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'paral_atom',tread0,'INT')
158 6959 : if(tread0==1) dtsets(idtset)%paral_atom=intarr(1)
159 :
160 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'paral_rf',tread0,'INT')
161 6959 : if (tread0==1.and.any(optdriver==[RUNL_RESPFN, RUNL_NONLINEAR])) dtsets(idtset)%paral_rf=intarr(1)
162 :
163 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'npimage',tread(2),'INT')
164 6959 : if(tread(2)==1) dtsets(idtset)%npimage=intarr(1)
165 :
166 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nppert',tread(3),'INT')
167 6959 : if (tread(3)==1.and.optdriver==RUNL_RESPFN) dtsets(idtset)%nppert=intarr(1)
168 :
169 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'np_spkpt',tread(4),'INT')
170 6959 : if(tread(4)==1)then
171 254 : dtsets(idtset)%np_spkpt=intarr(1)
172 : else
173 : ! npkpt is obsolete, but still read
174 6705 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'npkpt',tread(4),'INT')
175 6705 : if(tread(4)==1)then
176 72 : dtsets(idtset)%np_spkpt=intarr(1)
177 : endif
178 : endif
179 :
180 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'npspinor',tread(5),'INT')
181 6959 : if(tread(5)==1) dtsets(idtset)%npspinor=intarr(1)
182 :
183 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'npfft',tread(6),'INT')
184 6959 : if(tread(6)==1) dtsets(idtset)%npfft=intarr(1)
185 :
186 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'npband',tread(7),'INT')
187 6959 : if(tread(7)==1) dtsets(idtset)%npband=intarr(1)
188 :
189 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'bandpp',tread(8),'INT')
190 6959 : if(tread(8)==1) dtsets(idtset)%bandpp=intarr(1)
191 :
192 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'use_slk',tread(9),'INT')
193 6959 : if(tread(9)==1) dtsets(idtset)%use_slk=intarr(1)
194 :
195 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'np_slk',tread(10),'INT')
196 6959 : if(tread(10)==1) dtsets(idtset)%np_slk=intarr(1)
197 :
198 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slk_rankpp',tread(12),'INT')
199 6959 : if(tread(12)==1) dtsets(idtset)%slk_rankpp=intarr(1)
200 :
201 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pw_unbal_thresh',tread0,'DPR')
202 6959 : if(tread0==1) dtsets(idtset)%pw_unbal_thresh=dprarr(1)
203 6959 : mpi_enregs(idtset)%pw_unbal_thresh=dtsets(idtset)%pw_unbal_thresh
204 :
205 6959 : call intagm(dprarr,intarr,jdtset,marr,12,string(1:lenstr),'gpu_devices',tread0,'INT')
206 6959 : if(tread0==1) dtsets(idtset)%gpu_devices(1:12)=intarr(1:12)
207 :
208 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_nfft_blocks',tread0,'INT')
209 6959 : if(tread0==1) dtsets(idtset)%gpu_nfft_blocks=intarr(1)
210 :
211 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_kokkos_nthrd',tread0,'INT')
212 6959 : if(tread0==1) dtsets(idtset)%gpu_kokkos_nthrd=intarr(1)
213 :
214 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_linalg_limit',tread(11),'INT')
215 6959 : if(tread(11)==1) dtsets(idtset)%gpu_linalg_limit=intarr(1)
216 :
217 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_nl_distrib',tread0,'INT')
218 6959 : if(tread0==1) dtsets(idtset)%gpu_nl_distrib=intarr(1)
219 :
220 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_nl_splitsize',tread0,'INT')
221 6959 : if(tread0==1) dtsets(idtset)%gpu_nl_splitsize=intarr(1)
222 :
223 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gpu_thread_limit',tread0,'INT')
224 6959 : if(tread0==1) dtsets(idtset)%gpu_thread_limit=intarr(1)
225 :
226 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nphf',tread0,'INT')
227 6959 : if(tread0==1) dtsets(idtset)%nphf=intarr(1)
228 :
229 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'autoparal',tread0,'INT')
230 6959 : if(tread0==1) dtsets(idtset)%autoparal=intarr(1)
231 :
232 : ! Read paral_kgb and disable it if not supported in optdriver
233 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'paral_kgb',tread(1),'INT')
234 6959 : paral_kgb_read=(tread(1)==1)
235 6959 : if (paral_kgb_read) dtsets(idtset)%paral_kgb=intarr(1)
236 : if (xmpi_paral==0.and.dtsets(idtset)%paral_kgb==1) then
237 : dtsets(idtset)%paral_kgb=0
238 : write(msg, '(5a)' ) &
239 : 'When ABINIT is compiled without MPI flag,',ch10,&
240 : 'setting paral_kgb/=0 is useless. paral_kgb has been reset to 0.',ch10,&
241 : 'Action: modify compilation option or paral_kgb in the input file.'
242 : ABI_WARNING(msg)
243 : end if
244 6959 : if (ALL(optdriver /= [RUNL_GSTATE, RUNL_GWLS, RUNL_RTTDDFT]) .and. dtsets(idtset)%paral_kgb/=0) then
245 20 : dtsets(idtset)%paral_kgb=0
246 : write(msg, '(a,i0,a)') &
247 20 : "paral_kgb != 0 is not available in optdriver ",optdriver,". Setting paral_kgb to 0"
248 20 : ABI_COMMENT(msg)
249 : end if
250 :
251 6959 : wfoptalg_read=.false.
252 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wfoptalg',tread0,'INT')
253 6959 : if(tread0==1) then
254 353 : dtsets(idtset)%wfoptalg=intarr(1)
255 353 : wfoptalg_read=.true.
256 : else
257 6606 : if (dtsets(idtset)%usepaw==0) dtsets(idtset)%wfoptalg=0
258 6606 : if (dtsets(idtset)%usepaw/=0) dtsets(idtset)%wfoptalg=10
259 6606 : if ((optdriver==RUNL_GSTATE.or.optdriver==RUNL_GWLS).and.dtsets(idtset)%paral_kgb/=0) dtsets(idtset)%wfoptalg=114
260 : end if
261 :
262 : ! Dump the list of irreducible perturbations and exit.
263 6959 : if (dtsets(idtset)%paral_rf==-1.and.optdriver/=RUNL_NONLINEAR) then
264 0 : call dtsets(idtset)%get_npert_rbz(nband_rbz, nkpt_rbz, npert)
265 0 : ABI_FREE(nband_rbz)
266 0 : ABI_FREE(nkpt_rbz)
267 0 : iexit = iexit + 1
268 : end if
269 :
270 : ! From total number of procs, compute all possible distributions
271 : ! Ignore exit flag if GW/EPH calculations because autoparal section is performed in screening/sigma/bethe_salpeter/eph
272 46349 : if (any(optdriver == [RUNL_SCREENING, RUNL_SIGMA, RUNL_BSE, RUNL_EPH, RUNL_GWR, RUNL_NONLINEAR])) then
273 515 : iexit = 0
274 : else
275 6444 : call finddistrproc(dtsets,filnam,idtset,iexit,mband_upper,mpi_enregs(idtset),ndtset_alloc,tread)
276 : end if
277 :
278 6959 : call initmpi_img(dtsets(idtset),mpi_enregs(idtset),-1)
279 6959 : nproc=mpi_enregs(idtset)%nproc_cell
280 :
281 : ! Set paral_kgb to 1 when band-fft parallelism is activated
282 6959 : if (ANY(optdriver == [RUNL_GSTATE, RUNL_GWLS, RUNL_RTTDDFT])) then
283 5699 : if (mpi_enregs(idtset)%nproc_cell>1) then
284 994 : if (dtsets(idtset)%npband>1.or.dtsets(idtset)%npfft>1) then
285 626 : if (.not.paral_kgb_read) dtsets(idtset)%paral_kgb=1
286 : end if
287 : end if
288 : end if
289 :
290 6959 : if ((optdriver/=RUNL_GSTATE.and.optdriver/=RUNL_GWLS.and.optdriver/=RUNL_RTTDDFT).and. &
291 : & (dtsets(idtset)%np_spkpt/=1 .or.dtsets(idtset)%npband/=1.or.dtsets(idtset)%npfft/=1.or. &
292 : & dtsets(idtset)%npspinor/=1.or.dtsets(idtset)%bandpp/=1)) then
293 : !& .or.(dtsets(idtset)%iscf<0)) then
294 0 : dtsets(idtset)%np_spkpt=1 ; dtsets(idtset)%npspinor=1 ; dtsets(idtset)%npfft=1
295 0 : dtsets(idtset)%npband=1; dtsets(idtset)%nphf=1
296 0 : dtsets(idtset)%paral_kgb=0
297 0 : if(optdriver/=RUNL_RESPFN) dtsets(idtset)%bandpp=1
298 0 : dtsets(idtset)%wfoptalg=0
299 0 : ABI_COMMENT('For non ground state calculations, set bandpp, npfft, npband, npspinor, np_spkpt and nphf to 1')
300 : end if
301 :
302 : ! Take into account a possible change of paral_kgb (change of the default algorithm)
303 6959 : if (.not.wfoptalg_read) then
304 6606 : if (dtsets(idtset)%usepaw==0) dtsets(idtset)%wfoptalg=0
305 6606 : if (dtsets(idtset)%usepaw/=0) dtsets(idtset)%wfoptalg=10
306 6606 : if ((optdriver==RUNL_GSTATE.or.optdriver==RUNL_GWLS).and.dtsets(idtset)%paral_kgb/=0) dtsets(idtset)%wfoptalg=114
307 : #ifndef HAVE_DFTI
308 : if (mod(dtsets(idtset)%wfoptalg,10)==4) then
309 : do iikpt=1,dtsets(idtset)%nkpt
310 : if (any(abs(dtsets(idtset)%kpt(:,iikpt))>tol8)) dtsets(idtset)%istwfk(iikpt)=1
311 : end do
312 : end if
313 : #else
314 6606 : if (mod(dtsets(idtset)%wfoptalg,10)==4.and.dtsets(idtset)%wfoptalg/=114) then
315 0 : do iikpt=1,dtsets(idtset)%nkpt
316 0 : if (any(abs(dtsets(idtset)%kpt(:,iikpt))>tol8)) dtsets(idtset)%istwfk(iikpt)=1
317 : end do
318 : end if
319 : #endif
320 : end if
321 :
322 6959 : if (dtsets(idtset)%wfoptalg==114.or.dtsets(idtset)%wfoptalg==14.or.dtsets(idtset)%wfoptalg==4) then !if LOBPCG
323 485 : if (dtsets(idtset)%autoparal==0) then
324 : !Nband might have different values for different kpoint, but not bandpp.
325 : !In this case, we just use the largest nband (mband_upper), and the input will probably fail
326 : !at the bandpp check later on
327 436 : if(tread(8)==1) then
328 : ! check if nblock_lobpcg is read from the input, if so error msg
329 34 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nblock_lobpcg',tread0,'INT')
330 34 : if (tread0==1) then
331 0 : write(msg,'(3a)') 'Both bandpp and nblock_lobpcg are defined for the same dataset, this is confusing.',ch10,&
332 0 : 'Change the input to keep only nblock_lobpcg (preferably) or bandpp.'
333 0 : ABI_ERROR(msg)
334 : end if
335 34 : if (mod(mband_upper,dtsets(idtset)%bandpp*dtsets(idtset)%npband)==0) then
336 34 : dtsets(idtset)%nblock_lobpcg=mband_upper/(dtsets(idtset)%bandpp*dtsets(idtset)%npband)
337 : else
338 0 : write(msg,'(5a)') 'mband_upper( =max_{kpt}(nband) ) should be a mutltiple of npband*bandpp.',ch10,&
339 0 : 'Change nband, npband or bandpp in the input.',ch10,&
340 0 : 'A simpler solution is to use nblock_lobpcg instead of bandpp.'
341 0 : ABI_ERROR(msg)
342 : end if
343 : else
344 402 : if (mod(mband_upper,dtsets(idtset)%nblock_lobpcg*dtsets(idtset)%npband)==0) then
345 402 : dtsets(idtset)%bandpp=mband_upper/(dtsets(idtset)%nblock_lobpcg*dtsets(idtset)%npband)
346 : else
347 0 : write(msg,'(3a)') 'mband_upper( =max_{kpt}(nband) ) should be a mutltiple of nblock_lobpcg*npband.',ch10,&
348 0 : 'Change nband, npband or nblock_lobpcg in the input.'
349 0 : ABI_ERROR(msg)
350 : end if
351 : end if
352 : else ! autoparal /= 0
353 : ! check if nblock_lobpcg is read from the input, if so error msg
354 49 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nblock_lobpcg',tread0,'INT')
355 49 : if (tread0==1) then
356 0 : write(msg,'(3a)') 'When using autoparal, nblock_lobpcg is automatically set, so it cannot be in the input.',ch10,&
357 0 : 'Change the input to keep only autoparal or nblock_lobpcg.'
358 0 : ABI_ERROR(msg)
359 : end if
360 49 : dtsets(idtset)%nblock_lobpcg=mband_upper/(dtsets(idtset)%bandpp*dtsets(idtset)%npband)
361 : end if
362 : end if
363 :
364 : ! Warning when using different number of bands for different kpoints (occopt=2)
365 6959 : if ( dtsets(idtset)%occopt==2 .and. dtsets(idtset)%nkpt>1 .and. &
366 : & ((dtsets(idtset)%bandpp > 1) .or. (dtsets(idtset)%npband > 1)) ) then
367 : write(msg, '(4a)' ) &
368 : 'When working with blocks of bands (bandpp>1 or npband>1),'&
369 0 : &' the number of bands should be the same for every kpoints.',ch10,&
370 0 : &' The run will most probably fail on a other check. If it does not fail, ignore this message.'
371 0 : ABI_WARNING(msg)
372 : end if
373 :
374 6959 : dtsets(idtset)%densfor_pred=2
375 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'densfor_pred',tread0,'INT')
376 6959 : if(tread0==1) then
377 67 : dtsets(idtset)%densfor_pred=intarr(1)
378 : else
379 6892 : if (dtsets(idtset)%paral_kgb==1) dtsets(idtset)%densfor_pred=6
380 : end if
381 : if((dtsets(idtset)%iscf==5.or.dtsets(idtset)%iscf==6) &
382 6959 : .and. dtsets(idtset)%ionmov==4 .and. dtsets(idtset)%densfor_pred/=3 )then
383 1 : dtsets(idtset)%densfor_pred=3
384 : write(msg, '(a,a,a)' )&
385 1 : 'When ionmov==4 and iscf==5 or 6, densfor_pred must be 3.',ch10,&
386 2 : 'Set densfor_pred to 3.'
387 1 : ABI_COMMENT(msg)
388 : end if
389 :
390 : #ifdef HAVE_LOTF
391 : ! LOTF need densfor_pred=2
392 : if(dtsets(idtset)%ionmov==23) dtsets(idtset)%densfor_pred=2
393 : #endif
394 :
395 6959 : if (usepaw==0) then
396 5778 : dtsets(idtset)%ortalg=2
397 : else
398 1181 : dtsets(idtset)%ortalg=-2
399 : end if
400 6959 : ortalg_read=.false.
401 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ortalg',tread0,'INT')
402 6959 : if(tread0==1) then
403 17 : dtsets(idtset)%ortalg=intarr(1)
404 17 : ortalg_read=.true.
405 6942 : else if (dtsets(idtset)%wfoptalg>=10 .and. dtsets(idtset)%ortalg>0) then
406 205 : dtsets(idtset)%ortalg=-dtsets(idtset)%ortalg
407 : end if
408 :
409 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iomode',tread0,'INT')
410 6959 : if(tread0==1) then
411 186 : dtsets(idtset)%iomode=intarr(1)
412 : else
413 : if ((xmpi_mpiio==1).and.(dtsets(idtset)%paral_kgb==1)) dtsets(idtset)%iomode=IO_MODE_MPI
414 : #ifdef HAVE_NETCDF_DEFAULT
415 6773 : dtsets(idtset)%iomode=IO_MODE_ETSF
416 : #endif
417 : end if
418 :
419 6959 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawmixdg',tread0,'INT')
420 6959 : if(tread0==1) then
421 51 : dtsets(idtset)%pawmixdg=intarr(1)
422 6908 : else if (dtsets(idtset)%npfft>1.and.usepaw==1) then
423 136 : dtsets(idtset)%pawmixdg=1
424 : end if
425 :
426 : ! Cycle if the processor is not used
427 6959 : if (mpi_enregs(idtset)%me<0.or.iexit>0) then
428 1 : ABI_FREE(intarr)
429 1 : ABI_FREE(dprarr)
430 1 : cycle
431 : end if
432 :
433 6958 : response=0
434 : if (dtsets(idtset)%rfddk/=0 .or. dtsets(idtset)%rf2_dkdk/=0 .or. dtsets(idtset)%rf2_dkde/=0 .or. &
435 6958 : & dtsets(idtset)%rfelfd/=0 .or. dtsets(idtset)%rfphon/=0 .or. dtsets(idtset)%rfstrs/=0 .or. &
436 718 : & dtsets(idtset)%rfmagn/=0) response=1
437 :
438 : ! If no MPI, set all npxxx variables to 1
439 6958 : if (nproc==1) then
440 5732 : dtsets(idtset)%np_spkpt = 1 ; dtsets(idtset)%npband = 1
441 5732 : dtsets(idtset)%npfft = 1 ; dtsets(idtset)%npspinor = 1
442 5732 : dtsets(idtset)%nphf = 1
443 : end if
444 :
445 : ! --IF CUDA AND RECURSION:ONLY BAND PARALLELISATION
446 6958 : if(dtsets(idtset)%tfkinfunc==2 .and. nproc/=1)then
447 4 : dtsets(idtset)%npband = dtsets(idtset)%npband*dtsets(idtset)%np_spkpt*dtsets(idtset)%npspinor*dtsets(idtset)%npfft
448 4 : dtsets(idtset)%bandpp=1
449 4 : dtsets(idtset)%np_spkpt = 1
450 4 : dtsets(idtset)%npfft = 1
451 4 : dtsets(idtset)%npspinor = 1
452 : write(msg, '(5a,i6,a)' )&
453 4 : 'If the recursion scheme is used (tfkinfunc==2),',ch10,&
454 4 : 'only the band parallelisation is active, we set:',ch10,&
455 8 : 'npfft= 1, np_spkpt= 1, npband=',dtsets(idtset)%npband,' .'
456 4 : ABI_WARNING(msg)
457 : end if
458 :
459 6958 : if (dtsets(idtset)%npspinor>=2.and.dtsets(idtset)%nspinor==1) then
460 0 : dtsets(idtset)%npspinor=1
461 0 : dtsets(idtset)%npfft=2*dtsets(idtset)%npfft
462 : write(msg,'(3a)')&
463 0 : 'npspinor is bigger than nspinor !',ch10,&
464 0 : 'We set npspinor to 1 ; we set npfft to 2*npfft'
465 0 : ABI_WARNING(msg)
466 : end if
467 :
468 : ! Some checks on parallelization data
469 6958 : if(dtsets(idtset)%paral_kgb < 0 ) then
470 : cycle
471 6958 : else if(dtsets(idtset)%paral_kgb/=0.and.(dtsets(idtset)%bandpp/=1.or.dtsets(idtset)%npband/=1.or.&
472 : dtsets(idtset)%npfft/=1.or.dtsets(idtset)%np_spkpt/=1.or.dtsets(idtset)%npspinor/=1))then
473 624 : if(dtsets(idtset)%np_spkpt*dtsets(idtset)%npfft*dtsets(idtset)%npband*dtsets(idtset)%npspinor > nproc )then
474 : write(msg,'(7a)')&
475 0 : 'The product of np_spkpt, npfft, npband and npspinor is bigger than the number of processors.',ch10,&
476 0 : 'The user-defined values of np_spkpt, npfft, npband or npspinor will be modified,',ch10,&
477 0 : 'in order to bring this product below nproc .',ch10,&
478 0 : 'At present, only a very simple algorithm is used ...'
479 0 : ABI_WARNING(msg)
480 :
481 0 : if(dtsets(idtset)%np_spkpt*dtsets(idtset)%npband*dtsets(idtset)%npspinor <= nproc) then
482 0 : dtsets(idtset)%npfft=1
483 0 : ABI_WARNING('Set npfft to 1')
484 0 : else if(dtsets(idtset)%np_spkpt*dtsets(idtset)%npspinor <= nproc)then
485 0 : dtsets(idtset)%npfft=1
486 0 : dtsets(idtset)%npband=1
487 0 : ABI_WARNING('Set npfft and npband to 1')
488 0 : else if(dtsets(idtset)%np_spkpt <= nproc)then
489 0 : dtsets(idtset)%npfft=1
490 0 : dtsets(idtset)%npband=1
491 0 : dtsets(idtset)%npspinor=1
492 0 : ABI_WARNING('Set npfft ,npband and npspinor to 1')
493 : else
494 0 : dtsets(idtset)%npfft=1
495 0 : dtsets(idtset)%npband=1
496 0 : dtsets(idtset)%np_spkpt=1
497 0 : dtsets(idtset)%npspinor=1
498 0 : ABI_WARNING('Set npfft, npband, nspinor and np_spkpt to 1')
499 : end if
500 624 : else if(dtsets(idtset)%np_spkpt*dtsets(idtset)%npfft*dtsets(idtset)%npband*dtsets(idtset)%npspinor < nproc)then
501 : write(msg,'(a,5i6,4a)')&
502 0 : 'np_spkpt,npfft,npband,npspinor,nproc=',&
503 0 : & dtsets(idtset)%np_spkpt,dtsets(idtset)%npfft,dtsets(idtset)%npband,dtsets(idtset)%npspinor,nproc,ch10,&
504 0 : 'The number of processors must not be greater than npfft*npband*np_spkpt*npspinor ',ch10,&
505 0 : 'when npfft or np_spkpt or npband or npspinor are chosen manually in the input file.'
506 0 : ABI_ERROR(msg)
507 : end if
508 : end if
509 :
510 : ! LOBPCG and ChebFi need paral_kgb=1 in parallel
511 626 : if ((dtsets(idtset)%npband*dtsets(idtset)%npfft>1).and. &
512 6958 : & (mod(dtsets(idtset)%wfoptalg,10)==1.or.mod(dtsets(idtset)%wfoptalg,10)==4)) then
513 578 : dtsets(idtset)%paral_kgb=1
514 : end if
515 :
516 : ! Check size of Scalapack communicator
517 : #ifdef HAVE_LINALG_ELPA
518 : if(dtsets(idtset)%paral_kgb>0.and.dtsets(idtset)%np_slk>0) then
519 : np=min(dtsets(idtset)%np_slk,dtsets(idtset)%npband*dtsets(idtset)%npfft*dtsets(idtset)%npspinor)
520 : irow=int(sqrt(float(np)))
521 : do while(mod(np,irow)/=0)
522 : irow=irow-1
523 : end do
524 : icol=nproc/irow
525 : if (icol>mband_lower) then
526 : do while(icol>mband_lower)
527 : icol=icol-1
528 : do while(mod(np,icol)/=0)
529 : icol=icol-1
530 : end do
531 : end do
532 : dtsets(idtset)%np_slk=icol
533 : write(msg,'(5a,i6,a)')&
534 : 'The number of band*fft*spinor processors was not consistent with',ch10,&
535 : 'the size of communicator used for ELPA library (np_slk).',ch10,&
536 : 'np_slk value has been adjusted to ',dtsets(idtset)%np_slk,'.'
537 : ABI_COMMENT(msg)
538 : end if
539 : end if
540 : #endif
541 :
542 : !Additional check in case of a parallelized Hartree-Fock calculation
543 : ! %usefock == option to perform Fock exchange calculation
544 : ! %nphf == number of processors for Fock exchange calculation
545 6958 : if ((dtsets(idtset)%usefock==1).and.(dtsets(idtset)%nphf/=1)) then
546 :
547 0 : if ((dtsets(idtset)%nphf<0).or.(dtsets(idtset)%nphf==0)) then
548 0 : ABI_ERROR('The value of variable nphf should be a non negative integer.')
549 : end if
550 0 : if (dtsets(idtset)%paral_kgb/=0) then
551 0 : ABI_ERROR('Option paral_kgb should be turned off (value 0) for a parallelized Hartree-Fock calculation.')
552 : end if
553 0 : if (response/=0) then
554 0 : ABI_ERROR('A response function calculation is not yet possible with a parallelized Hartree-Fock calculation.')
555 : end if
556 0 : if (dtsets(idtset)%npspinor>1) then
557 0 : ABI_ERROR('The parallelism on spinors is not supported by a parallelized Hartree-Fock calculation.')
558 : end if
559 0 : if (dtsets(idtset)%np_spkpt*dtsets(idtset)%nphf > nproc )then
560 0 : write(msg,'(a,3(a,i0))') ch10,&
561 0 : 'The product of variables np_spkpt and nphf is bigger than the number of processors: np_spkpt= ',&
562 0 : dtsets(idtset)%np_spkpt,' nphf= ',dtsets(idtset)%nphf ,' and nproc= ', nproc
563 0 : ABI_ERROR(msg)
564 : end if
565 : end if ! Fock
566 :
567 : !When using chebfi, the number of blocks is equal to the number of processors
568 : if((dtsets(idtset)%wfoptalg == 1) .or. (dtsets(idtset)%wfoptalg == 111) .or. &
569 6958 : & (dtsets(idtset)%wfoptalg == 2) .or. (dtsets(idtset)%wfoptalg == 112)) then
570 : !Nband might have different values for different kpoint, but not bandpp.
571 : !In this case, we just use the largest nband (mband_upper), and the input will probably fail
572 : !at the bandpp check later on
573 208 : dtsets(idtset)%bandpp = mband_upper / dtsets(idtset)%npband
574 208 : if(tread(8)==1) then
575 : write(msg, '(a,i8,3a)' ) &
576 0 : 'bandpp has been internally set to ',dtsets(idtset)%bandpp,'.',ch10,&
577 : 'Indeed, there is no need to specify bandpp in the input when using chebfi (wfoptalg=1,111)&
578 0 : or spectrum slicing (wfoptalg=2,112).'
579 0 : ABI_COMMENT(msg)
580 : end if
581 : end if
582 :
583 : !Check parallelization in case of RTTDDFT
584 : !In particular ensure that bandpp = nband / npband
585 6958 : if (optdriver == RUNL_RTTDDFT) then
586 50 : dtsets(idtset)%bandpp = mband_upper / dtsets(idtset)%npband
587 50 : if ( tread(8) == 1 ) then
588 0 : write(msg, '(a,a)') 'Setting bandpp is useless in RT-TDDFT because it is automatically set to nband/npband.', ch10
589 0 : ABI_WARNING(msg)
590 : end if
591 50 : if (dtsets(idtset)%npfft/=1) then
592 0 : dtsets(idtset)%npfft=1
593 0 : write(msg, '(a,a)') 'RT-TDDFT is not compatible with FFT-parallelization. Remove npfft or set it to 1.', ch10
594 0 : ABI_ERROR(msg)
595 : end if
596 50 : if (dtsets(idtset)%npspinor/=1) then
597 0 : dtsets(idtset)%npspinor=1
598 0 : write(msg, '(a,a)') 'RT-TDDFT is not compatible with spinor parallelization. Remove npspinor or set it to 1.', ch10
599 0 : ABI_ERROR(msg)
600 : end if
601 50 : if (dtsets(idtset)%nphf/=1) then
602 0 : dtsets(idtset)%nphf=1
603 0 : write(msg, '(a,a)') 'RT-TDDFT is not compatible with HF parallelization. Remove nphf or set it to 1.', ch10
604 0 : ABI_ERROR(msg)
605 : end if
606 : end if
607 :
608 : ! Set mpi_enreg
609 6958 : mpi_enregs(idtset)%paral_kgb=dtsets(idtset)%paral_kgb
610 6958 : if(dtsets(idtset)%paral_kgb/=0)then
611 624 : mpi_enregs(idtset)%nproc_spkpt=dtsets(idtset)%np_spkpt
612 624 : mpi_enregs(idtset)%nproc_fft=dtsets(idtset)%npfft
613 624 : mpi_enregs(idtset)%nproc_band=dtsets(idtset)%npband
614 624 : mpi_enregs(idtset)%nproc_spinor=min(dtsets(idtset)%npspinor,dtsets(idtset)%nspinor)
615 624 : mpi_enregs(idtset)%bandpp=dtsets(idtset)%bandpp
616 : ! Additional setting in case of hybrid functional calculation => not yet tested (CMartins)
617 : ! if (dtsets(idtset)%usefock==1) then
618 : ! mpi_enregs(idtset)%nproc_hf = dtsets(idtset)%nphf
619 : ! if (dtsets(idtset)%nphf>1) mpi_enregs(idtset)%paral_hf=1
620 : ! end if
621 : else
622 6334 : mpi_enregs(idtset)%bandpp = dtsets(idtset)%bandpp
623 : ! Additional setting in case of a Fock exchange of PBE0 calculation
624 6334 : if (dtsets(idtset)%usefock==1) then
625 146 : if (dtsets(idtset)%nphf>1) mpi_enregs(idtset)%paral_hf=1
626 146 : mpi_enregs(idtset)%nproc_hf = dtsets(idtset)%nphf
627 146 : if (dtsets(idtset)%np_spkpt/=1) then
628 0 : mpi_enregs(idtset)%nproc_spkpt = dtsets(idtset)%np_spkpt
629 : else
630 146 : mpi_enregs(idtset)%nproc_spkpt = mpi_enregs(idtset)%nproc_cell/mpi_enregs(idtset)%nproc_hf
631 : end if
632 : else
633 6188 : mpi_enregs(idtset)%nproc_spkpt = mpi_enregs(idtset)%nproc_cell
634 : end if
635 : end if
636 :
637 6958 : if(dtsets(idtset)%paral_kgb>=0) then
638 :
639 : ! Compute processor distribution over perturbations
640 6958 : mpi_enregs(idtset)%paral_pert=dtsets(idtset)%paral_rf
641 6958 : if (mpi_enregs(idtset)%paral_pert==1) then
642 26 : dtsets(idtset)%nppert=max(1,dtsets(idtset)%nppert)
643 26 : if(dtsets(idtset)%nppert>mpi_enregs(idtset)%nproc) then
644 0 : ABI_ERROR('The number of processors must not be smaller than nppert !')
645 : end if
646 26 : call initmpi_pert(dtsets(idtset),mpi_enregs(idtset))
647 26 : mpi_enregs(idtset)%nproc_spkpt = mpi_enregs(idtset)%nproc_cell
648 26 : nproc=mpi_enregs(idtset)%nproc_cell
649 : end if
650 : ! Cycle if the processor is not used
651 6958 : if (mpi_enregs(idtset)%me<0) then
652 0 : ABI_FREE(intarr)
653 0 : ABI_FREE(dprarr)
654 0 : cycle
655 : end if
656 :
657 : ! Compute processor distribution over kpt (and eventually band-fft)
658 6958 : call initmpi_grid(mpi_enregs(idtset))
659 6958 : if(dtsets(idtset)%usewvl==1) mpi_enregs(idtset)%comm_fft=mpi_enregs(idtset)%comm_cell
660 :
661 : ! Initialize tabs used for k/spin parallelism (with sequential-type values)
662 34790 : ABI_MALLOC(mpi_enregs(idtset)%proc_distrb,(nkpt,mband_upper,nsppol))
663 20874 : ABI_MALLOC(mpi_enregs(idtset)%my_kpttab,(nkpt))
664 1053202 : mpi_enregs(idtset)%proc_distrb(:,:,:)=0
665 233568 : mpi_enregs(idtset)%my_kpttab(:)=(/(ii,ii=1,nkpt)/)
666 20874 : mpi_enregs(idtset)%my_isppoltab(:)=1;if (dtsets(idtset)%nsppol==1) mpi_enregs(idtset)%my_isppoltab(2)=0
667 :
668 : ! HF or hybrid calculation : initialization of the array distrb_hf
669 6958 : if (dtsets(idtset)%usefock==1) then
670 584 : ABI_MALLOC(mpi_enregs(idtset)%distrb_hf,(dtsets(idtset)%nkpthf,dtsets(idtset)%nbandhf,1))
671 : ! The dimension of distrb_hf are given by %nkpthf and %nbandhf.
672 : ! We assume that there will be no dependence in spinpol for all the occupied states.
673 12676 : mpi_enregs(idtset)%distrb_hf=0
674 : end if
675 :
676 : ! Define k-points distribution (determine who I am)
677 : ! Note that nkpt_me may differ from processor to processor
678 : ! This fact will NOT be taken into account when
679 : ! the memory needs will be evaluated in the subroutine memory.
680 : ! Also, the reduction of k points due to symmetry in RF calculations
681 : ! is NOT taken into account. This should be changed later ...
682 6958 : nkpt_me=nkpt
683 6958 : mband_mem=0
684 6958 : if(xmpi_paral==1 .and. dtsets(idtset)%usewvl == 0) then
685 6958 : nkpt_me=0
686 6958 : if(response==0 .or. (response==1 .and. dtsets(idtset)%efmas==1))then
687 6257 : mpi_enregs(idtset)%paralbd=0
688 6257 : call distrb2(mband_upper,mband_mem,dtsets(idtset)%nband,nkpt,nproc,nsppol,mpi_enregs(idtset))
689 64024 : do iikpt=1,nkpt
690 57767 : if(.not.(proc_distrb_cycle(mpi_enregs(idtset)%proc_distrb,iikpt,1,1,-1,mpi_enregs(idtset)%me_kpt)))&
691 56032 : & nkpt_me=nkpt_me+1
692 : end do ! ikpt=1,nkpt
693 : ! HF or hybrid calculation : define the occupied states distribution (in array distrb_hf)
694 6257 : if (dtsets(idtset)%usefock==1) then
695 146 : call distrb2_hf(dtsets(idtset)%nbandhf,dtsets(idtset)%nkpthf,nproc,nsppol,mpi_enregs(idtset))
696 : end if
697 : else ! response==1
698 : ! TODO: check or remove the following comment which seems outdated
699 : ! Wrongly assumes that the number of elements of the
700 : ! k-point sets of the two spin polarizations is the maximal
701 : ! value of one of these k-point sets ...
702 : ! This is to be corrected when RF is implemented
703 : ! for spin-polarized case.
704 : ! ENDTODO
705 701 : mpi_enregs(idtset)%paralbd=1
706 : ! nproc=mpi_enregs(idtset)%nproc_cell*mpi_enregs(idtset)%nproc_pert
707 701 : call distrb2(mband_upper,mband_mem,dtsets(idtset)%nband,nkpt,nproc,nsppol,mpi_enregs(idtset))
708 1463 : do isppol=1,nsppol
709 762 : nspink=0
710 57104 : do iikpt=1,nkpt
711 119704 : do iband=1,dtsets(idtset)%nband(iikpt+(isppol-1)*nkpt)
712 118942 : if(mpi_enregs(idtset)%proc_distrb(iikpt,iband,isppol)==mpi_enregs(idtset)%me_cell)then
713 44662 : nspink=nspink+1
714 44662 : exit
715 : end if
716 : end do ! iband
717 : end do ! iikpt
718 1463 : if(nspink>nkpt_me)nkpt_me=nspink
719 : end do ! isppol
720 : ! Is nband present in input file or automatically estimated ?
721 : tnband=0
722 701 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nband',tnband,'INT')
723 : ! If the number of bands was estimated, there might be a side effect
724 : ! when the definitive number of bands is known. k points
725 : ! might be attributed to different processors than the present
726 : ! proc_distrb describes. At most, the number of k points could increase by 1 ...
727 701 : if(tnband==0)nkpt_me=nkpt_me+1
728 : ! In any case, the maximal number of k points is nkpt
729 701 : if(nkpt_me>nkpt)nkpt_me=nkpt
730 :
731 : ! mband_mem
732 2103 : ABI_MALLOC (mybands, (mband_upper))
733 701 : mband_mem = 0
734 1463 : do isppol=1,nsppol
735 57805 : do iikpt=1,nkpt
736 395169 : mybands = 0
737 395149 : do iband=1,dtsets(idtset)%nband(iikpt+(isppol-1)*nkpt)
738 395149 : if(mpi_enregs(idtset)%proc_distrb(iikpt,iband,isppol)==mpi_enregs(idtset)%me_band)then
739 276135 : mybands(iband)=1
740 : end if
741 : end do ! iband
742 395931 : mband_mem = max(mband_mem, sum(mybands))
743 : end do ! iikpt
744 : end do ! isppol
745 701 : ABI_FREE (mybands)
746 : end if ! response case
747 : end if
748 : end if
749 6958 : if (mband_mem == 0) mband_mem = mband_upper
750 6958 : dtsets(idtset)%mband_mem = mband_mem
751 :
752 : ! Take care of mkmems. Use the generic name -mkmem- for mkmem as well as mkqmem
753 : ! and mk1mem.
754 6958 : nm_mkmem(1)='mkmem '
755 6958 : nm_mkmem(2)='mkqmem'
756 6958 : nm_mkmem(3)='mk1mem'
757 :
758 27832 : do ii=1,3
759 :
760 : ! Read in mkmem here if it is in the input file
761 : ! TODO: mkmem is not supported any longer. These variables can be removed.
762 20874 : if(ii==1)then
763 6958 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'mkmem',tread0,'INT')
764 13916 : else if(ii==2)then
765 6958 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'mkqmem',tread0,'INT')
766 : else if(ii==3)then
767 6958 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'mk1mem',tread0,'INT')
768 : end if
769 :
770 :
771 : ! Note that mkmem is used as a dummy variable, representing mkmem as well
772 : ! as mkqmem, and mk1mem.
773 20874 : if(tread0==1) then
774 0 : mkmem=intarr(1)
775 0 : if (mkmem<0) then
776 : ! mkmem is unreasonable; must be zero or positive
777 : write(msg, '(4a,i0,4a)')&
778 0 : nm_mkmem(ii),' must be positive or null but ',nm_mkmem(ii),' =',mkmem,ch10,&
779 0 : 'Use default ',nm_mkmem(ii),' = nkpt .'
780 0 : ABI_WARNING(msg)
781 0 : mkmem=nkpt
782 : end if
783 :
784 : else
785 :
786 : ! mkmem was not set in the input file so default to incore solution
787 : !write(msg,'(6a)') &
788 : !'mpi_setup: ',nm_mkmem(ii),' undefined in the input file.','Use default ',nm_mkmem(ii),' = nkpt'
789 : !call wrtout(std_out, msg)
790 20874 : mkmem=nkpt
791 : end if
792 :
793 : ! Check whether nkpt distributed on the processors <= mkmem;
794 : ! if so then may run entirely in core,
795 : ! avoiding i/o to disk for wavefunctions and kg data.
796 : ! mkmem/=0 to avoid i/o; mkmem==0 to use disk i/o for nkpt>=1.
797 20874 : if (nkpt_me<=mkmem .and. mkmem/=0 ) then
798 : write(msg, '(a,i0,a,a,a,i0,a)' ) &
799 20874 : ' mpi_setup: With nkpt_me=',nkpt_me,' and ',nm_mkmem(ii),' = ',mkmem,', ground state wf handled in core.'
800 20874 : if (prtvol > 0) call wrtout(std_out,msg)
801 20874 : if(nkpt_me<mkmem .and. nkpt_me/=0)then
802 1890 : write(msg,'(3a)')' Resetting ',nm_mkmem(ii),' to nkpt_me to save memory space.'
803 1890 : mkmem=nkpt_me
804 1890 : if (prtvol > 0) call wrtout(std_out,msg)
805 : end if
806 0 : else if(mkmem/=0)then
807 : write(msg, '(a,i0,3a,i0,5a)' ) &
808 0 : ' mpi_setup: With nkpt_me=',nkpt_me,'and ',nm_mkmem(ii),' = ',mkmem,&
809 0 : ' ground state wf require disk i/o.',ch10,&
810 0 : ' Resetting ',nm_mkmem(ii),' to zero to save memory space.'
811 0 : mkmem=0
812 0 : if (prtvol > 0) call wrtout(std_out,msg)
813 : end if
814 20874 : if(dtsets(idtset)%usewvl == 0 .or. dtsets(idtset)%usepaw==1)then
815 20874 : if(ii==1)dtsets(idtset)%mkmem=mkmem
816 : end if
817 20874 : if(ii==2)dtsets(idtset)%mkqmem=mkmem
818 20874 : if(ii==3)dtsets(idtset)%mk1mem=mkmem
819 :
820 27832 : if(dtsets(idtset)%usewvl == 1 .and. dtsets(idtset)%usepaw==1 )then
821 0 : if(dtsets(idtset)%mkmem .ne. dtsets(idtset)%nkpt) then
822 0 : ABI_ERROR("mkmem is not allowed for WVL+PAW")
823 : end if
824 : end if
825 :
826 : end do ! End the loop on the three possibilities mkmem, mkqmem, mk1mem.
827 :
828 6958 : if(dtsets(idtset)%paral_kgb==1) mpi_enregs(idtset)%paralbd=0
829 :
830 : ! Check if some MPI processes are empty (MBPT codes uses a complete different MPI algorithm)
831 39899 : do_check = all(optdriver /= [RUNL_SCREENING, RUNL_SIGMA, RUNL_BSE, RUNL_EPH, RUNL_GWR])
832 6958 : if (dtsets(idtset)%usewvl == 0 .and. do_check) then
833 6471 : if (.not.mpi_distrib_is_ok(mpi_enregs(idtset),mband_upper,&
834 : dtsets(idtset)%nkpt,dtsets(idtset)%mkmem,nsppol,msg=msg)) then
835 74 : write(msg,'(5a)') trim(msg),ch10,&
836 74 : 'YOU ARE STRONGLY ADVISED TO ACTIVATE AUTOMATIC PARALLELIZATION!',ch10,&
837 148 : 'USE "AUTOPARAL=1" IN THE INPUT FILE.'
838 74 : ABI_WARNING(msg)
839 : end if
840 : end if
841 :
842 : ! call mpi_setup1(dtsets(idtset),jdtset,lenstr,mband_upper,mpi_enregs(idtset),string)
843 : ! Printing of processor distribution
844 : ! MPIWF : here, set up the complete ngfft, containing the information
845 : ! for the parallelisation of the FFT
846 6958 : call abi_io_redirect(new_io_comm=mpi_enregs(idtset)%comm_world)
847 6958 : call libpaw_write_comm_set(mpi_enregs(idtset)%comm_world)
848 :
849 : ! Default values for sequential case
850 6958 : paral_fft=0; nproc_fft=1; me_fft=0
851 :
852 6958 : if(dtsets(idtset)%usewvl == 0)then
853 6958 : if(optdriver==RUNL_GSTATE.or.optdriver==RUNL_GWLS) then
854 5648 : paral_fft=1 ! parallelisation over FFT
855 5648 : if (mpi_enregs(idtset)%nproc_cell>0) then
856 5648 : if(mpi_enregs(idtset)%paral_kgb == 1) then
857 :
858 592 : if((dtsets(idtset)%gpu_option/=ABI_GPU_DISABLED).and.(mpi_enregs(idtset)%nproc_fft/=1))then
859 : write(msg,'(3a,i0)') &
860 0 : 'When the use of GPU is on, the number of FFT processors, npfft, must be 1',ch10,&
861 0 : 'However, npfft=',mpi_enregs(idtset)%nproc_fft
862 0 : ABI_ERROR(msg)
863 : end if
864 :
865 592 : if((dtsets(idtset)%cprj_in_memory/=0)) then
866 122 : if (mpi_enregs(idtset)%nproc_spinor/=1) then
867 : write(msg,'(3a,i0)') &
868 0 : 'If cprj_in_memory/=0, the number of processors for spinors, npspinor, must be 1',ch10,&
869 0 : 'However, npspinor=',mpi_enregs(idtset)%nproc_spinor
870 0 : ABI_ERROR(msg)
871 : end if
872 122 : if(mpi_enregs(idtset)%nproc_fft/=1)then
873 : write(msg,'(3a,i0)') &
874 0 : 'If cprj_in_memory/=0, the number of FFT processors, npfft, must be 1',ch10,&
875 0 : 'However, npfft=',mpi_enregs(idtset)%nproc_fft
876 0 : ABI_ERROR(msg)
877 : end if
878 : end if
879 :
880 592 : if(modulo(dtsets(idtset)%ngfft(2),mpi_enregs(idtset)%nproc_fft)/=0)then
881 : write(msg,'(3a,i0,a,i0)') &
882 0 : 'The number of FFT processors, npfft, should be a multiple of ngfft(2).',ch10,&
883 0 : 'However, npfft=',mpi_enregs(idtset)%nproc_fft,' and ngfft(2)=',dtsets(idtset)%ngfft(2)
884 0 : ABI_BUG(msg)
885 : end if
886 :
887 4032 : do iikpt=1,nkpt*nsppol
888 3440 : iikpt_modulo = modulo(iikpt,nkpt)+1
889 4032 : if ((dtsets(idtset)%istwfk(iikpt_modulo)==2)) then !.and.(dtsets(idtset)%ngfft(7)==401)) then
890 171 : if ((mpi_enregs(idtset)%bandpp==0).or. &
891 : ((mpi_enregs(idtset)%bandpp/=1).and.(modulo(mpi_enregs(idtset)%bandpp,2)/=0))) then
892 : write(msg,'(3a,i0)') &
893 0 : 'The number bandpp should be 1 or a multiple of 2',ch10,&
894 0 : 'However, bandpp=',mpi_enregs(idtset)%bandpp
895 0 : ABI_BUG(msg)
896 : end if
897 171 : if(modulo(dtsets(idtset)%nband(iikpt),mpi_enregs(idtset)%nproc_band*mpi_enregs(idtset)%bandpp)/=0)then
898 : write(msg,'(3a,i0,a,i0)') &
899 0 : 'The number of bands for the k-point, nband_k, should be a multiple of nproc_band*bandpp.',ch10,&
900 0 : 'However, nband_k=',dtsets(idtset)%nband(iikpt),' and nproc_band*bandpp=', &
901 0 : mpi_enregs(idtset)%nproc_band* mpi_enregs(idtset)%bandpp
902 0 : ABI_BUG(msg)
903 : end if
904 : else if ((dtsets(idtset)%istwfk(iikpt_modulo)==2) .and. (dtsets(idtset)%ngfft(7)==400)) then
905 : ABI_BUG('The fftalg=400 with istwfk=2 is not valid')
906 : else
907 3269 : if(modulo(dtsets(idtset)%nband(iikpt),mpi_enregs(idtset)%nproc_band*mpi_enregs(idtset)%bandpp)/=0)then
908 : write(msg,'(3a,i0,a,i0)') &
909 0 : 'The number of band for the k-point, nband_k, should be a multiple of nproc_band*bandpp.',ch10,&
910 0 : 'However, nband_k=',dtsets(idtset)%nband(iikpt),' and nproc_band*bandpp=', &
911 0 : mpi_enregs(idtset)%nproc_band* mpi_enregs(idtset)%bandpp
912 0 : ABI_BUG(msg)
913 : end if
914 3269 : if ((mpi_enregs(idtset)%bandpp==0)) then
915 : write(msg,'(a,i0,2a,i0,2a,i0)')&
916 0 : 'The number bandpp should not be 0 with fftalg=',dtsets(idtset)%ngfft(7),ch10,&
917 0 : 'and istwfk=',dtsets(idtset)%istwfk(iikpt_modulo),ch10,&
918 0 : 'However, bandpp=',mpi_enregs(idtset)%bandpp
919 0 : ABI_BUG(msg)
920 : end if
921 : end if
922 : end do
923 :
924 : if (xmpi_paral==1) then
925 592 : if(modulo(nkpt*nsppol,mpi_enregs(idtset)%nproc_spkpt)/=0)then
926 : write(msg,'(3a,i0,a,i0)') &
927 16 : 'The number of KPT processors, np_spkpt, should be a multiple of nkpt*nsppol.',ch10,&
928 32 : 'However, np_spkpt=',mpi_enregs(idtset)%nproc_spkpt,' and nkpt*nsppol=',nkpt*nsppol
929 16 : ABI_WARNING(msg)
930 : end if
931 : end if
932 : else
933 51183 : do iikpt=1,nkpt*nsppol
934 46127 : iikpt_modulo = modulo(iikpt,nkpt)+1
935 51183 : if(modulo(dtsets(idtset)%nband(iikpt),mpi_enregs(idtset)%nproc_band*mpi_enregs(idtset)%bandpp)/=0)then
936 : write(msg,'(3a,i0,a,i0)') &
937 0 : 'The number of band for the k-point, nband_k, should be a multiple of npband*bandpp.',ch10,&
938 0 : 'However, nband_k=',dtsets(idtset)%nband(iikpt),' and npband*bandpp=', &
939 0 : mpi_enregs(idtset)%nproc_band* mpi_enregs(idtset)%bandpp
940 0 : ABI_BUG(msg)
941 : end if
942 : end do
943 : end if
944 : end if
945 5648 : nproc_fft=mpi_enregs(idtset)%nproc_fft
946 5648 : me_fft=mpi_enregs(idtset)%me_fft
947 : end if
948 : end if
949 :
950 : ! Compute mgfft,mpw,nfft for this data set (it is dependent of mpi_enreg)
951 20874 : ABI_MALLOC(istwfk,(nkpt))
952 20874 : ABI_MALLOC(kpt_with_shift,(3,nkpt))
953 :
954 : ! Set the default value of fftalg for given npfft but allow the user to override it.
955 : ! Warning: If you need to change npfft, **DO IT** before this point so that here we get the correct fftalg
956 6958 : forbid_threads=(nthreads>1.and.dtsets(idtset)%bandpp>1)
957 6958 : dtsets(idtset)%ngfft(7) = fftalg_for_npfft(dtsets(idtset)%npfft,forbid_threads=forbid_threads)
958 6958 : dtsets(idtset)%ngfftdg(7) = fftalg_for_npfft(dtsets(idtset)%npfft,forbid_threads=forbid_threads)
959 :
960 : ! For RT-TDDFT make sure that we use the thread-safe version of FFT
961 : ! in case of Goedecker's FFT with more than one thread
962 6958 : if (optdriver==RUNL_RTTDDFT) then
963 50 : if (dtsets(idtset)%ngfft(7)/100==FFT_SG .and. nthreads>1) then
964 0 : write(msg,'(3a)') 'fftalg=1XX is not thread-safe, so it cannot be used with nthreads>1',ch10,&
965 0 : 'thus switching fftalg to a thread-safe version.'
966 0 : ABI_WARNING(msg)
967 0 : dtsets(idtset)%ngfft(7) = 401
968 0 : dtsets(idtset)%ngfftdg(7) = 401
969 : end if
970 : end if
971 :
972 6958 : fftalg_read=.false.
973 6958 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'fftalg',tread0,'INT')
974 :
975 6958 : if (tread0==1) then
976 261 : dtsets(idtset)%ngfft(7)=intarr(1)
977 261 : if (usepaw==1) dtsets(idtset)%ngfftdg(7)=intarr(1)
978 : fftalg_read=.true.
979 : end if
980 :
981 6958 : ecut =dtsets(idtset)%ecut
982 6958 : dilatmx =dtsets(idtset)%dilatmx
983 132202 : ngfft(:) =dtsets(idtset)%ngfft(:)
984 120263 : istwfk(:)=dtsets(idtset)%istwfk(1:nkpt)
985 6958 : nsym =dtsets(idtset)%nsym
986 :
987 6958 : nqpt=dtsets(idtset)%nqpt
988 9097 : qphon(:)=zero;if(nqpt/=0) qphon(:)=dtsets(idtset)%qptn(:)
989 :
990 20874 : ABI_MALLOC(symrel,(3,3,nsym))
991 2199265 : symrel(:,:,1:nsym)=dtsets(idtset)%symrel(:,:,1:nsym)
992 6958 : ecut_eff=ecut*dilatmx**2
993 :
994 6958 : if (usepaw==1) call wrtout(std_out,'getng is called for the coarse grid:')
995 6958 : kpt=k0; if (response==1.and.usepaw==1) kpt=qphon ! this is temporary
996 :
997 : call getng(dtsets(idtset)%boxcutmin,dtsets(idtset)%chksymtnons,ecut_eff,gmet,kpt,me_fft,mgfft,nfft,&
998 : & ngfft,nproc_fft,nsym,paral_fft,symrel,dtsets(idtset)%tnons,&
999 6958 : & gpu_option=dtsets(idtset)%gpu_option)
1000 :
1001 132202 : dtsets(idtset)%ngfft(:)=ngfft(:)
1002 6958 : dtsets(idtset)%mgfft=mgfft
1003 6958 : dtsets(idtset)%nfft=nfft
1004 460178 : kpt_with_shift(:,:)=dtsets(idtset)%kpt(:,1:nkpt)/dtsets(idtset)%kptnrm
1005 :
1006 6958 : exchn2n3d=dtsets(idtset)%exchn2n3d
1007 6958 : nproc_fft=ngfft(10) ; me_fft=ngfft(11)
1008 6958 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=mod(fftalg,10)
1009 :
1010 : ! Initialize tables for MPI-FFT.
1011 6958 : call mpi_enregs(idtset)%distribfft%init('c',mpi_enregs(idtset)%nproc_fft,ngfft(2),ngfft(3))
1012 :
1013 6958 : if(response/=0)then
1014 : ! This value of mpw is used in the first part of respfn.f
1015 718 : call getmpw(ecut_eff,exchn2n3d,gmet,istwfk,kpt_with_shift,mpi_enregs(idtset),mpw_k,nkpt)
1016 : end if
1017 6958 : if(nqpt/=0)then
1018 50414 : kpt_with_shift(1,:)=kpt_with_shift(1,:)+qphon(1)
1019 50414 : kpt_with_shift(2,:)=kpt_with_shift(2,:)+qphon(2)
1020 50414 : kpt_with_shift(3,:)=kpt_with_shift(3,:)+qphon(3)
1021 : end if
1022 6958 : if (dtsets(idtset)%usewvl == 0) then
1023 6958 : if (dtsets(idtset)%neb_cell_algo==NEB_CELL_ALGO_NONE) then
1024 6957 : call getmpw(ecut_eff,exchn2n3d,gmet,istwfk,kpt_with_shift,mpi_enregs(idtset),mpw,nkpt)
1025 : else
1026 1 : max_mpw=0
1027 8 : do ii=1,dtsets(idtset)%nimage
1028 7 : call mkrdim(dtsets(idtset)%acell_orig(1:3,ii),dtsets(idtset)%rprim_orig(1:3,1:3,ii),rprimd)
1029 7 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1030 7 : call getmpw(ecut_eff,exchn2n3d,gmet,istwfk,kpt_with_shift,mpi_enregs(idtset),mpw,nkpt)
1031 8 : if (mpw>max_mpw) max_mpw=mpw
1032 : end do
1033 1 : mpw=max_mpw
1034 : end if
1035 :
1036 : ! Allocate tables for parallel IO of the wavefunctions.
1037 624 : if( xmpi_mpiio==1 .and. mpi_enregs(idtset)%paral_kgb == 1 .and. &
1038 6958 : & any(dtsets(idtset)%iomode == [IO_MODE_MPI, IO_MODE_ETSF])) then
1039 1968 : ABI_MALLOC(mpi_enregs(idtset)%my_kgtab,(mpw,dtsets(idtset)%mkmem))
1040 : end if
1041 : else
1042 0 : mpw = 0
1043 : end if
1044 :
1045 : ! The dimensioning, in the RF case, should be done only with mpw,
1046 : ! but mpw is used in the first part of respfn.f, and should at least
1047 : ! be equal to mpw_k . The chosen way to code is not optimal, only convenient :
1048 : ! it leads to a small waste of memory.
1049 6958 : if(response/=0 .and. mpw_k>mpw)mpw=mpw_k
1050 132202 : dtsets(idtset)%ngfft(:)=ngfft(:)
1051 :
1052 : ! Initialize ngfftc to the initial guess for the coarse mesh
1053 27832 : ngfftc(:) = 2
1054 :
1055 : ! In case of PAW, compute fine FFT parameters
1056 6958 : if (usepaw==1) then
1057 1180 : ecutdg_eff=dtsets(idtset)%pawecutdg*dtsets(idtset)%dilatmx**2
1058 22420 : ngfftdg(:)=dtsets(idtset)%ngfftdg(:)
1059 1180 : call wrtout(std_out,'getng is called for the fine grid:')
1060 : ! Start with the coarse mesh as an initial guess for the fine mesh
1061 : ! This ensures that the fine mesh will not be any coarser than the coarse mesh in each dimension
1062 4720 : ngfftc(:) = ngfft(1:3)
1063 1180 : kpt=k0; if (response==1.and.usepaw==1) kpt=qphon ! this is temporary
1064 :
1065 : call getng(dtsets(idtset)%bxctmindg,dtsets(idtset)%chksymtnons,&
1066 : & ecutdg_eff,gmet,kpt,me_fft,mgfftdg,&
1067 : & nfftdg,ngfftdg,nproc_fft,nsym,paral_fft,symrel,dtsets(idtset)%tnons,ngfftc,&
1068 1180 : & gpu_option=dtsets(idtset)%gpu_option)
1069 :
1070 22420 : dtsets(idtset)%ngfftdg(:)=ngfftdg(:)
1071 1180 : dtsets(idtset)%mgfftdg=mgfftdg
1072 1180 : dtsets(idtset)%nfftdg=nfftdg
1073 : ! Compute fft distribution for fine grid
1074 1180 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=mod(fftalg,10)
1075 1180 : call mpi_enregs(idtset)%distribfft%init('f', mpi_enregs(idtset)%nproc_fft,ngfftdg(2),ngfftdg(3))
1076 : end if
1077 :
1078 6958 : dtsets(idtset)%mpw=mpw
1079 6958 : ABI_FREE(symrel)
1080 6958 : ABI_FREE(istwfk)
1081 6958 : ABI_FREE(kpt_with_shift)
1082 6958 : ABI_FREE(intarr)
1083 6958 : ABI_FREE(dprarr)
1084 :
1085 : ! Initialize data for the parallelization over atomic sites (PAW)
1086 6958 : if (dtsets(idtset)%natom==1) dtsets(idtset)%paral_atom=0
1087 6958 : if (dtsets(idtset)%usepaw==0) dtsets(idtset)%paral_atom=0
1088 6958 : if (dtsets(idtset)%usewvl/=0) dtsets(idtset)%paral_atom=0
1089 6958 : if (dtsets(idtset)%usedmft==1) dtsets(idtset)%paral_atom=0
1090 6958 : if (dtsets(idtset)%usedmft==10) dtsets(idtset)%paral_atom=0
1091 6958 : if (optdriver/=RUNL_GSTATE.and.optdriver/=RUNL_RESPFN.and.optdriver/=RUNL_GWLS) dtsets(idtset)%paral_atom=0
1092 6958 : if (dtsets(idtset)%macro_uj/=0) dtsets(idtset)%paral_atom=0
1093 :
1094 6958 : call initmpi_atom(dtsets(idtset),mpi_enregs(idtset))
1095 :
1096 : ! In case of the use of a GPU (Cuda), some defaults can change
1097 : ! according to a threshold on matrix sizes
1098 6958 : if (dtsets(idtset)%gpu_option==ABI_GPU_LEGACY) then
1099 0 : if (optdriver==RUNL_GSTATE.or.optdriver==RUNL_GWLS) then
1100 0 : vectsize=dtsets(idtset)%mpw*dtsets(idtset)%nspinor/dtsets(idtset)%npspinor
1101 0 : if (all(dtsets(idtset)%istwfk(:)==2)) vectsize=2*vectsize
1102 0 : blocksize=dtsets(idtset)%npband*dtsets(idtset)%bandpp
1103 0 : if (dtsets(idtset)%paral_kgb==0) blocksize=dtsets(idtset)%npfft
1104 0 : if ((vectsize*blocksize**2)>=dtsets(idtset)%gpu_linalg_limit) then
1105 0 : if (.not.wfoptalg_read) then
1106 0 : dtsets(idtset)%wfoptalg=14
1107 0 : if (tread(8)==1) then ! bandpp read, so set nblock_lobpcg
1108 0 : dtsets(idtset)%nblock_lobpcg=mband_upper/(dtsets(idtset)%bandpp*dtsets(idtset)%npband)
1109 : else ! bandpp set by nblock_lobpcg
1110 0 : dtsets(idtset)%bandpp=mband_upper/(dtsets(idtset)%nblock_lobpcg*dtsets(idtset)%npband)
1111 : end if
1112 0 : if (.not.fftalg_read) then
1113 0 : forbid_threads=(nthreads>1.and.dtsets(idtset)%bandpp>1)
1114 0 : dtsets(idtset)%ngfft(7) = fftalg_for_npfft(dtsets(idtset)%npfft,forbid_threads=forbid_threads)
1115 0 : if (usepaw==1) dtsets(idtset)%ngfftdg(7) = fftalg_for_npfft(dtsets(idtset)%npfft,forbid_threads=forbid_threads)
1116 : end if
1117 0 : if (.not.ortalg_read) dtsets(idtset)%ortalg=-abs(dtsets(idtset)%ortalg)
1118 : end if
1119 : end if
1120 : end if
1121 : end if
1122 :
1123 : ! initialize data for the parallelization for WVL:
1124 15344 : if(dtsets(idtset)%usewvl==1) then
1125 0 : mpi_enregs(idtset)%comm_wvl=mpi_enregs(idtset)%comm_cell
1126 0 : mpi_enregs(idtset)%nproc_wvl=xmpi_comm_size(mpi_enregs(idtset)%comm_wvl)
1127 0 : mpi_enregs(idtset)%me_wvl=xmpi_comm_rank(mpi_enregs(idtset)%comm_wvl)
1128 : end if
1129 :
1130 : end do
1131 :
1132 : !This is not a very clean exit in case of paral_kgb<0
1133 1428 : if (iexit/=0)then
1134 1 : call wrtout(std_out, "Stopping now!")
1135 1 : call abi_abort('PERS', exit_status=0, print_config=.false.)
1136 : end if
1137 :
1138 : DBG_EXIT("COLL")
1139 :
1140 2855 : end subroutine mpi_setup
1141 : !!***
1142 :
1143 : !!****f* ABINIT/finddistrproc
1144 : !! NAME
1145 : !! finddistrproc
1146 : !!
1147 : !! FUNCTION
1148 : !! Given a total number of processors, find a suitable distribution
1149 : !! that fill all the different levels of parallelization
1150 : !! (npimage, nppert, np_spkpt, npspinor, npband, npfft, bandpp)
1151 : !! Also determine parameters of parallel Linear Algebra routines
1152 : !! (use_slk, np_slk, gpu_linalg_limit)
1153 : !!
1154 : !! INPUTS
1155 : !! dtsets(0:ndtset_alloc)=<type datafiles_type>contains all input variables,
1156 : !! for all datasets; at this stage only datasets with index lower than
1157 : !! idtset are already initialized
1158 : !! filnam(5)=character strings giving file names
1159 : !! idtset=number of the current dataset
1160 : !! mpi_enreg=information about MPI parallelization
1161 : !! mband=maximum number of bands.
1162 : !! ndtset_alloc=number of datasets, corrected for allocation of at least one data set
1163 : !! tread(11)=flags indicating whether parallel input parameters were read from input file
1164 : !! tread(1) : paral_kgb tread(6) : npfft
1165 : !! tread(2) : npimage tread(7) : npband
1166 : !! tread(3) : nppert tread(8) : bandpp
1167 : !! tread(4) : np_spkpt tread(9) : use_slk
1168 : !! tread(5) : nspinor tread(10): np_slk
1169 : !! tread(11) : gpu_linalg_limit
1170 : !!
1171 : !! SIDE EFFECTS
1172 : !! iexit= if incremented, an exit is required
1173 : !! dtset%paral_kgb= flag for band-fft parallelism
1174 : !! dtset%npimage = number of processors for parallelisation over image
1175 : !! dtset%nppert = number of processors for parallelisation over perturbations
1176 : !! dtset%npspinor = number of processors for parallelisation on spinor components
1177 : !! dtset%np_spkpt = number of processors for parallelisation on spin / k points
1178 : !! dtset%npfft = number of processors for parallelisation on fft grid
1179 : !! dtset%npband = number of processors for parallelisation on bands
1180 : !! dtset%nphf = number of processors for parallelisation on occupied states for fock exchange
1181 : !! dtset%bandpp = internal parameter for lobpcg parallelisation algorithm
1182 : !! dtset%use_slk = flag for ScalaPAck use
1183 : !! dtset%np_slk = number of processors used in ScaLapack routines
1184 : !! dtset%gpu_linalg_limit=threshold activating Linear Algebra on GPU
1185 : !!
1186 : !! SOURCE
1187 :
1188 6444 : subroutine finddistrproc(dtsets,filnam,idtset,iexit,mband,mpi_enreg,ndtset_alloc,tread)
1189 :
1190 : !Arguments ------------------------------------
1191 : !scalars
1192 : integer,intent(in) :: idtset,mband,ndtset_alloc
1193 : integer,intent(inout) :: iexit
1194 : type(dataset_type),intent(inout),target :: dtsets(0:ndtset_alloc)
1195 : type(MPI_type),intent(inout) :: mpi_enreg
1196 : !arrays
1197 : integer,intent(in) :: tread(11)
1198 : character(len=fnlen),intent(in) :: filnam(5)
1199 :
1200 : !Local variables-------------------------------
1201 : !scalars
1202 : !128 should be a reasonable maximum for npfft (scaling is very poor for npfft>20)
1203 : integer,parameter :: ALGO_NOT_SET=-1, ALGO_DEFAULT_PAR=2
1204 : integer,parameter :: ALGO_CG=0, ALGO_LOBPCG_OLD=1, ALGO_LOBPCG_NEW=2, ALGO_CHEBFI=3, ALGO_CHEBFI_NEW=4
1205 : integer,parameter :: NPFMAX=128,BLOCKSIZE_MAX=3000,MAXBAND_PRINT=10
1206 : integer,parameter :: MAXCOUNT=250,MAXPRINT=10,MAXBENCH=25,MAXABIPY=25,NPF_CUTOFF=20
1207 : real(dp),parameter :: relative_nband_range=0.025
1208 : integer :: wf_algo,wf_algo_global,bpp,bpp_max,bpp_min,optdriver,autoparal,nblocks,blocksize
1209 : integer :: npi_max,npi_min,npc,npc_max,npc_min
1210 : integer :: np_sk,np_sk_max,np_sk_min,npp_max,npp_min
1211 : integer :: nps,nps_max,nps_min,npf,npf_max,npf_min
1212 : integer :: npb,npb_max,npb_min,max_ncpus,ount,paral_kgb
1213 : integer :: work_size,nks_per_proc,tot_ncpus
1214 : integer :: ib1,ib2,ibest,icount,ii,imin,jj,kk,mcount,mcount_eff,mpw
1215 : integer :: n2,n3,ncell_eff,ncount,nimage_eff,nkpt_eff,npert_eff
1216 : integer :: nproc,nproc1,nprocmin,np_slk,nthreads,use_linalg_gpu,omp_ncpus
1217 : logical :: dtset_found,file_found,first_bpp,iam_master
1218 : logical :: with_image,with_pert,with_kpt,with_spinor,with_fft,with_band,with_bandpp,with_thread
1219 : real(dp):: acc_c,acc_k,acc_kgb,acc_kgb_0,acc_s,ecut_eff,eff,ucvol,weight0
1220 : character(len=10) :: suffix
1221 : character(len=20) :: strg
1222 : character(len=500) :: msg,msgttl
1223 : character(len=fnlen) :: filden
1224 6444 : type(hdr_type) :: hdr0
1225 : !arrays
1226 : integer :: idum(1),idum3(3),ngmax(3),ngmin(3)
1227 6444 : integer,allocatable :: nband_best(:),isort(:),jdtset_(:)
1228 6444 : integer,allocatable :: my_algo(:),my_distp(:,:),nproc_best(:)
1229 6444 : integer,pointer :: nkpt_rbz(:)
1230 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1231 6444 : real(dp),allocatable :: weight(:)
1232 6444 : real(dp),pointer :: nband_rbz(:,:)
1233 : type(dataset_type),pointer :: dtset
1234 :
1235 : !******************************************************************
1236 :
1237 : DBG_ENTER("COLL")
1238 :
1239 : !Select current dataset
1240 6444 : dtset => dtsets(idtset)
1241 :
1242 : !Is automatic parallelization activated?
1243 6444 : autoparal = dtset%autoparal
1244 6444 : if (autoparal==0) return
1245 :
1246 : !Is it available
1247 100 : if ((dtset%usefock==1).AND.(dtset%nphf/=1)) then
1248 0 : ABI_ERROR("autoparal>0 not available for Hartree-Fock or hybrid XC calculations!")
1249 : end if
1250 100 : if ((autoparal>1).and.dtset%wfoptalg/=4.and.dtset%wfoptalg/=14) then
1251 0 : ABI_ERROR("autoparal>1 only available for the old LOBPCG algorithm (wfoptalg=4/14)!")
1252 : end if
1253 :
1254 : ! Unit number used for outputting the autoparal sections
1255 100 : ount = ab_out
1256 :
1257 : ! From the documentation:
1258 : !
1259 : ! If autoparal > 1 and max_ncpus is greater than 0, ABINIT analyzes the
1260 : ! efficiency of the process distribution for each possible number of processors
1261 : ! from 2 to max_ncpus. After having printed out the efficiency, the code stops.
1262 :
1263 : ! Handy local variables
1264 100 : iam_master = (mpi_enreg%me==0)
1265 100 : optdriver = dtset%optdriver
1266 100 : max_ncpus = dtset%max_ncpus ; if (dtset%paral_kgb<0) max_ncpus=abs(dtset%paral_kgb)
1267 100 : nthreads=xomp_get_max_threads()
1268 100 : nproc=mpi_enreg%nproc
1269 100 : if (max_ncpus>0) nproc = dtset%max_ncpus/nthreads
1270 : if (xmpi_paral==0.and.max_ncpus<=0) nproc=1
1271 :
1272 100 : nprocmin=2
1273 100 : if (xmpi_paral==1.and.max_ncpus<=0) nprocmin=max(2,nproc-100)
1274 100 : if (max_ncpus>0.and.autoparal/=0) nprocmin=1
1275 :
1276 100 : wf_algo_global=ALGO_NOT_SET
1277 100 : if (dtset%wfoptalg==0.and.tread(1)==1) wf_algo_global=ALGO_CG
1278 100 : if (dtset%wfoptalg==4.or.dtset%wfoptalg==14) wf_algo_global=ALGO_LOBPCG_OLD
1279 100 : if (dtset%wfoptalg==114) wf_algo_global=ALGO_LOBPCG_NEW
1280 100 : if (dtset%wfoptalg==1) wf_algo_global=ALGO_CHEBFI
1281 100 : if (dtset%wfoptalg==111) wf_algo_global=ALGO_CHEBFI_NEW
1282 :
1283 : ! Some peculiar cases (with direct exit)
1284 : ! MG: What is the meaning of max_ncpus < 0. This is not documented!
1285 100 : if (max_ncpus<=0) then
1286 99 : if (nproc==1.and.max_ncpus<=0) then
1287 5 : if (tread(1)==0.or.xmpi_paral==0) dtset%paral_kgb= 0
1288 5 : if (tread(2)==0.or.xmpi_paral==0) dtset%npimage = 1
1289 5 : if (tread(3)==0.or.xmpi_paral==0) dtset%nppert = 1
1290 5 : if (tread(4)==0.or.xmpi_paral==0) dtset%npspinor = 1
1291 5 : if (tread(5)==0.or.xmpi_paral==0) dtset%np_spkpt = 1
1292 5 : if (tread(6)==0.or.xmpi_paral==0) dtset%npfft = 1
1293 5 : if (tread(7)==0.or.xmpi_paral==0) dtset%npband = 1
1294 5 : if (tread(8)==0.or.xmpi_paral==0) dtset%bandpp = 1
1295 5 : if (tread(9)==0.or.xmpi_paral==0) dtset%use_slk = 0
1296 5 : if (tread(10)==0.or.xmpi_paral==0) dtset%np_slk = 1000000
1297 5 : return
1298 : end if
1299 94 : if ((optdriver/=RUNL_GSTATE.and. optdriver/=RUNL_RESPFN.and. optdriver/=RUNL_GWLS).or. &
1300 : (optdriver==RUNL_GSTATE.and.dtset%usewvl==1)) then
1301 0 : dtset%paral_kgb= 0
1302 0 : dtset%npimage = max(1,dtset%npimage)
1303 0 : dtset%nppert = max(1,dtset%nppert)
1304 0 : dtset%npspinor = max(1,dtset%npspinor)
1305 0 : dtset%np_spkpt = max(1,dtset%np_spkpt)
1306 0 : dtset%npfft = max(1,dtset%npfft)
1307 0 : dtset%npband = max(1,dtset%npband)
1308 0 : dtset%bandpp = max(1,dtset%bandpp)
1309 0 : return
1310 : end if
1311 : end if
1312 :
1313 : ! Need the metric tensor
1314 95 : call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd)
1315 95 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1316 :
1317 : ! Determine some quantities related to plane waves
1318 : ! - Crude estimation of the number of PW
1319 : ! - Number of G vectors in each direction
1320 95 : mpw=0;ngmin=0;ngmax=0
1321 95 : if (optdriver==RUNL_GSTATE) then
1322 75 : ecut_eff = dtset%ecut*dtset%dilatmx**2
1323 75 : mpw = nint(ucvol*((two*ecut_eff)**1.5_dp)/(six*pi**2)) ! Crude estimation
1324 79 : if (all(dtset%istwfk(1:dtset%nkpt)>1)) mpw=mpw/2+1
1325 75 : call kpgcount(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,dtset%kpt,ngmax,ngmin,dtset%nkpt)
1326 75 : write(msg,'(a,i0)') ' getmpw sequential formula gave: ',mpw
1327 75 : call wrtout(std_out,msg)
1328 : end if
1329 :
1330 : ! Parallelization over images
1331 95 : npi_min=1;npi_max=1;nimage_eff=1
1332 95 : if (optdriver==RUNL_GSTATE) then
1333 75 : nimage_eff=dtset%ndynimage
1334 75 : if (dtset%ntimimage<=1) nimage_eff=dtset%nimage
1335 75 : npi_min=max(1,dtset%npimage)
1336 75 : npi_max=min(nproc,nimage_eff)
1337 75 : if (tread(2)==1) npi_max=dtset%npimage
1338 : end if
1339 :
1340 : !Parallelization over k-points and spin components (GS)
1341 95 : np_sk_min=1;np_sk_max=1;nkpt_eff=0
1342 95 : if (optdriver==RUNL_GSTATE) then
1343 75 : nkpt_eff=dtset%nkpt*dtset%nsppol
1344 75 : np_sk_min=max(1,dtset%np_spkpt)
1345 75 : np_sk_max=min(nproc,nkpt_eff)
1346 75 : if (tread(4)==1) np_sk_max=dtset%np_spkpt
1347 : end if
1348 :
1349 : !Parallelization over perturbations, k-points and spin components (DFPT)
1350 95 : npp_min=1;npp_max=1;npert_eff=1
1351 95 : if (any(optdriver == [RUNL_RESPFN, RUNL_LONGWAVE])) then
1352 20 : if (dtset%paral_rf==1) then
1353 20 : call dtset%get_npert_rbz(nband_rbz, nkpt_rbz, npert_eff)
1354 76 : do jj=1,npert_eff
1355 10928 : ii=dtset%nsppol*nkpt_rbz(jj)*maxval(nband_rbz(:,jj))
1356 76 : nkpt_eff=max(nkpt_eff,ii)
1357 : end do
1358 20 : npp_min=max(1,dtset%nppert)
1359 20 : npp_max=min(nproc,npert_eff)
1360 20 : if (tread(3)==1) then
1361 20 : npp_max=dtset%nppert
1362 20 : if (npp_max>npert_eff) then
1363 4 : npp_min=npert_eff;npp_max=npert_eff
1364 4 : ABI_WARNING('nppert is bigger than npert; we set nppert=npert')
1365 : end if
1366 : end if
1367 20 : np_sk_min=1
1368 20 : np_sk_max=min(nproc,nkpt_eff)
1369 20 : ABI_FREE(nkpt_rbz)
1370 20 : ABI_FREE(nband_rbz)
1371 : else
1372 0 : nkpt_eff=nproc
1373 0 : np_sk_min=nproc-5
1374 0 : np_sk_max=nproc
1375 : end if
1376 : end if
1377 :
1378 : !Parallelization over spinorial components
1379 95 : nps_min=1;nps_max=1
1380 95 : if (optdriver==RUNL_GSTATE) then
1381 75 : nps_min=max(1,dtset%npspinor)
1382 75 : nps_max=min(nproc,dtset%nspinor)
1383 75 : if (tread(5)==1) nps_max=dtset%npspinor
1384 : end if
1385 :
1386 : !KGB Parallelization
1387 :
1388 95 : npf_min=1;npf_max=1
1389 95 : npb_min=1;npb_max=1
1390 95 : bpp_min=1;bpp_max=1
1391 95 : n2=0;n3=0
1392 95 : if (optdriver==RUNL_GSTATE) then
1393 :
1394 : ! >> FFT level
1395 75 : npf_min=max(1,dtset%npfft)
1396 75 : npf_min=min(npf_min,ngmin(2))
1397 75 : npf_max=min(nproc,NPFMAX)
1398 75 : if (tread(6)==1) then
1399 0 : npf_max=dtset%npfft
1400 0 : if (npf_max>ngmin(2)) then
1401 : write(msg,'(3a)') &
1402 0 : "Value of npfft given in input file is too high for the FFT grid!",ch10,&
1403 0 : "Action: decrease npfft or increase FFT grid (ecut, ngfft, ...)."
1404 0 : ABI_ERROR(msg)
1405 : end if
1406 : end if
1407 75 : npf_max=min(npf_max,ngmin(2))
1408 : ! Deactivate MPI FFT parallelism for GPU
1409 75 : if (dtset%gpu_option/=ABI_GPU_DISABLED) then
1410 0 : npf_min=1;npf_max=1
1411 : end if
1412 : !Deactivate MPI FFT parallelism for GPU
1413 75 : if (tread(1)==1.and.dtset%paral_kgb==0) then
1414 18 : npf_min=1;npf_max=1
1415 : end if
1416 : !Deactivate MPI FFT parallelism for multi-threaded LOBPCG / CHEBFI
1417 75 : if ((wf_algo_global==ALGO_LOBPCG_NEW.or.wf_algo_global==ALGO_CHEBFI.or.wf_algo_global==ALGO_CHEBFI_NEW).and.nthreads>1) then
1418 0 : npf_min=1;npf_max=1
1419 : end if
1420 :
1421 : ! Number of FFT procs has to be a multiple of FFT grid sizes
1422 : ! In case of a restart from a density file, it has to be
1423 : ! compatible with the FFT grid used for the density
1424 75 : n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
1425 75 : if (n2==0.and.n3==0) then
1426 : if (dtset%getden/=0.or.dtset%irdden/=0.or.&
1427 39 : & dtset%getkden/=0.or.dtset%irdkden/=0.or.dtset%iscf<0) then
1428 6 : dtset_found=.false.;file_found=.false.
1429 : !1-Try to find ngfft from previous dataset
1430 6 : if (dtset%getden/=0.or.dtset%getkden/=0) then
1431 30 : do ii=1,ndtset_alloc
1432 24 : jj=dtset%getden;if (jj==0) jj=dtset%getkden
1433 24 : if (jj<0) jj=dtset%jdtset+jj
1434 30 : if (dtsets(ii)%jdtset==jj) then
1435 6 : dtset_found=.true.
1436 6 : n2=dtsets(ii)%ngfftdg(2);n3=dtsets(ii)%ngfftdg(3)
1437 : end if
1438 : end do
1439 : end if
1440 : !2-If not found, try to extract ngfft from density file
1441 6 : if (.not.dtset_found) then
1442 : !Retrieve file name
1443 0 : if (dtset%getden/=0.or.dtset%irdden/=0) then
1444 0 : suffix='_DEN';if (dtset%nimage>1) suffix='_IMG1_DEN'
1445 0 : else if (dtset%getkden/=0.or.dtset%irdkden/=0) then
1446 0 : suffix='_KDEN';if (dtset%nimage>1) suffix='_IMG1_KDEN'
1447 : end if
1448 0 : ABI_MALLOC(jdtset_,(0:ndtset_alloc))
1449 0 : jdtset_=0;if(ndtset_alloc/=0) jdtset_(0:ndtset_alloc)=dtsets(0:ndtset_alloc)%jdtset
1450 0 : call mkfilename(filnam,filden,dtset%getden,idtset,dtset%irdden,jdtset_,ndtset_alloc,suffix,'den',ii)
1451 0 : ABI_FREE(jdtset_)
1452 : !Retrieve ngfft from file header
1453 0 : idum3=0
1454 0 : if (mpi_enreg%me==0) then
1455 0 : inquire(file=trim(filden),exist=file_found)
1456 0 : if (file_found) then
1457 0 : call hdr0%from_fname(filden,ii,xmpi_comm_self)
1458 0 : idum3(1:2)=hdr0%ngfft(2:3);if (file_found) idum3(3)=1
1459 0 : call hdr0%free()
1460 0 : ABI_WARNING("Cannot find filden "//filden)
1461 : end if
1462 : end if
1463 0 : call xmpi_bcast(idum3,0,mpi_enreg%comm_world,ii)
1464 0 : n2=idum3(1);n3=idum3(2);file_found=(idum3(3)/=0)
1465 : end if
1466 : end if
1467 : end if
1468 :
1469 : ! >> Band level
1470 75 : npb_min=max(1,dtset%npband)
1471 75 : npb_max=min(nproc,mband)
1472 75 : if (tread(7)==1) npb_max=dtset%npband
1473 75 : if (tread(1)==1.and.dtset%paral_kgb==0) then
1474 18 : npb_min=1;npb_max=1
1475 : end if
1476 :
1477 : ! >> banddp level
1478 75 : if (tread(8)==1) then
1479 0 : bpp_min = dtset%bandpp
1480 : else
1481 75 : bpp_min = 1
1482 : end if
1483 75 : bpp_max=mband
1484 75 : if (wf_algo_global==ALGO_LOBPCG_OLD) bpp_max=max(4,nint(mband/10.)) ! reasonable bandpp max
1485 75 : if (tread(8)==1) bpp_max=dtset%bandpp
1486 75 : if (wf_algo_global==ALGO_CHEBFI) bpp_min=1 ! bandpp not used with ChebFi
1487 0 : if (wf_algo_global==ALGO_CHEBFI) bpp_max=1
1488 75 : if (wf_algo_global==ALGO_CHEBFI_NEW) bpp_min=1 ! bandpp not used with ChebFi
1489 0 : if (wf_algo_global==ALGO_CHEBFI_NEW) bpp_max=1 ! bandpp not used with ChebFi
1490 :
1491 : end if ! RUNL_GSTATE
1492 :
1493 : !Disable KGB parallelisation in some cases:
1494 : ! - no GS
1495 : ! - paral_kgb=0 present in input file
1496 : ! - nstep=0
1497 : ! - Hartree-Fock or hybrid calculation (for now on)
1498 : if ( (optdriver/=RUNL_GSTATE).or.(dtset%paral_kgb==0.and.tread(1)==1).or. &
1499 95 : (dtset%nstep==0).or.(dtset%usefock==1)) then
1500 38 : nps_min=1; nps_max=1
1501 38 : npf_min=1; npf_max=1
1502 38 : npb_min=1; npb_max=1
1503 38 : bpp_min=1; bpp_max=1
1504 : end if
1505 :
1506 : ! Which levels of parallelism do we have?
1507 95 : with_image =(npi_min/=1.or.npi_max/=1)
1508 95 : with_pert =(npp_min/=1.or.npp_max/=1)
1509 95 : with_kpt =(np_sk_min/=1.or.np_sk_max/=1)
1510 95 : with_spinor=(nps_min/=1.or.nps_max/=1)
1511 95 : with_fft =(npf_min/=1.or.npf_max/=1)
1512 95 : with_band =(npb_min/=1.or.npb_max/=1)
1513 95 : with_bandpp=(bpp_min/=1.or.bpp_max/=1)
1514 95 : with_thread=(nthreads>1)
1515 :
1516 : !Allocate lists
1517 95 : ABI_MALLOC(my_distp,(10,MAXCOUNT))
1518 95 : ABI_MALLOC(weight,(MAXCOUNT))
1519 95 : ABI_MALLOC(my_algo,(MAXCOUNT))
1520 213845 : my_distp(1:7,:)=0;weight(:)=zero
1521 23845 : my_distp(8,:)=dtset%use_slk
1522 23845 : my_distp(9,:)=dtset%np_slk
1523 23845 : my_distp(10,:)=dtset%gpu_linalg_limit
1524 23845 : my_algo(:)=wf_algo_global
1525 95 : icount=0;imin=1
1526 :
1527 : !Cells= images or perturbations
1528 : npc_min=1;npc_max=1;ncell_eff=1
1529 : if (optdriver==RUNL_GSTATE) then
1530 95 : ncell_eff=nimage_eff;npc_min=npi_min;npc_max=npi_max
1531 : end if
1532 95 : if (any(optdriver == [RUNL_RESPFN, RUNL_LONGWAVE])) then
1533 : ncell_eff=npert_eff;npc_min=npp_min;npc_max=npp_max
1534 : end if
1535 :
1536 : !Loop over all possibilities
1537 : !Computation of weight~"estimated acceleration"
1538 : !================================================================
1539 :
1540 : !Cells= images or perturbations
1541 95 : npc_min=1;npc_max=1;ncell_eff=1
1542 95 : if (optdriver==RUNL_GSTATE) then
1543 75 : ncell_eff=nimage_eff;npc_min=npi_min;npc_max=npi_max
1544 : end if
1545 95 : if (any(optdriver == [RUNL_RESPFN, RUNL_LONGWAVE])) then
1546 20 : ncell_eff=npert_eff;npc_min=npp_min;npc_max=npp_max
1547 : end if
1548 :
1549 : !>>>>> CELLS
1550 232 : do npc=npc_min,npc_max
1551 137 : acc_c=one;if (npc>1) acc_c=0.99_dp*speedup_fdp(ncell_eff,npc)
1552 :
1553 : ! >>>>> K-POINTS
1554 669 : do np_sk=np_sk_min,np_sk_max
1555 : ! -> for DFPT runs, impose that nsppol divides np_sk
1556 437 : if (any(optdriver == [RUNL_RESPFN, RUNL_LONGWAVE]) .and. modulo(np_sk,dtset%nsppol)>0.and.np_sk>1) cycle
1557 437 : acc_k=one;if (np_sk>1) acc_k=0.96_dp*speedup_fdp(nkpt_eff,np_sk)
1558 :
1559 : ! >>>>> SPINORS
1560 1011 : do nps=nps_min,nps_max
1561 437 : acc_s=one;if (nps>1) acc_s=0.85_dp*speedup_fdp(dtset%nspinor,nps)
1562 :
1563 : ! >>>>> FFT
1564 1736 : do npf=npf_min,npf_max
1565 : ! -> npf should divide ngfft if set (if unset, ngfft=0 so the modulo test is ok)
1566 862 : if((modulo(n2,npf)>0).or.(modulo(n3,npf)>0)) cycle
1567 : ! -> npf should be only divisible by 2, 3 or 5
1568 774 : ii=npf
1569 1021 : do while (modulo(ii,2)==0)
1570 247 : ii=ii/2
1571 : end do
1572 915 : do while (modulo(ii,3)==0)
1573 141 : ii=ii/3
1574 : end do
1575 778 : do while (modulo(ii,5)==0)
1576 4 : ii=ii/5
1577 : end do
1578 774 : if(ii/=1) cycle
1579 :
1580 : ! Change algo if npfft>1
1581 766 : wf_algo=wf_algo_global
1582 766 : if (optdriver==RUNL_GSTATE.and.npf>1.and. wf_algo_global==ALGO_NOT_SET) wf_algo=ALGO_DEFAULT_PAR
1583 :
1584 : ! FFT parallelism not compatible with multithreading
1585 498 : if (wf_algo==ALGO_LOBPCG_NEW.or.wf_algo==ALGO_CHEBFI.or.wf_algo==ALGO_CHEBFI_NEW) then
1586 346 : if (nthreads>1.and.npf>1) cycle
1587 : end if
1588 :
1589 : ! >>>>> BANDS
1590 4187 : do npb=npb_min,npb_max
1591 2984 : nproc1=npc*np_sk*nps*npf*npb
1592 2984 : if (nproc1<nprocmin) cycle
1593 2906 : if (nproc1>nproc) cycle
1594 819 : if (modulo(mband,npb)>0) cycle
1595 :
1596 : ! Change algo if npband>1
1597 682 : if (optdriver==RUNL_GSTATE.and.npb>1.and. wf_algo_global==ALGO_NOT_SET) wf_algo=ALGO_DEFAULT_PAR
1598 :
1599 : ! Base speedup
1600 682 : acc_kgb_0=one;if (npb*npf*nthreads>1) acc_kgb_0=0.7_dp*speedup_fdp(mpw,(npb*npf*nthreads))
1601 :
1602 682 : if (npb*npf>4.and.wf_algo==ALGO_LOBPCG_OLD) then
1603 : ! Promote npb=npf
1604 0 : acc_kgb_0=acc_kgb_0*min((one*npf)/(one*npb),(one*npb)/(one*npf))
1605 : ! Promote npf<=20
1606 0 : if (npf>20)then
1607 : acc_kgb_0=acc_kgb_0* &
1608 : & 0.2_dp+(one-0.2_dp)*(sin((pi*(npf-NPF_CUTOFF))/(one*(NPFMAX-NPF_CUTOFF))) &
1609 0 : & /((pi*(npf-NPF_CUTOFF))/(one*(NPFMAX-NPF_CUTOFF))))**2
1610 : end if
1611 : end if
1612 :
1613 682 : first_bpp=.true.
1614 53476 : do bpp=bpp_min,bpp_max
1615 :
1616 51932 : if (wf_algo==ALGO_LOBPCG_NEW) then
1617 51088 : blocksize=npb*bpp;nblocks=mband/blocksize
1618 51088 : if (modulo(bpp,nthreads)>0) cycle
1619 51088 : if ((bpp>1).and.(modulo(bpp,2)>0)) cycle
1620 26014 : if (modulo(mband,npb*bpp)>0) cycle
1621 844 : else if (wf_algo==ALGO_LOBPCG_OLD) then
1622 0 : blocksize=npb*bpp;nblocks=mband/blocksize
1623 0 : if (modulo(mband/npb,bpp)>0) cycle
1624 0 : if ((bpp>1).and.(modulo(bpp,2)>0)) cycle
1625 0 : if (one*npb*bpp >max(1.,mband/3.).and.(mband>30)) cycle
1626 0 : if (npb*npf<=4.and.(.not.first_bpp)) cycle
1627 844 : else if (wf_algo==ALGO_CHEBFI .or. wf_algo==ALGO_CHEBFI_NEW) then
1628 : !Nothing
1629 : else
1630 844 : if (bpp/=1.or.npb/=1) cycle
1631 : end if
1632 :
1633 2229 : first_bpp=.false.
1634 :
1635 0 : acc_kgb=acc_kgb_0
1636 : ! OLD LOBPCG: promote bpp*npb>mband/3
1637 2229 : if (wf_algo==ALGO_LOBPCG_OLD) then
1638 0 : if (npb*npf>4.and.mband>30) acc_kgb=acc_kgb*(one-(three*bpp*npb)/(one*mband))
1639 : end if
1640 : ! NEW LOBPCG: promote minimal number of blocks
1641 : ! promote block size <= BLOCKSIZE_MAX
1642 2229 : if (wf_algo==ALGO_LOBPCG_NEW) then
1643 2017 : acc_kgb=acc_kgb*(one-0.9_dp*dble(nblocks-1)/dble(mband-1))
1644 2017 : if (blocksize>BLOCKSIZE_MAX) acc_kgb=acc_kgb*max(0.1_dp,one-dble(blocksize)/dble(10*BLOCKSIZE_MAX))
1645 2017 : if (nthreads==1) then
1646 : ! Promote npband vs bandpp & npfft
1647 2017 : if (blocksize>1) acc_kgb=acc_kgb*(0.1_dp*bpp+0.9_dp-blocksize)/(one-blocksize)
1648 2017 : if (npb*npf>4.and.mband>100) acc_kgb=acc_kgb*(one-0.8_dp*((three*bpp*npb)/(one*mband)-one)**2)
1649 2017 : tot_ncpus=max(npb,npf);if (tot_ncpus==2) tot_ncpus=0
1650 2017 : acc_kgb=acc_kgb*(one-0.8_dp*((dble(npb)/dble(npf))-2_dp)**2/(tot_ncpus-2_dp)**2)
1651 2017 : eff=max(npf,20);acc_kgb=acc_kgb*(one-0.8_dp*min(one,(eff-20)**2))
1652 : end if
1653 : end if
1654 :
1655 : ! CHEBFI: promote npfft=npband and nband>=npfft
1656 2229 : if (wf_algo==ALGO_CHEBFI .or. wf_algo==ALGO_CHEBFI_NEW) then
1657 0 : if (npf>1) then
1658 0 : if (npb>npf) then
1659 0 : acc_kgb=acc_kgb*(one-0.8_dp*0.25_dp*((dble(npb)/dble(npf))-one)**2/(nproc1-one)**2)
1660 : else
1661 0 : acc_kgb=acc_kgb*(one-0.8_dp*nproc1**2*((dble(npb)/dble(npf))-one)**2/(nproc1-one)**2)
1662 : end if
1663 : end if
1664 : end if
1665 :
1666 : ! Resulting "weight"
1667 : ! weight0=acc_c*acc_k*acc_s*acc_kgb
1668 2229 : weight0=nproc1*(acc_c+acc_k+acc_s+acc_kgb)/(npc+np_sk+nps+(npf*npb))
1669 :
1670 : ! Store data
1671 2229 : icount=icount+1
1672 5213 : if (icount<=MAXCOUNT) then
1673 1706 : my_algo(icount)=merge(ALGO_CG,wf_algo,wf_algo==ALGO_NOT_SET)
1674 13648 : my_distp(1:7,icount)=(/npc,np_sk,nps,npf,npb,bpp,nproc1/)
1675 1706 : weight(icount)=weight0
1676 1706 : if (weight0<weight(imin)) imin=icount
1677 : else
1678 523 : if (weight0>weight(imin)) then
1679 393 : my_algo(imin)=merge(ALGO_CG,wf_algo,wf_algo==ALGO_NOT_SET)
1680 3144 : my_distp(1:7,imin)=(/npc,np_sk,nps,npf,npb,bpp,nproc1/)
1681 393 : weight(imin)=weight0
1682 99036 : idum=minloc(weight);imin=idum(1)
1683 : end if
1684 : end if
1685 :
1686 : end do ! bpp
1687 : end do ! npb
1688 : end do ! npf
1689 : end do ! nps
1690 : end do ! np_sk
1691 : end do ! npc
1692 :
1693 : !Compute number of selected distributions
1694 95 : mcount_eff=icount
1695 95 : mcount=min(mcount_eff,MAXCOUNT)
1696 :
1697 : !Stop if no solution found
1698 95 : if (mcount==0) then
1699 : ! Override here the 0 default value changed in indefo1
1700 0 : dtset%npimage = max(1,dtset%npimage)
1701 0 : dtset%nppert = max(1,dtset%nppert)
1702 0 : dtset%np_spkpt = max(1,dtset%np_spkpt)
1703 0 : dtset%npspinor = max(1,dtset%npspinor)
1704 0 : dtset%npfft = max(1,dtset%npfft)
1705 0 : dtset%npband = max(1,dtset%npband)
1706 0 : dtset%bandpp = max(1,dtset%bandpp)
1707 : write(msg,'(a,i0,2a,i0,a)') &
1708 0 : 'Your input dataset does not let Abinit find an appropriate process distribution with nCPUs=',nproc*nthreads,ch10, &
1709 0 : 'Try to comment all the np* vars and set max_ncpus=',nthreads*nproc,' to have advice on process distribution.'
1710 0 : ABI_WARNING(msg)
1711 0 : if (max_ncpus>0) call wrtout(ab_out,msg, do_flush=.True.)
1712 0 : iexit=iexit+1
1713 : end if
1714 :
1715 : !Sort data by increasing weight
1716 95 : if (mcount>0) then
1717 285 : ABI_MALLOC(isort,(mcount))
1718 5308 : isort=(/(ii,ii=1,mcount)/)
1719 95 : call sort_dp(mcount,weight,isort,tol6)
1720 95 : ncount=min(mcount,MAXPRINT)
1721 : end if
1722 :
1723 : !Deduce a global value for paral_kgb
1724 95 : paral_kgb=dtset%paral_kgb
1725 95 : if (tread(1)==0) then
1726 52 : if (any(my_algo(:)/=ALGO_CG)) paral_kgb=1
1727 : end if
1728 :
1729 : ! ======================================
1730 : ! Print output for abipy in Yaml format
1731 : ! ======================================
1732 :
1733 : ! Please DO NOT CHANGE this part without contacting gmatteo first
1734 : ! since ANY CHANGE can easily break the interface with AbiPy.
1735 95 : if (iam_master .and. max_ncpus > 0.and. (mcount>0 .or. wf_algo_global == ALGO_CG)) then
1736 1 : write(ount,'(2a)')ch10,"--- !Autoparal"
1737 1 : if (optdriver==RUNL_GSTATE .and. paral_kgb == 0) then
1738 0 : write(ount,"(a)")"# Autoparal section for GS run (band-by-band CG method)"
1739 1 : else if (optdriver==RUNL_GSTATE) then
1740 1 : write(ount,'(a)')'# Autoparal section for GS calculations with paral_kgb 1'
1741 0 : else if (optdriver==RUNL_RESPFN) then
1742 0 : write(ount,'(a)')'# Autoparal section for DFPT calculations'
1743 0 : else if (optdriver==RUNL_LONGWAVE) then
1744 0 : write(ount,'(a)')'# Autoparal section for LONGWAVE calculations'
1745 : else
1746 0 : ABI_ERROR(sjoin('Unsupported optdriver:', itoa(optdriver)))
1747 : end if
1748 1 : write(ount,"(a)") "info:"
1749 1 : write(ount,"(a,i0)")" autoparal: ",autoparal
1750 1 : write(ount,"(a,i0)")" paral_kgb: ",paral_kgb
1751 1 : write(ount,"(a,i0)")" max_ncpus: ",max_ncpus
1752 1 : write(ount,"(a,i0)")" nspinor: ",dtset%nspinor
1753 1 : write(ount,"(a,i0)")" nsppol: ",dtset%nsppol
1754 1 : write(ount,"(a,i0)")" nkpt: ",dtset%nkpt
1755 1 : write(ount,"(a,i0)")" mband: ",mband
1756 1 : write(ount,"(a)")"configurations:"
1757 :
1758 1 : if (optdriver==RUNL_GSTATE.and.paral_kgb==0) then
1759 0 : work_size = dtset%nkpt * dtset%nsppol
1760 0 : do ii=1,max_ncpus
1761 0 : if (ii > work_size) cycle
1762 0 : do omp_ncpus=1,nthreads
1763 0 : nks_per_proc = work_size / ii
1764 0 : nks_per_proc = nks_per_proc + MOD(work_size, ii)
1765 0 : eff = (one * work_size) / (ii * nks_per_proc)
1766 0 : write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
1767 0 : write(ount,"(a,i0)")" mpi_ncpus: ",ii
1768 0 : write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus
1769 0 : write(ount,"(a,f12.9)")" efficiency: ",eff
1770 : !write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
1771 : end do
1772 : end do
1773 :
1774 1 : else if (optdriver==RUNL_GSTATE) then
1775 1 : omp_ncpus=nthreads
1776 11 : do jj=mcount,mcount-min(ncount,MAXABIPY)+1,-1
1777 10 : ii=isort(jj)
1778 10 : tot_ncpus = my_distp(7,ii)
1779 10 : eff = weight(jj) / tot_ncpus
1780 10 : write(ount,'(a,i0)')' - tot_ncpus: ',tot_ncpus
1781 10 : write(ount,'(a,i0)')' mpi_ncpus: ',tot_ncpus
1782 10 : write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus
1783 10 : write(ount,'(a,f12.9)')' efficiency: ',eff
1784 : !write(ount,'(a,f12.2)')' mem_per_cpu: ',mempercpu_mb
1785 10 : write(ount,'(a)' )' vars: {'
1786 10 : write(ount,'(a,i0,a)')' npimage: ',my_distp(1,ii),','
1787 : ! Keep on using legacy npkpt instead of np_spkpt to maintain compatibility with AbiPy
1788 10 : write(ount,'(a,i0,a)')' npkpt: ',my_distp(2,ii),','
1789 : !write(ount,'(a,i0,a)')' np_spkpt: ',my_distp(2,ii),','
1790 10 : write(ount,'(a,i0,a)')' npspinor: ',my_distp(3,ii),','
1791 10 : write(ount,'(a,i0,a)')' npfft: ', my_distp(4,ii),','
1792 10 : write(ount,'(a,i0,a)')' npband: ',my_distp(5,ii),','
1793 10 : write(ount,'(a,i0,a)')' bandpp: ',my_distp(6,ii),','
1794 11 : write(ount,'(a)') ' }'
1795 : end do
1796 :
1797 0 : else if (any(optdriver == [RUNL_RESPFN, RUNL_LONGWAVE])) then
1798 0 : do jj=mcount,mcount-min(ncount,MAXABIPY)+1,-1
1799 0 : ii=isort(jj)
1800 0 : tot_ncpus = my_distp(7,ii)
1801 0 : eff = weight(jj) / tot_ncpus
1802 0 : write(ount,'(a,i0)')' - tot_ncpus: ',tot_ncpus
1803 0 : write(ount,'(a,i0)')' mpi_ncpus: ',tot_ncpus
1804 : !write(ount,'(a,i0)')' omp_ncpus: ',omp_ncpus !OMP not supported (yet)
1805 0 : write(ount,'(a,f12.9)')' efficiency: ',eff
1806 : !write(ount,'(a,f12.2)')' mem_per_cpu: ',mempercpu_mb
1807 0 : write(ount,'(a)' )' vars: {'
1808 0 : write(ount,'(a,i0,a)')' nppert: ', my_distp(1,ii),','
1809 : ! Keep on using legacy npkpt instead of np_spkpt to maintain compatibility with AbiPy
1810 0 : write(ount,'(a,i0,a)')' npkpt: ', my_distp(2,ii),','
1811 : !write(ount,'(a,i0,a)')' np_spkpt: ', my_distp(2,ii),','
1812 0 : write(ount,'(a)') ' }'
1813 : end do
1814 : end if
1815 1 : write(ount,'(a)')"..."
1816 : end if
1817 :
1818 : !Print out tab with selected choices
1819 95 : if (mcount>0.and.iam_master) then
1820 27 : if (nthreads==1) then
1821 27 : write(msg,'(a,1x,100("="),2a,i0,2a)') ch10,ch10,&
1822 54 : & ' Searching for all possible proc distributions for this input with #CPUs<=',nthreads*nproc,':',ch10
1823 : else
1824 0 : write(msg,'(a,1x,100("="),2a,i0,a,i0,2a)') ch10,ch10,&
1825 0 : & ' Searching for all possible proc distributions for this input with #CPUs<=',nthreads*nproc,&
1826 0 : & ' and ',nthreads,' openMP threads:',ch10
1827 : end if
1828 27 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1829 : !Titles of columns
1830 27 : msgttl='~'
1831 27 : if (with_image) msgttl=trim(msgttl)//'~~~~~~~~~~~'
1832 27 : if (with_pert) msgttl=trim(msgttl)//'~~~~~~~~~~~'
1833 27 : msgttl=trim(msgttl)//'~~~~~~~~~~~~~' ! kpt
1834 27 : if (with_spinor) msgttl=trim(msgttl)//'~~~~~~~~~~'
1835 27 : if (with_fft) msgttl=trim(msgttl)//'~~~~~~~~~~~~~'
1836 27 : if (with_band) msgttl=trim(msgttl)//'~~~~~~~~~~~~~'
1837 27 : if (with_bandpp) msgttl=trim(msgttl)//'~~~~~~~~~~~~~'
1838 27 : if (with_thread) msgttl=trim(msgttl)//'~~~~~~~~~~'
1839 27 : msgttl=trim(msgttl)//'~~~~~~~~~~~~~' ! nproc
1840 27 : if (with_thread) msgttl=trim(msgttl)//'~~~~~~~~~~~~~'
1841 27 : msgttl=trim(msgttl)//'~~~~~~~~~~~' ! CPUs
1842 27 : msgttl=' '//trim(msgttl)
1843 27 : call wrtout(std_out,msgttl);if(max_ncpus>0) call wrtout(ab_out,msgttl)
1844 27 : msg='|'
1845 27 : if (with_image) msg=trim(msg)//' npimage|'
1846 27 : if (with_pert) msg=trim(msg)//' nppert|'
1847 27 : msg=trim(msg)//' np_spkpt|'
1848 27 : if (with_spinor) msg=trim(msg)//' npspinor|'
1849 27 : if (with_fft) msg=trim(msg)//' npfft|'
1850 27 : if (with_band) msg=trim(msg)//' npband|'
1851 27 : if (with_bandpp) msg=trim(msg)//' bandpp|'
1852 27 : if (with_thread) msg=trim(msg)//' #Threads|'
1853 27 : msg=trim(msg)//' #MPI(proc)|'
1854 27 : if (with_thread) msg=trim(msg)//' #CPUs|'
1855 27 : msg=trim(msg)//' WEIGHT|'
1856 27 : msg=' '//trim(msg)
1857 27 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1858 27 : msg='|'
1859 27 : write(strg,'(i4,a,i4,a)') npi_min,'<<',npi_max,'|';if (with_image) msg=trim(msg)//trim(strg)
1860 27 : write(strg,'(i4,a,i4,a)') npp_min,'<<',npp_max,'|';if (with_pert) msg=trim(msg)//trim(strg)
1861 27 : write(strg,'(i5,a,i5,a)') np_sk_min,'<<',np_sk_max,'|'; msg=trim(msg)//trim(strg)
1862 27 : write(strg,'(i5,a,i2,a)') nps_min,'<<',nps_max,'|';if (with_spinor) msg=trim(msg)//trim(strg)
1863 27 : write(strg,'(i5,a,i5,a)') npf_min,'<<',npf_max,'|';if (with_fft) msg=trim(msg)//trim(strg)
1864 27 : write(strg,'(i5,a,i5,a)') npb_min,'<<',npb_max,'|';if (with_band) msg=trim(msg)//trim(strg)
1865 27 : write(strg,'(i5,a,i5,a)') bpp_min,'<<',bpp_max,'|';if (with_bandpp) msg=trim(msg)//trim(strg)
1866 27 : write(strg,'(i9,a)' ) nthreads ,'|';if (with_thread) msg=trim(msg)//trim(strg)
1867 27 : write(strg,'(i5,a,i5,a)') 1 ,'<<',nproc ,'|'; msg=trim(msg)//trim(strg)
1868 27 : write(strg,'(i4,a,i6,a)') nthreads,'<<',nthreads*nproc,'|';if (with_thread) msg=trim(msg)//trim(strg)
1869 27 : write(strg,'(a,i6,a)') ' <=',nthreads*nproc,'|'; msg=trim(msg)//trim(strg)
1870 27 : msg=' '//trim(msg)
1871 27 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1872 27 : call wrtout(std_out,msgttl);if(max_ncpus>0) call wrtout(ab_out,msgttl)
1873 : !Loop over selected choices
1874 217 : do jj=mcount,mcount-ncount+1,-1
1875 190 : ii=isort(jj)
1876 190 : msg='|'
1877 190 : write(strg,'(i10,a)') my_distp(1,ii),'|';if (with_image) msg=trim(msg)//trim(strg)
1878 190 : write(strg,'(i10,a)') my_distp(1,ii),'|';if (with_pert) msg=trim(msg)//trim(strg)
1879 190 : write(strg,'(i12,a)') my_distp(2,ii),'|'; msg=trim(msg)//trim(strg)
1880 190 : write(strg,'(i9,a)') my_distp(3,ii),'|';if (with_spinor) msg=trim(msg)//trim(strg)
1881 190 : write(strg,'(i12,a)') my_distp(4,ii),'|';if (with_fft) msg=trim(msg)//trim(strg)
1882 190 : write(strg,'(i12,a)') my_distp(5,ii),'|';if (with_band) msg=trim(msg)//trim(strg)
1883 190 : write(strg,'(i12,a)') my_distp(6,ii),'|';if (with_bandpp) msg=trim(msg)//trim(strg)
1884 190 : write(strg,'(i9,a)') nthreads ,'|';if (with_thread) msg=trim(msg)//trim(strg)
1885 190 : write(strg,'(i12,a)') my_distp(7,ii),'|'; msg=trim(msg)//trim(strg)
1886 190 : write(strg,'(i12,a)') nthreads*my_distp(7,ii),'|';if (with_thread) msg=trim(msg)//trim(strg)
1887 190 : write(strg,'(f10.3,a)') weight(jj) ,'|'; msg=trim(msg)//trim(strg)
1888 190 : msg=' '//trim(msg)
1889 217 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1890 : end do
1891 : !End of tab
1892 27 : call wrtout(std_out,msgttl);if(max_ncpus>0) call wrtout(ab_out,msgttl)
1893 27 : write(msg,'(a,i6,a,i6,a)')' Only the best possible choices for nproc are printed...'
1894 27 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1895 : end if ! mcount>0
1896 :
1897 : !Determine an optimal number of bands
1898 235 : if (optdriver==RUNL_GSTATE.and. &
1899 : & (any(my_algo(1:mcount)==ALGO_LOBPCG_OLD.or. &
1900 : & my_algo(1:mcount)==ALGO_LOBPCG_NEW.or. &
1901 : & my_algo(1:mcount)==ALGO_CHEBFI.or. &
1902 : & my_algo(1:mcount)==ALGO_CHEBFI_NEW))) then
1903 57 : if (mcount>0) then
1904 57 : icount=isort(mcount)
1905 57 : npc=my_distp(1,icount);np_sk=my_distp(2,icount)
1906 57 : nps=my_distp(3,icount);npf=my_distp(4,icount)
1907 : else
1908 0 : npc=1;if (with_image ) npc=npi_min
1909 0 : np_sk=1;if (with_kpt ) np_sk=np_sk_min
1910 0 : nps=1;if (with_spinor) nps=nps_min
1911 0 : npf=1;if (with_fft ) npf=npf_min
1912 : end if
1913 57 : nproc1=npc*np_sk*nps*npf
1914 57 : msg=ch10//' >>> Possible (best) choices for the number of bands (nband) are:'
1915 57 : if (with_image.or.with_kpt.or.with_spinor.or.with_fft) msg=trim(msg)//ch10//' with:'
1916 57 : write(strg,'(a,i0)') ' npimage=' ,npc;if (with_image) msg=trim(msg)//trim(strg)
1917 57 : write(strg,'(a,i0)') ' np_spkpt=' ,np_sk;if (with_kpt) msg=trim(msg)//trim(strg)
1918 57 : write(strg,'(a,i0)') ' npspinor=',nps;if (with_spinor) msg=trim(msg)//trim(strg)
1919 57 : write(strg,'(a,i0)') ' npfft=' ,npf;if (with_fft) msg=trim(msg)//trim(strg)
1920 57 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1921 57 : ib1=mband-int(mband*relative_nband_range);if (my_algo(icount)==ALGO_CHEBFI .or. my_algo(icount)==ALGO_CHEBFI_NEW) ib1=mband
1922 57 : ib2=mband+int(mband*relative_nband_range)
1923 171 : ABI_MALLOC(nproc_best,(1+ib2-ib1))
1924 114 : ABI_MALLOC(nband_best,(1+ib2-ib1))
1925 146 : nproc_best(:)=1
1926 381 : nband_best=(/(ii,ii=ib1,ib2)/)
1927 57 : bpp=merge(1,nthreads,my_algo(icount)==ALGO_CHEBFI .or. my_algo(icount)==ALGO_CHEBFI_NEW)
1928 146 : do ii=ib1,ib2
1929 778 : do jj=1,nproc/nproc1
1930 : ibest=1
1931 5272 : do kk=1,jj
1932 4640 : if (mod(jj,kk)/=0) cycle
1933 2434 : if (mod(ii,kk*bpp)==0) ibest=max(ibest,kk)
1934 : end do
1935 721 : nproc_best(1+ii-ib1)=max(nproc_best(1+ii-ib1),ibest)
1936 : end do
1937 : end do
1938 57 : call sort_int(1+ib2-ib1,nproc_best,nband_best)
1939 57 : kk=-1
1940 125 : do ii=1+ib2-ib1,max(ib2-ib1-MAXBAND_PRINT,1),-1
1941 68 : write(msg,'(3(a,i6),a,i3,a,i5,a)') ' nband=',nband_best(ii),' using ',nproc1*nproc_best(ii)*nthreads,&
1942 136 : & ' CPUs =',nproc1*nproc_best(ii),' MPI x',nthreads,' threads (npband=',nproc_best(ii),')'
1943 68 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1944 125 : if (nband_best(ii)==mband) kk=nproc_best(ii)
1945 : end do
1946 146 : if (kk==maxval(nproc_best(:))) then
1947 : if (my_algo(icount)/=ALGO_CHEBFI .or. my_algo(icount)/=ALGO_CHEBFI_NEW) then
1948 57 : write(msg,'(a,i6,a)') ' >>> The present nband value (',mband,') seems to be the best choice!'
1949 : end if
1950 57 : if (my_algo(icount)==ALGO_CHEBFI .or. my_algo(icount)/=ALGO_CHEBFI_NEW) then
1951 57 : write(msg,'(a,i6,a)') ' >>> The present nband value (',mband,') seems to be a good choice!'
1952 : end if
1953 57 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1954 : end if
1955 57 : ABI_FREE(nproc_best)
1956 57 : ABI_FREE(nband_best)
1957 : end if
1958 :
1959 1801 : if (optdriver==RUNL_GSTATE.and.(any(my_algo(1:mcount)==ALGO_CHEBFI .or. my_algo(1:mcount)==ALGO_CHEBFI_NEW))) then
1960 : write(msg,'(5a)') &
1961 0 : & ' >>> Note that with the "Chebyshev Filtering" algorithm, it is often',ch10,&
1962 0 : & ' better to increase the number of bands (10% more or a few tens more).',ch10,&
1963 0 : & ' Advice: increase nband and put nbdbuf input variable to (nband_new-nband_old).'
1964 0 : call wrtout(std_out,msg);if(max_ncpus>0) call wrtout(ab_out,msg)
1965 : end if
1966 :
1967 : !Refinement of the process distribution by mean of a LinAlg routines benchmarking
1968 95 : if (mcount>0.and.optdriver==RUNL_GSTATE.and.autoparal/=1) then
1969 0 : icount=isort(mcount)
1970 0 : if (autoparal/=3) then
1971 0 : if (autoparal==2) then
1972 0 : write(msg,'(5a,9(a10,a1))') ch10, &
1973 0 : & ' Values below have been tested with respect to Linear Algebra performance;',ch10,&
1974 0 : & ' Weights below are corrected according:',ch10,&
1975 0 : & 'npimage','|','np_spkpt' ,'|','npspinor' ,'|','npfft' ,'|','npband','|',' bandpp ' ,'|',&
1976 0 : & 'nproc' ,'|','weight','|','new weight','|'
1977 : else
1978 0 : write(msg,'(5a,11(a10,a1))') ch10, &
1979 0 : & ' Values below have been tested with respect to Linear Algebra performance;',ch10,&
1980 0 : & ' Weights below are corrected according:',ch10,&
1981 0 : & 'npimage','|','np_spkpt' ,'|','npspinor' ,'|','npfft' ,'|','npband','|',' bandpp ' ,'|',&
1982 0 : & 'nproc' ,'|','weight','|','new weight','|','best npslk','|','linalggpu' ,'|'
1983 : end if
1984 0 : call wrtout(std_out,msg);if (max_ncpus > 0) call wrtout(ab_out,msg)
1985 : end if
1986 0 : acc_k=zero
1987 0 : ncount=min(MAXBENCH,mcount);if (autoparal==3) ncount=1
1988 0 : do jj=mcount,mcount-ncount+1,-1
1989 0 : ii=isort(jj)
1990 0 : npf=my_distp(4,ii);npb=my_distp(5,ii);bpp=my_distp(6,ii)
1991 0 : if ((npb*npf*bpp>1).and.(npf*npb<=mpi_enreg%nproc)) then
1992 0 : use_linalg_gpu=dtset%gpu_option
1993 0 : call compute_kgb_indicator(acc_kgb,bpp,xmpi_world,mband,mpw,npb,npf,np_slk,use_linalg_gpu)
1994 0 : if (autoparal/=2) then
1995 0 : my_distp(9,ii)=np_slk
1996 0 : if (np_slk>0) my_distp(8,ii)=1
1997 : ! * gpu_linalg_limit:
1998 : ! No use of GPU: htgspw_01.outuge value ~2 *vectsize*blocksize**2 tested
1999 : ! Use of GPU: tiny value ~0.5*vectsize*blocksize**2 tested
2000 0 : my_distp(10,ii)=2*dtset%mpw*(npb*bpp)**2/npf
2001 0 : if (use_linalg_gpu/=ABI_GPU_DISABLED) my_distp(10,ii)=my_distp(10,ii)/4
2002 : end if
2003 0 : if (abs(acc_k)<=tol12) acc_k=acc_kgb ! Ref value : the first one computed
2004 : ! * Weight (corrected by 10% of the computed ratio)
2005 0 : weight0=weight(jj)*(one + 0.1_dp*acc_k/acc_kgb)
2006 0 : if (autoparal==2) then
2007 : write(msg, '(7(i10,a1),f9.2,a2,f9.5,a2)') &
2008 0 : & my_distp(1,ii),'|',my_distp(2,ii),'|',my_distp(3,ii),'|',my_distp(4,ii),'|',&
2009 0 : & my_distp(5,ii),'|',my_distp(6,ii),'|',my_distp(7,ii),'|',weight(jj),'=>', weight0,' |'
2010 0 : else if (autoparal==3) then
2011 0 : write(msg,'(a,5(a,i3))') ch10,' For npband=',npb,', npfft=',npf,' and bandpp=',bpp, &
2012 0 : & ', compute_kgb_indicator recommends you to set np_slk=',my_distp(9,ii),&
2013 0 : & ' and use_linalg_gpu=',use_linalg_gpu
2014 : else
2015 : write(msg, '(7(i10,a1),f9.2,a2,f9.5,a2,2(i10,a1))') &
2016 0 : & my_distp(1,ii),'|',my_distp(2,ii),'|',my_distp(3,ii),'|',my_distp(4,ii),'|',&
2017 0 : & my_distp(5,ii),'|',my_distp(6,ii),'|',my_distp(7,ii),'|',weight(jj),'=>', weight0,' |',&
2018 0 : & my_distp(9,ii),'|',use_linalg_gpu,'|'
2019 : end if
2020 0 : call wrtout(std_out,msg);if (max_ncpus>0) call wrtout(ab_out,msg)
2021 : ! We store the best value in weight(mcount) and keep icount
2022 0 : if (weight0 > weight(mcount)) then
2023 0 : icount=ii;weight(mcount)=weight0
2024 : end if
2025 : end if
2026 : end do
2027 : end if
2028 :
2029 : !Store new process distribution
2030 95 : if (mcount>0.and.max_ncpus<=0) then
2031 94 : icount=isort(mcount)
2032 94 : nproc1=my_distp(7,icount)
2033 : ! Work load distribution
2034 94 : if (optdriver==RUNL_GSTATE) then
2035 74 : dtset%npimage= my_distp(1,icount)
2036 74 : dtset%nppert = 1
2037 74 : dtset%np_spkpt = my_distp(2,icount)
2038 : end if
2039 94 : if (optdriver==RUNL_RESPFN) then
2040 20 : dtset%npimage= 1
2041 20 : dtset%nppert = my_distp(1,icount)
2042 20 : dtset%np_spkpt = 1
2043 : end if
2044 94 : dtset%npspinor = my_distp(3,icount)
2045 94 : dtset%npfft = my_distp(4,icount)
2046 94 : dtset%npband = my_distp(5,icount)
2047 94 : dtset%bandpp = my_distp(6,icount)
2048 138 : if (tread(1)==0) dtset%paral_kgb= merge(0,1,my_algo(icount)==ALGO_CG)
2049 : ! The following lines are mandatory : the DFT+DMFT must use ALL the
2050 : ! available procs specified by the user. So nproc1=nproc.
2051 : ! Works only if paral_kgb is not activated??
2052 94 : if (dtset%usedmft/=0.and.optdriver==RUNL_GSTATE) then
2053 0 : if (dtset%paral_kgb==0) then
2054 0 : dtset%npspinor = 1 ; dtset%npfft = 1
2055 0 : dtset%npband = 1 ; dtset%bandpp = 1
2056 0 : dtset%npimage = 1
2057 : end if
2058 0 : nproc1 = nproc
2059 : end if
2060 94 : if (dtset%npband*dtset%npfft*dtset%bandpp>1) dtset%paral_kgb=1
2061 : ! LinAlg parameters: we change values only if they are not present in input file
2062 94 : if (dtset%paral_kgb==1) then
2063 48 : if (tread(9)==0) dtset%use_slk=my_distp(8,icount)
2064 48 : if (tread(10)==0) dtset%np_slk=my_distp(9,icount)
2065 48 : if (tread(11)==0) dtset%gpu_linalg_limit=my_distp(10,icount)
2066 : end if
2067 : ! New definition of "world" MPI communicator
2068 94 : if (optdriver==RUNL_RESPFN.and.dtset%paral_rf==1) then
2069 20 : nproc1=max(nproc1,dtset%nsppol*dtset%nkpt) ! Take into account the code in respfn.F90
2070 20 : nproc1=min(nproc1,nproc)
2071 20 : nproc1=(nproc1/dtset%nppert)*dtset%nppert
2072 : end if
2073 94 : call initmpi_world(mpi_enreg,nproc1)
2074 : end if
2075 :
2076 : !Final advice in case max_ncpus > 0
2077 95 : if (max_ncpus>0.and.mcount>0) then
2078 1 : write(msg,'(6a)') ch10,&
2079 1 : ' Launch a parallel version of ABINIT with a distribution of processors among the above list,',ch10,&
2080 1 : ' and the associated input variables (np_spkpt, npband, npfft, bandpp, etc.).',ch10,&
2081 2 : ' The higher weight should be better.'
2082 1 : call wrtout(std_out,msg);if (max_ncpus>0) call wrtout(ab_out,msg)
2083 : end if
2084 :
2085 95 : if (mcount>0) then
2086 95 : ABI_FREE(isort)
2087 : end if
2088 95 : ABI_FREE(my_distp)
2089 95 : ABI_FREE(my_algo)
2090 95 : ABI_FREE(weight)
2091 :
2092 : !Final line
2093 95 : write(msg,'(a,100("="),2a)') " ",ch10,ch10
2094 95 : call wrtout(std_out,msg);if (max_ncpus>0) call wrtout(ab_out,msg)
2095 :
2096 : !max_ncpus requires a stop
2097 12983 : if (max_ncpus>0) then
2098 1 : iexit = iexit + 1 ! will stop in the parent.
2099 : end if
2100 :
2101 : DBG_EXIT("COLL")
2102 :
2103 : contains
2104 :
2105 815 : real(dp) pure function speedup_fdp(nn, mm)
2106 : ! Expected linear speedup for a nn-sized problem and mm processes
2107 : integer,intent(in) :: nn, mm
2108 1280 : speedup_fdp = (one*nn) / (one* ((nn / mm) + merge(0, 1, mod(nn, mm) == 0)))
2109 : end function speedup_fdp
2110 :
2111 : end subroutine finddistrproc
2112 : !!***
2113 :
2114 : !!****f* ABINIT/compute_kgb_indicator
2115 : !! NAME
2116 : !! compute_kgb_indicator
2117 : !!
2118 : !! FUNCTION
2119 : !! Only for "KGB" parallelism (LOBPCG algorithm for Ground-state):
2120 : !! Give an indicator of performance for a given distribution of processors
2121 : !! (npband, npfft and bandpp).
2122 : !! Determine best choice of parameters for Scalapack and/or Magma Linear Algebra routines.
2123 : !!
2124 : !! INPUTS
2125 : !! bandpp=internal lobpcg optimization variable
2126 : !! glb_comm=communicator for global MPI communications
2127 : !! mband=maximum number of bands.
2128 : !! mband=maximum number of plane waves
2129 : !! npband=number of processor 'band'
2130 : !! npfft = number of processor 'fft'
2131 : !! use_linalg_gpu=indicate if we also test the gpu linear algebra (compatible only with the legacy 2013 GPU code)
2132 : !!
2133 : !! OUTPUT
2134 : !! acc_kgb = indicator of performance
2135 : !! npslk = number of process to used in communicators
2136 : !!
2137 : !! SIDE EFFECTS
2138 : !! This routine can be used to find an indicator in order to refine automatic process distribution.
2139 : !! This indicator is returned in acc_kgb
2140 : !! This routine can be used to find the optimal values of np_slk parameter (ScaLapack)
2141 : !! and wheter or not we should use Magma for Linear Algebra in lobpcgwf
2142 : !!
2143 : !! SOURCE
2144 :
2145 0 : subroutine compute_kgb_indicator(acc_kgb,bandpp,glb_comm,mband,mpw,npband,npfft,npslk,use_linalg_gpu)
2146 :
2147 : use m_abi_linalg
2148 :
2149 : !Arguments ------------------------------------
2150 : !scalars
2151 : integer,intent(in) :: bandpp,glb_comm,mband,mpw,npband,npfft
2152 : integer,intent(inout) :: npslk,use_linalg_gpu
2153 : real(dp),intent(inout) :: acc_kgb
2154 :
2155 : !Local variables-------------------------------
2156 : !scalars
2157 : integer,parameter :: max_number_of_npslk=10,max_number_of_iter=10
2158 : integer :: blocksize,bigorder,ierr,ii,islk,islk1,iter,jj,keep_gpu
2159 : integer :: kgb_comm,my_rank,np_slk,np_slk_max,np_slk_best,nranks
2160 : integer :: use_lapack_gpu,use_slk,vectsize,wfoptalg
2161 : real(dp) :: min_eigen,min_ortho,time_xeigen,time_xortho
2162 : character(len=500) :: msg
2163 : !arrays
2164 0 : integer,allocatable :: ranks(:),val_npslk(:)
2165 0 : real(dp),allocatable :: eigen(:),grama(:,:),gramb(:,:)
2166 0 : complex(dp),allocatable :: blockvectorbx(:,:),blockvectorx(:,:),sqgram(:,:)
2167 : !******************************************************************
2168 :
2169 : DBG_ENTER("COLL")
2170 :
2171 : #ifdef DEBUG_MODE
2172 : write(msg,'(a,3i3)') 'compute_kgb_indicator : (bpp,npb,npf) = ', bandpp, npband, npfft
2173 : call wrtout(std_out,msg,'PERS')
2174 : #endif
2175 :
2176 : !Create local communicator for test
2177 : if (xmpi_paral==1) then
2178 0 : nranks=npfft*npband
2179 0 : ABI_MALLOC(ranks,(nranks))
2180 0 : ranks=(/((my_rank-1),my_rank=1,nranks)/)
2181 0 : kgb_comm=xmpi_subcomm(glb_comm,nranks,ranks,my_rank_in_group=my_rank)
2182 0 : ABI_FREE(ranks)
2183 : else
2184 : kgb_comm=xmpi_comm_self
2185 : my_rank=0
2186 : end if
2187 :
2188 : !Only for process in the new subgroup
2189 0 : if (my_rank/=xmpi_undefined) then
2190 :
2191 : ! We enforce vectsize >=blocksize (This is not true in lobpcg but
2192 : ! these are rare cases and this simplify the matrix constructions below...)
2193 0 : blocksize=npband*bandpp
2194 0 : vectsize=max(1+mpw/(npband*npfft),blocksize)
2195 0 : bigorder=3*blocksize
2196 :
2197 0 : ABI_MALLOC(blockvectorx,(vectsize,blocksize))
2198 0 : ABI_MALLOC(blockvectorbx,(vectsize,blocksize))
2199 0 : ABI_MALLOC(sqgram,(blocksize,blocksize))
2200 0 : ABI_MALLOC(grama,(2*bigorder,bigorder))
2201 0 : ABI_MALLOC(gramb,(2*bigorder,bigorder))
2202 0 : ABI_MALLOC(eigen,(bigorder))
2203 0 : ABI_MALLOC(val_npslk,(max_number_of_npslk)) ! not too much values tested
2204 :
2205 0 : min_eigen=greatest_real
2206 0 : min_ortho=greatest_real
2207 0 : np_slk_best=-1 ; np_slk_max=0
2208 : #ifdef HAVE_LINALG_SCALAPACK
2209 : np_slk_max=min(max_number_of_npslk,npband*npfft)
2210 : #endif
2211 :
2212 : ! Preselect a range of available np_slk values
2213 0 : val_npslk(1:)=0 ; val_npslk(2)=1
2214 : do islk=3,np_slk_max
2215 : np_slk=val_npslk(islk-1)*2
2216 : do while ((modulo(npband*npfft,np_slk)>0).and.(np_slk<(npband*npfft)))
2217 : np_slk=np_slk+1
2218 : end do
2219 : if(np_slk>(npband*npfft).or.np_slk>mband) exit
2220 : val_npslk(islk)=np_slk
2221 : end do
2222 0 : np_slk_max=islk-1
2223 :
2224 : ! Loop over np_slk values
2225 0 : islk1=1
2226 : #ifdef HAVE_LINALG_MAGMA
2227 : if (use_linalg_gpu==ABI_GPU_LEGACY) islk1=0
2228 : #endif
2229 0 : do islk=islk1,np_slk_max
2230 :
2231 0 : time_xortho=zero ; time_xeigen=zero
2232 :
2233 0 : use_slk=0
2234 : if (islk==0) then
2235 : ! This is the test for the GPU
2236 : use_lapack_gpu=1 ; np_slk=0
2237 : else
2238 0 : use_lapack_gpu=0 ; np_slk=val_npslk(islk)
2239 0 : if (np_slk>0) use_slk=1
2240 : end if
2241 :
2242 : ! Initialize linalg parameters for this np_slk value
2243 : ! For the first np_slk value, everything is initialized
2244 : ! For the following np_slk values, only Scalapack parameters are updated
2245 0 : wfoptalg=14 ! Simulate use of LOBPCG
2246 : call abi_linalg_init(bigorder,RUNL_GSTATE,wfoptalg,1,&
2247 0 : & use_lapack_gpu,use_slk,np_slk,kgb_comm)
2248 :
2249 : ! We could do mband/blocksize iter as in lobpcg but it's too long
2250 0 : do iter=1,max_number_of_iter
2251 :
2252 : ! Build matrixes
2253 0 : do ii=1,vectsize
2254 0 : do jj=1,blocksize
2255 0 : if (ii>jj) then
2256 0 : blockvectorx(ii,jj) =czero
2257 0 : blockvectorbx(ii,jj)=czero
2258 : else
2259 0 : blockvectorx(ii,jj) =cone
2260 0 : blockvectorbx(ii,jj)=cone
2261 : end if
2262 : end do
2263 : end do
2264 0 : grama=zero;gramb=zero
2265 0 : do jj=1,bigorder
2266 0 : do ii=jj,bigorder
2267 0 : if (ii==jj) then
2268 0 : grama(2*ii-1,jj)=one
2269 0 : gramb(2*ii-1,jj)=one
2270 : else
2271 0 : grama(2*ii-1:2*ii,jj)=one
2272 0 : grama(2*jj-1,ii)= one
2273 0 : grama(2*jj ,ii)=-one
2274 : end if
2275 : end do
2276 : end do
2277 :
2278 : ! Call to abi_xorthonormalize
2279 0 : time_xortho=time_xortho-abi_wtime()
2280 0 : call abi_xorthonormalize(blockvectorx,blockvectorbx,blocksize,kgb_comm,sqgram,vectsize)
2281 0 : time_xortho = time_xortho + abi_wtime()
2282 :
2283 : ! Call to abi_xhegv
2284 0 : time_xeigen=time_xeigen-abi_wtime()
2285 : call abi_xhegv(1,'v','u',bigorder,grama,bigorder,gramb,bigorder,eigen,&
2286 0 : & x_cplx=2,use_slk=use_slk,use_gpu_magma=use_lapack_gpu)
2287 0 : time_xeigen=time_xeigen+abi_wtime()
2288 :
2289 : end do ! iter
2290 :
2291 : ! Finalize linalg parameters for this np_slk value
2292 : ! For the last np_slk value, everything is finalized
2293 : ! For the previous np_slk values, only Scalapack parameters are updated
2294 0 : call abi_linalg_finalize(use_lapack_gpu)
2295 :
2296 0 : time_xortho= time_xortho*mband/blocksize
2297 0 : time_xeigen= time_xeigen*mband/blocksize
2298 : if (time_xortho<min_ortho) min_ortho=time_xortho
2299 0 : if (time_xeigen<min_eigen) then
2300 0 : min_eigen=time_xeigen
2301 0 : np_slk_best=np_slk
2302 0 : keep_gpu=use_lapack_gpu
2303 : end if
2304 :
2305 : end do ! np_slk
2306 :
2307 : #ifdef DEBUG_MODE
2308 : write(msg,'(2(a,es15.3),a,i3)') ' In the best case, xortho took ',min_ortho,&
2309 : ' and xeigen took ',min_eigen,' for np_slk=',np_slk_best
2310 : call wrtout(std_out,msg,'PERS')
2311 : #endif
2312 :
2313 : ! Final values to be sent to others process
2314 0 : acc_kgb=min_ortho+four*min_eigen
2315 0 : npslk=max(np_slk_best,1)
2316 0 : use_linalg_gpu=keep_gpu
2317 :
2318 0 : ABI_FREE(blockvectorx)
2319 0 : ABI_FREE(blockvectorbx)
2320 0 : ABI_FREE(sqgram)
2321 0 : ABI_FREE(grama)
2322 0 : ABI_FREE(gramb)
2323 0 : ABI_FREE(eigen)
2324 0 : ABI_FREE(val_npslk)
2325 :
2326 : end if ! my_rank in group
2327 :
2328 : !Free local MPI communicator
2329 0 : call xmpi_comm_free(kgb_comm)
2330 :
2331 : !Broadcast of results to be sure every process has them
2332 0 : call xmpi_bcast(acc_kgb,0,glb_comm,ierr)
2333 0 : call xmpi_bcast(npslk,0,glb_comm,ierr)
2334 0 : call xmpi_bcast(use_linalg_gpu,0,glb_comm,ierr)
2335 :
2336 : #ifndef DEBUG_MODE
2337 : ABI_UNUSED(msg)
2338 : #endif
2339 :
2340 : DBG_EXIT("COLL")
2341 :
2342 0 : end subroutine compute_kgb_indicator
2343 : !!***
2344 :
2345 : end module m_mpi_setup
2346 : !!***
|
